text:
# Types:

################################################################
# musicxml yes-no type
struct YN
    yn::String
    bool::Bool

    function YN(bool::Bool)
        if bool
            yn = "yes"
        else
            yn = "no"
        end
        new(yn, bool)
    end

    function YN(yn::String)
        if yn == "yes"
            bool = true
        elseif yn == "no"
            bool = false
        end
        new(yn, bool)
    end

    function YN(yn::String, bool::Bool)
        if (yn == "yes" && bool == true) || (yn == "no" && bool == false)
            new(yn, bool)
        else
            error("invalid YN definition")
        end
    end
end
YN(::Nothing) = nothing
################################################################
# Common function checks
midi16(x) = 1<= x <= 16
midi128(x) = 1 <= x <= 128
midi16384(x) = 1 <= x <= 16384
percent(x) = 0 <= x <=100
rot180(x) = -180 <= x <= 180
positive(x) = x>0
################################################################
"""
    ScoreInstrument

# Arguments
```julia
- id::String = "P1-I1", att"id"
- name::String = "Piano", "instrument-name"
- abbreviation::UN{String} = nothing, "instrument-abbreviation"
- sound::UN{String} = nothing, "instrument-sound"
- # ensemble::UN{Int64} = nothing, empty"~", positive
- # solo::UN{Int64} = nothing, empty"~"
- # VST::VST, "virtual-instrument"
```

The score-instrument type represents a single instrument within a score-part.

If you don't give scoreinstrument in ScorePart, it will automatically generate one.

As with the score-part type, each score-instrument has a required id, a name, and an optional abbreviation. A score-instrument type is also required if the score specifies MIDI 1.0 channels, banks, or programs. An initial midi-instrument assignment can also be made here. MusicXML software should be able to automatically assign reasonable channels and instruments without these elements in simple cases, such as where part names match General MIDI instrument names.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-score-instrument.htm)

# Example
```julia
ScoreInstrument(name = "Violin", id = "P1-I1")
```
"""
@aml mutable struct ScoreInstrument "score-instrument"
    id::String = "P1-I1", att"id"
    name::String = "Piano", "instrument-name"
    abbreviation::UN{String} = nothing, "instrument-abbreviation"
    sound::UN{String} = nothing, "instrument-sound"
    # ensemble::UN{Int64} = nothing, empty"~", positive
    # solo::UN{Int64} = nothing, empty"~"
    # VST::VST, "virtual-instrument"
end
################################################################
"""
    MidiDevice

# Arguments
```julia
- id::String = "P1-I1", att"id"
- port::String = "1", att"port"
```
The midi-device type corresponds to the DeviceName meta event in Standard MIDI Files. Unlike the DeviceName meta event, there can be multiple midi-device elements per MusicXML part starting in MusicXML 3.0.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-midi-device.htm)

# Example
```julia
MidiDevice(port = "1", id = "P1-I1")
```
"""
@aml mutable struct MidiDevice "midi-device"
    id::String = "P1-I1", att"id"
    port::String = "1", att"port"
end
################################################################
"""
    MidiInstrument

# Arguments
```julia
- id::String = "P1-I1", att"id"
- channel::Int64 = 1, "midi-channel",  midi16
- name::UN{String} = nothing, "midi-name"
- bank::UN{Int64} = nothing, "midi-bank", midi16384
- program::Int64 = 1, "midi-program", midi128
- unpitched::UN{Int64} = nothing, "midi-unpitched", midi16
- volume::Float64 = 100.0, "volume", percent
- pan::Float64 = 0, "pan", rot180
- elevation::UN{Float64} = nothing, "elevation", rot180
```

MidiInstrument type holds information about the sound of a midi instrument.

See https://en.wikipedia.org/wiki/General_MIDI#Program_change_events for different instrument programs in General Midi.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-midi-instrument.htm)

Pan: -90 is hard left, 0 is center, -180 is behind

# Example
```julia
MidiInstrument(channel= 1, program =1, volume = 100, pan =0, id = "P1-I1")
```
"""
@aml mutable struct MidiInstrument "midi-instrument"
    id::String = "P1-I1", att"id"
    channel::Int64 = 1, "midi-channel",  midi16
    name::UN{String} = nothing, "midi-name"
    bank::UN{Int64} = nothing, "midi-bank", midi16384
    program::Int64 = 1, "midi-program", midi128
    unpitched::UN{Int64} = nothing, "midi-unpitched", midi16
    volume::Float64 = 100.0, "volume", percent
    pan::Float64 = 0, "pan", rot180
    elevation::UN{Float64} = nothing, "elevation", rot180
end
################################################################
"""
    ScorePart

# Arguments
```julia
- id::String = "P1", att"id"
- # identification
- name::String = "Piano", "part-name"
- nameDisplay::UN{String} = name, "part-name-display"
- abbreviation::UN{String} = nothing, "part-abbreviation"
- abbreviationDisplay::UN{String} = abbreviation, "part-abbreviation-display"
- scoreinstrument::UN{ScoreInstrument} = ScoreInstrument(name = name, abbreviation = - abbreviation, id = id *"-I1"), "score-instrument"
- mididevice::UN{MidiDevice} = nothing, "midi-device"
- midiinstrument::MidiInstrument = MidiInstrument(), "midi-instrument"
```

Holds information about one ScorePart in a score

Each MusicXML part corresponds to a track in a Standard MIDI Format 1 file. The score-instrument elements are used when there are multiple instruments per track. The midi-device element is used to make a MIDI device or port assignment for the given track or specific MIDI instruments. Initial midi-instrument assignments may be made here as well.

scoreinstrument: See [`ScoreInstrument`](@ref) doc
mididevice: See [`MidiDevice`](@ref) doc
midiinstrument: See [`MidiInstrument`](@ref) doc

[More info](https://usermanuals.musicxml.com/MusicXML/Content/CT-MusicXML-score-part.htm)

# Examples
```julia
ScorePart(name = "Piano",midiinstrument = MidiInstrument(), id = "P1")
```
"""
@aml mutable struct ScorePart "score-part"
    id::String = "P1", att"id"
    # identification
    name::String = "Piano", "part-name"
    nameDisplay::UN{String} = name, "part-name-display"
    abbreviation::UN{String} = nothing, "part-abbreviation"
    abbreviationDisplay::UN{String} = abbreviation, "part-abbreviation-display"
    scoreinstrument::UN{ScoreInstrument} = ScoreInstrument(name = name, abbreviation = abbreviation, id = id *"-I1"), "score-instrument"
    mididevice::UN{MidiDevice} = nothing, "midi-device"
    midiinstrument::MidiInstrument = MidiInstrument(), "midi-instrument"
end
################################################################
"""
    PartList

# Arguments
```julia
- # TODO partgroup
- scoreparts::Vector{ScorePart} = [ScorePart()], "score-part"
```

Holds scoreparts and partgroup.

See [`ScorePart`](@ref) doc

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-part-list.htm)

# Example
```
PartList([
    ScorePart(name = "Piano 1", midiinstrument = MidiInstrument(), id = "P1"),
    ScorePart(name = "Piano 2", midiinstrument = MidiInstrument(), id = "P2"),
])
```
"""
@aml mutable struct PartList "part-list"
    # TODO partgroup
    scoreparts::Vector{ScorePart} = [ScorePart()], "score-part"
end
################################################################
"""
    Key

# Arguments
```julia
- # cancel
- fifths::Int8 = 0, "~"
- mode::UN{String} = nothing, "~", modeCheck
- # key-octave
```
A type to hold key information for a measure in musicxml file.

The key element represents a key signature. Both traditional and non-traditional key signatures are supported. The optional number attribute refers to staff numbers. If absent, the key signature applies to all staves in the part.

fifth: number of flats or sharps in a traditional key signature. Negative numbers are used for flats and positive numbers for sharps, reflecting the key's placement within the circle of fifths

mode:  major, minor, dorian, phrygian, lydian, mixolydian, aeolian, ionian, locrian, none

[More info](https://usermanuals.musicxml.com/MusicXML/Content/CT-MusicXML-key.htm)

# Example
```julia
Key(fifths = 0, mode = "major")
```
"""
@aml mutable struct Key "key"
    # cancel
    fifths::Int8 = 0, "~"
    mode::UN{String} = nothing, "~", modeCheck
    # key-octave
end

modeCheck(x) = in(x, ["major", "minor", "dorian", "phrygian", "lydian", "mixolydian", "aeolian", "ionian", "locrian", "none"])
################################################################
"""
    Clef

# Arguments
```julia
- sign::String = "G", "sign"
- line::Int16 = 2, "line"
- octave::UN{Int64} = nothing, "clef-octave-change"
```

A type to hold clef information for a measure in musicxml file.

Clefs are represented by a combination of sign, line, and clef-octave-change elements. Clefs appear at the start of each system unless the print-object attribute has been set to "no" or the additional attribute has been set to "yes".

sign: The sign element represents the clef symbol: G (treble), F (bass), C, percussion, TAB, jianpu, none. [More info](https://usermanuals.musicxml.com/MusicXML/Content/ST-MusicXML-clef-sign.htm)

line: Line numbers are counted from the bottom of the staff. Standard values are 2 for the G sign (treble clef), 4 for the F sign (bass clef), 3 for the C sign (alto clef) and 5 for TAB (on a 6-line staff).

octave: The clef-octave-change element is used for transposing clefs. A treble clef for tenors would have a value of -1.

number: For multi-clef instrument, give numbers for each clef.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/CT-MusicXML-clef.htm)

# Example
```julia
Clef(sign = "TAB")
```
"""
@aml mutable struct Clef "clef"
    sign::String = "G", "~"
    line::Int16 = 2, "~"
    octave::UN{Int64} = nothing, "clef-octave-change"
    number::UN{Int64} = nothing, att"number"
end

# Standard values
function Clef(sign::String) # ; needs KeywordDispatch
    if sign == "G"
        line = 2
    elseif sign == "F"
        line = 4
    elseif sign == "C"
        line = 3
    elseif sign == "TAB"
        line = 5
    end
    return Clef(sign = sign, line = line, octave = nothing)
end
################################################################
"""
    Transpose

# Arguments
```julia
- diatonic::Int8 = 0, "~"
- chromatic::Int8 = 0, "~"
- octaveChange::UN{Int8} = nothing, "octave-change"
- double::UN{Bool} = nothing, "~"
```

A type to hold transpose information for a measure in musicxml file.

If the part is being encoded for a transposing instrument in written vs. concert pitch, the transposition must be encoded in the transpose element using the transpose type.

diatonic: The diatonic element specifies the number of pitch steps needed to go from written to sounding pitch. This allows for correct spelling of enharmonic transpositions.

chromatic: The chromatic element represents the number of semitones needed to get from written to sounding pitch. This value does not include octave-change values; the values for both elements need to be added to the written pitch to get the correct sounding pitc

octaveChange: The octave-change element indicates how many octaves to add to get from written pitch to sounding pitch.

double: If the double element is present, it indicates that the music is doubled one octave down from what is currently written (as is the case for mixed cello / bass parts in orchestral literature).

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-transpose.htm)

# Example
```
Transpose() # default values
```
"""
@aml mutable struct Transpose "transpose"
    diatonic::Int8 = 0, "~"
    chromatic::Int8 = 0, "~"
    octaveChange::UN{Int8} = nothing, "octave-change"
    double::UN{Bool} = nothing, "~"
end
################################################################
"""
    Time

# Arguments
```julia
- beats::Int8 = 4, "~"
- beatType::Int8 = 4, "beat-type"
- # interchangeable
```

Time signatures are represented by the beats element for the numerator and the beat-type element for the denominator.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-time.htm)

```julia
Time(beats=6, beattype = 8)
```
"""
@aml mutable struct Time "time"
    beats::Int8 = 4, "~"
    beattype::Int8 = 4, "beat-type"
    # interchangeable
end
################################################################
"""
    Attributes

# Arguments
```julia
- time::Time = Time(), "~"
- divisions::Int16 = 1, "~"
- clef::UN{Vector{Clef}} = nothing, "~"
- staves::UN{UInt16} = nothing, "~"
- key::Key = Key(), "~"
- transpose::UN{Transpose} = nothing, "~"
- instruments::UN{UInt16} = nothing, "~"
```

A type to hold the data for the attributes of a musicxml measure

The attributes element contains musical information that typically changes on measure boundaries. This includes key and time signatures, clefs, transpositions, and staving. When attributes are changed mid-measure, it affects the music in score order, not in MusicXML document order.

key: See [`Key`](@ref) doc

divisions: The divisions element indicates how many divisions per quarter note are used to indicate a note's duration. `Note.duration/Attributes.divisions` gives the actual duration of a note. Should be given based on the shortest note in the measure. If you use 16th notes give `4`, and give duration for the notes as 1. For example, if `divisions = 2`, and we have a `Note` with `duration = 1` this is an eighth note duration. The default value is `4`, which means each `duration = 1` means a 16th note.

Duration and divisions are used directly for generating sound output, so they must be chosen to take tuplets into account. Using a divisions element lets us use just one number to represent a duration for each note in the score, while retaining the full power of a fractional representation. If maximum compatibility with Standard MIDI 1.0 files is important, do not have the divisions value exceed 16383.

time: See [`Time`](@ref) doc

staves: The staves element is used if there is more than one staff represented in the given part (e.g., 2 staves for typical piano parts). If absent, a value of 1 is assumed. Staves are ordered from top to bottom in a part in numerical order, with staff 1 above staff 2.

instruments: The instruments element is only used if more than one instrument is represented in the part (e.g., oboe I and II where they play together most of the time). If absent, a value of 1 is assumed.

clef: See [`Clef`](@ref) doc

[More info](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-attributes.htm)
"""
@aml mutable struct Attributes "attributes"
    time::Time = Time(), "~"
    divisions::Int16 = 4, "~"
    clef::UN{Vector{Clef}} = nothing, "~"
    staves::UN{UInt16} = nothing, "~"
    key::Key = Key(), "~"
    transpose::UN{Transpose} = nothing, "~"
    instruments::UN{UInt16} = nothing, "~"
