text:
using Phonetics
using Test

#Tests for correctness of Caverphone
@test caverphone("Lee") ==  "l11111"
@test caverphone("Thompson") ==  "tmpsn1"

##Tests for correctness of Soundex
@test soundex("Euler") == "e460"
@test soundex("Ellery") == "e460"
@test soundex("Gauss") == "g200"
@test soundex("Ghosh") == "g200"
@test soundex("Hilbert") ==  "h416"
@test soundex("Heilbronn") ==  "h416"
@test soundex("Kant") ==  "k530"
@test soundex("Knuth") ==  "k530"
@test soundex("Ladd") == "l300"
@test soundex("Lloyd") == "l300"
@test soundex("Christina") ==  "c623"
@test soundex("Kristina") ==  "k623"

#Tests for correctness of Metaphone
@test metaphone("School") ==  "skl"
@test metaphone("Shubert") ==  "xbrt"
@test metaphone("Bonner") ==  "bnr"
@test metaphone("Baymore") ==  "bmr"
@test metaphone("Smith") ==  "sm0"
@test metaphone("Saneed") ==  "snt"
@test metaphone("Aardvark") ==  "artf"
@test metaphone("persuade") ==  "prst"
@test metaphone("pressed") ==  "prst"

#NB: reference seemed confused about these two examples.
@test metaphone("Van Hoesen") ==  "fnhs"
@test metaphone("Vincenzo") ==  "fnsn"

#Tests for correctness of Phonex
@test phonex("Peter") ==  "b360"
@test phonex("Pete") ==  "b300"
@test phonex("Pedro") ==  "b360"
@test phonex("Stephen") ==  "s315"
@test phonex("Steve") ==  "s310"
@test phonex("Smythe") ==  "s530"
@test phonex("Smith") ==  "s530"
@test phonex("Gail") ==  "g400"
@test phonex("Gayle") ==  "g400"
@test phonex("Christina") ==  "c623"
@test phonex("Kristina") ==  "c623"

#Test for correctnes of Phonix
@test phonix("Peter") ==  "p300"
@test phonix("Pete") ==  "p300"
@test phonix("Pedro") ==  "p360"
@test phonix("Stephen") ==  "s375"
@test phonix("Steve") ==  "s370"
@test phonix("Smith") ==  "s530"
@test phonix("Smythe") ==  "s530"
@test phonix("Gail") ==  "g400"
@test phonix("Gayle") ==  "g400"
@test phonix("Christina") ==  "k683"
@test phonix("Kristina") ==  "k683"

#Test for correctness of NYSIIS
@test nysiis("Peter") ==  "patar"
@test nysiis("Pete") ==  "pat"
@test nysiis("Pedro") ==  "padr"
@test nysiis("Stephen") ==  "stafan"
@test nysiis("Steve") ==  "staf"
@test nysiis("Smith") ==  "snat"
@test nysiis("Smythe") ==  "snyt"
@test nysiis("Gail") ==  "gal"
@test nysiis("Gayle") ==  "gayl"
@test nysiis("Christina") ==  "chrast"
@test nysiis("Kristina") ==  "crasta"
@test nysiis("Lawson") ==  "lasan"
@test nysiis("Greene") ==  "gran"
@test nysiis("Lynch") ==  "lync"
@test nysiis("Wheeler") ==  "whalar"
@test nysiis("Mitchell") ==  "matcal"
@test nysiis("Anderson") ==  "andars"
@test nysiis("Truman") ==  "tranan"
@test nysiis("Jellyfish") ==  "jalyf"

#Test for correctness of Double-Metaphone
@test double_metaphone("Peter") == "ptr"
@test double_metaphone("Pete") ==  "pt"
@test double_metaphone("Pedro") ==  "ptr"
@test double_metaphone("Stephen") ==  "stfn"
@test double_metaphone("Steve") ==  "stf"
@test double_metaphone("Smith") ==  ["sm0","xmt"]
@test double_metaphone("Smythe") ==  ["sm0","xmt"]
@test double_metaphone("Gail") ==  "kl"
@test double_metaphone("Gayle") ==  "kl"
@test double_metaphone("Christina") ==  "krstn"
@test double_metaphone("Kristina") ==  "krstn"
@test double_metaphone("Whithers") ==  ["a0rs","atrs"]
@test double_metaphone("Wasserman") ==  ["asrmn","fsrmn"]
@test double_metaphone("Vasserman") ==  "fsrmn"
@test double_metaphone("Arnoff") ==  "arnf"
@test double_metaphone("Arnow") ==  ["arn","arnf"]
@test double_metaphone("Filipowicz") == ["flpts", "flpfx"]
@test double_metaphone("Wrainwright") ==  "rnrt"
@test double_metaphone("Hugh") ==  "h"
@test double_metaphone("Loughbridge") ==  "lfbrj"
@test double_metaphone("Edge") ==  "aj"
@test double_metaphone("Edgar") ==  "atkr"
@test double_metaphone("Pepsi") ==  "pps"
@test double_metaphone("Phillipa") ==  "flp"
@test double_metaphone("Xavier") ==  ["sf","sfr"]
@test double_metaphone("Hochmeier") ==  "hkmr"
@test double_metaphone("Caeser") ==  "ssr"
@test double_metaphone("Chianti") ==  "knt"
@test double_metaphone("Michael") ==  ["mkl","mxl"]
@test double_metaphone("Chemise") ==  "kms"
@test double_metaphone("McHugh") ==  "mkh"
@test double_metaphone("Czerny") ==  ["srn","xrn"]
@test double_metaphone("Focaccia") ==  ["fxx","fkx"]
@test double_metaphone("Bellaccio") ==  "blx"
@test double_metaphone("Bacchus") ==  "bkhs"
@test double_metaphone("Accident") ==  "akstnt"
@test double_metaphone("Succeed") ==  "skst"
@test double_metaphone("Belacqua") ==  "blk"
@test double_metaphone("Francia") ==  ["frns","frnx"]
@test double_metaphone("Mac Gregor") ==  "mkrkr"
@test double_metaphone("Mac Caffrey") ==  "mkfr"
@test double_metaphone("McLaughin") ==  "mklfn"
@test double_metaphone("Tagliaro") ==  ["tklr", "tlr"]
@test double_metaphone("Gelman") ==  ["klmn", "jlmn"]
@test double_metaphone("Knowles") ==  ["nls","nfls"]
@test double_metaphone("Jose") ==  "hs"
@test double_metaphone("Yankelofich") ==  "anklfx"
@test double_metaphone("Jankelowicz") ==  ["jnklts","anklfx"]
@test double_metaphone("Cabrillo") ==  ["kbrl","kbr"]

#Test for correctness of MRA
@test match_rating_encode("Smith") ==  "smth"
@test match_rating_encode("Smythe") ==  "smyth"
@test match_rating_encode("Christina") ==  "chrstn"
@test match_rating_encode("Kristina") ==  "krstn"
@test match_rating_encode("Catherine") ==  "cthrn"
@test match_rating_encode("Kathryn") ==  "kthryn"
@test match_rating_encode("Byrne") ==  "byrn"
@test match_rating_encode("Boern") ==  "brn"

@test match_rating("Byrne","Boern") == 5
@test meets_match_rating("Byrne","Boern")
@test match_rating("Smith","Smythe") == 5
@test meets_match_rating("Smith","Smythe")
@test match_rating("Catherine","Kathryn") ==  4
@test meets_match_rating("Catherine","Kathryn")

#Test for correctness of Fuzzy Soundex
@test fuzzy_soundex("Kristen") == "k6935"
@test fuzzy_soundex("Krissy") == "k6900"
@test fuzzy_soundex("Christen") == "k6935"

#Test some examples of code_similarity
@test code_similarity("Kristina", "Christina") == 1
@test code_similarity("Kristina", "Kristian") == 0.6
@test code_similarity("Kristian", "Kristina", phonix) == 1

#Test some examples of code_match
helpstrings = ["Halp", "Elf", "Hulk", "Heelp","Half", "Abba", "Any"]
@test code_match("Help", helpstrings) == ["Halp", "Heelp", "Half"]
@test code_match("Help", helpstrings, phonex) == ["Halp", "Elf", "Heelp", "Half", "Abba"]
@test code_match("Help", helpstrings, metaphone) == ["Halp", "Heelp"]
@test code_match("Help", helpstrings, double_metaphone) == ["Halp", "Heelp"]
@test code_match("Help", helpstrings, metaphone, 0.5) == ["Halp", "Hulk", "Heelp", "Half"]
@test code_match("Help", helpstrings, phonex, 0.1) == helpstrings

#Test some examples of code_cluster
strings = ["Sing", "Sink", "Song", "Sunk", "Sinking", "Singing", "Single"]
ccs1 = code_cluster(strings, phonix, 0.7, 0.9, false)
@test ccs1[1] == ["Sing","Sink","Song","Sunk","Sinking","Singing"]
@test ccs1[2] == ["Sinking","Singing"]
@test ccs1[3] == ["Single"]
@test length(ccs1) == 3
ccs2 = code_cluster(strings, phonix, 1, 1, false)
@test ccs2[1] == ["Sing","Sink","Song","Sunk"]
@test ccs2[2] == ["Sinking","Singing"]
@test ccs2[3] == ["Single"]
@test length(ccs2) == 3

#Test some examples for editex
@test editex("Hello", "Hullo") == 1
@test editex("Hellophant", "Hullo") == 10
@test editex("ello", "Hola") == 4
@test editex("lo", "Hola") == 6
@test editex("Singings", "Cingingz") == 2
@test editex("Singings", "Singjings") == 1
@test editex("Singings", "Singlings") == 2

#Test examples of syllable counts
@test syllable_count("syllable") == 3
@test syllable_count("syllables") == 3
@test syllable_count("Mary had a little lamb") == 7
@test syllable_count("calibre") == 3
@test syllable_count("sabre-toothed") == 4
@test syllable_count("are") == 1
@test syllable_count("area") == 2
@test syllable_count("underlies") == 3

#Test for correctness of Roger Root
@test roger_root("Johnson") == "32020"
@test roger_root("Williams") == "45300"
@test roger_root("Smith") == "00310"
@test roger_root("Jones") == "32000"
@test roger_root("Brown") == "09420"
@test roger_root("Davis") == "01800"
@test roger_root("Jackson") == "37020"
@test roger_root("Wilson") == "45020"
@test roger_root("Lee") == "05000"
@test roger_root("Thomas") == "01300"

@test spoken_length("Mary had a little lamb") ≈ 2.52064687529
@test spoken_length("Mary hatte ein kleines lamm", "de") ≈ 2.1091149688610002


text:
cd(expanduser("~/dev/julia/"))
using Pkg
Pkg.activate("wiktionary")
## using FilingForest
## using ParserAlchemy
using OrgParser
using ParserAlchemy.Tokens
using WikitextParser
using Test
#include(joinpath(dirname(pathof(FilingForest)),"wiktionary.jl"))
import WikitextParser: wiki_freelink, wiki_external_link, wiki_template

@testset "wikitext links" begin
    @test tokenize(wiki_freelink, "[[Test|1]]") ==
        Token(Symbol("Test"), "1")

    @test tokenize(wiki_freelink, "[[Test]]") ==
        Token(Symbol("wiktionary.de"), "Test")

    @test tokenize(wiki_freelink, "[[Texas|Lone Star State]]") ==
        Token(Symbol("Texas"), "Lone Star State")
    # @test tokenize(wiki_freelink, "[[Texas#Racism|Lone Star State]]") ==
    # (target="Texas",section="Racism",title="Lone Star State")
    tokenize(wiki_freelink, "[[File:wiki.png|thumb|Wikipedia logo]]")

    @test tokenize(wiki_external_link,
                   "[http://www.bictxt.de/Schule]") ==
                       Token("http://www.bictxt.de/Schule","")
    ## todo: bare urls