end
################################################################
"""
    Grace

# Arguments
```julia
steal_time_previous::UN{Float64} = nothing, "steal-time-previous", percent
steal_time_following::UN{Float64} = nothing, "steal-time-following", percent
# TODO make-time "make-time"
slash::UN{YN} = nothing, "~"
```

steal_time_previous: The steal-time-previous attribute indicates the percentage of time to steal from the previous note for the grace note.

steal_time_following: The steal-time-following attribute indicates the percentage of time to steal from the following note for the grace note, as for appoggiaturas.

slash: The slash attribute for a grace note is yes for slashed eighth notes.

[More info](https://usermanuals.musicxml.com/MusicXML/Content/CT-MusicXML-grace.htm)
```
"""
@aml mutable struct Grace empty"grace"
    steal_time_previous::UN{Float64} = nothing, "steal-time-previous", percent
    steal_time_following::UN{Float64} = nothing, "steal-time-following", percent
    # TODO make-time "make-time"
    slash::UN{YN} = nothing, "~"
end

################################################################
"""
    Pitch

# Arguments
```julia
- step::String = "C", "~"
- alter::UN{Int16} = nothing, "~"
- octave::Int8 = 1, "~"
```

Holds musicxml pitch data. MusicXML pitch data is represented as a combination of the step of the diatonic scale, the chromatic alteration, and the octave.

Use step, alter, octave = pitch2xml(midipitch) and  midipitch = xml2pitch(step, alter, octave)
for conversions between midi pitch and musicxml pitch

![Step Alter Octave on Staff](../pitchesonstaff.png)
![Pitch on Guitar](../pitchesonguitar.png)
![Pitch on Full Keyboard](../fullpiano.gif)

step: The step type represents a step of the diatonic scale, represented using the English letters A through G.

alter: The alter element represents chromatic alteration in number of semitones (e.g., -1 for flat, 1 for sharp). Decimal values like 0.5 (quarter tone sharp) are used for microtones. The octave element is represented by the numbers 0 to 9, where 4 indicates the octave started by middle C.  In the first example below, notice an accidental element is used for the third note, rather than the alter element, because the pitch is not altered from the the pitch designated to that staff position by the key signature.

octave: Octaves are represented by the numbers 0 to 9, where 4 indicates the octave started by middle C.

[MoreInfo](https://usermanuals.musicxml.com/MusicXML/Content/EL-MusicXML-pitch.htm)

"""
@aml mutable struct Pitch "pitch"
    step::String = "C", "~"
    alter::UN{Int16} = nothing, "~"
    octave::Int8 = 1, "~"
end

function Pitch(midipitch::UInt8)
    step, alter, octave = pitch2xml(midipitch)
    return Pitch(step = step, alter = alter, octave = octave)
end

################################################################
"""
    Rest

# Arguments
```julia
measure::UN{YN} = nothing, att"~"
dispStep::UN{String} = nothing, "display-step"
dispOctave::UN{Int8} = nothing, "display-octave"
```

The rest element indicates notated rests or silences. Rest elements are usually empty, but placement on the staff can be specified using display-step and display-octave elements. If the measure attribute is set to yes, this indicates this is a complete measure rest.

The display-step-octave group contains the sequence of elements used by both the rest and unpitched elements. This group is used to place rests and unpitched elements on the staff without implying that these elements have pitch. Positioning follows the current clef. If percussion clef is used, the display-step and display-octave elements are interpreted as if in treble clef, with a G in octave 4 on line 2. If not present, the note is placed on the middle line of the staff, generally used for a one-line staff.

"""
@aml mutable struct Rest empty"rest"
    measure::UN{YN} = nothing, att"~"
    dispStep::UN{String} = nothing, "display-step"
    dispOctave::UN{Int8} = nothing, "display-octave"
end
################################################################
"""
    Unpitched

# Arguments
```julia
measure::UN{YN} = nothing, att"~"
dispStep::UN{String} = nothing, "display-step"
dispOctave::UN{Int8} = nothing, "display-octave"
```

The unpitched type represents musical elements that are notated on the staff but lack definite pitch, such as unpitched percussion and speaking voice.
"""
@aml mutable struct Unpitched empty"unpitched"
    measure::UN{YN} = nothing, att"~"
    dispStep::UN{String} = nothing, "display-step"
    dispOctave::UN{Int8} = nothing, "display-octave"
end

################################################################
"""
    Chord

If become present shows a chord
"""
@aml mutable struct Chord empty"chord"
end

################################################################
"""
    Note

# Arguments
```julia
- grace::UN{Grace} = nothing, "~"
- pitch::UN{Pitch} = nothing, "~"
- rest::UN{Rest} = nothing, "~"
- unpitched::UN{Unpitched} = nothing, "~"
- duration::UInt, "~"
- chord::UN{Chord} = nothing, "~"
- # voice
- type::UN{String} = nothing, "~"
- accidental::UN{String} = nothing, "~"
- tie::UN{String} = nothing, "~" # start, stop, nothing TODO
```

Notes are the most common type of MusicXML data. The MusicXML format keeps the MuseData distinction between elements used for sound information and elements used for notation information (e.g., tie is used for sound, tied for notation). Thus grace notes do not have a duration element. Cue notes have a duration element, as do forward elements, but no tie elements. Having these two types of information available can make interchange considerably easier, as some programs handle one type of information much more readily than the other.

grace: See [`Grace`](@ref) doc

pitch: See [`Pitch`](@ref) doc

duration : See [`MIDI.Note`](https://juliamusic.github.io/JuliaMusic_documentation.jl/latest/midi/notes/) doc

type: Type indicates the graphic note type, Valid values (from shortest to longest) are 1024th, 512th, 256th, 128th, 64th, 32nd, 16th, eighth, quarter, half, whole, breve, long, and maxima. The size attribute indicates full, cue, or large size, with full the default for regular notes and cue the default for cue and grace notes.

accidental: The accidental type represents actual notated accidentals. Editorial and cautionary indications are indicated by attributes. Values for these attributes are "no" if not present. Specific graphic display such as parentheses, brackets, and size are controlled by the level-display attribute group. Empty accidental objects are not allowed. If no accidental is desired, it should be omitted. sharp, flat, natural, double sharp, double flat, parenthesized accidental

tie:

[More info](https://usermanuals.musicxml.com/MusicXML/Content/CT-MusicXML-note.htm)
"""
@aml mutable struct Note "note"
    grace::UN{Grace} = nothing, "~"
    pitch::UN{Pitch} = nothing, "~"
    rest::UN{Rest} = nothing, "~"
    unpitched::UN{Unpitched} = nothing, "~"
    duration::UN{UInt} = grace === nothing ? nothing : 1, "~"
    chord::UN{Chord} = nothing, "~"
    # voice
    type::UN{String} = nothing, "~"
    accidental::UN{String} = nothing, "~"
    tie::UN{String} = nothing, "~" # start, stop, nothing TODO
end
################################################################
"""
    Measure

# Arguments
```julia
- attributes::UN{Attributes} = nothing, "~"
- notes::Vector{Note}, "note"
```

A type to hold the data for a musicxml measure

attributes: See [`Attributes`](@ref) doc
notes: See [`Note`](@ref) doc

"""
@aml mutable struct Measure "measure"
    attributes::UN{Attributes} = nothing, "~"
    notes::Vector{Note}, "note"
end
################################################################
"""
    Part

# Arguments
```julia
- id::String = "P1", att"~"
- measures::Vector{Measure}, "measure"
```

A type to hold the data for a part in musicxml file.

measures: See [`Measure`](@ref) doc

"""
@aml mutable struct Part "part"
    id::String = "P1", att"~"
    measures::Vector{Measure}, "measure"
end
################################################################
"""
    ScorePartwise

# Arguments
```julia
- # TODO identification
- # TODO defaults
- partlist::PartList, "part-list"
- parts::Vector{Part}, "part"
```

A type to hold the data for a musicxml file.
"""
@aml mutable struct ScorePartwise doc"score-partwise"
    # TODO identification
    # TODO defaults
    partlist::PartList, "part-list"
    parts::Vector{Part}, "part"
end


text:
using ShapML
using Documenter

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

makedocs(
    sitename = "ShapML",
    authors = "Nickalus Redell",
    doctest = true,
    format = Documenter.HTML(
        highlights = ["yaml"]
        ),
    pages = [
        "Introduction" => "index.md",
        "Vignettes" => Any["Stochastic vs. TreeSHAP" => "vignettes/consistency.md"],
        "Functions" => "functions/functions.md"
    ]
)

deploydocs(
    repo = "github.com/nredell/ShapML.jl",
    devbranch="master",
)

ext: jl

max_stars_repo_name: nredell/shapML

max_issues_repo_name: nredell/shapML

max_forks_repo_name: nredell/shapML


text:
const JNI_VERSION_1_1 =  convert(Cint, 0x00010001)
const JNI_VERSION_1_2 =  convert(Cint, 0x00010002)
const JNI_VERSION_1_4 =  convert(Cint, 0x00010004)
const JNI_VERSION_1_6 =  convert(Cint, 0x00010006)
const JNI_VERSION_1_8 =  convert(Cint, 0x00010008)

const JNI_TRUE = convert(Cchar, 1)
const JNI_FALSE = convert(Cchar, 0)

# Return Values
const JNI_OK           = convert(Cint, 0)               #/* success */
const JNI_ERR          = convert(Cint, -1)              #/* unknown error */
const JNI_EDETACHED    = convert(Cint, -2)              #/* thread detached from the VM */
const JNI_EVERSION     = convert(Cint, -3)              #/* JNI version error */
const JNI_ENOMEM       = convert(Cint, -4)              #/* not enough memory */
const JNI_EEXIST       = convert(Cint, -5)              #/* VM already created */
const JNI_EINVAL       = convert(Cint, -6)              #/* invalid arguments */

function javahome_winreg()
    try
        keypath = "SOFTWARE\\JavaSoft\\Java Runtime Environment"
        value = querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "CurrentVersion")
        keypath *= "\\"*value
        return querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "JavaHome")
    catch
        try
            keypath = "SOFTWARE\\JavaSoft\\Java Development Kit"
            value = querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "CurrentVersion")
            keypath *= "\\"*value
            return querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "JavaHome")
        catch
            error("Cannot find an installation of Java in the Windows Registry.  Please install Java from https://java.com/en/download/manual.jsp and subsequently reload JavaCall.")
        end
    end
end

@static is_unix() ? (global const libname = "libjvm") : (global const libname = "jvm")

function findjvm()
    javahomes = Any[]
    libpaths = Any[]

    if haskey(ENV,"JAVA_HOME")
        push!(javahomes,ENV["JAVA_HOME"])
    else
        @static is_windows() ? ENV["JAVA_HOME"] = javahome_winreg() : nothing
        @static is_windows() ? push!(javahomes, ENV["JAVA_HOME"]) : nothing
    end

    if isfile("/usr/libexec/java_home")
        push!(javahomes,chomp(readstring(`/usr/libexec/java_home`)))
    end

    if isdir("/usr/lib/jvm/default-java/")
        push!(javahomes, "/usr/lib/jvm/default-java/")
    end

    push!(libpaths,pwd())
    for n in javahomes
        @static if is_windows()
            push!(libpaths, joinpath(n, "bin", "server"))
            push!(libpaths, joinpath(n, "jre", "bin", "server"))
            push!(libpaths, joinpath(n, "bin", "client"))
        end
        @static if is_linux()
             if Sys.WORD_SIZE==64
                push!(libpaths, joinpath(n, "jre", "lib", "amd64", "server"))
			elseif Sys.WORD_SIZE==32
                push!(libpaths, joinpath(n, "jre", "lib", "i386", "server"))
             end
         end
        push!(libpaths, joinpath(n, "jre", "lib", "server"))
    end

    ext = @static is_windows()?"dll":(@static is_apple()?"dylib":"so")
    ext = "."*ext

    try
        for n in libpaths
            libpath = joinpath(n,libname*ext);
            if isfile(libpath)
                if is_windows()
                    bindir = dirname(dirname(libpath))
                    m = filter(x -> ismatch(r"msvcr(?:.*).dll",x), readdir(bindir))
                    Libdl.dlopen(joinpath(bindir,m[1]))
                end
                global libjvm = Libdl.dlopen(libpath)
                println("Loaded $libpath")
                return
            end
        end
    end

    errorMsg =
        [
         "Cannot find java library $libname$ext\n",
         "Search Path:"
         ];
    for path in libpaths
        push!(errorMsg,"\n   $path")
    end
    error(reduce(*,errorMsg));
end

immutable JavaVMOption
    optionString::Ptr{UInt8}
    extraInfo::Ptr{Void}
end

immutable JavaVMInitArgs
    version::Cint
    nOptions::Cint
    options::Ptr{JavaVMOption}
    ignoreUnrecognized::Cchar
end

@static is_unix() ? (const sep = ":") : nothing
@static is_windows() ? (const sep = ";") : nothing
cp=Array(String, 0)
opts=Array(String, 0)
addClassPath(s::String) = isloaded()?warn("JVM already initialised. This call has no effect"): push!(cp, s)
addOpts(s::String) = isloaded()?warn("JVM already initialised. This call has no effect"): push!(opts, s)

function init()
    if isempty(cp)
        init(opts)
    else
        init(vcat(opts, reduce((x,y)->string(x,sep,y),"-Djava.class.path=$(cp[1])",cp[2:end])))
    end
end

isloaded() = isdefined(JavaCall, :jnifunc) && isdefined(JavaCall, :penv) && penv != C_NULL

assertloaded() = isloaded()?nothing:error("JVM not initialised. Please run init()")
assertnotloaded() = isloaded()?error("JVM already initialised"):nothing

# Pointer to pointer to pointer to pointer alert! Hurrah for unsafe load
function init{T<:AbstractString}(opts::Array{T, 1})
    assertnotloaded()
    opt = Array(JavaVMOption, length(opts))
    for i in 1:length(opts)
        opt[i]=JavaVMOption(pointer(opts[i]), C_NULL)
    end
    ppjvm=Array(Ptr{JavaVM},1)
    ppenv=Array(Ptr{JNIEnv},1)
    vm_args = JavaVMInitArgs(JNI_VERSION_1_6, convert(Cint, length(opts)), convert(Ptr{JavaVMOption}, pointer(opt)), JNI_TRUE)

    res = ccall(create, Cint, (Ptr{Ptr{JavaVM}}, Ptr{Ptr{JNIEnv}}, Ptr{JavaVMInitArgs}), ppjvm, ppenv, &vm_args)
    if res < 0; error("Unable to initialise Java VM: $(res)"); end
    global penv = ppenv[1]
    global pjvm = ppjvm[1]
    jnienv=unsafe_load(penv)
    jvm = unsafe_load(pjvm)
    global jvmfunc = unsafe_load(jvm.JNIInvokeInterface_)
    global jnifunc = unsafe_load(jnienv.JNINativeInterface_) #The JNI Function table
    return
end

function destroy()
    if (!isdefined(JavaCall, :penv) || penv == C_NULL) ; error("Called destroy without initialising Java VM"); end
    res = ccall(jvmfunc.DestroyJavaVM, Cint, (Ptr{Void},), pjvm)
    if res < 0; error("Unable to destroy Java VM"); end
    global penv=C_NULL; global pjvm=C_NULL;
end


text:
using ITensors, Test

a = [-0.1, -0.12]
@test ITensors.truncate!(a) == (0., 0.)
@test length(a) == 1
a = [0.1, 0.01, 1e-13]
@test ITensors.truncate!(a,absoluteCutoff=true,cutoff=1e-5) == (1e-13, (0.01 + 1e-13)/2)
@test length(a) == 2

i = Index(2,"i")
j = Index(2,"j")
A = randomITensor(i,j)
@test_throws ArgumentError factorize(A, i, dir="fakedir")


text:
# generalise this into coupler-specific function

struct ConservationCheck{A}
    ρe_tot_atmos::A
    ρe_tot_slab::A
end

function check_conservation(cs, atmos_sim, slab_sim)
    atmos_field = atmos_sim.integrator.u.thermodynamics.ρe_tot
    slab_field = get_slab_energy(slab_sim)

    ρe_tot_atmos = sum(atmos_field)
    ρe_tot_slab = sum(slab_field)

    push!(cs.ρe_tot_atmos, ρe_tot_atmos)
    push!(cs.ρe_tot_slab, ρe_tot_slab)
end

function get_slab_energy(slab_sim)
    ρe_tot = slab_sim.params.ρ .* slab_sim.params.c .* slab_sim.integrator.u.T_sfc .* slab_sim.params.h
end

# struct ConservationCheck{F, E, I, S}
#     fields::F
#     exception::E
#     interval::I
#     show::S
# end

max_stars_repo_name: CliMA/ClimaCoupler


text:
########################################################################################
# OpenGL Arrays
"""
Transform julia datatypes to opengl enum type
"""
julia2glenum(x::Type{T}) where {T <: FixedPoint} = julia2glenum(FixedPointNumbers.rawtype(x))
# julia2glenum(x::Type{OffsetInteger{O, T}}) where {O, T} = julia2glenum(T)
julia2glenum(::Type{GLubyte})  = GL_UNSIGNED_BYTE
julia2glenum(::Type{GLbyte})   = GL_BYTE
julia2glenum(::Type{GLuint})   = GL_UNSIGNED_INT
julia2glenum(::Type{GLushort}) = GL_UNSIGNED_SHORT
julia2glenum(::Type{GLshort})  = GL_SHORT
julia2glenum(::Type{GLint})    = GL_INT
julia2glenum(::Type{GLfloat})  = GL_FLOAT
julia2glenum(::Type{GLdouble}) = GL_DOUBLE
julia2glenum(::Type{Float16})  = GL_HALF_FLOAT
function julia2glenum(::Type{T}) where T
    glasserteltype(T)
    julia2glenum(eltype(T))
end

gl_convert(a::T) where {T <: NATIVE_TYPES} = a
gl_convert(::Type{T}, a::NATIVE_TYPES; kw_args...) where {T <: NATIVE_TYPES} = a

gl_promote(x::Type{T}) where {T <: Integer} = Cint
gl_promote(x::Type{Union{Int16, Int8}}) = x

gl_promote(x::Type{T}) where {T <: Unsigned} = Cuint
gl_promote(x::Type{Union{UInt16, UInt8}}) = x

gl_promote(x::Type{T}) where {T <: AbstractFloat} = Float32
gl_promote(x::Type{Float16}) = x

gl_promote(x::Type{T}) where {T <: Normed} = N0f32
gl_promote(x::Type{N0f16}) = x
gl_promote(x::Type{N0f8}) = x

gl_promote(x::Type{Bool}) = GLboolean

function gl_promote(x::Type{T}) where T
    glasserteltype(T)
    similar_type(T, gl_promote(eltype(T)))
end

gl_convert(x::T) where {T <: Number} = gl_promote(T)(x)
gl_convert(s::Nothing) = s

isa_gl_struct(x::Array) = false
function isa_gl_struct(x::T) where T
    !isleaftype(T) && return false
    if T <: Tuple
        return false
    end
    fnames = fieldnames(T)
    !isempty(fnames) && all(name -> isleaftype(fieldtype(T, name)) && isbits(getfield(x, name)), fnames)
end
function gl_convert_struct(x::T, uniform_name::Symbol) where T
    if isa_gl_struct(x)
        return Dict{Symbol, Any}(map(fieldnames(x)) do name
            (Symbol("$uniform_name.$name") => gl_convert(getfield(x, name)))
        end)
    else
        @error "can't convert $x to a OpenGL type. Make sure all fields are of a concrete type and isbits(FieldType)-->true"
    end
end

function gl_convert(::Type{T}, a::Array{X, N}; kw_args...) where {T <: GPUArray, X, N}
    T(map(gl_promote(X), a); kw_args...)
end
function gl_convert(::Type{T}, a::Vector{Array{X, 2}}; kw_args...) where {T <: Texture, X}
    T(a; kw_args...)
end

gl_convert(f::Function, a) = f(a)


text:
using Printf

println("hello world :-) !")

max_stars_repo_name: mathbiol/julia

max_issues_repo_name: mathbiol/julia

max_forks_repo_name: mathbiol/julia


text:
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

import Cbc.CbcCInterface
const CbcCI = CbcCInterface

using SparseArrays

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::CbcCI.CbcModel
    silent::Bool
    # Cache the params so they can be reset on `empty!`.
    params::Dict{String, String}
    # Cache the objective constant (if there is one).
    objective_constant::Float64

    """
        Optimizer(; kwargs...)

    Create a new Cbc Optimizer.
    """
    function Optimizer(; kwargs...)
        model = CbcCI.CbcModel()
        optimizer = new(model, false, Dict{String, String}(), 0.0)
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawParameter(key), value)
        end
        return optimizer
    end
end

function MOI.supports(optimizer::Optimizer, param::MOI.RawParameter)
    # FIXME find a way to check if `param.name` is recognized by Cbc
    #       a list of parameters can be obtained by running `cbc`
    #       and issuing the command `???` to the CLI.
    #       We could make a `Set` out of this list but it would be
    #       better to have a Cbc `hasParameter` to the Cbc C API.
    return true
end
function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    if !MOI.supports(optimizer, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    name = string(param.name)
    optimizer.params[name] = string(value)
    if !(optimizer.silent && name == "logLevel")
        CbcCI.setParameter(optimizer.inner, name, string(value))
    end
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    # TODO: This gives a poor error message if the name of the parameter is invalid.
    return optimizer.params[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    if value != optimizer.silent
        if value
            CbcCI.setParameter(optimizer.inner, "logLevel", "0")
        else
            log_level = get(optimizer.params, "logLevel", "1")
            CbcCI.setParameter(optimizer.inner, "logLevel", log_level)
        end
    end
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, value)
    if value === nothing
        delete!(optimizer.params, "seconds")
        CbcCI.setParameter(optimizer.inner, "seconds", "??")
    else
        MOI.set(optimizer, MOI.RawParameter("seconds"), value)
    end
end
function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
    value = get(optimizer.params, "seconds", nothing)
    if value === nothing
        return value
    else
        return parse(Float64, value)
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "COIN Branch-and-Cut (Cbc)"

function MOI.empty!(model::Optimizer)
    model.inner = CbcCI.CbcModel()
    model.objective_constant = 0.0
    for (name, value) in model.params
        CbcCI.setParameter(model.inner, name, value)
    end
    if model.silent
        CbcCI.setParameter(model.inner, "logLevel", "0")
    end
    return
end

function MOI.is_empty(model::Optimizer)
    return CbcCI.getNumCols(model.inner) == 0 &&
        CbcCI.getNumRows(model.inner) == 0
end

mutable struct CbcModelFormat
    num_rows::Int
    num_cols::Int
    # (row_idx, col_idx, values) are the sparse elements of the constraint matrix.
    row_idx::Vector{Int}
    col_idx::Vector{Int}
    values::Vector{Float64}
    # Constraint bounds.
    row_lb::Vector{Float64}
    row_ub::Vector{Float64}
    # Variable bounds.
    col_lb::Vector{Float64}
    col_ub::Vector{Float64}
    # Columns that are binary or integer
    binary::Vector{Int}
    integer::Vector{Int}
    # SOS1 constraints
    sos1_weights::Vector{Vector{Float64}}
    sos1_indices::Vector{Vector{Int32}}
    # SOS2 constraints
    sos2_weights::Vector{Vector{Float64}}
    sos2_indices::Vector{Vector{Int32}}
    # Objective coefficients.
    obj::Vector{Float64}
    objective_constant::Float64

    function CbcModelFormat(num_rows::Integer, num_cols::Integer)
        # An `InexactError` might occur if `num_rows` or `num_cols` is too
        # large, e.g., if `num_cols isa Int64` and is larger than 2^31 on
        # 32-bit hardware.
        new(
            num_rows,
            num_cols,
            Int[],  # row_idx
            Int[],  # col_idx
            Float64[],  # values
            fill(-Inf, num_rows),  # row_lb
            fill(Inf, num_rows),  # row_ub
            fill(-Inf, num_cols),  # col_lb
            fill(Inf, num_cols),  # col_ub
            Int[],  # binary
            Int[],  # integer
            Vector{Vector{Float64}}(), # SOS1 weights
            Vector{Vector{Int32}}(),     # SOS1 column indices
            Vector{Vector{Float64}}(), # SOS2 weights
            Vector{Vector{Int32}}(),     # SOS2 column indices
            fill(0.0, num_cols),  # obj
            0.0  # objective_constant
        )
    end
end

###
### SingleVariable-in-{EqualTo,LessThan,GreaterThan,Interval,ZeroOne,Integer}
###

function column_value(map::MOIU.IndexMap, index::MOI.VariableIndex)
    return map.varmap[index].value
end

function column_value(map::MOIU.IndexMap, func::MOI.SingleVariable)
    return column_value(map, func.variable)
end

function MOI.supports_constraint(
        ::Optimizer, ::Type{MOI.SingleVariable}, ::Type{<:Union{
            MOI.EqualTo{Float64}, MOI.LessThan{Float64},
            MOI.GreaterThan{Float64}, MOI.Interval{Float64},
            MOI.ZeroOne, MOI.Integer}})
    return true
end

MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.SOS1{Float64}}) = true

MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.SOS2{Float64}}) = true

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, set::MOI.EqualTo)
    column = column_value(mapping, func)
    model.col_lb[column] = set.value
    model.col_ub[column] = set.value
    return
end

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, set::MOI.LessThan)
    column = column_value(mapping, func)
    model.col_ub[column] = set.upper
    return
end

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, set::MOI.GreaterThan)
    column = column_value(mapping, func)
    model.col_lb[column] = set.lower
    return
end

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, set::MOI.Interval)
    column = column_value(mapping, func)
    model.col_lb[column] = set.lower
    model.col_ub[column] = set.upper
    return
end

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, ::MOI.ZeroOne)
    push!(model.binary, column_value(mapping, func))
    return
end

function load_constraint(
        ::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.SingleVariable, ::MOI.Integer)
    push!(model.integer, column_value(mapping, func))
    return
end

###
### ScalarAffineFunction-in-{EqualTo, LessThan, GreaterThan, Interval}
###

function MOI.supports_constraint(
        ::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{<:Union{
            MOI.EqualTo{Float64}, MOI.LessThan{Float64},
            MOI.GreaterThan{Float64}, MOI.Interval{Float64}}})
    return true
end

function add_terms(
        model::CbcModelFormat, mapping::MOIU.IndexMap,
        index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S},
        func::MOI.ScalarAffineFunction{Float64}) where {S}
    for term in func.terms
        push!(model.row_idx, mapping.conmap[index].value)
        push!(model.col_idx, column_value(mapping, term.variable_index))
        push!(model.values, term.coefficient)
    end
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.ScalarAffineFunction,
        set::MOI.EqualTo)
    add_terms(model, mapping, index, func)
    row = mapping.conmap[index].value
    model.row_lb[row] = set.value - func.constant
    model.row_ub[row] = set.value - func.constant
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.ScalarAffineFunction,
        set::MOI.GreaterThan)
    add_terms(model, mapping, index, func)
    row = mapping.conmap[index].value
    model.row_lb[row] = set.lower - func.constant
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.ScalarAffineFunction,
        set::MOI.LessThan)
    add_terms(model, mapping, index, func)
    row = mapping.conmap[index].value
    model.row_ub[row] = set.upper - func.constant
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.ScalarAffineFunction,
        set::MOI.Interval)
    add_terms(model, mapping, index, func)
    row = mapping.conmap[index].value
    model.row_ub[row] = set.upper - func.constant
    model.row_lb[row] = set.lower - func.constant
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.VectorOfVariables,
        set::MOI.SOS1{Float64})
    var_indices = Int32[v.value for v in func.variables]
    push!(model.sos1_weights, set.weights)
    push!(model.sos1_indices, var_indices)
    return
end

function load_constraint(
        index::MOI.ConstraintIndex, model::CbcModelFormat,
        mapping::MOIU.IndexMap, func::MOI.VectorOfVariables,
        set::MOI.SOS2{Float64})
    var_indices = [v.value for v in func.variables]
    push!(model.sos2_weights, set.weights)
    push!(model.sos2_indices, var_indices)
    return
end

###
### ObjectiveSense
###

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    if sense == MOI.MAX_SENSE
        CbcCI.setObjSense(model.inner, -1)
    else ## Other senses are set as minimization (cbc default)
        CbcCI.setObjSense(model.inner, 1)
    end
end

###
### ObjectiveFunction{ScalarAffineFunction}
###

function MOI.supports(
        ::Optimizer, ::MOI.ObjectiveFunction{
            <:Union{MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable}})
    return true
end

function load_objective(model::CbcModelFormat, mapping::MOIU.IndexMap,
                        func::MOI.ScalarAffineFunction)
    # We need to increment values of objective function with += to handle
    # cases like $x_1 + x_2 + x_1$. This is safe because objective function
    # is initialized with zeros in the constructor and `load_objective` only
    # gets called once.
    for term in func.terms
        column = column_value(mapping, term.variable_index)
        model.obj[column] += term.coefficient
    end
    model.objective_constant = func.constant
    return
end

function load_objective(model::CbcModelFormat, mapping::MOIU.IndexMap,
                        func::MOI.SingleVariable)
    column = column_value(mapping, func)
    model.obj[column] = 1.0
    model.objective_constant = 0.0
    return
end

###
### This main copy_to function.
###

"""
    create_constraint_indices(src::MOI.ModelLike, mapping::MOIU.IndexMap)

Create a new set of constraint indices. Importantly:

 - The `.value` field of each `ConstraintIndex{ScalarAffineFunction, S}` is its
   row in the constraint matrix.
 - The `.value` field of each `ConstraintIndex{SingleVariable, S}` is the
   column of the associated variable in the constraint matrix.
"""
function create_constraint_indices(src::MOI.ModelLike, mapping::MOIU.IndexMap)
    num_rows = 0
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        for index in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
            if F == MOI.SingleVariable
                c_func = MOI.get(src, MOI.ConstraintFunction(), index)
                column = mapping.varmap[c_func.variable].value
                mapping.conmap[index] = MOI.ConstraintIndex{F, S}(column)
            else
                num_rows += 1
                mapping.conmap[index] = MOI.ConstraintIndex{F, S}(num_rows)
            end
        end
    end
    return num_rows
end

"""
    create_variable_indices(src::MOI.ModelLike, mapping::MOIU.IndexMap)

Create a new set of variable indices. Importantly, the `.value` field of each
`VariableIndex` is its column in the constraint matrix.
"""
function create_variable_indices(src::MOI.ModelLike, mapping::MOIU.IndexMap)
    var_index = MOI.get(src, MOI.ListOfVariableIndices())
    for (column, variable) in enumerate(var_index)
        mapping.varmap[variable] = MOI.VariableIndex(column)
    end
    return length(var_index)
end

function MOI.copy_to(cbc_dest::Optimizer, src::MOI.ModelLike;
                     copy_names=false)
    # Begin by checking that we support all the constraints. This means that
    # we don't have to do this later.
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        if !(MOI.supports_constraint(cbc_dest, F, S))
            throw(MOI.UnsupportedConstraint{F, S}(
                "Cbc.Optimizer does not support constraints of type $F -in- $S."
            ))
        end
    end

    # A mapping between the indices of `src` and the ones we create.
    mapping = MOIU.IndexMap()

    num_cols = create_variable_indices(src, mapping)
    num_rows = create_constraint_indices(src, mapping)

    # Create a new temporary storage instance.
    tmp_model = CbcModelFormat(num_rows, num_cols)

    # Copy the constraints out of `src` into `tmp_model`.
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        for index in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
            load_constraint(
                index,
                tmp_model,
                mapping,
                MOI.get(src, MOI.ConstraintFunction(), index),
                MOI.get(src,  MOI.ConstraintSet(), index)
            )
        end
    end

    # Since Cbc doesn't have an explicit binary variable, we need to add [0, 1]
    # bounds and make it integer (which is done at the end of this function).
    for column in tmp_model.binary
        tmp_model.col_lb[column] = max(tmp_model.col_lb[column], 0.0)
        tmp_model.col_ub[column] = min(tmp_model.col_ub[column], 1.0)
    end

    # Copy the objective function.
    objective_function_type = MOI.get(src, MOI.ObjectiveFunctionType())
    if !MOI.supports(cbc_dest, MOI.ObjectiveFunction{objective_function_type}())
        error("Objective function type $(objective_function_type) not supported.")
    end
    objective_function = MOI.get(
        src, MOI.ObjectiveFunction{objective_function_type}())
    load_objective(tmp_model, mapping, objective_function)

    sense = MOI.get(src, MOI.ObjectiveSense())
    MOI.set(cbc_dest, MOI.ObjectiveSense(), sense)

    # Load the problem into Cbc.
    CbcCI.loadProblem(
        cbc_dest.inner,
        sparse(tmp_model.row_idx, tmp_model.col_idx, tmp_model.values,
               tmp_model.num_rows, tmp_model.num_cols),
        tmp_model.col_lb, tmp_model.col_ub,
        tmp_model.obj,
        tmp_model.row_lb, tmp_model.row_ub
    )

    cbc_dest.objective_constant = tmp_model.objective_constant

    # Set the integer variables.
    for column in tmp_model.integer
        CbcCI.setInteger(cbc_dest.inner, column - 1)
    end

    # Set the binary variables.
    for column in tmp_model.binary
        CbcCI.setInteger(cbc_dest.inner, column - 1)
    end

    # SOS constraints
    # type 1
    for i in eachindex(tmp_model.sos1_weights)
        CbcCI.addSOS(cbc_dest.inner, 1, Cint[1, length(tmp_model.sos1_weights[i]) + 1], tmp_model.sos1_indices[i], tmp_model.sos1_weights[i], 1)
    end

    # type 2
    for i in eachindex(tmp_model.sos2_weights)
        CbcCI.addSOS(cbc_dest.inner, 1, Cint[1, length(tmp_model.sos2_weights[i]) + 1], tmp_model.sos2_indices[i], tmp_model.sos2_weights[i], 2)
    end

    return MOIU.IndexMap(mapping.varmap, mapping.conmap)
end

###
### Optimize and post-optimize functions
###

function MOI.optimize!(model::Optimizer)
    CbcCI.solve(model.inner)
    return
end

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return CbcCI.getNumCols(model.inner)
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    return CbcCI.getBestPossibleObjValue(model.inner) + model.objective_constant
end

function MOI.get(model::Optimizer, ::MOI.NodeCount)
    return CbcCI.getNodeCount(model.inner)
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    return CbcCI.getObjValue(model.inner) + model.objective_constant
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal,
                 ref::MOI.VariableIndex)
    primal_solution = CbcCI.unsafe_getColSolution(model.inner)
    return primal_solution[ref.value]
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal,
                 indices::Vector{MOI.VariableIndex})
    primal_solution = CbcCI.unsafe_getColSolution(model.inner)
    return [primal_solution[index.value] for index in indices]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any})
    primal_solution = CbcCI.unsafe_getRowActivity(model.inner)
    return primal_solution[index.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 index::MOI.ConstraintIndex{MOI.SingleVariable, <:Any})
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(index.value))
end

function MOI.get(model::Optimizer, object::MOI.ResultCount)
    if CbcCI.isProvenInfeasible(model.inner) ||
            CbcCI.isContinuousUnbounded(model.inner) ||
            CbcCI.isAbandoned(model.inner) ||
            CbcCI.getObjValue(model.inner) >= 1e300
        return 0
    end
    return 1
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    sense = CbcCI.getObjSense(model.inner)
    if sense == 1
        return MOI.MIN_SENSE
    elseif sense == -1
        return MOI.MAX_SENSE
    else
        error("Objective sense $(sense) not recognized.")
    end
end

struct Status <: MOI.AbstractModelAttribute end
# If we don't define this method, the `CachingOptimizer` will redirect the call
# to the cache instead of the optimizer.
MOI.is_set_by_optimize(::Status) = true

struct SecondaryStatus <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::SecondaryStatus) = true