end
import WikitextParser: wikitext, wikitextParser

import WikitextParser: expand_numbers
@testset "number references" begin
    @test tokenize(expand_numbers,"[1]") == ["1"]
    @test tokenize(expand_numbers,"[1-3]") == ["1","2","3"]
    @test tokenize(expand_numbers,"[1-3,10]") == ["1","2","3","10"]
end

tokenize(wikitext,
         "<
ref name=
\"
Georges
\"
>
{{Ref-Georges}}
<
/ref
>")

t = ":„Die Monatsnamen werden stark flektiert. Besonderheiten: Das Dativ-e tritt nicht mehr auf: ''im Januar'' (nicht: ''im Januare''). Auch der Genitiv ist im heutigen Sprachgebrauch häufig ohne Endung: ''des Januar[s],'' […]; […], ''des 12. Januar[s].'' […] Allgemein gilt: Die endungs- und artikellose Form des Monatsnamens steht vor allem dann, wenn ein Substantiv vorangeht: ''Anfang Januar,'' […]. Auch wenn die Monatsnamen als Apposition bei Monat stehen, bleiben sie ohne Flexionsendung: ''des Monats Januar.''“"
using BenchmarkTools
##@btime tokenize(wikitextParser, t)

@testset "wikitext templates" begin
    @test tokenize(wikitextParser, "{{test}}")[1] == Line(Token[Token(:whitespace,"")],LineContent[Template("test")])
    x = tokenize(wiki_template(), "{{Ü-Tabelle|Ü-links|Ü-rechts}}")
end

x= tokenize(rep(wikitext),
            "token [http://bictxt.de] {{IPA|\n[[Test]]}}")

tokenize(rep(wikitext),
         "token [http://bictxt.de] {{Pl.}} ")

import WikitextParser: template_inner

tokenize(template_inner, "-\n")
tokenize(wikitext, " -") |> dump

tokenize(wikitext, "([[Reich]]e)")

## todo:
tokenize(wikitext, "[[Kategorie][ Afrikaans]]")
# while !isempty(inner_templates.els)
#     pop!(inner_templates.els)
# end
# for t in  inner_template_names
#     push!(inner_templates.els, t)
# end

## this hack makes the parser recursive
ParserAlchemy.result_type(wikitext)
import WikitextParser: wiki_template
ParserAlchemy.result_type(wiki_template(nothing))

## @test tokenize(heading(2), "== test ==\n") ==
##    (level =2, title=[l"test", delim" "])
##    (indent=[Token(:headline, "2")], tokens=[l"test", delim" "])

import WikitextParser: template_inner
x = tokenize(template_inner,
":[1] Königtum, (königliche) Herrschaft, Regentschaft, Monarchie
:[2] Königreich
")
x[1] |> dump
x = tokenize(seq(
            wiki_template("Bedeutungen"),emptyline,
            template_inner;
            transform = 3),
"{{Bedeutungen}}
:[1] Königtum, (königliche) Herrschaft, Regentschaft, Monarchie
:[2] Königreich
")

import WikitextParser: wiki_lines
tokenize(alternate(wiki_lines, emptyline),
         ":{{IPA}} {{Lautschrift|haˈloː}}\n")

tokenize(
    seq(
            wiki_template("Aussprache"),emptyline,
            alternate(wiki_lines, emptyline);
            transform = 3, partial=true),
"{{Aussprache}}
:{{IPA}} {{Lautschrift|haˈloː}}
:{{Hörbeispiele}} {{Audio|De-Hallo.ogg}}
:{{Reime}} {{Reim|oː|Deutsch}}
")

tokenize(seq(
            wiki_template("Aussprache"),emptyline,
            alternate(wiki_lines, emptyline);
            transform = 3),
         "{{Aussprache}}
:{{IPA}} {{Lautschrift|ˈspiːʃiːz}}, {{Pl.}} {{Lautschrift|ˈspiːʃiːz}}
:{{Hörbeispiele}} {{Audio|En-us-species.ogg|spr=us}}")

x= tokenize(seq(
            wiki_template("Beispiele"),emptyline,
            alternate(wiki_lines, emptyline);
    transform = 3),
            "{{Beispiele}}
:[1] {{Beispiele fehlen|spr=la}}
")

import WikitextParser: heading
tokenize(heading(2), "== species ==")
tokenize(heading(2),
         "== species ({{Sprache|Englisch}}) ==")

tokenize(heading(3),
         "=== {{Wortart|Substantiv|Deutsch}}, {{f}} ===")

tokenize(wiki_freelink, "[[File:Backhandspring.ogv|miniatur|[9] vom [[turnen|Turnen]] [[kennen]] wir die ''Figur'' [[Flickflack]]]]")

familia = """
{{Siehe auch|[[família]]}}
== familia ({{Sprache|Latein}}) ==
=== {{Wortart|Substantiv|Latein}}, {{f}} ===

{{Latein Substantiv Übersicht|familia|famili|ae}}

{{Worttrennung}}
:fa·mi·lia, {{Pl.}} fa·mi·li·ae

{{Grammatische Merkmale}}
:[[Genitiv]] Singular nach [[pater]], [[mater]], [[filius]], [[filia]]: ''familiās''<ref name="Georges">{{Ref-Georges}}</ref>

{{Bedeutungen}}
:[1] das [[Haus]], der [[Hausstand]] (mit Sklaven)
:[2] das [[Geschlecht]], die [[Familie]]
:[3] die [[Dienerschaft]]

{{Herkunft}}
:zu dem Adjektiv ''{{Ü|la|famulus}}''<ref name="Georges" /> „dienend“; ursprüngliche Bedeutung: die Gesamtheit der unter einem ''[[dominus]]'' stehenden Sklaven, besonders als Angehörige und Teil der Familie, dann übertragen die ganze Hausgenossenschaft (Freie und Sklaven)

{{Beispiele}}
:[1] {{Beispiele fehlen|spr=la}}

==== {{Übersetzungen}} ====
{{Ü-Tabelle|Ü-links=
:{{Übersetzungen umleiten|1|Hausstand|}} {{m}}, {{Übersetzungen umleiten||Haus|}} {{m}}
:{{Übersetzungen umleiten|2|Geschlecht|}}, {{Übersetzungen umleiten||Familie|}}
:{{Übersetzungen umleiten|3|Dienerschaft|}}
|Ü-rechts=
}}

{{Referenzen}}
:[1–3] {{Wikipedia|spr=la|familia}}
:[1–3] {{Ref-Georges}}
:[1–3] {{Lit-Stowasser: Lateinisch-deutsches Schulwörterbuch|A=1994}}, Seite 202, Eintrag „familia“

{{Quellen}}
"""

r=tokenize(wikitextParser, familia)
v = tokenize(wiktionary_defs,r, delta=3)
v[1].defs[1][2]
v = wiki_meaning(v[1])[1][2]

using WikitextParser: number_line, wikitextParser
tokenize(alt(wiki_lines...),":[1] ''kurz für:'' [[Schraubenmutter]], welche das Gegenstück zu der [[Schraube]] bildet\n")

mutter = read("WikitextParser/test/mutter.wikitext",String);
r=tokenize(wikitextParser, mutter)
v = tokenize(wiktionary_defs,r, delta=3)
m = wiki_meaning(v[1])
typeof(m[1][2])
using StructArrays
using IndexedTables
table(StructArray( m[1][2] ))
v[1].defs[1][2]
 wiki_meaning(v[1]) |> typeof
v

ext: jl


text:
"""
    FluxSummary

A struct used to store summary information about the solution
of a constraint based analysis result.
"""
struct FluxSummary
    biomass_fluxes::OrderedDict{String,Float64}
    import_fluxes::OrderedDict{String,Float64}
    export_fluxes::OrderedDict{String,Float64}
    unbounded_fluxes::OrderedDict{String,Float64}
end

"""
    FluxSummary()

A default empty constructor for `FluxSummary`.
"""
function FluxSummary()
    FluxSummary(
        OrderedDict{String,Float64}(),
        OrderedDict{String,Float64}(),
        OrderedDict{String,Float64}(),
        OrderedDict{String,Float64}(),
    )
end

"""
    flux_summary(flux_result::Dict{String, Float64};
        exclude_exchanges = false,
        exchange_prefixes = _constants.exchange_prefixes,
        biomass_strings = _constants.biomass_strings,
        exclude_biomass = false,
        small_flux_bound = 1.0/_constants.default_reaction_bound^2,
        large_flux_bound = _constants.default_reaction_bound,
        keep_unbounded = false,
    )::FluxSummary

Summarize a dictionary of fluxes into small, useful representation of the most
important information contained. Useful for pretty-printing and quickly
exploring the results. Internally this function uses
[`looks_like_biomass_reaction`](@ref) and
[`looks_like_exchange_reaction`](@ref). The corresponding keyword arguments
passed to these functions. Use this if your model has non-standard ids for
reactions. Fluxes smaller than `small_flux_bound` are not stored, while fluxes
larger than `large_flux_bound` are only stored if `keep_unbounded` is `true`.

# Example
```
julia> sol = flux_dict(flux_balance_analysis(model, Tulip.Optimizer))
julia> fr = flux_summary(sol)
Biomass:
  BIOMASS_Ecoli_core_w_GAM: 0.8739
Import:
  EX_o2_e:     -21.7995
  EX_glc__D_e: -10.0
  EX_nh4_e:    -4.7653
  EX_pi_e:     -3.2149
Export:
  EX_h_e:      17.5309
  EX_co2_e:    22.8098
  EX_h2o_e:    29.1758
```
"""
function flux_summary(
    flux_result::Maybe{Dict{String,Float64}};
    exclude_exchanges = false,
    exchange_prefixes = _constants.exchange_prefixes,
    biomass_strings = _constants.biomass_strings,
    exclude_biomass = false,
    small_flux_bound = 1.0 / _constants.default_reaction_bound^2,
    large_flux_bound = _constants.default_reaction_bound,
    keep_unbounded = false,
)
    isnothing(flux_result) && return FluxSummary()

    rxn_ids = collect(keys(flux_result))
    ex_rxns = filter(
        x -> looks_like_exchange_reaction(
            x,
            exclude_biomass = exclude_biomass,
            biomass_strings = biomass_strings,
            exchange_prefixes = exchange_prefixes,
        ),
        rxn_ids,
    )
    bmasses = filter(
        x -> looks_like_biomass_reaction(
            x;
            exclude_exchanges = exclude_exchanges,
            exchange_prefixes = exchange_prefixes,
            biomass_strings = biomass_strings,
        ),
        rxn_ids,
    )

    ex_fluxes = [flux_result[k] for k in ex_rxns]
    bmass_fluxes = [flux_result[k] for k in bmasses]

    idx_srt_fluxes = sortperm(ex_fluxes)
    import_fluxes = [
        idx for
        idx in idx_srt_fluxes if -large_flux_bound < ex_fluxes[idx] <= -small_flux_bound
    ]
    export_fluxes = [
        idx for
        idx in idx_srt_fluxes if small_flux_bound < ex_fluxes[idx] <= large_flux_bound
    ]

    if keep_unbounded
        lower_unbounded =
            [idx for idx in idx_srt_fluxes if ex_fluxes[idx] <= -large_flux_bound]
        upper_unbounded =
            [idx for idx in idx_srt_fluxes if ex_fluxes[idx] >= large_flux_bound]
        return FluxSummary(
            OrderedDict(k => v for (k, v) in zip(bmasses, bmass_fluxes)),
            OrderedDict(
                k => v for (k, v) in zip(ex_rxns[import_fluxes], ex_fluxes[import_fluxes])
            ),
            OrderedDict(
                k => v for (k, v) in zip(ex_rxns[export_fluxes], ex_fluxes[export_fluxes])
            ),
            OrderedDict(
                k => v for (k, v) in zip(
                    [ex_rxns[lower_unbounded]; ex_rxns[upper_unbounded]],
                    [ex_fluxes[lower_unbounded]; ex_fluxes[upper_unbounded]],
                )
            ),
        )
    else
        return FluxSummary(
            OrderedDict(k => v for (k, v) in zip(bmasses, bmass_fluxes)),
            OrderedDict(
                k => v for (k, v) in zip(ex_rxns[import_fluxes], ex_fluxes[import_fluxes])
            ),
            OrderedDict(
                k => v for (k, v) in zip(ex_rxns[export_fluxes], ex_fluxes[export_fluxes])
            ),
            OrderedDict{String,Float64}(),
        )
    end
end


text:
using Oceananigans.BoundaryConditions: FieldBoundaryConditions, regularize_field_boundary_conditions

#####
##### Velocity fields tuples
#####

"""
    VelocityFields(arch, grid, user_bcs = NamedTuple())

Return a `NamedTuple` with fields `u`, `v`, `w` initialized on the architecture `arch`
and `grid`. Boundary conditions `bcs` may be specified via a named tuple of
`FieldBoundaryCondition`s.
"""
function VelocityFields(arch, grid, user_bcs = NamedTuple())

    template = FieldBoundaryConditions()

    default_bcs = (
        u = regularize_field_boundary_conditions(template, grid, :u),
        v = regularize_field_boundary_conditions(template, grid, :v),
        w = regularize_field_boundary_conditions(template, grid, :w)
    )

    bcs = merge(default_bcs, user_bcs)

    u = XFaceField(arch, grid, bcs.u)
    v = YFaceField(arch, grid, bcs.v)
    w = ZFaceField(arch, grid, bcs.w)

    return (u=u, v=v, w=w)
end

#####
##### Tracer fields tuples
#####

"""
    TracerFields(tracer_names, arch, grid, user_bcs)

Return a `NamedTuple` with tracer fields specified by `tracer_names` initialized as
`CenterField`s on the architecture `arch` and `grid`. Boundary conditions `user_bcs`
may be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function TracerFields(tracer_names, arch, grid, user_bcs)
    default_bcs = NamedTuple(name => FieldBoundaryConditions(grid, (Center, Center, Center)) for name in tracer_names)
    bcs = merge(default_bcs, user_bcs) # provided bcs overwrite defaults
    return NamedTuple(c => CenterField(arch, grid, bcs[c]) for c in tracer_names)
end

"""
    TracerFields(tracer_names, arch, grid; kwargs...)

Return a `NamedTuple` with tracer fields specified by `tracer_names` initialized as
`CenterField`s on the architecture `arch` and `grid`. Fields may be passed via optional
keyword arguments `kwargs` for each field.

This function is used by `OutputWriters.Checkpointer` and `TendencyFields`.
```
"""
TracerFields(tracer_names, arch, grid; kwargs...) =
    NamedTuple(c => c ∈ keys(kwargs) ? kwargs[c] : CenterField(arch, grid) for c in tracer_names)

# 'Nothing', or empty tracer fields
TracerFields(::Union{Tuple{}, Nothing}, arch, grid, bcs) = NamedTuple()

"Shortcut constructor for empty tracer fields."
TracerFields(::NamedTuple{(), Tuple{}}, arch, grid, bcs) = NamedTuple()

#####
##### Pressure fields tuples
#####

"""
    PressureFields(arch, grid, bcs::NamedTuple)

Return a `NamedTuple` with pressure fields `pHY′` and `pNHS` initialized as
`CenterField`s on the architecture `arch` and `grid`.  Boundary conditions `bcs`
may be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function PressureFields(arch, grid, bcs=NamedTuple())

    default_pressure_boundary_conditions =
        (pHY′ = FieldBoundaryConditions(grid, (Center, Center, Center)),
         pNHS = FieldBoundaryConditions(grid, (Center, Center, Center)))

    bcs = merge(default_pressure_boundary_conditions, bcs)

    pHY′ = CenterField(arch, grid, bcs.pHY′)
    pNHS = CenterField(arch, grid, bcs.pNHS)

    return (pHY′=pHY′, pNHS=pNHS)
end

#=
function PressureFields(arch, grid::AbstractGrid{<:Any, <:Any, <:Any, <:Flat}, bcs=NamedTuple())
    default_pressure_boundary_conditions =
        (pHY′ = FieldBoundaryConditions(grid, (Center, Center, Center)),
         pNHS = FieldBoundaryConditions(grid, (Center, Center, Center)))

    bcs = merge(default_pressure_boundary_conditions, bcs)
    pNHS = CenterField(arch, grid, bcs.pNHS)

    return (; pNHS=pNHS)
end
=#

"""
    TendencyFields(arch, grid, tracer_names;
                   u = XFaceField(arch, grid),
                   v = YFaceField(arch, grid),
                   w = ZFaceField(arch, grid),
                   kwargs...)

Return a `NamedTuple` with tendencies for all solution fields (velocity fields and
tracer fields), initialized on the architecture `arch` and `grid`. Optional `kwargs`
can be specified to assign data arrays to each tendency field.
"""
function TendencyFields(arch, grid, tracer_names;
                        u = XFaceField(arch, grid),
                        v = YFaceField(arch, grid),
                        w = ZFaceField(arch, grid),
                        kwargs...)

    velocities = (u=u, v=v, w=w)

    tracers = TracerFields(tracer_names, arch, grid; kwargs...)

    return merge(velocities, tracers)
end

#####
##### Helper functions for NonhydrostaticModel constructor
#####

VelocityFields(::Nothing, arch, grid, bcs) = VelocityFields(arch, grid, bcs)
PressureFields(::Nothing, arch, grid, bcs) = PressureFields(arch, grid, bcs)

"""
    VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, arch, grid, bcs)

Return a `NamedTuple` of velocity fields, overwriting boundary conditions
in `proposed_velocities` with corresponding fields in the `NamedTuple` `bcs`.
"""
function VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, arch, grid, bcs)

    validate_field_tuple_grid("velocities", proposed_velocities, grid)

    u = XFaceField(arch, grid, bcs.u, proposed_velocities.u.data)
    v = YFaceField(arch, grid, bcs.v, proposed_velocities.v.data)
    w = ZFaceField(arch, grid, bcs.w, proposed_velocities.w.data)

    return (u=u, v=v, w=w)
end

"""
    TracerFields(proposed_tracers::NamedTuple, arch, grid, bcs)

Return a `NamedTuple` of tracers, overwriting boundary conditions
in `proposed_tracers` with corresponding fields in the `NamedTuple` `bcs`.
"""
function TracerFields(proposed_tracers::NamedTuple, arch, grid, bcs)

    validate_field_tuple_grid("tracers", proposed_tracers, grid)

    tracer_names = propertynames(proposed_tracers)
    tracer_fields = Tuple(CenterField(arch, grid, bcs[c], proposed_tracers[c].data) for c in tracer_names)

    return NamedTuple{tracer_names}(tracer_fields)
end

"""
    PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, arch, grid, bcs)

Return a `NamedTuple` of pressure fields with, overwriting boundary conditions
in `proposed_tracer_fields` with corresponding fields in the `NamedTuple` `bcs`.
"""
function PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, arch, grid, bcs)
    validate_field_tuple_grid("pressures", proposed_pressures, grid)

    pHY′ = CenterField(arch, grid, bcs.pHY′, proposed_pressures.pHY′.data)
    pNHS = CenterField(arch, grid, bcs.pNHS, proposed_pressures.pNHS.data)

    return (pHY′=pHY′, pNHS=pNHS)
end


text:
#-----------------------------------------------------------------------# IndexedTable
Tables.istable(::Type{<:IndexedTable}) = true

Tables.materializer(t::IndexedTable) = table

Tables.columnaccess(::Type{<:IndexedTable}) = true
Tables.columns(t::IndexedTable) = columns(t)

Tables.rowaccess(::Type{<:IndexedTable}) = true
Tables.rows(t::IndexedTable) = rows(t)

Tables.schema(t::IndexedTable) = Tables.Schema(eltype(t))

# table(x; copy=false, kw...) = table(Tables.columntable(x); copy=copy, kw...)


text:
"""
    create_agents(graph, config)

Create agents from a graph with specified configuration.

# Arguments
- `graph`: A network as created by create_network
- `config`: Config object as provided by Config

See also: [Agent](@ref), [create_network](@ref), [generate_opinion](@ref), [generate_inclin_interact](@ref), [generate_check_regularity](@ref)
"""
function create_agents(
    graph::AbstractGraph,
    config::Config
)

    agent_list = Array{Agent, 1}(undef, length(vertices(graph)))

    for agent_idx in 1:length(agent_list)
        agent_list[agent_idx] = Agent(
            agent_idx,
            generate_opinion(),
            generate_inclin_interact(),
            generate_check_regularity(),
            generate_desired_input_count(
                config.agent_props.mean_desired_input_count
            )
        )
    end

    return agent_list

end

"""
    create_agents(config)

Create agents with specified configuration.

# Arguments
- `config`: Config object as provided by Config

See also: [Agent](@ref), [create_network](@ref), [generate_opinion](@ref), [generate_inclin_interact](@ref), [generate_check_regularity](@ref)
"""
function create_agents(
    config::Config
)

    agent_list = Array{Agent, 1}(undef, config.network.agent_count)

    for agent_idx in 1:length(agent_list)
        agent_list[agent_idx] = Agent(
            agent_idx,
            generate_opinion(),
            generate_inclin_interact(),
            generate_check_regularity(),
            generate_desired_input_count(
                config.agent_props.mean_desired_input_count
            )
        )
    end

    return agent_list

end

# suppress output of include()
;


text:
using Test
using LowRankModels, Random, SparseArrays

# Check that sparse algorithm converges to the same solution
# from the same initial conditions for simple pca.
m,n,k = 100,100,3
A = randn(m,k)*randn(k,n)
loss = QuadLoss()
r = ZeroReg()

glrm_1 = GLRM(A,loss,r,r,k) # solve with prox algorithm
glrm_2 = deepcopy(glrm_1)   # solve with sparse prox algorithm

X1,Y1,ch1 = fit!(glrm_1,ProxGradParams())
X2,Y2,ch2 = fit!(glrm_2,SparseProxGradParams())

@test X1 ≈ X2
@test Y1 ≈ Y2

# Check that the sparsity pattern in the data is correctly identified.
A = sprand(m,n,0.5)

glrm = GLRM(A,loss,r,r,k) # create glrm from sparse matrix

@test length(glrm.observed_features) == m
@test length(glrm.observed_examples) == n

for i = 1:m
	for j = 1:n
		if j in glrm.observed_features[i]
			@test A[i,j] != 0.0
		else
			@test A[i,j] == 0.0
		end
	end
end

for j = 1:n
	for i = 1:m
		if i in glrm.observed_examples[j]
			@test A[i,j] != 0.0
		else
			@test A[i,j] == 0.0
		end
	end
end


text:
function master_formulation(rd::RoutingData, type::FlowFormulation)
    m = _create_model(rd)
    rm = basic_routing_model_unitary(m, rd, type) # Includes flow-conservation constraints.

    @variable(m, mu >= 0)
    @objective(m, Min, mu)

    @variable(m, dual_alpha[e in edges(rd), d in demands(rd), v in vertices(rd)])
    @variable(m, dual_beta[e in edges(rd), e2 in edges(rd)] >= 0)

    # Robust uncertainty sets (one per edge).
    if rd.model_robust_reformulation_traffic_matrices
        dual_constraints = Dict{Edge{Int}, Dict{Edge{Int}, ConstraintRef}}(
            e => Dict{Edge{Int}, ConstraintRef}() for e in edges(rd)
        ) # Edge -> demand -> constraint reference.
    end
    for e in edges(rd)
        for d in demands(rd)
            # Dual constraint from the flow variables.
            for e2 in edges(rd)
                @constraint(m, dual_beta[e, e2] + dual_alpha[e, d, dst(e2)] - dual_alpha[e, d, src(e2)] >= 0)
            end

            # Dual constraint from the demand variables.
            c = @constraint(m, dual_alpha[e, d, src(d)] - dual_alpha[e, d, dst(d)] >= rm.routing[d, e])
            if rd.model_robust_reformulation_traffic_matrices
                dual_constraints[e][d] = c
            end
        end

        # Relate the main decision variables to the uncertainty sets.
        lhs = sum(dual_beta[e, e2] * capacity(rd, e2) for e2 in edges(rd))
        @constraint(m, lhs / capacity(rd, e) <= mu)
    end

    return RoutingModel(rd, m, UnitaryFlows, rm.routing, mu=mu, dual_alpha=dual_alpha, dual_beta=dual_beta)
end


text:
"""
    FEMMBaseModule

Module for comments/base operations on interiors and boundaries of domains.
"""
module FEMMBaseModule

import LinearAlgebra: mul!, Transpose
my_At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
import SparseArrays: sparse
import LinearAlgebra: norm
using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, manifdim, nodesperelem, subset, map2parametric, inparametric, centroidparametric, bfun
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianmdim, Jacobianvolume
import FinEtools.CSysModule: CSys
import FinEtools.FieldModule: ndofs, nents, gatherdofnums!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.ElementalFieldModule: ElementalField, nelems
import FinEtools.AssemblyModule: SysvecAssemblerBase, startassembly!, assemble!, makematrix!, SysvecAssembler, makevector!
import FinEtools.ForceIntensityModule: ForceIntensity, updateforce!
import FinEtools.MatrixUtilityModule: locjac!
import FinEtools.BoxModule: initbox!, boundingbox, inflatebox!
import FinEtools.MeshModificationModule: nodepartitioning, compactnodes, renumberconn!
import FinEtools.MeshSelectionModule: selectelem, vselect, findunconnnodes, connectednodes
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler

"""
    FEMMAbstractBase

Abstract base type for all finite element models.
"""
abstract type FEMMAbstractBase; end

"""
    FEMMBase

Class for base finite element modeling machine.
"""
mutable struct FEMMBase{S<:FESet, F<:Function} <: FEMMAbstractBase
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
end

"""
    FEMMBase(integdomain::IntegDomain{S, F}) where {S<:FESet, F<:Function}

Construct with the default orientation matrix (identity).
"""
function FEMMBase(integdomain::IntegDomain{S, F}) where {S<:FESet, F<:Function}
    return FEMMBase(integdomain, CSys(manifdim(integdomain.fes)))
end

"""
    associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})

Associate geometry field with the FEMM.

There may be operations that could benefit from pre-computations
that involve a geometry field. If so, associating the geometry
field gives the FEMM a chance to save on repeated computations.

Geometry field is normally passed into any routine that evaluates some
forms (integrals) over the mesh.  Whenever the geometry passed into a
routine is not consistent with the one for which `associategeometry!()`
was called before, `associategeometry!()` needs to be called with
the new geometry field.
"""
function associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})
    return self # default is no-op
end

"""
    inspectintegpoints(self::FEMM, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F,  idat, quantity=:Cauchy; context...) where {FEMM<:FEMMAbstractBase, T<:Number, F<:Function}

Inspect integration points.
"""
function inspectintegpoints(self::FEMM, geom::NodalField{FFlt}, felist::FIntVec, inspector::F,  idat, quantity=:Cauchy; context...) where {FEMM<:FEMMAbstractBase, T<:Number, F<:Function}
    return idat # default is no-op
end

"""
    integratefieldfunction(self::FEMMAbstractBase,
        geom::NodalField{FFlt},  afield::FL, fh::F,  initial::R;
        m::FInt=-1) where {T<:Number, FL<:NodalField{T}, R, F<:Function}

Integrate a nodal-field function over the discrete manifold.

`afield` = NODAL field to be supply the values
`fh` = function taking position and the field value as arguments, returning value of type `R`.

Returns value of type `R`, which is initialized by `initial`.
"""
function integratefieldfunction(self::FEMMAbstractBase, geom::NodalField{FFlt},  afield::FL, fh::F,  initial::R;
    m::FInt=-1) where {T<:Number, FL<:NodalField{T}, R, F<:Function}
    fes = self.integdomain.fes  # finite elements
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(afield); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes per element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain);
    a = fill(zero(T), nne,ndn); # array of field DOFS-- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    val = fill(zero(T), 1, ndn); # field value at the point -- used as a buffer
    J = fill(zero(FFlt), sdim,mdim); # Jacobian matrix -- used as a buffer
    if m >= 0
        # Either the manifold dimension was supplied
    else
        m=mdim;# ...Or it is implied
    end
    result = initial;           # initial value for the result
    for i=1:count(fes) #Now loop over all fes in the block
        gathervalues_asmat!(afield, a, fes.conn[i]);# retrieve element dofs
        for j = 1:npts #Loop over all integration points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            my_At_mul_B!(val, Ns[j], a);# Field value at the quadrature point
            Jac = Jacobianmdim(self.integdomain, J, loc, fes.conn[i],  Ns[j], m);
            result = result + fh(loc,val)*Jac*w[j];
        end
    end
    return result
end

"""
    integratefieldfunction(self::FEMMAbstractBase,
        geom::NodalField{FFlt},  afield::FL, fh::F, initial::R;
        m::FInt=-1) where {T<:Number, FL<:ElementalField{T}, R, F<:Function}

Integrate a elemental-field function over the discrete manifold.

`afield` = ELEMENTAL field to be supply the values
`fh` = function taking position and the field value as arguments, returning value of type `R`.

Returns value of type `R`, which is initialized by `initial`.
"""
function integratefieldfunction(self::FEMMAbstractBase, geom::NodalField{FFlt},  afield::FL, fh::F,  initial::R;
    m::FInt=-1) where {T<:Number, FL<:ElementalField{T}, R, F<:Function}
    fes = self.integdomain.fes  # finite elements
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(afield); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes per element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain);
    a = fill(zero(FFlt), nne,ndn); # array of field DOFS-- used as a buffer
    loc = fill(zero(FFlt), 1,sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim,mdim); # Jacobian matrix -- used as a buffer
    if m >= 0
        # Either the manifold dimension was supplied
    else
        m = mdim;# ...Or it is implied
    end
    result = initial;           # initial value for the result
    for i=1:count(fes) #Now loop over all fes in the block
        gathervalues_asmat!(afield, a, [i]);# retrieve element dofs
        for j = 1:npts #Loop over all integration points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianmdim(self.integdomain, J, loc, fes.conn[i],  Ns[j], m);
            result = result + fh(loc, a)*Jac*w[j];
        end
    end
    return result
end

"""
    integratefunction(self::FEMMAbstractBase,
        geom::NodalField{FFlt}, fh::F, m::FInt = -1) where {F<:Function}

Integrate a function over the discrete manifold.

Integrate some scalar function over the geometric cells. The function takes a
single argument, the position vector.

When the scalar function returns just +1 (such as `(x) ->  1.0`), the result
measures the volume (number of points, length, area, 3-D volume, according to
the manifold dimension). When the function returns the mass density, the
method measures the mass, when the function returns the x-coordinate equal
measure the static moment with respect to the y- axis, and so on.

# Example:
Compute the volume of the mesh and then its center of gravity:
```
V = integratefunction(femm, geom, (x) ->  1.0)
Sx = integratefunction(femm, geom, (x) ->  x[1])
Sy = integratefunction(femm, geom, (x) ->  x[2])
Sz = integratefunction(femm, geom, (x) ->  x[3])
CG = vec([Sx Sy Sz]/V)
```
Compute the moments of inertia:
```
% Now compute the moments of inertia
Ixx = measure(femm,geom,@(x)(norm(x-CG)^2*eye(3)-(x-CG)'*(x-CG)))
mass = V*rhos;
Inertia = I*rhos;
```
"""
function integratefunction(self::FEMMAbstractBase, geom::NodalField{FFlt}, fh::F, m::FInt = -1) where {F<:Function}
    fes = self.integdomain.fes
    if m < 0
        m = manifdim(fes);  # native  manifold dimension
    end
    # Constants
    nfes = count(fes); # number of finite elements in the set
    nne = nodesperelem(fes); # number of nodes per element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain);
    loc = fill(zero(FFlt), 1,sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim,mdim); # Jacobian matrix -- used as a buffer
    result = 0.0;# Initialize the result
    for i = 1:count(fes)  # Now loop over all fes in the set
        for j=1:npts #Loop over all integration points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianmdim(self.integdomain, J, loc, fes.conn[i],  Ns[j], m);
            result = result + fh(vec(loc))*Jac*w[j];
        end
    end
    return result
end

"""
    transferfield!(ff::F, fensf::FENodeSet, fesf::FESet,
        fc::F, fensc::FENodeSet, fesc::FESet, tolerance::FFlt
        )  where {T<:Number, F<:NodalField{T}}

Transfer a nodal field from a coarse mesh to a finer one.
- `ff` = the fine-mesh field (modified and also returned)
- `fensf` = finite element node set for the fine-mesh
- `fc` = the coarse-mesh field
- `fensc` = finite element node set for the fine-mesh,
- `fesc` = finite element set for the coarse mesh
- `geometricaltolerance` = tolerance in physical space for searches of the
  adjacent nodes
- `parametrictolerance` = tolerance in parametric space for for check whether
  node is inside an element
"""
function transferfield!(ff::F, fensf::FENodeSet, fesf::FESet, fc::F, fensc::FENodeSet, fesc::FESet, geometricaltolerance::FFlt; parametrictolerance::FFlt = 0.01)  where {T<:Number, F<:NodalField{T}}
    fill!(ff.values, Inf) # the "infinity" value indicates a missed node
    @assert count(fensf) == nents(ff)
    parametrictol = 0.01
    nodebox = initbox!([], vec(fensc.xyz[1, :]))
    # Find out how many partitions of the nodes on the fine mesh we should use
    npartitions = max(2, Int(round(count(fensf)/1000)))
    # Partition the nodes of the fine mesh
    npf = nodepartitioning(fensf, npartitions)
    partitionnumbers = unique(npf)
    npartitions = length(partitionnumbers)
    # Go through all the partitions
    for p = 1:npartitions
        pnl = findall(x -> x == partitionnumbers[p], npf) # subset of fine-mesh nodes
        # Find the bounding box
        subbox = boundingbox(fensf.xyz[pnl, :])
        tol = 2*geometricaltolerance # increase the box a bit
        # Construct a sub mesh of the coarse mesh that covers the nodes from this partition
        sublist = selectelem(fensc, fesc, overlappingbox = subbox, inflate = tol)
        if !isempty(sublist) # there are some finite elements to work with
            fescsub = subset(fesc, sublist)
            connected = findunconnnodes(fensc, fescsub);
            fenscsub, newnumber = compactnodes(fensc, connected); # nodes of the sub mesh
            fescsub = renumberconn!(fescsub, newnumber); # elements of the sub mesh
            present = findall(x -> x > 0, newnumber)
            fcsub  =  NodalField(fc.values[present, :]) # reduce the coarse-mesh field to the sub mesh
            # Now we can find the values at the nodes of the subset of the fine mesh
            # working only with the sub mesh of the coarse mesh
            for i in pnl # for all nodes in the subset
                nl = vselect(fenscsub.xyz; nearestto = fensf.xyz[i, :])
                # For each node in the fine field try to find one in  the coarse field
                if !isempty(nl) && norm(fensf.xyz[i, :] - fenscsub.xyz[nl[1], :]) < geometricaltolerance
                    # These nodes correspond in the  refined  and coarse mesh
                    ff.values[i, :] = fcsub.values[nl[1], :]
                else
                    # Obviously, some nodes in the fine mesh are not located "at" the
                    # location of a coarse-mesh node. Then we have to search which
                    # element they fall into.
                    nodebox = initbox!(nodebox, vec(fensf.xyz[i, :]))
                    nodebox = inflatebox!(nodebox, geometricaltolerance)
                    el = selectelem(fenscsub, fescsub; overlappingbox = nodebox)
                    for e = el
                        c = [k for k in fescsub.conn[e]] #c = view(fescsub.conn, e, :)
                        pc, success = map2parametric(fescsub, fenscsub.xyz[c, :],
                            vec(fensf.xyz[i, :]); tolerance = 0.000001, maxiter =7)
                        @assert success # this shouldn't be tripped; normally we succeed
                        if inparametric(fescsub, pc; tolerance = parametrictolerance) # coarse mesh element encloses the node
                            N = bfun(fescsub,  pc)
                            ff.values[i, :] = transpose(N) * fcsub.values[c, :]
                            break
                        end
                    end
                end
            end # for i in pnl # for all nodes in the subset
        end # if !isempty(sublist)
    end # for p = 1:npartitions

    # Check that we haven't missed any node connected to some finite elements.
    cnl = connectednodes(fesf)
    for i = cnl
        if any(v -> v == Inf, ff.values[i, :])
            println("fensf.xyz[$i, :] = $(fensf.xyz[i, :])")
            error("Missed node in transfer")
        end
    end

    return ff
end

# Old, non-performing version. It can get very expensive for large meshes.
# function transferfield!(ff::F, fensf::FENodeSet, fesf::FESet,
#     fc::F, fensc::FENodeSet, fesc::FESet, tolerance::FFlt
#     )  where {T<:Number, F<:NodalField{T}}
#     @assert count(fensf) == nents(ff)
#     nodebox = initbox!([], vec(fensc.xyz[1, :]))
#     for i = 1:count(fensf)
#         nl = vselect(fensc.xyz; nearestto = fensf.xyz[i, :])
#         # For each node in the fine field try to find one in  the coarse field
#         if !isempty(nl) && norm(fensf.xyz[i, :] - fensc.xyz[nl[1], :]) < tolerance
#             # These nodes correspond in the  refined  and coarse mesh
#             ff.values[i, :] = fc.values[nl[1], :]
#         else
#             # Obviously, some nodes in the fine mesh are not located "at" the
#             # location of a coarse-mesh node. Then we have to search which
#             # element they fall into.
#             nodebox = initbox!(nodebox, vec(fensf.xyz[i, :]))
#             nodebox = inflatebox!(nodebox, tolerance)
#             el = selectelem(fensc, fesc; overlappingbox = nodebox)
#             foundone = false
#             for e = el
#                 c = view(fesc.conn, e, :)
#                 pc, success = map2parametric(fesc, fensc.xyz[c, :],
#                 vec(fensf.xyz[i, :]); Tolerance = 0.000001, maxiter =7)
#                 @assert success # this shouldn't be tripped; normally we succeed
#                 if  inparametric(fesc, pc; tolerance = 0.01) # coarse mesh element encloses the node
#                     N = bfun(fesc,  pc)
#                     ff.values[i, :] = transpose(N) * fc.values[c, :]
#                     foundone = true
#                     break
#                 end
#             end
#             @assert foundone
#         end
#     end
#     return ff
# end

"""
    transferfield!(ff::F, fensf::FENodeSet, fesf::FESet, fc::F,
        fensc::FENodeSet, fesc::FESet, geometricaltolerance::FFlt;
        parametrictolerance::FFlt = 0.01 )  where {T<:Number,
        F<:ElementalField{T}}

Transfer an elemental field from a coarse mesh to a finer one.
- `ff` = the fine-mesh field (modified and also returned)
- `fensf` = finite element node set for the fine-mesh
- `fc` = the coarse-mesh field
- `fensc` = finite element node set for the fine-mesh,
- `fesc` = finite element set for the coarse mesh
- `tolerance` = tolerance in physical space for searches of the adjacent nodes
"""
function transferfield!(ff::F, fensf::FENodeSet, fesf::FESet, fc::F, fensc::FENodeSet, fesc::FESet, geometricaltolerance::FFlt; parametrictolerance::FFlt = 0.01)  where {T<:Number, F<:ElementalField{T}}
    @assert count(fesf) == nents(ff)
    nodebox = initbox!([], vec(fensc.xyz[1, :]))
    centroidpc = centroidparametric(fesf)
    N = bfun(fesf, centroidpc)
    NT = transpose(N)
    for i = 1:count(fesf) # For all finite elements in the fine mesh
        c = [k for k in fesf.conn[i]]
        centroid = NT * fensf.xyz[c, :]
        nodebox = initbox!(nodebox, vec(centroid))
        nodebox = inflatebox!(nodebox, geometricaltolerance)
        el = selectelem(fensc, fesc; overlappingbox = nodebox)
        foundone = false
        for e = el
            c = [k for k in fesc.conn[e]]
            pc, success = map2parametric(fesc, fensc.xyz[c, :],
                vec(centroid); tolerance = 0.000001, maxiter =9)
            # if !success
            # println("pc = $(pc)")
            # N1 = bfun(fesf, pc)
            # p = transpose(N1) * fensc.xyz[view(fesc.conn, e, :), :]
            # println("p = $(p)")
            # println("centroid = $(centroid)")
            # end
            # @assert success # this shouldn't be tripped; normally we succeed
            if success &&  inparametric(fesc, pc; tolerance = 0.001) # coarse mesh element encloses the centroid
                ff.values[i, :] = fc.values[e, :]
                foundone = true
                break
            end
        end
        @assert foundone
    end
    return ff
end

"""
    distribloads(self::FEMM, assembler::A, geom::NodalField{FFlt}, P::NodalField{T},
      fi::ForceIntensity,
      m::FInt) where {FEMM<:FEMMAbstractBase, T<:Number, A<:SysvecAssemblerBase}

Compute the distributed-load vector.

`fi`=force intensity object
`m`= manifold dimension, 1= curve, 2= surface, 3= volume
"""
function distribloads(self::FEMM, assembler::A,
    geom::NodalField{FFlt}, P::NodalField{T}, fi::ForceIntensity, m::FInt) where {FEMM<:FEMMAbstractBase, T<:Number, A<:SysvecAssemblerBase}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(P); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes per element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes);     # manifold dimension of the element
    Cedim = ndn*nne;       # dimension of the element matrix/vector
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain);
    # Prepare some buffers:
    dofnums = fill(zero(FInt), Cedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), (1,sdim)); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), (sdim, mdim)) # Jac. matrix -- used as a buffer
    Fe = fill(zero(T), (Cedim, ))
    startassembly!(assembler, P.nfreedofs);
    for i = 1:nfes # Loop over elements
        fill!(Fe, 0.0);
        for j = 1:npts
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianmdim(self.integdomain, J, loc, fes.conn[i],  Ns[j], m);
            force = updateforce!(fi, loc, J, fes.label[i]); # retrieve the applied load
            Factor::FFlt = (Jac * w[j]);
            NkxF::FFlt = 0.0
            rx::FInt=1;
            for kx = 1:nne # all the nodes
                NkxF = Ns[j][kx] * Factor
                for mx = 1:ndn   # all the degrees of freedom
                    Fe[rx] = Fe[rx] + NkxF * force[mx];
                    rx = rx+1;    # next component of the vector
                end
            end
        end
        gatherdofnums!(P, dofnums, fes.conn[i]);
        assemble!(assembler, Fe, dofnums);
    end
    F = makevector!(assembler);
    return F
end

function distribloads(self::FEMM, geom::NodalField{FFlt}, P::NodalField{T}, fi::ForceIntensity, m::FInt) where {FEMM<:FEMMAbstractBase, T<:Number}
    assembler = SysvecAssembler(0.0*P.values[1])#T(0.0))
    return distribloads(self, assembler, geom, P, fi, m)
end

"""
    connectionmatrix(self::FEMM, nnodes::FInt) where {FEMM<:FEMMAbstractBase}

Compute the connection matrix.

The matrix has a nonzero in all the rows and columns which correspond to nodes
connected by some finite element.
"""
function connectionmatrix(self::FEMM, nnodes::FInt) where {FEMM<:FEMMAbstractBase}
    fes = self.integdomain.fes
    nfes = length(fes.conn)
    nconns = nodesperelem(fes)
    N = nfes*nconns*nconns
    rb = FInt[]; sizehint!(rb, N)
    cb = FInt[]; sizehint!(cb, N)
    vb = ones(FInt, N);
    @inbounds for  j = 1:nfes
        @inbounds for  k = 1:nconns
            append!(rb, fes.conn[j])
            @inbounds for  m = 1:nconns
                push!(cb, fes.conn[j][k])
            end
        end
    end
    return sparse(rb, cb, vb, nnodes, nnodes)
end

struct InverseDistanceInspectorData
    component::Vector{FInt}
    d::Vector{FFlt} # nodesperelem(integdomain.fes)
    sum_inv_dist::Vector{FFlt} # nnodes(geom)
    sum_quant_inv_dist::Array{FFlt, 2} # nnodes(geom) x length(component)
end

# This is an inverse-distance interpolation inspector.
function _idi_inspector(idat, elnum, conn, xe,  out,  xq)
    # xe = coordinates of the nodes of the element
    # xq = coordinate of the quadrature point
    mindn = Inf
    for jjj = 1:size(xe,1)
        idat.d[jjj]  = sum((vec(xe[jjj, :])-vec(xq)).^2);
        if (idat.d[jjj] > 0.0)
            mindn = min(mindn, idat.d[jjj])
        end
    end
    mindn = mindn/1.0e9;
    for jjj = 1:length(idat.d)
        invdjjj  = 1.0/(idat.d[jjj]+mindn);
        quant =  out[idat.component]
        for kkk = 1:length(quant)
            idat.sum_quant_inv_dist[conn[jjj], kkk] += invdjjj*quant[kkk];
        end
        idat.sum_inv_dist[conn[jjj]] += invdjjj;
    end
    return idat
end

struct AveragingInspectorData
    component::Vector{FInt}
    d::Vector{FFlt} # nodesperelem(fes)
    ncontrib::Vector{FInt} # nnodes(geom)
    sum_quant::Array{FFlt, 2} # nnodes(geom) x length(component)
end

# This is a simple nodal averaging inspector.
# The quadrature point is assumed contribute only to the node that is nearest to
# it: The idea is that the inspected element  will produce outputs  at the
# locations of the nodes.
function _avg_inspector(idat, elnum, conn, xe,  out,  xq)
    # xe = coordinates of the nodes of the element
    # xq = coordinate of the quadrature point
    for jjj = 1:size(xe, 1)
        idat.d[jjj] = sum((vec(xe[jjj, :])-vec(xq)).^2);
    end
    # Find  the node nearest to the quadrature point
    (minval, ix) = findmin(idat.d)
    quant =  out[idat.component]
    for kkk = 1:length(quant)
        idat.sum_quant[conn[ix], kkk] += quant[kkk];
    end
    idat.ncontrib[conn[ix]] += 1
    return idat
end

"""
    fieldfromintegpoints(self::FEMM,
      geom::NodalField{FFlt},  u::NodalField{T},
      dT::NodalField{FFlt},  quantity::Symbol,  component::FInt;
      context...) where {FEMM<:FEMMAbstractBase, T<:Number}

Construct nodal field from integration points.

Input arguments
`geom`     - reference geometry field
`u`        - displacement field
`dT`       - temperature difference field
`quantity`   - this is what you would assign to the 'quantity' argument
           of the material update!() method.
`component`- component of the 'quantity' array: see the material update()
           method.
Keyword arguments
`nodevalmethod` = `:invdistance` (the default) or `:averaging`;
`reportat` = at which point should the  element quantities be reported?
    This argument is interpreted inside the `inspectintegpoints()` method.
Output argument
 - the new field that can be used to map values to colors and so on
"""
function fieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    dT::NodalField{FFlt},  quantity::Symbol,  component::FIntVec;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    fes = self.integdomain.fes
    # Constants
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    nodevalmethod = :invdistance
    reportat = :default
    for apair in pairs(context)
        sy, val = apair
        if sy == :nodevalmethod
            nodevalmethod = val
        end
        if sy == :reportat
            reportat = val
        end
    end
    if nodevalmethod == :averaging
        # Container of intermediate results
        idat = AveragingInspectorData(
        component,
        zeros(FFlt, nne),
        zeros(Int, nnodes(geom)),
        zeros(FFlt, nnodes(geom), length(component))
        );
        # Loop over cells to interpolate to nodes
        idat = inspectintegpoints(self,  geom,  u,  dT,
            collect(1:count(fes)),  _avg_inspector,  idat, quantity;
            context...);
        # The data for the field to be constructed is initialized
        nvals = zeros(FFlt, nnodes(geom), length(component));
        # compute the data array
        for kkk = 1:size(nvals, 1)
            for j = 1:length(component)
                if (idat.ncontrib[kkk] > 0)
                    nvals[kkk, j] = idat.sum_quant[kkk, j]/idat.ncontrib[kkk];
                end
            end
        end
    else # inverse distance
        @assert (reportat == :default) ||
            (reportat == :meanonly) "Inverse-distance interpolation requires :meanonly"
        # Container of intermediate results
        idat = InverseDistanceInspectorData(
        component,
        zeros(FFlt, nne),
        zeros(FFlt, nnodes(geom)),
        zeros(FFlt, nnodes(geom), length(component))
        );
        # Loop over cells to interpolate to nodes
        idat = inspectintegpoints(self,  geom,  u,  dT,
            collect(1:count(fes)),  _idi_inspector,  idat, quantity;
            context...);
        # The data for the field to be constructed is initialized
        nvals = zeros(FFlt, nnodes(geom), length(component));
        # compute the data array
        for kkk = 1:size(nvals, 1)
            for j = 1:length(component)
                if (idat.sum_inv_dist[kkk] > 0.0)
                    nvals[kkk, j] = idat.sum_quant_inv_dist[kkk, j]/idat.sum_inv_dist[kkk];
                end
            end
        end
    end
    # Make the field
    return NodalField(nvals);
end

function fieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    dT::NodalField{FFlt},  quantity::Symbol,  component::FInt;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    return fieldfromintegpoints(self, geom, u, dT, quantity, [component];
        context...)
end

function fieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    quantity::Symbol,  component::FIntVec;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    dT = NodalField(zeros(FFlt, nnodes(geom), 1)) # zero difference in temperature
    return fieldfromintegpoints(self, geom, u, dT, quantity, component; context...)
end

function fieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    quantity::Symbol,  component::FInt;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    dT = NodalField(zeros(FFlt, nnodes(geom), 1)) # zero difference in temperature
    return fieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
end

struct MeanValueInspectorData
    n_quant::Vector{FInt}
    sum_quant_value::Array{FFlt, 2}
end

"""
    elemfieldfromintegpoints(self::FEMM,
      geom::NodalField{FFlt},  u::NodalField{T},
      dT::NodalField{FFlt},  quantity::Symbol,  component::FInt;
      context...) where {FEMM<:FEMMAbstractBase, T<:Number}

Construct nodal field from integration points.

Input arguments
`geom`     - reference geometry field
`u`        - displacement field
`dT`       - temperature difference field
`quantity`   - this is what you would assign to the 'quantity' argument
           of the material update!() method.
`component`- component of the 'quantity' array: see the material update()
           method.
Output argument
 - the new field that can be used to map values to colors and so on
"""
function elemfieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    dT::NodalField{FFlt},  quantity::Symbol,  component::FIntVec;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    fes = self.integdomain.fes
    # Constants
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    # Container of intermediate results
    idat = MeanValueInspectorData(
    zeros(FInt, count(fes)),
    zeros(FFlt, count(fes), length(component))
    );
    # This is an mean-value interpolation inspector. The mean of the
    # quadrature-point quantities is reported per element.
    function mv_inspector(idat, elnum, conn, xe,  out,  xq)
        # xe = coordinates of the nodes of the element
        # xq = coordinate of the quadrature point
        idat.n_quant[elnum] += 1
        quant = out[component]
        for kkk = 1:length(quant)
            idat.sum_quant_value[elnum, kkk] += quant[kkk];
        end
        return idat
    end
    # Loop over cells to interpolate to nodes
    idat = inspectintegpoints(self,  geom,  u,  dT,
        collect(1:count(fes)), mv_inspector,  idat, quantity;
        context...);
    # The data for the field to be constructed is initialized
    evals = zeros(FFlt, count(fes), length(component));
    # compute the data array
    for j = 1:size(evals, 1)
        for kkk = 1:size(evals, 2)
            evals[j, kkk] = idat.sum_quant_value[j, kkk]/idat.n_quant[j];
        end
    end
    # Make the field
    return ElementalField(evals);
end

function elemfieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    dT::NodalField{FFlt},  quantity::Symbol,  component::FInt;
    context...) where {FEMM<:FEMMAbstractBase, T<:Number}
    return elemfieldfromintegpoints(self,
        geom, u, dT,  quantity,  [component]; context...)
end

function elemfieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    quantity::Symbol,  component::FInt; context...) where {FEMM<:FEMMAbstractBase,
    T<:Number}
    dT = NodalField(zeros(FFlt, nnodes(geom), 1)) # zero difference in temperature
    return elemfieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
end

function elemfieldfromintegpoints(self::FEMM,
    geom::NodalField{FFlt},  u::NodalField{T},
    quantity::Symbol,  component::FIntVec; context...) where {FEMM<:FEMMAbstractBase,
    T<:Number}
    dT = NodalField(zeros(FFlt, nnodes(geom), 1)) # zero difference in temperature
    return elemfieldfromintegpoints(self, geom, u, dT, quantity, component; context...)
end

function  buffers(self::FEMM, geom::NodalField{FFlt}, afield::NodalField{T}) where {FEMM<:FEMMAbstractBase, T}
    # Constants
    fes = self.integdomain.fes
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(afield); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Kedim = ndn*nne;      # dimension of the element matrix
    elmat = fill(zero(FFlt), Kedim, Kedim); # buffer
    dofnums = fill(zero(FInt), Kedim); # buffer
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    gradN = fill(zero(FFlt), nne, mdim); # buffer
    return dofnums, loc, J, gradN, elmat
end

"""
    innerproduct(self::FEMMHeatDiff,
      assembler::A, geom::NodalField{FFlt},
      temp::NodalField{FFlt}) where {A<:SysmatAssemblerBase}

Compute the inner-product (Gram) matrix.
"""
function innerproduct(self::FEMM, assembler::A, geom::NodalField{FFlt}, afield::NodalField{T}) where {FEMM<:FEMMAbstractBase, A<:SysmatAssemblerBase, T}
    fes = self.integdomain.fes
    dofnums, loc, J, gradN, elmat = buffers(self, geom, afield)
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    NexpTNexp = FFltMat[];# basis f. matrix -- buffer
    ndn = ndofs(afield)
    Indn = [i==j ? one(FFlt) : zero(FFlt) for i=1:ndn, j=1:ndn] # "identity"
    for j = 1:npts # This quantity is the same for all quadrature points
        Nexp = fill(zero(FFlt), ndn, size(elmat,1))
        for l1 = 1:nodesperelem(fes)
            Nexp[1:ndn, (l1-1)*ndn+1:(l1)*ndn] = Indn * Ns[j][l1];
        end
        push!(NexpTNexp, Nexp'*Nexp);
    end
    startassembly!(assembler,  size(elmat,1),  size(elmat,2),  count(fes), afield.nfreedofs,  afield.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            thefactor::FFlt =(Jac*w[j]);
            elmat .+= NexpTNexp[j]*thefactor
        end # Loop over quadrature points
        gatherdofnums!(afield,  dofnums,  fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler,  elmat,  dofnums,  dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function innerproduct(self::FEMM, geom::NodalField{FFlt}, afield::NodalField{T}) where {FEMM<:FEMMAbstractBase, T}
    assembler = SysmatAssemblerSparseSymm();
    return innerproduct(self, assembler, geom, afield);
end

end


text:
# global packages -
using DifferentialEquations
using DiffEqFinancial
using Statistics
using Distributions
using KernelDensity
using Random
using Dates
using DataFrames
using Reexport
@reexport using PooksoftBase

# include my code ...
include("./base/Types.jl")
include("./base/Compute.jl")
include("./base/Utility.jl")
include("./base/Return.jl")
include("./base/Factory.jl")
include("./base/Technical.jl")
include("./base/Lattice.jl")


text:
## AbstractTensor ##
abstract type AbstractTensor end
## Methods to be implemented (see MemTensor{T} for an example
##      size
#       data
#       permutedims
#       getindex
#       zeroed initialization without intermediate array (defeats the purpose of e.g. distributed memory)

import Base.size
import Base.permutedims
import Base.getindex

## TODO: Implement zeroed initialization for all AbstractTensor

"""
    PackedTensor{T} <: AbstractTensor where T <: Number

Symmetry packed tesor. Not functional."""
struct PackedTensor{T} <: AbstractTensor where T <: Number
    symmetries::Dict{Tuple,String}
    data::Array{T,1}
    ndim::Int8
    dims::Array{UInt16,1}
end

## start GeneratedTensor ##
"""
    GeneratedTensor{T} <: AbstractTensor where T <: Number

Tensor for which elements are computed on the fly by `generator`.
"""
struct GeneratedTensor{T} <: AbstractTensor where T <: Number
    generator::Function
    data::Dict{String,Any}
    ndim::Int8
    dims::Array{UInt16,1}
end

function getindex(g::GeneratedTensor,I::Vararg{Int,N}) where N
    try
        g.generator(g,I...)
    catch
        error("Incorrect number of arguments passed to getindex(::GeneratedTensor) 😦")
    end
end
## end GeneratedTensor ##

## start DiskTensor ##
struct DiskTensor{T} <: AbstractTensor where T <: Number
    fname::String
    buf::Array{T}
    ndim::Int8
    dims::Array{UInt32,1}
end
## end DiskTensor ##

struct DistributedTensor{T} <: AbstractTensor where T <: Number
    data::DistributedArrays.DArray{T}
    ndim::Int8
    dims::Array{UInt32,1}
end

## start MemTensor ##
"""
    MemTensor{T} <: AbstractTensor where T <: Number

Tensor object held entirely in memory. Thin wrapper around a standard Julia array,
but tagged as a MemTensor to fit into Fermi's dispatch system.
"""
struct MemTensor{T} <: AbstractTensor where T <: Number
    data::Array{T}
    ndim::Int8
    dims::Array{UInt16,1}
end

function MemTensor(data::Array)
    @assert eltype(data) <: Number "MemTensor was passed an array of non-numeric values."
    MemTensor{eltype(data)}(data)
end
function MemTensor{T}(data::Array{T}) where T <: Number
    ndim = convert(Int8,length(size(data)))
    dims = zeros(UInt16,ndim)
    for (idx,val) in enumerate(size(data))
        dims[idx] = convert(UInt16,val)
    end
    MemTensor{T}(data,ndim,dims)
end

function size(A::MemTensor)
    size(A.data)
end
function data(A::MemTensor)
    A.data
end

function data(A::Array)
    A
end

function permutedims(A::MemTensor,tup)
    A.data .= permutedims(A.data,tup)
end

function getindex(A::MemTensor,I...)
    return A.data[I...]
end
## end MemTensor ##


text:
using ImageView, Images
#

cloud = load("/home/travis/Documents/Mael Engine/sprites/cloud.png")
cloud2 = load("/home/travis/Documents/Mael Engine/sprites/cloud2.png")
LeftPlayer = load("/home/travis/Documents/Mael Engine/sprites/leftplayer.png")
MiddlePlayer = load("/home/travis/Documents/Mael Engine/sprites/middleplayer.png")
img = fill(RGBA(0.65625,0.796875,1,1), 500, 800)

# ==============================================================================
# This draws a line to a picture!
# ==============================================================================
function line(image, x1, y1, x2, y2, color)
        dx, dy = x2-x1,  y2-y1

        for x in x1:x2
              y = y1 + dy * (x - x1) / dx
              image[ round(Int32, y), round(Int32, x) ] = color
        end
end
# ==============================================================================
# This draws a rectrangle to a picture!
# ==============================================================================
function rect(image, left, top, width, height, color)
    for y in 1:height, x in 1:width
        pixel(image, top+y, left+x, color )
    end
end
# ==============================================================================
# This draws a rectrangle to a picture! 1="centered" 2= topCenter
# ==============================================================================
function drawSprite(image, sprite, x, y, position=1)
    # width, height = size(image)
    SpHeight, SpWidth = size(sprite)
    if position==1
        cx,cy = round(Int32, SpWidth/2), round(Int32, SpHeight/2)
    elseif  position==2
        cx,cy = round(Int32, SpWidth/2), SpHeight
    elseif position==3
        cx,cy = round(Int32, SpWidth/2), 0
    end
    left, top = x-cy, y-cx

    for y in 1:SpWidth,  x in 1:SpHeight
        pixel(image, left+x, top+y, sprite[ x, y ] )
    end
end
# ==============================================================================
# This draws a pixel with Alpha chanel
# ==============================================================================
function pixel(image, X, Y, color )
    a = alpha(color)
    width, height = size(image)
    if X>0 && X < width && Y>0 && Y < height
        image[ X,Y ] = image[ X,Y ]*(1-a) + (color*a)
    end
end

rect(img, 200, 400, 150, 25, RGBA(0.65,0,0.5, 0.5) )
drawSprite(img, LeftPlayer, 400, 274, 2)
drawSprite(img, MiddlePlayer, 401, 477, 2)
rect(img, 0, 0, 50, 50, RGBA(0.65,0,0.5, 0.5) )

rect(img, 300, 410, 150, 25, RGBA(0.72,0.65,0.5, 0.8) )
# line(img, 100, 200, 250, 275, RGB(0,0.3,0.5))
drawSprite(img, cloud, 100, 100)backup
drawSprite(img, cloud, 100, 600)

drawSprite(img, cloud2, 100, 400)
drawSprite(img, cloud, 150, 200)
imshow(img)

#  type Player
#  	sprite::Sprite
#  	velocity::Vector2f
#  end


text:
module Aligner

export align,
       levenshtein,
       generic_edit_distance,
       Alignment,
       null

type Alignment
    source_to_target
    target_to_source
    source
    target
    mat
end

null(::Type{Char}) = '-'
null{S<:AbstractString}(::Type{S}) = S("N")
null{S<:AbstractString}(::Type{SubString{S}}) = S("N")

cost_func(cost) = (x, y) -> cost
equality(cost) = (x, y) -> x == y ? 0 : cost

function generic_distance_matrix{A,B}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::Function,
    del_cost::Function,
    sub_cost::Function
)
    seq1len = length(seq1)
    seq2len = length(seq2)
    len = max(seq1len, seq2len)

    d_type = typeof(promote(map(f->f(seq1[1], seq2[1]),
        [ins_cost, del_cost, sub_cost])...)[1])
    d = zeros(d_type, seq1len+1, seq2len+1)
    d[:,1] = 0:seq1len
    d[1,:] = 0:seq2len

    for i = 2:seq1len+1, j = 2:seq2len+1
        ins = j < 3 ? 1 : ins_cost(seq1[i-1], seq2[j-2])
        del = i < 3 ? 1 : del_cost(seq1[i-2], seq2[j-1])
        # ins = ins_cost(seq1[i-1], seq2[j-1])
        # del = del_cost(seq1[i-1], seq2[j-1])
        sub = sub_cost(seq1[i-1], seq2[j-1])
        d[i,j] = min(d[i-1,j]+ins, d[i,j-1]+del, d[i-1,j-1]+sub)
    end
    d
end

function generic_distance_matrix{A,B,C<:Real}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::C,
    del_cost::C,
    sub_cost::C
)
    generic_distance_matrix(seq1, seq2, cost_func(ins_cost), cost_func(del_cost), equality(sub_cost))
end

function generic_distance_matrix{A,B,C<:Real}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::C,
    del_cost::C,
    sub_cost::Function
)
    generic_distance_matrix(seq1, seq2, cost_func(ins_cost), cost_func(del_cost), sub_cost)
end

function generic_distance_matrix{S<:AbstractString,T<:AbstractString}(
    seq1::S, seq2::T, ins_cost, del_cost, sub_cost)
    generic_distance_matrix(collect(seq1), collect(seq2),
        ins_cost, del_cost, sub_cost)
end

# levenshtein edit distance in which costs are (ins, del, sub) == (1, 1, 2)
levenshtein(seq1, seq2) = generic_distance_matrix(seq1, seq2, 1, 1, 2)[end, end]

generic_edit_distance(seq1, seq2, ins, del, sub) = generic_distance_matrix(seq1, seq2, ins, del, sub)[end, end]

function align{A,B}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::Function,
    del_cost::Function,
    sub_cost::Function
)
    mat = generic_distance_matrix(seq1, seq2, ins_cost, del_cost, sub_cost)
    res1 = similar(seq1, 0)
    res2 = similar(seq2, 0)
    idx1 = Int[]
    idx2 = Int[]
    i = size(mat, 1); j = size(mat, 2)
    while true
        i <= 1 && j <= 1 && break
        if i > 1 && j > 1 && mat[i,j] == mat[i-1,j-1] + sub_cost(seq1[i-1], seq2[j-1])
            i -= 1; j -= 1
            push!(idx1, j)
            push!(idx2, i)
            push!(res1, seq1[i])
            push!(res2, seq2[j])
        elseif i > 1 && mat[i,j] == mat[i-1,j] + (j == size(mat, 2) ? 1 : ins_cost(seq1[i-1], seq2[j]) )
            i -= 1
            push!(idx1, 0)
            push!(idx2, i)
            push!(res1, seq1[i])
            push!(res2, null(B))
        elseif j > 1 && mat[i,j] == mat[i,j-1] + (i == size(mat, 1) ? 1 : del_cost(seq1[i], seq2[j-1]) )
            j -= 1
            push!(idx1, j)
            push!(idx2, 0)
            push!(res1, null(A))
            push!(res2, seq2[j])
        end
    end
    Alignment(reverse!(idx1), reverse!(idx2), reverse!(res1), reverse!(res2), mat)
end

function align{A,B,C<:Real}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::C,
    del_cost::C,
    sub_cost::C
)
    align(seq1, seq2, cost_func(ins_cost), cost_func(del_cost), equality(sub_cost))
end

function align{A,B,C<:Real}(
    seq1::Vector{A},
    seq2::Vector{B},
    ins_cost::C,
    del_cost::C,
    sub_cost::Function
)
    align(seq1, seq2, cost_func(ins_cost), cost_func(del_cost), sub_cost)
end

# for debugging
function printmat(m)
    mat = map(x-> @sprintf("%1.2f", x), m)
    for i = 1:size(m, 1)
        for j = 1:size(m, 2)
            print(mat[i,j], " ")
        end
        println()
    end
end

end # module


text:
"""
```ϵ_hole_layers(x, y, ps, interstice = 0.5, hole = 0.75)```
return the permittivity of a unit-cell which consists of air holes in silica

Arguments:

- ps : widths of the air holes (need to be unit-less)
- refractive_indexes : optional argument with refractive indexes of background, hole and substrate. For reference simulation: set refractive indexes to ones(3)*eps_substrate
- interstice (optional): number of length_scale in between holes acounting for the permittivity of the substrate
- hole (optional): number of length_scale in the holes

Returns:

- geometry : a complex array with the epsilon data of the unit-cell
"""
function ϵ_hole_layers(x, y, ps; refractive_indexes=zeros(3), interstice = 0.5, hole = 0.75)
    @assert length(x)> 2 # makes sure δ is defined

    nx, ny = length(x), length(y)
    δ = x[2] - x[1]
    Ly_pml = y[end] - y[1] + δ
    Lx = x[end] - x[1] + δ

    @assert all(ps .> δ) # makes sure pixel-averaging handles all cases
    @assert all(ps .<= Lx) # makes sure that the holes are not bigger than the period

    # material properties of the unit-cell
    if refractive_indexes == zeros(3)
        refractive_index_background = 1.0
        refractive_index_hole = 1.0
        refractive_index_substrate = 1.45
    else
        refractive_index_background, refractive_index_hole,
        refractive_index_substrate = refractive_indexes
    end

    eps_background, eps_hole, eps_substrate =
    refractive_index_background^2, refractive_index_hole^2,
    refractive_index_substrate^2

    geometry = ones(ComplexF64, ny, nx) * eps_background

    index_top_substrate = floor(Int64, Ly_pml * 0.35/ δ) # 80% of the domain is substrate

    # substrate
    geometry[index_top_substrate:end, :] .= eps_substrate

    # handles case for sub-pixel averaging
    if x[nx÷2] == 0
        w_offset=1/2
    else
        w_offset=0
    end

    # holes
    number_holes = length(ps)
    n_inter_hole = floor(Int64, interstice / refractive_index_substrate / δ)
    n_hole_height = floor(Int64, hole / δ)

    @assert index_top_substrate + number_holes * (n_inter_hole + n_hole_height) < ny
    # makes sure that Ly is big enough for the holes (possibly in PML)

    for it_holes = 1:number_holes
        half_width = ps[it_holes]/2δ - w_offset
        n_half_width = floor(Int64, half_width)
        weight_eps_hole = half_width - n_half_width

        # inside holes
        n_start = floor(Int64, (nx - 2*n_half_width)/2 - w_offset)
        geometry[index_top_substrate + it_holes * n_inter_hole + (it_holes-1) *
        n_hole_height + 1:index_top_substrate + it_holes *
        (n_inter_hole + n_hole_height),
        n_start + 1: n_start + floor(Int64, 2*(n_half_width + w_offset)) + 1] .=
        eps_hole

        # pixel averaging
        # left
        geometry[index_top_substrate + it_holes * n_inter_hole + (it_holes-1) *
        n_hole_height + 1:index_top_substrate + it_holes *
        (n_inter_hole + n_hole_height), n_start] .=
        weight_eps_hole * eps_hole + (1 - weight_eps_hole) * eps_substrate
        # right
        geometry[index_top_substrate + it_holes * n_inter_hole + (it_holes-1) *
        n_hole_height + 1:index_top_substrate + it_holes *
        (n_inter_hole + n_hole_height), end-n_start+1] .=
        weight_eps_hole * eps_hole + (1 - weight_eps_hole) * eps_substrate
    end

    return geometry
end

"""
```function ϵ_pillar_function(x, y, ps; refractive_indexes=zeros(3))```
return the permittivity of a unit-cell which consists of a pillar on top of a substrate

Arguments:

- ps : width and height of the pillar (adimensional)
- refractive_indexes : optional argument with refractive indexes of background, hole and substrate. For reference simulation: set refractive indexes to ones(3)*eps_substrate

Returns:

- geometry : a complex array with the epsilon data of the unit-cell
"""
function ϵ_pillar_function(x, y, ps; refractive_indexes=zeros(3))
    @assert length(x)> 2 # makes sure δ is defined
    @assert length(ps) == 2
    width = ps[1]
    height = ps[2]

    nx, ny = length(x), length(y)
    δ = x[2] - x[1]
    Ly_pml = y[end] - y[1] + δ
    Lx = x[end] - x[1] + δ
    p = width/Lx

    @assert all(ps .> δ) # makes sure pixel-averaging handles all cases
    @assert (width <= Lx) # makes sure that the pillar width are not bigger than the period
    @assert (height <= Ly_pml) # makes sure that height fits in domain (TODO: need to be careful of PML)

    # material properties of the unit-cell
    if refractive_indexes == zeros(3)
        refractive_index_background = 1.0
        refractive_index_pillar = 2.0
        refractive_index_substrate = 1.4
    else
        refractive_index_background, refractive_index_pillar,
        refractive_index_substrate = refractive_indexes
    end

    eps_background, eps_pillar, eps_substrate =
    refractive_index_background^2, refractive_index_pillar^2,
    refractive_index_substrate^2

    geometry = ones(ComplexF64, ny, nx) * eps_background

    index_top_substrate = floor(Int64, Ly_pml * 0.5/ δ) # 80% of the domain is substrate

    # substrate
    geometry[index_top_substrate:end, :] .= eps_substrate

    # handles case for sub-pixel averaging
    if x[nx÷2] == 0
        w_offset=1/2
    else
        w_offset=0
    end

    # pillar
    n_pillar_height = floor(Int64, height / δ)

    half_width = p/2δ - w_offset
    n_half_width = floor(Int64, half_width)
    weight_eps_hole = half_width - n_half_width

    # inside pillar
    n_start = floor(Int64, (nx - 2*n_half_width)/2 - w_offset)
    geometry[index_top_substrate-n_pillar_height:index_top_substrate-1,
    n_start + 1: n_start + floor(Int64, 2*(n_half_width + w_offset)) + 1] .=
    eps_pillar

    # pixel averaging
    # left
    geometry[index_top_substrate-n_pillar_height:index_top_substrate-1, n_start] .=
    weight_eps_hole * eps_pillar + (1 - weight_eps_hole) * eps_background
    # right
    geometry[index_top_substrate-n_pillar_height:index_top_substrate-1, end-n_start+1] .=
    weight_eps_hole * eps_pillar + (1 - weight_eps_hole) * eps_background

    return geometry
end


text:
export PEBlock, projection_analysis
"""
    PEBlock(UG, n_reg, n_b) -> ChainBlock

phase estimation circuit.

    * `UG`: the input unitary matrix.
    * `n_reg`: the number of bits to store phases,
    * `n_b`: the number of bits to store vector.
"""
function PEBlock(UG::GeneralMatrixBlock, n_reg::Int, n_b::Int)
    nbit = n_b + n_reg
    # Apply Hadamard Gate.
    hs = repeat(nbit, H, 1:n_reg)

    # Construct a control circuit.
    control_circuit = chain(nbit)
    for i = 1:n_reg
        push!(control_circuit, control(nbit, (i,), (n_reg+1:nbit...,)=>UG))
        if i != n_reg
            UG = matblock(mat(UG) * mat(UG))
        end
    end

    # Inverse QFT Block.
    iqft = subroutine(nbit, QFTBlock{n_reg}()',[1:n_reg...,])
    chain(hs, control_circuit, iqft)
end

"""
    projection_analysis(evec::Matrix, reg::ArrayReg) -> Tuple

Analyse using state projection.
It returns a tuple of (most probable configuration, the overlap matrix, the relative probability for this configuration)
"""
function projection_analysis(evec::Matrix, reg::ArrayReg)
    overlap = evec'*state(reg)
    amp_relative = Float64[]
    bs = Int[]

    for b in basis(overlap)
        mc = argmax(view(overlap, b+1, :) .|> abs)-1
        push!(amp_relative, abs2(overlap[b+1, mc+1])/sum(overlap[b+1, :] .|> abs2))
        push!(bs, mc)
    end
    bs, overlap, amp_relative
end


text:
using Distributions, Random, DataFrames, GLM

export generate_synthetic_data, RegressionProblem, ols
function generate_synthetic_data(population_size, start_yr, end_yr)

    Random.seed!(1234)

    years = range(start_yr, end_yr+1)

    flu_mu = population_size*rand(Uniform(0.25, 0.75))
    flu_std = population_size*rand(Uniform(0.1, 0.2))

    flu_distribution = Normal(flu_mu, flu_std)
    flu_patients = [abs(rand(flu_distribution)) for i in start_yr:end_yr+1]

    beta_0 = population_size*rand(Uniform(0.1, 0.25))
    beta_1 = rand(Uniform(0.5, 0.75))

    error_distribution = Normal(0.1, 2.0)
    errors = [abs(rand(error_distribution)) for i in start_yr:end_yr]

    vaccines_produced = [0.0 for i in start_yr:end_yr]

    for i in range(2, stop=length(years)-1, step=1)
        vaccines_produced[i] = ceil(beta_0 + beta_1*flu_patients[i-1] + errors[i-1])
    end

    df = DataFrame(year=Int64[], flu_patients=Float64[], vaccines_produced=Union{Float64, Missing}[])

    push!(df, [years[1], flu_patients[1], missing])

    for i in range(2, length(years)-1)
        push!(df, [years[i], flu_patients[i], vaccines_produced[i]])
    end

    return df
end

function ols(df, num_flu_patients_from_sim, year_to_predict)

    model =  lm(@formula(vaccines_produced ~ flu_patients), df[2:length(df.year), [:year, :flu_patients, :vaccines_produced]])
    targetDF = DataFrame(year=year_to_predict, flu_patients=num_flu_patients_from_sim, vaccines_produced=missing)
    predicted_num_vaccines = predict(model, targetDF)
    println("GLM Model:")
    println(model)
    println("Predicted number of vaccines based on simulated number of flu patients for year ", year_to_predict, " = ", ceil(predicted_num_vaccines[1]))
    return predicted_num_vaccines

end

struct RegressionProblem{F,T,U,Y} <: AbstractModel
    formula::F
    β::T
    X::U
    y::Y
end

function solve(m::RegressionProblem{F,T,U,Missing}) where {F,T, U}
    return predict(m.β, m.X)
end

function solve(m::RegressionProblem{F, Missing, U,V}) where {F,U,V}
    return lm(m.formula, m.X)
end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
using SolidStateDetectors
using Test

T = Float32
max_refinements = 0

@testset "Package SolidStateDetectors" begin

    @testset "Simulate example detector: Inverted Coax" begin
        sim = Simulation(SSD_examples[:InvertedCoax])
        simulate!(sim, max_refinements = max_refinements)
        evt = SSD.Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 40e-3 )]))
        simulate!(evt, sim)
        signalsum::T = 0
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
    @testset "Simulate example detector: Coax" begin
        sim = Simulation(SSD_examples[:Coax])
        simulate!(sim, max_refinements = max_refinements)
        evt = SSD.Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 20e-3 )]))
        simulate!(evt, sim)
        signalsum::T = 0
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
    @testset "Simulate example detector: BEGe" begin
        sim = Simulation(SSD_examples[:BEGe])
        simulate!(sim, max_refinements = max_refinements)
        evt = SSD.Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 20e-3 )]))
        simulate!(evt, sim)
        signalsum::T = 0
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
    @testset "Simulate example detector: CGD" begin
        sim = Simulation(SSD_examples[:CGD])
        simulate!(sim, max_refinements = max_refinements)
        evt = SSD.Event([CartesianPoint{T}(5e-3, 5e-3, 5e-3)])
        simulate!(evt, sim)
        signalsum::T = 0
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end

end


text:
#regex escape from https://github.com/JuliaLang/julia/pull/29643
"convert string to an exact match regex."
function _regex_escape(s::AbstractString)
    res = replace(s, r"([()[\]{}?*+\-|^\$\\.&~#\s=!<>|:])" => s"\\\1")
    Regex(replace(res, "\0" => "\\0"))
end

"isolate given glossary in a word"
isolate_gloss(word::String, gloss::String)::Vector{String} = isolate_gloss(word, _regex_escape(gloss))
function isolate_gloss(word::String, gloss::Regex)
    if occursin(gloss, word) && !occursin(Regex("^"*gloss.pattern*"\$"), word)
        splits = split(word, gloss)
        matched = collect(eachmatch(gloss, word))
        res = String[]
        foreach(zip(splits, matched)) do (s, m)
            s != "" && push!(res, intern(s))
            m.match != "" && push!(res, intern(m.match))
        end
        splits[end] != "" && push!(res, intern(splits[end]))
    else
        res = String[intern(word)]
    end
    res
end

isolate_gloss(word::String, glosses::Vector{Union{Regex, String}}) =
    isolate_gloss(word,
                  map(x->typeof(x)===String ? _regex_escape(x) : x,
                      glosses)
                  )

isolate_gloss(word::String, glosses::Vector{String}) =
    isolate_gloss(word, map(_regex_escape, glosses))

function isolate_gloss(word::String, glosses::Vector{Regex})
    word = String[word]
    for gloss ∈ glosses
        len = length(word)
        for i ∈ 1:len
            wp = popfirst!(word)
            foreach(isolate_gloss(wp, gloss)) do nwp
                push!(word, nwp)
            end
        end
    end
    word
end

function isolate_gloss(xs::Vector{String}, glosses::Vector{Regex})
  ys = Vector{String}()
  for x in xs
    append!(ys, isolate_gloss(x, glosses))
  end
  return ys
end

struct Glossary
  glossaries::Vector{Regex}
end

Glossary(gs::Glossary) = gs
Glossary(gs::Vector) = Glossary(map(g->typeof(g) <: Regex ? g : _regex_escape(g), gs))

(g::Glossary)(x) = isolate_gloss(x, g.glossaries)

isgloss(g::Glossary, x) = any(gloss->(match(gloss, x) !== nothing), g.glossaries)