MOI.get(model::Optimizer, ::Status) = CbcCI.status(model.inner)
MOI.get(model::Optimizer, ::SecondaryStatus) = CbcCI.secondaryStatus(model.inner)

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    # FIXME We should do something more helpful here
    return string("status = ", MOI.get(model, Status()),
                  ", secondaryStatus = ", MOI.get(model, SecondaryStatus()))
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if CbcCI.isProvenInfeasible(model.inner)
        return MOI.INFEASIBLE
    elseif CbcCI.isContinuousUnbounded(model.inner)
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif CbcCI.isNodeLimitReached(model.inner)
        return MOI.NODE_LIMIT
    elseif CbcCI.isSecondsLimitReached(model.inner)
        return MOI.TIME_LIMIT
    elseif CbcCI.isSolutionLimitReached(model.inner)
        return MOI.SOLUTION_LIMIT
    elseif CbcCI.isProvenOptimal(model.inner) ||
            CbcCI.isInitialSolveProvenOptimal(model.inner) ||
            MOI.get(model, MOI.ResultCount()) == 1
        return MOI.OPTIMAL
    elseif CbcCI.isAbandoned(model.inner)
        return MOI.INTERRUPTED
    elseif CbcCI.optimizeNotCalled(model.inner)
        return MOI.OPTIMIZE_NOT_CALLED
    else
        error("Internal error: Unrecognized solution status: ",
              MOI.get(model, MOI.RawStatusString()))
    end
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if CbcCI.isProvenOptimal(model.inner) ||
            CbcCI.isInitialSolveProvenOptimal(model.inner)
        return MOI.FEASIBLE_POINT
    elseif CbcCI.isProvenInfeasible(model.inner)
        return MOI.INFEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.get(::Optimizer, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end


text:
module NonlinearSolversTests

using Gridap.Algebra
using Gridap.Algebra: NonlinearOperatorMock

op = NonlinearOperatorMock()

ls = LUSolver()
tol = 1.e-10
maxiters = 20
nls = NewtonRaphsonSolver(ls,tol,maxiters)

x0 = zero_initial_guess(op)
x = [1.0, 3.0]
test_nonlinear_solver(nls,op,x0,x)

x0 = [2.1,2.9]
x = [2.0, 3.0]
test_nonlinear_solver(nls,op,x0,x)

end # module


text:
module NRLTB

using ForwardDiff

using JuLIP.Potentials: ZeroSitePotential
using SKTB: TBModel, SKHamiltonian, NONORTHOGONAL,
                  norbitals, ChemicalPotential
import SKTB: hop!, overlap!, onsite!, hop, overlap, onsite_grad!, GammaPoint
import JuLIP.Potentials: cutoff

export NRLTBModel, NRLHamiltonian

const BOHR = 0.52917721092::Float64  # atomic unit of length 1 Bohr = 0.52917721092 Å
const half_eV = 13.60569301::Float64     # conversion from Ha to eV:  1 eV = 27.21138602 Ha

# Norbital: 4 if(s, p) : s, px, py, pz
#           9 if(s, p, d) : s, px, py, pz, dxy, dyz, dzx, dx^2-y^2, d3z^2-r^2
# Nbond: 4 if (s, p) : ssσ, spσ, ppσ, ppπ,
#        10 if (s, p ,d) : ssσ, spσ, ppσ, ppπ, sdσ, pdσ, pdπ, ddσ, ddπ, ddδ

"""
`NRLHamiltonian `: specifies the parameters for the NRLTB hamiltonian
"""
mutable struct NRLHamiltonian{NORB, FCUT} <: SKHamiltonian{NONORTHOGONAL, NORB}
    Norbital::Int64
    Nbond::Int64
# cutoff parameters
    Rc::Float64
    lc::Float64
    fcut::FCUT
# onsite
    λ::Float64
    a::Vector{Float64}
    b::Vector{Float64}
    c::Vector{Float64}
    d::Vector{Float64}
# hopping H
    e::Vector{Float64}
    f::Vector{Float64}
    g::Vector{Float64}
    h::Vector{Float64}
# hopping M
    p::Vector{Float64}
    q::Vector{Float64}
    r::Vector{Float64}
    s::Vector{Float64}
    t::Vector{Float64}
end

cutoff(H::NRLHamiltonian) = H.Rc

# ================= NRL CUTOFF ===========
# (Modified??) cutoff function in the NRL ansatz for Hamiltonian matrix elements
#    r  : variable
#    Rc : cutoff radius
#    lc : cutoff weight
#    Mc : we changed NRL's 5.0 to 10.0, but I changed it back to 5.0 (CO)
# NB: this cutoff is CRAP - this cut-off is not even differentiable; we should
#     at least force-shift it?!?!?
# TODO: talk to Noam about this.

# TODO: probably we shoudn't vectorise cutoff_NRL_*** here but rather do it
#       at the call site where it is known whether r is a scalar or a vector

_nrlcut_(r, Rc, lc, Mc) = 1.0 / (1.0 + exp( (r-Rc) / lc + Mc ))

cutoff_NRL_original(r, Rc, lc, Mc=5.0) =
   _nrlcut_.(r, Rc, lc, Mc) .* (r .<= Rc)

cutoff_NRL_Eshift(r, Rc, lc, Mc = 5.0) =
   (_nrlcut_.(r, Rc, lc, Mc) .- 1.0/(1.0+exp(Mc))) .* (r .<= Rc)

cutoff_NRL_Fshift(r, Rc, lc, Mc=5.0) =
   (_nrlcut_.(r, Rc, lc, Mc) .- 1.0/(1.0+exp(Mc))
         .+ (exp(Mc)/lc)/(1.0+exp(Mc))^2 * (r .- Rc)) .* (r .<= Rc)

# default
cutoff_NRL = cutoff_NRL_Fshift

# Loading routines; and spec of minimal set of parameters
include("NRLTB_data.jl")

"""
`NRLTBModel(species, fs; kwargs...)`: constructs the NRL tight binding model.

* `species`: symbol or string specifying the chemical symbol
* `fs` : chemical potential

### Keyword Parameters

* `orbitals`: either `:sp` or `:spd` (the default option makes a predetermined choice)
* `hfd` = 1e-6 : finite difference step for computing hessians
* `cutoff`: NRLTB has one of the most awful cutoff multipliers in the history
of interatomic potentials: use `cutoff=:original` for the original cutoff
and enjoy the benefits of a discontinuous potential. Other options:
`cutoff = :energyshift` and `cutoff = :forceshift`.
"""
NRLTBModel(species, fs::ChemicalPotential;
           orbitals=default_orbitals(species), bzquad=GammaPoint(), hfd=1e-6,
            cutoff = :forceshift) =
   TBModel(NRLHamiltonian(species, orbitals=orbitals, cutoff=cutoff),
           ZeroSitePotential(), fs, bzquad, hfd)

# ================= HOPPING INTEGRALS =====================

nrl_hop(H::NRLHamiltonian, r, i) = (H.e[i] + (H.f[i] + H.g[i] * r) * r) * exp( - H.h[i]^2 * r)

hop(H::NRLHamiltonian, r, i) =
   half_eV * nrl_hop(H, r/BOHR, i) * H.fcut(r/BOHR, H.Rc, H.lc)

# ================= OVERLAP INTEGRALS  =====================

nrl_olap(H, r, i) = (H.p[i] + (H.q[i] + (H.r[i] + H.s[i] * r) * r) * r) * exp(-H.t[i]^2 * r)

overlap(H::NRLHamiltonian, r::Real, i::Integer) =
      nrl_olap(H, r/BOHR, i) * H.fcut(r/BOHR, H.Rc, H.lc)

# on-site overlap block
function overlap!(H::NRLHamiltonian, M_nn)
   fill!(M_nn, 0.0)
   for i = 1:norbitals(H); M_nn[i,i] = 1.0; end
   return M_nn
end

# we don't need a derivative of overlap! (onsite) since it is constant

# ================= IMPLEMENTATION OF ONSITE FUNCTION =====================

# Preliminaries: Pseudo electron density on site l : ρ_l
# r, R : distances and displacements of the neighboring atoms
# elem : NRLParams
# OUTPUT
# ρ    : return the pseudo density on site n = 1, ... , length(atm)
# note that the NRL pseudo density has ignored the self-distance
pseudoDensity(H::NRLHamiltonian, r::AbstractVector) =
   sum( exp.(- H.λ^2 * r) .* H.fcut(r, H.Rc, H.lc) )

# auxiliary functions for computing the onsite terms
nrl_os(H::NRLHamiltonian, ρ, i) =
   H.a[i] + H.b[i] * ρ^(2/3) + H.c[i] * ρ^(4/3) + H.d[i] * ρ^2

nrl_os_d(H::NRLHamiltonian, ρ, i) =
   H.b[i] * (2/3) * ρ^(-1/3) + H.c[i] * (4/3) * ρ^(1/3) + H.d[i] * 2 * ρ

function onsite!(H::NRLHamiltonian, r, _, H_nn)
   ρ = pseudoDensity(H, r / BOHR)
   fill!(H_nn, 0.0)
   for i = 1:norbitals(H)
      H_nn[i,i] = half_eV * nrl_os(H, ρ, i)
   end
   return H_nn
end

function onsite_grad!(H::NRLHamiltonian, r, R, dH_nn)
   ρ = pseudoDensity(H, r / BOHR)
   ∇ρ = ForwardDiff.gradient( r_ -> pseudoDensity(H, r_ / BOHR), r )
   fill!(dH_nn, 0.0)
   for i = 1:norbitals(H), a = 1:3, j = 1:length(r)
      dH_nn[a,i,i,j] = half_eV * nrl_os_d(H, ρ, i) * ∇ρ[j] * R[j][a] / r[j]
   end
   return dH_nn
end

end


text:
using LinearAlgebra, SparseArrays, Random, SharedArrays
import Dagger: chunks, DArray, domainchunks, treereduce_nd

@testset "treereduce_nd" begin
    xs = rand(1:10, 8,8,8)
    concats = [(x...)->cat(x..., dims=n) for n in 1:3]
    @test treereduce_nd(concats, xs) == xs
    @test treereduce_nd(reverse(concats), xs) == permutedims(xs, [3,2,1])
end

@testset "DArray constructor" begin
    x = compute(rand(Blocks(2,2), 3,3))
    @test collect(x) == DArray{Float64, 2}(x.domain, x.subdomains, x.chunks) |> collect
end

@testset "rand" begin
    function test_rand(X)
        X1 = compute(X)
        X2 = collect(X1)

        @test isa(X, Dagger.ArrayOp)
        @test isa(X1, Dagger.DArray)
        @test X2 |> size == (100, 100)
        @test all(X2 .>= 0.0)
        @test size(chunks(X1)) == (10, 10)
        @test domain(X1) == ArrayDomain(1:100, 1:100)
        @test domainchunks(X1) |> size == (10, 10)
        @test domainchunks(X1) == partition(Blocks(10, 10), ArrayDomain(1:100, 1:100))
        @test collect(X1) == collect(X1)
    end
    X = rand(Blocks(10, 10), 100, 100)
    test_rand(X)
    Xsp = sprand(Blocks(10, 10), 100, 100, 0.1)
    test_rand(Xsp)
    R = rand(Blocks(10), 20)
    r = collect(R)
    @test r[1:10] != r[11:20]
end
@testset "sum" begin
    X = ones(Blocks(10, 10), 100, 100)
    @test sum(X) == 10000
    Y = zeros(Blocks(10, 10), 100, 100)
    @test sum(Y) == 0
end

@testset "distributing an array" begin
    function test_dist(X)
        X1 = Distribute(Blocks(10, 20), X)
        @test X1 == X
        Xc = compute(X1)
        @test chunks(Xc) |> size == (10, 5)
        @test domainchunks(Xc) |> size == (10, 5)
        @test map(x->size(x) == (10, 20), domainchunks(Xc)) |> all
    end
    x = [1 2; 3 4]
    @test Distribute(Blocks(1,1), x) == x
    #test_dist(rand(100, 100))
    #test_dist(sprand(100, 100, 0.1))

    x = distribute(rand(10), 2)
    @test collect(distribute(x, 3)) == collect(x)
end

@testset "transpose" begin
    function test_transpose(X)
        x, y = size(X)
        X1 = Distribute(Blocks(10, 20), X)
        @test X1' == X'
        Xc = compute(X1')
        @test chunks(Xc) |> size == (div(y, 20), div(x,10))
        @test domainchunks(Xc) |> size == (div(y, 20), div(x, 10))
        @test map(x->size(x) == (20, 10), domainchunks(Xc)) |> all
    end
    test_transpose(rand(100, 100))
    test_transpose(rand(100, 120))
    test_transpose(sprand(100, 100, 0.1))
    test_transpose(sprand(100, 120, 0.1))
end

@testset "matrix-matrix multiply" begin
    function test_mul(X)
        tol = 1e-12
        X1 = Distribute(Blocks(10, 20), X)
        @test_throws DimensionMismatch compute(X1*X1)
        X2 = compute(X1'*X1)
        X3 = compute(X1*X1')
        @test norm(collect(X2) - X'X) < tol
        @test norm(collect(X3) - X*X') < tol
        @test chunks(X2) |> size == (2, 2)
        @test chunks(X3) |> size == (4, 4)
        @test map(x->size(x) == (20, 20), domainchunks(X2)) |> all
        @test map(x->size(x) == (10, 10), domainchunks(X3)) |> all
    end
    test_mul(rand(40, 40))

    x = rand(10,10)
    X = Distribute(Blocks(3,3), x)
    y = rand(10)
    @test norm(collect(X*y) - x*y) < 1e-13
end

@testset "matrix powers" begin
    x = compute(rand(Blocks(4,4), 16, 16))
    @test collect(x^1) == collect(x)
    @test collect(x^2) == collect(x*x)
    @test collect(x^3) == collect(x*x*x)
end

@testset "concat" begin
    m = rand(75,75)
    x = Distribute(Blocks(10,20), m)
    y = Distribute(Blocks(10,10), m)
    @test hcat(m,m) == collect(hcat(x,x))
    @test vcat(m,m) == collect(vcat(x,x))
    @test hcat(m,m) == collect(hcat(x,y))
    @test_throws DimensionMismatch compute(vcat(x,y))
end

@testset "scale" begin
    x = rand(10,10)
    X = Distribute(Blocks(3,3), x)
    y = rand(10)

    @test Diagonal(y)*x == collect(Diagonal(y)*X)
end

@testset "Getindex" begin
    function test_getindex(x)
        X = Distribute(Blocks(3,3), x)
        @test collect(X[3:8, 2:7]) == x[3:8, 2:7]
        ragged_idx = [1,2,9,7,6,2,4,5]
        @test collect(X[ragged_idx, 2:7]) == x[ragged_idx, 2:7]
        @test collect(X[ragged_idx, reverse(ragged_idx)]) == x[ragged_idx, reverse(ragged_idx)]
        ragged_idx = [1,2,9,7,6,2,4,5]
        @test collect(X[[2,7,10], :]) == x[[2,7,10], :]
        @test collect(X[[], ragged_idx]) == x[[], ragged_idx]
        @test collect(X[[], []]) == x[[], []]

        @testset "dimensionality reduction" begin
        # THESE NEED FIXING!!
            @test vec(collect(X[ragged_idx, 5])) == vec(x[ragged_idx, 5])
            @test vec(collect(X[5, ragged_idx])) == vec(x[5, ragged_idx])
            @test collect(X[5, 5]) == x[5,5]
        end
    end

    test_getindex(rand(10,10))
    test_getindex(sprand(10,10,0.5))

    y = rand(10, 10)
    xs = distribute(y, Blocks(2,2))
    for i=1:10, j=1:10
        @test compute(xs[i:j, j:i]) == y[i:j, j:i]
    end
end

@testset "cleanup" begin
    X = Distribute(Blocks(10,10), rand(10,10))
    @test collect(sin.(X)) == collect(sin.(X))
end

@testset "reducedim" begin
    x = rand(1:10, 10, 5)
    X = Distribute(Blocks(3,3), x)
    @test reduce(+, x, dims=1) == collect(reduce(+, X, dims=1))
    @test reduce(+, x, dims=2) == collect(reduce(+, X, dims=2))

    x = rand(1:10, 10, 5)
    X = Distribute(Blocks(10, 10), x)
    @test sum(x, dims=1) == collect(sum(X, dims=1))
    @test sum(x, dims=2) == collect(sum(X, dims=2))
end

@testset "setindex" begin
    x=rand(10,10)
    y=copy(x)
    y[3:8, 2:7] .= 1.0
    X = Distribute(Blocks(3,3), x)
    @test collect(setindex(X,1.0, 3:8, 2:7)) == y
    @test collect(X) == x
end

@testset "sort" begin
    x = shuffle(1:10)
    X = distribute(x, 4)
    @test collect(sort(X)) == sort(x)

    X = distribute(x, 10)
    @test collect(sort(X)) == sort(x)
    @test collect(sort(X, rev=true)) == sort(x, rev=true)

    X = distribute(x, 1)
    @test collect(sort(X)) == sort(x)
    @test collect(sort(X, rev=true)) == sort(x, rev=true)

    x = [("A",1), ("A",2), ("B",1)]
    y = distribute(x, 3)
    @test collect(sort(y)) == x

    x = ones(10)
    y = compute(Distribute(Blocks(3), x))
    #@test map(x->length(collect(x)), compute(sort(y)).chunks) == [3,3,3,1]
end

using MemPool

@testset "affinity" begin
    x = Dagger.tochunk([1:10;])
    aff = Dagger.affinity(x)
    @test aff[1] == Dagger.OSProc(myid())
    @test aff[2] == sizeof(Int)*10
    @test Dagger.tochunk(x) === x
    f = MemPool.FileRef("/tmp/d", aff[2])
    aff = Dagger.affinity(f)
    #@test length(aff) == 3
    @test (aff[1]).pid in procs()
    @test aff[2] == sizeof(Int)*10
end

@testset "show_plan" begin
    @test !isempty(Dagger.show_plan(Dagger.Thunk(()->10)))
end

#=
@testset "darray distributed refcount" begin
    D2 = remotecall_fetch(2, compute(Distribute(Blocks(10, 20), rand(40,40)))) do D
        D2 = D
    end
    @test size(collect(D2)) == (40,40)
    GC.gc()
    @test size(collect(D2)) == (40,40)
end
=#

@testset "sharedarray" begin
    A = SharedArray{Int}((1024,))
    B = SharedArray{Int}((1024,))
    C = Dagger.merge_sorted(Base.Order.Forward, A, B)
    @test length(C) === length(A) + length(B)
    @test typeof(C) === (Dagger.use_shared_array[] ? SharedArray{Int,1} : Array{Int,1})
end


text:
using Dates: DateTime, now, UTC
using CStructures: default_value
using Base: CFunction

export timespec_now, main_loop

"""
    noop(args...)

Dummy function - should never by actually called.
If it is called, it behaves as if the function was not supported.
"""
noop(args...) = UV_ENOTSUP
dummy = @cfunction $noop Cint ()

"""
    filter_ops(m::Module)

Return a vector of C-callable function pointers for functions with names
as defined in `FuseApi.FuseLowlevelOps`.

For functions, which are not exported by `module`, a null pointer is generated
to indicate to `fuselib3`, that this function is not supported.
"""
function filter_ops(mod::Module, fs, cbhandles)
    n = length(fieldnames(FuseLowlevelOps))
    resize!(cbhandles, n)
    fill!(cbhandles, dummy)
    res = fill(C_NULL, n)
    for (i, f) in enumerate(fieldnames(FuseLowlevelOps))
        if (cf = getp(mod, f, noop)) != noop
            cb = gen_callback(f, cf, fs)
            cbhandles[i] = cb
            res[i] = Base.unsafe_convert(Ptr{Cvoid}, Base.cconvert(Ptr{Cvoid}, cb))
        end
    end
    FuseLowlevelOps(res...)
end

"""
    getp(m::Module, f::Symbol, default::Function)

Return the function named `f`, if defined and exported in module `m`,
otherwise the default value.
"""
function getp(m, f::Symbol, default::Function)
    try
        if f in names(m) && isdefined(m, f)
            getproperty(m, f)
        else
            default
        end
    catch
        default
    end
end

"""
    gen_callback(name::Symbol, f::Function, fs::Any)

Call the function named `Symbol('G', name)` (for example `Gread`), which
generates a C-callable closure of type `CFunction`, dedicated to call back the Julia
function `f(fs, ...)` (for example `read(fs, req, ...)`).

This closure serves as a GC-handle for the function pointer, which is actually passed
as a callback to be called from a C-environment.
"""
function gen_callback(name::Symbol, c::Function, fs::Any)
    gf = Symbol('G', name)
    gfunction = getproperty(FuseApi, gf)
    gfunction(c, fs)
end

"""
    docall(f, req)

Provides common frame for all req-related C-callback functions
If the actual callback function throws or returns a value != 0
An error reply is returned.
"""
function docall(f::Function, req::FuseReq)
    error = Base.UV_ENOTSUP
    try
        if req.pointer === C_NULL
            ccall(:jl_breakpoint, Cvoid, (Any,), f)
            throw(ArgumentError("docall($f, $req) - req was zero!!!"))
        end
        error = f()
    catch
        error = Base.UV_EACCES
        rethrow()
    finally
        if error != 0 && req.pointer !== C_NULL
            fuse_reply_err(req, abs(error))
        end
        nothing
    end
    nothing
end

"""
    docall(f)

Frame for `init` and `destroy`
"""
function docall(f::Function)
    try
        f()
    finally
        nothing
    end
    nothing
end

"""
    create_args(cmd::String, args::AbstractVector{String})

Return CStructGuarded{FuseCmdlineArgs} object to pass arguments to C.
"""
function create_args(cmd::String, arg::AbstractVector{String})
    data = Cserialize(FuseCmdlineArgs, (argc = length(arg) + 1, argv = [cmd, arg..., nothing]))
    CStructGuarded{FuseCmdlineArgs}(data)
end

"""
    main_loop(cmdline_args::Vector{String}, fs::Module)

Set up and start a Fuse session. Use commandline arguments `cmdline_args`
and the filesystem implementation defined in module `fs`.

Module `fs` must define and export all supported callback functions, the names
of which are specified in `FuseLowlevelOps`.
"""
function main_loop(args::AbstractVector{String}, mod::Module, user_data=nothing)
    user_data_r = Ref(user_data)
    cbhandles = CFunction[]
    # it is essential, that user_data and cbhandles are protected from GC
    # during the active time of _main_loop.
    GC.@preserve user_data_r cbhandles _main_loop(args, mod, user_data, cbhandles)
end
function _main_loop(args, mod, user_data, cbhandles)
        fargs = create_args("command", args)
        opts = fuse_parse_cmdline(fargs)
        callbacks = filter_ops(mod, user_data, cbhandles)
        se = fuse_session_new(fargs, callbacks)
        se == C_NULL && throw(ArgumentError("fuse_session_new failed $rc"))
        rc = fuse_session_mount(se, opts.mountpoint)
        rc != 0 && throw(ArgumentError("fuse_session_mount failed $rc"))

        cfg = FuseLoopConfig(opts.clone_fd, opts.max_idle_threads)
        rc = opts.singlethread ? fuse_session_loop(se) : fuse_session_loop_mt(se, cfg)
        rc != 0 && throw(ArgumentError("fuse_session_loop failed $rc"))
        fuse_session_unmount(se)
        fuse_session_destroy(se)
end

"""
    fuse_parse_cmdline(args::CStructAccess{FuseCmdlineArgs})

Return FuseCmdlineOpts object with selected options and mountpoint from the cmdline.
"""
function fuse_parse_cmdline(args::CStructAccess{FuseCmdlineArgs})
    opts = Cserialize(LFuseCmdlineOpts, ())
    popts = pointer(opts)
    ccall((:fuse_parse_cmdline, :libfuse3), Cint, (Ptr{FuseCmdlineArgs}, Ptr{UInt8}), args, popts)
    default_value(FuseCmdlineOpts, CStruct{LFuseCmdlineOpts}(opts))
end

function log(args...)
    flush(stderr)
    println(stderr, args...)
end

# Base.unsafe_convert(::Type{Ptr{FuseReq}}, cs::FuseReq) where T = Ptr{FuseReq}(cs.pointer)
function Base.unsafe_convert(::Type{Ptr{T}}, s::SubArray{T,1}) where T
    Ptr{T}(pointer(s.parent)) + first(s.indices[1]) - 1
end

Base.convert(::Type{Timespec}, t::CStruct{Timespec}) = Timespec(t.seconds, t.nanoseconds)

"""
    timespec_now()

Current time in `Timespec` format (seconds and nanoseconds)
"""
function timespec_now()
    a = (now(UTC) - DateTime(1970, 1, 1)).value
    s, n = fldmod(a, 10^3)
    n *= 10^6
    Timespec(s, n)
end


text:
# ---------------------------------------------------------------
# Read all the NCA5 files and reformat them to a file
# thats usable throughout the report. All heights in mm
#
# For each scenario, the following files are made:
#   NCA5_scn_gmsl.nc: Global-mean sea level for each scenario
#   NCA5_scn_grid.nc: Local sea level on a 1° grid
#   NCA5_scn_tg.nc: Local sea level at PSMSL tide-gauge locations
# This script only has to be run once.
# ---------------------------------------------------------------
module ConvertNCA5ToGrid

using NetCDF
using DelimitedFiles
using Statistics
using NCDatasets
using Interpolations
using DelimitedFiles

function RunConvertNCA5ToGrid(settings)
    println("\nProcessing NCA5 scenarios...")
    # ConvertGMSL(settings)
    for scn ∈ 1:length(settings["NCA5_scenarios"])
        println(" "*settings["NCA5_scenarios"][scn]*"...")
        ConvertGMSL(scn,settings)
        ConvertRegional(scn,settings)
        println(" "*settings["NCA5_scenarios"][scn]*" done")

    end
    println("Processing NCA5 scenarios done\n")
    return nothing
end

function ConvertGMSL(scn,settings)
    # -------------------------------------------
    # This reads the Sep 14 data
    # GMSL is in different file
    # Read the truncated file with GMSL values
    # Save as netcdf
    # -------------------------------------------
    println("  Processing GMSL...")
    scn_gmsl   = Dict()
    scenario_target = [30,50,100,150,200]

    # Read file info
    fn = settings["dir_NCA5_raw"] * "gmsl/gmsl"*lpad(scenario_target[scn],3,"0")*"/AIS_gmsl"*lpad(scenario_target[scn],3,"0")*".nc"
    quantiles_raw = ncread(fn,"quantiles")
    idx_low = findfirst(quantiles_raw.>0.1699)
    idx_med = findfirst(quantiles_raw.>0.4999)
    idx_high = findfirst(quantiles_raw.>0.8299)
    scn_gmsl["quantiles"] = quantiles_raw[[idx_low,idx_med,idx_high]]
    scn_gmsl["years"] = ncread(fn,"years")

    # Read individual fields
    for process ∈ settings["processes"]
        if process != "verticallandmotion"
            println("   Reading "*process)
            fn = settings["dir_NCA5_raw"] * "gmsl/gmsl"*lpad(scenario_target[scn],3,"0")*"/"*process*"_gmsl"*lpad(scenario_target[scn],3,"0")*".nc"
            if (process == "landwaterstorage")
                scn_gmsl[process] = (ncread(fn,"sea_level_change",start=[1,2,1],count=[1,-1,-1])[:,:,[idx_low,idx_med,idx_high]])[1,:,:];
            else
                scn_gmsl[process]  = (ncread(fn,"sea_level_change",start=[1,1,1],count=[1,-1,-1])[:,:,[idx_low,idx_med,idx_high]])[1,:,:];
            end
        end
    end

    fn_gmsl_out = settings["dir_NCA5"] * "NCA5_"*settings["NCA5_scenarios"][scn]*"_gmsl.nc"
    fh = Dataset(fn_gmsl_out,"c")
    defDim(fh,"years", length(scn_gmsl["years"]))
    defDim(fh,"percentiles",3)
    defVar(fh,"years",trunc.(Int16,scn_gmsl["years"]),("years",),deflatelevel=6)
    defVar(fh,"percentiles",trunc.(Int16,[17,50,83]),("percentiles",),deflatelevel=6)
    for process ∈ settings["processes"]
        if process != "verticallandmotion"
            defVar(fh,process,trunc.(Int16,scn_gmsl[process]),("years","percentiles",),deflatelevel=6)
        end
    end
    close(fh)
    println("  Processing GMSL done")
end

# function ConvertGMSL(settings)
#     # -------------------------------------------
#     # GMSL is in different file
#     # Read the truncated file with GMSL values
#     # Save as netcdf
#     # Currently, I don't have the GMSL components
#     # -------------------------------------------
#     println("  Processing GMSL...")
#     fn_in = settings["dir_NCA5_raw"]*"NCA5_GMSL.csv"
#     data_raw = readdlm(fn_in,',',skipstart=1)
#     traw = [2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100, 2110, 2120,2130,2140, 2150, 2200, 2250, 2300]
#     gmsl = Dict()
#     gmsl["years"] = traw
#     for (idx,scn) ∈ enumerate(settings["NCA5_scenarios"])
#         gmsl[scn] = 10 .* convert.(Int,data_raw[3*(idx-1)+1:3*idx,7:end]')
#         @. gmsl[scn] = gmsl[scn][:,[2,1,3]]
#     end

#     # Save as netcdf
#     for (idx,scn) ∈ enumerate(settings["NCA5_scenarios"])
#         fn_gmsl_out = settings["dir_NCA5"] * "NCA5_"*scn*"_gmsl.nc"
#         fh = Dataset(fn_gmsl_out,"c")
#         defDim(fh,"years", length(gmsl["years"]))
#         defDim(fh,"percentiles",3)
#         defVar(fh,"years",trunc.(Int16,gmsl["years"]),("years",),deflatelevel=6)
#         defVar(fh,"percentiles",trunc.(Int16,[17,50,83]),("percentiles",),deflatelevel=6)
#         defVar(fh,"total",trunc.(Int16,gmsl[scn]),("years","percentiles",),deflatelevel=6)
#         close(fh)
#     end
#     println("  Processing GMSL done")
# end

function ConvertRegional(scn,settings)
    # -------------------------------------------
    # Read all processes for individual scenarios
    # Save lat/lon grid and tide-gauge values in
    # individual files
    # -------------------------------------------
    println("  Processing regional sea level...")
    scn_grid = Dict()
    scn_tg   = Dict()
    scenario_target = [30,50,100,150,200]

    # Read file info
    fn = settings["dir_NCA5_raw"] * "regional/gmsl"*lpad(scenario_target[1],3,"0")*"/AIS_gmsl"*lpad(scenario_target[1],3,"0")*".nc"
    quantiles_raw = ncread(fn,"quantiles")
    idx_low = findfirst(quantiles_raw.>0.1699)
    idx_med = findfirst(quantiles_raw.>0.4999)
    idx_high = findfirst(quantiles_raw.>0.8299)
    scn_grid["quantiles"] = quantiles_raw[[idx_low,idx_med,idx_high]]
    scn_tg["quantiles"] = quantiles_raw[[idx_low,idx_med,idx_high]]
    scn_grid["years"] = ncread(fn,"years")
    scn_tg["years"] = ncread(fn,"years")
    acc_loc_grid = ncread(fn,"locations") .> 4000
    acc_loc_tg = ncread(fn,"locations") .< 4000
    latmat = reshape(ncread(fn,"lat")[acc_loc_grid],(360,181))
    lonmat = reshape(ncread(fn,"lon")[acc_loc_grid],(360,181))
    @. lonmat = mod(lonmat,360)
    scn_grid["lon"] = lonmat[:,1]
    scn_grid["lat"] = reverse(latmat[1,:])

    scn_tg["loc"] = ncread(fn,"locations")[acc_loc_tg]
    scn_tg["lon"] = mod.(ncread(fn,"lon")[acc_loc_tg],360)
    scn_tg["lat"] = ncread(fn,"lat")[acc_loc_tg]

    # Read individual fields
    for process ∈ settings["processes"]
        println("   Reading "*process)
        fn = settings["dir_NCA5_raw"] * "regional/gmsl"*lpad(scenario_target[scn],3,"0")*"/"*process*"_gmsl"*lpad(scenario_target[scn],3,"0")*".nc"
        if (process == "landwaterstorage") || (process == "verticallandmotion")
            scn_grid[process] = reshape(ncread(fn,"sea_level_change",start=[findfirst(acc_loc_grid),2,1],count=[-1,-1,-1])[:,:,[idx_low,idx_med,idx_high]],(360,181,length(scn_grid["years"]),3));
            scn_tg[process] = ncread(fn,"sea_level_change",start=[1,2,1],count=[sum(acc_loc_tg),-1,-1])[:,:,[idx_low,idx_med,idx_high]];
        else
            scn_grid[process] = reshape(ncread(fn,"sea_level_change",start=[findfirst(acc_loc_grid),1,1],count=[-1,-1,-1])[:,:,[idx_low,idx_med,idx_high]],(360,181,length(scn_grid["years"]),3));
            scn_tg[process]  = ncread(fn,"sea_level_change",start=[1,1,1],count=[sum(acc_loc_tg),-1,-1])[:,:,[idx_low,idx_med,idx_high]];
        end
        scn_grid[process] = convert.(Float32,scn_grid[process]);
        scn_tg[process] = convert.(Float32,scn_tg[process]);
        reverse!(scn_grid[process],dims=2);
    end

    # Estimate climate-driven sea-level changes (ice, sterodynamic, land water storage, but not VLM)
    # Assume errors are independent, so they're added in quadrature
    println("   Summing climate-related fields...")

    scn_grid["total_climate"] = zeros(Float32,size(scn_grid["total"]))
    scn_tg["total_climate"] = zeros(Float32,size(scn_tg["total"]))
    for process ∈ settings["processes_clim"]
        @. scn_grid["total_climate"][:,:,:,2] += scn_grid[process][:,:,:,2]
        @. scn_grid["total_climate"][:,:,:,1] += (abs((scn_grid[process][:,:,:,1] -  scn_grid[process][:,:,:,2]))^2)
        @. scn_grid["total_climate"][:,:,:,3] += (abs((scn_grid[process][:,:,:,3] -  scn_grid[process][:,:,:,2]))^2)
        @. scn_tg["total_climate"][:,:,2] += scn_tg[process][:,:,2]
        @. scn_tg["total_climate"][:,:,1] += (abs((scn_tg[process][:,:,1] -  scn_tg[process][:,:,2]))^2)
        @. scn_tg["total_climate"][:,:,3] += (abs((scn_tg[process][:,:,3] -  scn_tg[process][:,:,2]))^2)
    end
    @. scn_grid["total_climate"][:,:,:,1] = scn_grid["total_climate"][:,:,:,2] - sqrt(@. scn_grid["total_climate"][:,:,:,1])
    @. scn_grid["total_climate"][:,:,:,3] = scn_grid["total_climate"][:,:,:,2] + sqrt(@. scn_grid["total_climate"][:,:,:,3])
    @. scn_tg["total_climate"][:,:,1] = scn_tg["total_climate"][:,:,2] - sqrt(@. scn_tg["total_climate"][:,:,1])
    @. scn_tg["total_climate"][:,:,3] = scn_tg["total_climate"][:,:,2] + sqrt(@. scn_tg["total_climate"][:,:,3])

    # Fix NaNs
    @. scn_grid["total_climate"][scn_grid["total_climate"]<-32000] = -32000
    @. scn_grid["total_climate"][scn_grid["total_climate"]>32000] = -32000
    @. scn_tg["total_climate"][scn_tg["total_climate"]<-32000] = -32000
    @. scn_tg["total_climate"][scn_tg["total_climate"]>32000] = -32000
    for process ∈ settings["processes"]
        @. scn_grid[process][scn_grid[process]<-32000] = -32000
        @. scn_grid[process][scn_grid[process]>32000] = -32000
        @. scn_tg[process][scn_tg[process]<-32000] = -32000
        @. scn_tg[process][scn_tg[process]>32000] = -32000
    end

    # Save scenario file
    println("   Saving")
    fn_grid_out = settings["dir_NCA5"] * "NCA5_"*settings["NCA5_scenarios"][scn]*"_grid.nc"
    fn_tg_out   = settings["dir_NCA5"] * "NCA5_"*settings["NCA5_scenarios"][scn]*"_tg.nc"

    fh = Dataset(fn_grid_out,"c")
    defDim(fh,"years", length(scn_grid["years"]))
    defDim(fh,"percentiles",3)
    defDim(fh,"lon", length(scn_grid["lon"]))
    defDim(fh,"lat", length(scn_grid["lat"]))
    defVar(fh,"lon",scn_grid["lon"],("lon",),deflatelevel=6)
    defVar(fh,"lat",scn_grid["lat"],("lat",),deflatelevel=6)
    defVar(fh,"years",trunc.(Int16,scn_grid["years"]),("years",),deflatelevel=6)
    defVar(fh,"percentiles",trunc.(Int16,[17,50,83]),("percentiles",),deflatelevel=6)
    for process ∈ settings["processes"]
        defVar(fh,process,trunc.(Int16,scn_grid[process]),("lon","lat","years","percentiles",),deflatelevel=6)
    end
    defVar(fh,"total_climate",trunc.(Int16,scn_grid["total_climate"]),("lon","lat","years","percentiles",),deflatelevel=6)
    close(fh)

    fh = Dataset(fn_tg_out,"c")
    defDim(fh,"years", length(scn_tg["years"]))
    defDim(fh,"percentiles",3)
    defDim(fh,"loc", length(scn_tg["loc"]))
    defVar(fh,"loc",scn_tg["loc"],("loc",),deflatelevel=6)
    defVar(fh,"lon",scn_tg["lon"],("loc",),deflatelevel=6)
    defVar(fh,"lat",scn_tg["lat"],("loc",),deflatelevel=6)
    defVar(fh,"years",trunc.(Int16,scn_tg["years"]),("years",),deflatelevel=6)
    defVar(fh,"percentiles",trunc.(Int16,[17,50,83]),("percentiles",),deflatelevel=6)
    for process ∈ settings["processes"]
        defVar(fh,process,trunc.(Int16,scn_tg[process]),("loc","years","percentiles",),deflatelevel=6)
    end
    defVar(fh,"total_climate",trunc.(Int16,scn_tg["total_climate"]),("loc","years","percentiles",),deflatelevel=6)
    close(fh)
    println("  Processing regional sea level done")
    return nothing
end

end


text:
#
# Plotting routines.
#

"""
$(TYPEDSIGNATURES)

Heuristic check if Plotter is PyPlot
"""
ispyplot(Plotter)=typeof(Plotter)==Module&& isdefined(Plotter,:Gcf)

"""
$(TYPEDSIGNATURES)

Heuristic check if Plotter is Plots
"""
isplots(Plotter)=typeof(Plotter)==Module&& isdefined(Plotter,:gr)

# Plot color scale for grid colors.
function frgb(Plotter,i,max;pastel=false)
    x=Float64(i-1)/Float64(max)
    if (x<0.5)
        r=1.0-2.0*x
        g=2.0*x
        b=0.0
    else
        r=0.0
        g=2.0-2.0*x
        b=2.0*x-1.0
    end
    if pastel
        r=0.5+0.5*r
        g=0.5+0.5*g
        b=0.5+0.5*b
    end
    if ispyplot(Plotter)
        return (r,g,b)
    end
    if isplots(Plotter)
        return Plotter.RGB(r,g,b)
    end
end

"""
$(SIGNATURES)

Find the circumcenter of a triangle.

Created from C source of Jonathan R Shewchuk <jrs@cs.cmu.edu>

Modified to return absolute coordinates.
"""
function tricircumcenter!(circumcenter,a,b,c)

    # Use coordinates relative to point `a' of the triangle.
    xba = b[1] - a[1]
    yba = b[2] - a[2]
    xca = c[1] - a[1]
    yca = c[2] - a[2]

    # Squares of lengths of the edges incident to `a'.
    balength = xba * xba + yba * yba
    calength = xca * xca + yca * yca

    # Calculate the denominator of the formulae.
    # if EXACT
    #    Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html
    #    to ensure a correctly signed (and reasonably accurate) result,
    #    avoiding any possibility of division by zero.
    #  denominator = 0.5 / orient2d((double*) b, (double*) c, (double*) a)

    # Take your chances with floating-point roundoff
    denominator = 0.5 / (xba * yca - yba * xca)

    # Calculate offset (from `a') of circumcenter.
    xcirca = (yca * balength - yba * calength) * denominator
    ycirca = (xba * calength - xca * balength) * denominator

# The result is returned both in terms of x-y coordinates and xi-eta
# coordinates, relative to the triangle's point `a' (that is, `a' is
# the origin of both coordinate systems).	 Hence, the x-y coordinates
# returned are NOT absolute; one must add the coordinates of `a' to
# find the absolute coordinates of the circumcircle.  However, this means
# that the result is frequently more accurate than would be possible if
# absolute coordinates were returned, due to limited floating-point
# precision.  In general, the circumradius can be computed much more
# accurately.

    circumcenter[1] = xcirca+a[1]
    circumcenter[2] = ycirca+a[2]

    return circumcenter
end

"""
$(TYPEDSIGNATURES)

The plot function is not exported, but kept for backward compatibility.
"""
plot(Plotter, tio; kwargs...)=plot_triangulateio(Plotter, tio; kwargs...)

"""
$(TYPEDSIGNATURES)

Plot contents of triangulateio structure.
The plot module (currently only PyPlot is possible, GLMakie may follow)
is passed as parameter.
This allows to keep the package free of heavy
plot package dependencies.
"""
function plot_triangulateio(Plotter,
                            tio::TriangulateIO;
                            voronoi=nothing,
                            aspect=1,
                            circumcircles=false
                            )

    if ispyplot(Plotter)
        PyPlot=Plotter
        ax = PyPlot.matplotlib.pyplot.gca()
        ax.set_aspect(aspect)
        if numberofpoints(tio)==0
            return
        end
        x=tio.pointlist[1,:]
        y=tio.pointlist[2,:]
        xminmax=extrema(x)
        yminmax=extrema(y)

        dx=xminmax[2]-xminmax[1]
        ax.set_xlim(xminmax[1]-dx/10,xminmax[2]+dx/10)
        dy=yminmax[2]-yminmax[1]
        ax.set_ylim(yminmax[1]-dy/10,yminmax[2]+dy/10)

        PyPlot.scatter(x,y, s=10,color="b")
        t=transpose(tio.trianglelist.-1)
        if size(tio.triangleattributelist,2)>0
            PyPlot.tripcolor(x,y,t,tio.triangleattributelist[1,:],cmap="Pastel2")
        end
        if numberoftriangles(tio)>0
            PyPlot.triplot(x,y,t,color="k")
            if circumcircles
                t=0:0.025*π:2π
                circle(x,y,r)=PyPlot.plot(x.+r.*cos.(t),y.+r.*sin.(t),color=(0.4,0.05,0.4), linewidth=0.3)
                cc=zeros(2)
                for itri=1:numberoftriangles(tio)

                    trinodes=tio.trianglelist[:,itri]
                    a=tio.pointlist[:,trinodes[1]]
                    b=tio.pointlist[:,trinodes[2]]
                    c=tio.pointlist[:,trinodes[3]]
                    tricircumcenter!(cc,a,b,c)
                    r=sqrt((cc[1]-a[1])^2+(cc[2]-a[2])^2)
                    circle(cc[1],cc[2],r)
                    PyPlot.scatter([cc[1]],[cc[2]], s=20,color=(0.4,0.05,0.4))
                end
            end
        end
        if numberofsegments(tio)>0
            lines=Any[]
            rgb=Any[]
            markermax=maximum(tio.segmentmarkerlist)
            # see https://gist.github.com/gizmaa/7214002
            for i=1:numberofsegments(tio)
                x1=tio.pointlist[1,tio.segmentlist[1,i]]
                y1=tio.pointlist[2,tio.segmentlist[1,i]]
                x2=tio.pointlist[1,tio.segmentlist[2,i]]
                y2=tio.pointlist[2,tio.segmentlist[2,i]]
                push!(lines,collect(zip([x1,x2],[y1,y2])))
                push!(rgb,frgb(PyPlot,tio.segmentmarkerlist[i],markermax))
            end
            ax.add_collection(PyPlot.matplotlib.collections.LineCollection(lines,colors=rgb,linewidth=3))
        end

        if numberoftriangles(tio)==0 && numberofregions(tio)>0
            x=tio.regionlist[1,:]
            y=tio.regionlist[2,:]
            r=tio.regionlist[3,:]
            PyPlot.scatter(x,y, s=20,c=r,cmap="Dark2")
        end

        if numberoftriangles(tio)==0 && numberofholes(tio)>0
            x=tio.holelist[1,:]
            y=tio.holelist[2,:]
            PyPlot.scatter(x,y, s=20,color="k")
        end

        if voronoi!=nothing && numberofedges(voronoi)>0
            bcx=sum(x)/length(x)
            bcy=sum(y)/length(y)
            wx=maximum(abs.(x.-bcx))
            wy=maximum(abs.(y.-bcy))
            ww=max(wx,wy)

            x=voronoi.pointlist[1,:]
            y=voronoi.pointlist[2,:]
            PyPlot.scatter(x,y, s=10,color="g")
            for i=1:numberofedges(voronoi)
                i1=voronoi.edgelist[1,i]
                i2=voronoi.edgelist[2,i]
                if i1>0 && i2>0
                    PyPlot.plot([voronoi.pointlist[1,i1],voronoi.pointlist[1,i2]],
                                 [voronoi.pointlist[2,i1],voronoi.pointlist[2,i2]],
                                 color="g")
                else
                    x0=voronoi.pointlist[1,i1]
                    y0=voronoi.pointlist[2,i1]
                    xn=voronoi.normlist[1,i]
                    yn=voronoi.normlist[2,i]
                    normscale=1.0e10
                    if x0+normscale*xn>bcx+ww
                        normscale=abs((bcx+ww-x0)/xn)
                    end
                    if x0+normscale*xn<bcx-ww
                        normscale=abs((bcx-ww-x0)/xn)
                    end
                    if y0+normscale*yn>bcy+ww
                        normscale=abs((bcy+ww-y0)/yn)
                    end
                    if y0+normscale*yn<bcy-ww
                        normscale=abs((bcy-ww-y0)/yn)
                    end
                    PyPlot.plot([x0, x0+normscale*xn], [y0,y0+normscale*yn],color="g")
                end

            end
        end
    end

    ### Currently not working
    if isplots(Plotter)
        #=
        if p==nothing
            p=Plotter.plot()
        end
        for icell=1:num_cells(grid)
            rgb=frgb(Plotter,grid.cellregions[icell],num_cellregions(grid),pastel=true)
            coord1=nodecoord(grid,cellnode(grid,1,icell))
            coord2=nodecoord(grid,cellnode(grid,2,icell))
            coord3=nodecoord(grid,cellnode(grid,3,icell))
            # https://github.com/JuliaPlotter/Plotter.jl/issues/605
            tri=Plotter.Shape([coord1[1],coord2[1], coord3[1]],[coord1[2],coord2[2],coord3[2]])
            Plotter.plot!(p,tri,color=rgb,label="")
        end
        for icell=1:num_cells(grid)
                coord1=nodecoord(grid,cellnode(grid,1,icell))
            coord2=nodecoord(grid,cellnode(grid,2,icell))
            coord3=nodecoord(grid,cellnode(grid,3,icell))
            Plotter.plot!(p, [coord1[1],coord2[1]],[coord1[2],coord2[2]]  ,linewidth=0.5,color=:black,label="")
            Plotter.plot!(p, [coord1[1],coord3[1]],[coord1[2],coord3[2]]  ,linewidth=0.5,color=:black,label="")
            Plotter.plot!(p, [coord2[1],coord3[1]],[coord2[2],coord3[2]]  ,linewidth=0.5,color=:black,label="")
        end
        for ibface=1:num_bfaces(grid)
            rgb=frgb(Plotter,grid.bfaceregions[ibface],num_bfaceregions(grid))
            coord1=nodecoord(grid,bfacenode(grid,1,ibface))
            coord2=nodecoord(grid,bfacenode(grid,2,ibface))
            Plotter.plot!(p,[coord1[1],coord2[1]],[coord1[2],coord2[2]]  ,linewidth=5,color=rgb,label="")
        end
        if show
            Plotter.gui(p)
        end
        return p
        =#
    end
end

"""
$(TYPEDSIGNATURES)

Plot  pair of triangulateio structures arranged
in two subplots. This is intendd for visualizing both input
and output data.
"""
function plot_in_out(Plotter, triin, triout;voronoi=nothing, circumcircles=false,title="")
    if ispyplot(Plotter)
        PyPlot=Plotter
        PyPlot.clf()
        PyPlot.subplots(nrows=1, ncols=2)
        PyPlot.suptitle(title)
        PyPlot.subplot(121)
        PyPlot.title("In")
        plot_triangulateio(PyPlot,triin)
        PyPlot.subplot(122)
        PyPlot.title("Out")
        plot_triangulateio(PyPlot,triout,voronoi=voronoi,circumcircles=circumcircles)
        PyPlot.tight_layout()
    end
end


text:
using TOML
using Configurations

the_julia_language = """

    Julia Programming Language

"""

head = """

    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=edge">
    <meta name="viewport" content="width=device-width, initial-scale=1">

    <title>Roger's packages</title>

    <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-giJF6kkoqNQ00vy+HMDP7azOuL0xtbfIcaT9wjKHr8RbDVddVHyTfAAsrekwKmP1" crossorigin="anonymous">
    <link rel="stylesheet" href="assets/main.css">
    <link rel="shortcut icon" type="image/x-icon" href="/images/favicon.jpeg">

"""

intro = """

    rogerluo.dev

        this page is made for convenience when someone (including myself)
        wants to check the status of packages I created/maintain.

        Most of the packages are written in $the_julia_language.

        My Julia software projects are usually organized as many small
        packages so that people can easily use different components
        for other purpose.

        The documentation of each project are maintained
        as a curated page under the project package. You can also find the
        link to them below.

"""

theme_script = """

"""

@option struct Package
    name::String
    description::String
    doc::Bool = true # has doc
    codecov::Bool = true # has codecov
    codecov_token::Maybe{String} = nothing
    gh_action::Bool = true # has github action
    user::String = "Roger-luo" # github user
end

@option struct Project
    name::String
    description::String
    user::Maybe{String} = nothing
    page::Maybe{String} = nothing
    packages::Maybe{Vector{Package}} = nothing
end

@option struct Info
    project::Vector{Project}
end

repo_link(user, name) = "https://github.com/$(user)/$(name).jl"
page_link(user, name) = "https://$(user).github.io/$(name).jl"

package_badges(pkg::Package) = package_badges(pkg.user, pkg.name, pkg.doc, pkg.gh_action, pkg.codecov, pkg.codecov_token)

function package_badges(user, name, doc::Bool, gh_action::Bool, codecov::Bool, codecov_token)
    repo = repo_link(user, name)
    page = page_link(user, name)
    stable_link = "$page/stable"
    dev_link = "$page/dev"
    gha_link = "$repo/actions"
    codecov_link = "https://codecov.io/gh/$(user)/$(name).jl"

    stable_img = "https://img.shields.io/badge/docs-stable-blue.svg"
    dev_img = "https://img.shields.io/badge/docs-dev-blue.svg"
    gha_img = "$repo/workflows/CI/badge.svg"
    codecov_img = "https://codecov.io/gh/$(user)/$(name).jl/branch/master/graph/badge.svg"
    if !isnothing(codecov_token)
        codecov_img *= "?token=g0FlmB4YLj"
    end

    badge(img, link, alt) = """

    """

    return """
    $(doc ? badge(stable_img, stable_link, "doc-stable") : "")
    $(doc ? badge(dev_img, dev_link, "doc-dev") : "")
    $(gh_action ? badge(gha_img, gha_link, "github-action") : "")
    $(codecov ? badge(codecov_img, codecov_link, "codecov") : "")
    """
end

function html(pkg::Package)
    """

                $(pkg.name)

$(package_badges(pkg))

$(pkg.description)

    """
end

function html(project::Project)
    name = "
$(project.name)
"
    description = """

        $(project.description)

    """

    page = isnothing(project.page) ? "" : """

        $(project.page)

    """

    if isnothing(project.packages)
        user = isnothing(project.user) ? "Roger-luo" : project.user
        """

            $name
            $page

                $(package_badges(user, project.name, true, true, true, nothing))

            $description

        """
    else
        packages = join([html(each) for each in project.packages], "\n")
        """

            $name
            $page
            $description

                $packages

        """
    end
end

function html(info::Info)
    join(map(html, info.project), "\n")
end

info = from_toml(Info, "packages.toml")::Info

index = """
<!DOCTYPE html>

$head

        $intro
        $(html(info))

        罗秀哲 Xiu-Zhe(Roger) Luo - me@rogerluo.dev  - References on request

    $theme_script

"""

if ispath("page")
    rm("page"; force=true, recursive=true)
else
    mkpath("page")
    mkpath(joinpath("page", "assets"))
end

cp("assets", joinpath("page", "assets"); force=true)

open(joinpath("page", "index.html"), "w+") do io
    write(io, index)
end

max_stars_repo_name: Roger-luo/packages

max_issues_repo_name: Roger-luo/packages

max_forks_repo_name: Roger-luo/packages