text:
# ......................................... FITS public sector .................................................................

# .................................................... fits_copy ...................................................

"""
    fits_copy(filenameA [, filenameB="" [; protect=true]])

Copy "filenameA" to "filenameB" (with mandatory ".fits" extension)
Key:
* `protect::Bool`: overwrite protection
#### Examples:
```
fits_copy("T01.fits")
  'T01.fits' was saved as 'T01 - Copy.fits'

fits_copy("T01.fits", "T01a.fits")
  FitsError: 'T01a.fits' in use (set ';protect=false' to lift overwrite protection)

fits_copy("T01.fits", "T01a.fits"; protect=false)
  'T01.fits' was saved as 'T01a.fits'
```
"""
function fits_copy(filenameA::String, filenameB::String=" "; protect=true)

    o = _fits_read_IO(filenameA)
    f = cast_FITS_name(filenameA)

    filenameB = filenameB == " " ? "$(f.name) - Copy.fits" : filenameB

    _validate_FITS_name(filenameB)
    _isavailable(filenameB, protect) || error("FitsError: '$filenameB' in use (set ';protect=false' to lift overwrite protection)")
    _fits_write_IO(o,filenameB)

    return println("'$filenameA' was saved as '$filenameB'")

end

# .................................................... fits_copy ...................................................

"""
    fits_combine(strFirst, strLast [; protect=true])

Copy "filenameFirst" to "filenameLast" (with mandatory ".fits" extension)

Key:
* `protect::Bool`: overwrite protection
#### Example:
```
fits_combine("T01.fits", "T22.fits")
  'T01-T22.fits': file created
```
"""
function fits_combine(filnamFirst::String, filnamLast::String; protect=true)

    Base.Filesystem.isfile(filnamFirst) || error("FitsError: '$filnamFirst': file not found")
    Base.Filesystem.isfile(filnamLast) ||  error("FitsError: '$filnamLast': file not found")

    filnamFirst = uppercase(filnamFirst)
    filnamLast = uppercase(filnamLast)

    nam = cast_FITS_name(filnamFirst)
    strPre = nam.prefix
    strNum = nam.numerator
    strExt = nam.extension
    valNum = parse(Int,strNum )
    numLeadingZeros = length(strNum) - length(string(valNum))

    nam2 = cast_FITS_name(filnamLast)
    strPre2 = nam2.prefix
    strNum2 = nam2.numerator
    strExt2 = nam2.extension
    valNum2 = parse(Int,strNum2 )
    numLeadingZeros2 = length(strNum2) - length(string(valNum2))

    if strPre ≠ strPre2
        error(strPre * " ≠ " * strPre2 * " (prefixes must be identical)")
    elseif strExt ≠ strExt2
        error(strExt * " ≠ " * strExt2 * " (file extensions must be identical)")
    elseif uppercase(strExt) ≠ ".FITS"
        error("file extension must be '.fits'")
    end

    numFiles = 1 + valNum2 - valNum
    f = fits_read(filnamFirst)
    dataFirst = f[1].dataobject.data  # read an image from disk
    t = typeof(f[1].dataobject.data[1,1,1])
    s = size(f[1].dataobject.data)

    dataStack =  Array{t,3}(undef, s[1], s[2] , numFiles)

    itr = valNum:valNum2
    filnamNext = filnamFirst
    for i ∈ itr
        l = length(filnamNext)
        filnamNext = strPre * "0"^numLeadingZeros * string(i) * ".fits"
        if l < length(filnamNext)
            numLeadingZeros = numLeadingZeros -1
            filnamNext = strPre * "0"^numLeadingZeros * string(i) * ".fits"
        end
        f = fits_read(filnamNext)
        dataNext = f[1].dataobject.data                # read an image from disk
        dataStack[:, :,i] = dataNext[:, :,1]
    end

    filnamOut = strPre * strNum * "-" * strPre * strNum2 * strExt

    fits_create(filnamOut, dataStack; protect)

    return println("'$filnamOut': file created")

end

# .................................................... fits_create ...................................................

"""
    fits_create(filename [, data [; protect=true]])

Create FITS file of given filename [, optional data block [, default overwrite protection]] and return Array of HDUs.
Key:
* `protect::Bool`: overwrite protection
#### Examples:
```
strExample = "minimal.fits"
fits_create(strExample;protect=false)

f = fits_read(strExample)
a = f[1].dataobject.data
b = f[1].header.keys
println(a);println(b)
  Any[]
  ["SIMPLE", "NAXIS", "EXTEND", "COMMENT", "END"]

strExample = "remove.fits"
data = [11,21,31,12,22,23,13,23,33]
data = reshape(data,(3,3,1))
fits_create(strExample, data; protect=false)

f = fits_read(strExample)
fits_info(f[1])

  File: remove.fits
  hdu: 1
  hdutype: PRIMARY
  DataType: Int64
  Datasize: (3, 3, 1)

  Metainformation:
  SIMPLE  =                    T / file does conform to FITS standard
  BITPIX  =                   64 / number of bits per data pixel
  NAXIS   =                    3 / number of data axes
  NAXIS1  =                    3 / length of data axis 1
  NAXIS2  =                    3 / length of data axis 2
  NAXIS3  =                    1 / length of data axis 3
  BZERO   =                  0.0 / offset data range to that of unsigned integer
  BSCALE  =                  1.0 / default scaling factor
  EXTEND  =                    T / FITS dataset may contain extensions
  COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg
  END

  3×3×1 Array{Int64, 3}:
  [:, :, 1] =
   11  12  13
   21  22  23
   31  23  33
```
"""
function fits_create(filename::String, data=[]; protect=true)

    strErr = "FitsError: '$filename': creation failed (filename in use - set ';protect=false' to overrule overwrite protection)"

    _validate_FITS_name(filename)
    _isavailable(filename, protect) || error(strErr)

    nhdu = 1
    hdutype = "PRIMARY"

       FITS_data = [_cast_data(i, hdutype, data) for i=1:nhdu]
    FITS_headers = [_cast_header(_PRIMARY_input(FITS_data[i]), i) for i=1:nhdu]

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return _fits_save(FITS)

end
# test ...
function fits_create()

    strExample = "minimal.fits"
    fits_create(strExample, protect=false)

    f = fits_read(strExample)
    a = f[1].header.keys[1]  == "SIMPLE"
    b = f[1].dataobject.data == Any[]
    c = get(Dict(f[1].header.dict),"SIMPLE",0)
    d = get(Dict(f[1].header.dict),"NAXIS",0) == 0;

    rm(strExample)

    test = .![a, b, c, d];

    return !convert(Bool,sum(test))

end

# .................................................... fits_info ...................................................
"""
    fits_info(hdu)

Print metafinformation and data of given `FITS_HDU`
#### Example:
```
strExample = "remove.fits"
data = [11,21,31,12,22,23,13,23,33]
data = reshape(data,(3,3,1))
fits_create(strExample, data; protect=false)

f = fits_read(strExample)
fits_info(f[1])

  File: remove.fits
  hdu: 1
  hdutype: PRIMARY
  DataType: Int64
  Datasize: (3, 3, 1)

  Metainformation:
  SIMPLE  =                    T / file does conform to FITS standard
  BITPIX  =                   64 / number of bits per data pixel
  NAXIS   =                    3 / number of data axes
  NAXIS1  =                    3 / length of data axis 1
  NAXIS2  =                    3 / length of data axis 2
  NAXIS3  =                    1 / length of data axis 3
  BZERO   =                  0.0 / offset data range to that of unsigned integer
  BSCALE  =                  1.0 / default scaling factor
  EXTEND  =                    T / FITS dataset may contain extensions
  COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg
  END

  3×3×1 Array{Int64, 3}:
  [:, :, 1] =
   11  12  13
   21  22  23
   31  23  33

```
"""
function fits_info(hdu::FITS_HDU)

    typeof(hdu) <: FITS_HDU || error("FitsWarning: FITS_HDU not found")

    info = [
        "\r\nFile: " * hdu.filename,
        "hdu: " * Base.string(hdu.hduindex),
        "hdutype: " * hdu.dataobject.hdutype,
        "DataType: " * Base.string(Base.eltype(hdu.dataobject.data)),
        "Datasize: " * Base.string(Base.size(hdu.dataobject.data)),
        "\r\nMetainformation:"
        ]

    records = hdu.header.records

    Base.append!(info, records)

    println(Base.join(info .* "\r\n"))

    return hdu.dataobject.data

end

# .................................................... fits_read ...................................................
"""
    fits_read(filename)

Read FITS file and return Array of `FITS_HDU`s
#### Example:
```
strExample = "minimal.fits"
fits_create(strExample;protect=false)

f = fits_read(strExample)
f[1].dataobject.data
  Any[]

rm(strExample); f = nothing
```
"""
function fits_read(filename::String)

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return FITS

end
# test ...
function fits_read()

    strExample = "minimal.fits"
    fits_create(strExample, protect=false)

    f = fits_read(strExample)
    a = f[1].header.keys[1]  == "SIMPLE"
    b = f[1].dataobject.data == Any[]
    c = get(Dict(f[1].header.dict),"SIMPLE",0)
    d = get(Dict(f[1].header.dict),"NAXIS",0) == 0;

    rm(strExample)

    test = .![a, b, c, d];

    return !convert(Bool,sum(test))

end

# .................................................... fits_extend ...................................................
"""
    fits_extend(filename, data_extend, hdutype="IMAGE")

Extend the FITS file of given filename with the data of `hdutype` from `data_extend`  and return Array of HDUs.
#### Examples:
```
strExample = "test_example.fits"
data = [0x0000043e, 0x0000040c, 0x0000041f]
fits_create(strExample, data, protect=false)

f = fits_read(strExample)
a = Float16[1.01E-6,2.0E-6,3.0E-6,4.0E-6,5.0E-6]
b = [0x0000043e, 0x0000040c, 0x0000041f, 0x0000042e, 0x0000042f]
c = [1.23,2.12,3.,4.,5.]
d = ['a','b','c','d','e']
e = ["a","bb","ccc","dddd","ABCeeaeeEEEEEEEEEEEE"]
data = [a,b,c,d,e]
fits_extend(strExample, data, "TABLE")

f = fits_read(strExample)
f[2].dataobject.data
  5-element Vector{String}:
   "1.0e-6 1086 1.23 a a                    "
   "2.0e-6 1036 2.12 b bb                   "
   "3.0e-6 1055 3.0  c ccc                  "
   "4.0e-6 1070 4.0  d dddd                 "
   "5.0e-6 1071 5.0  e ABCeeaeeEEEEEEEEEEEE "

rm(strExample); f = data = a = b = c = d = e = nothing
```
"""
function fits_extend(filename::String, data_extend, hdutype="IMAGE")

    hdutype == "IMAGE"    ? (records, data) = _IMAGE_input(data_extend)    :
    hdutype == "TABLE"    ? (records, data) = _TABLE_input(data_extend)    :
    hdutype == "BINTABLE" ? (records, data) = _BINTABLE_input(data_extend) : error("FitsError: unknown HDU type")

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    nhdu = nhdu + 1

    push!(FITS_headers, _cast_header(records, nhdu))              # update FITS_header object
    push!(FITS_data, _cast_data(nhdu, hdutype, data))             # update FITS_data object

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    _fits_save(FITS)

    return FITS

end
# test ...
function fits_extend()

    strExample = "test_example.fits"
    data = [0x0000043e, 0x0000040c, 0x0000041f]
    fits_create(strExample, data, protect=false)

    f = fits_read(strExample)
    a = Float16[1.01E-6,2.0E-6,3.0E-6,4.0E-6,5.0E-6]
    b = [0x0000043e, 0x0000040c, 0x0000041f, 0x0000042e, 0x0000042f]
    c = [1.23,2.12,3.,4.,5.]
    d = ['a','b','c','d','e']
    e = ["a","bb","ccc","dddd","ABCeeaeeEEEEEEEEEEEE"]
    data = [a,b,c,d,e]
    fits_extend(strExample, data, "TABLE")

    f = fits_read(strExample)
    a = f[1].header.keys[1]  == "SIMPLE"
    b = f[1].dataobject.data[1] == 0x0000043e
    c = f[2].header.keys[1]  == "XTENSION"
    d = f[2].dataobject.data[1] == "1.0e-6 1086 1.23 a a                    "
    e = get(Dict(f[2].header.dict),"NAXIS",0) == 2

    rm(strExample)

    test = .![a, b, c, d, e]

    return !convert(Bool,sum(test))

end

# .................................................... fits_add_key ...................................................

"""
    fits_add_key(filename, hduindex, key, value, comment)

Add a header record of given 'key, value and comment' to 'HDU[hduindex]' of file with name 'filename'
#### Example:
```
strExample="minimal.fits"
fits_create(strExample;protect=false)
fits_add_key(strExample, 1, "KEYNEW1", true, "FITS dataset may contain extension")

f = fits_read(strExample)
fits_info(f[1])

  File: minimal.fits
  hdu: 1
  hdutype: PRIMARY
  DataType: Any
  Datasize: (0,)

  Metainformation:
  SIMPLE  =                    T / file does conform to FITS standard
  NAXIS   =                    0 / number of data axes
  EXTEND  =                    T / FITS dataset may contain extensions
  COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg
  KEYNEW1 =                    T / FITS dataset may contain extension
  END

  Any[]
```
"""
function fits_add_key(filename::String, hduindex::Int, key::String, val::Any, com::String)

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    key = _format_key(key)

    h = FITS_headers[hduindex]
    Base.get(h.maps, key, 0) > 0 && return println("FitsError: '$key': key in use (use different name or edit key)")

    newrecords = _fits_new_records(key, val, com)

    Base.pop!(h.records)
   [Base.push!(h.records, newrecords[i]) for i ∈ eachindex(newrecords)]
    Base.push!(h.records, "END" * Base.repeat(" ",77))

    FITS_headers[hduindex] = _cast_header(h.records, hduindex)

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return _fits_save(FITS)

end
# test ...
function fits_add_key()

    strExample="minimal.fits"
    fits_create(strExample;protect=false)
    fits_add_key(strExample, 1, "KEYNEW1", true, "FITS dataset may contain extension")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW1",0)
    r = f[1].header.records;

    test = r[i] == "KEYNEW1 =                    T / FITS dataset may contain extension             "

    rm(strExample)

    return test

end

"""
    fits_edit_key(filename, hduindex, key, value, comment)

Edit a header record of given 'key, value and comment' to 'HDU[hduindex]' of file with name 'filename'
#### Example:
```
data = DateTime("2020-01-01T00:00:00.000")
strExample="minimal.fits"
fits_create(strExample;protect=false)
fits_add_key(strExample, 1, "KEYNEW1", true, "this is record 5")
fits_edit_key(strExample, 1, "KEYNEW1", data, "record 5 changed to a DateTime type")

f = fits_read(strExample)
fits_info(f[1])

  File: minimal.fits
  hdu: 1
  hdutype: PRIMARY
  DataType: Any
  Datasize: (0,)

  Metainformation:
  SIMPLE  =                    T / file does conform to FITS standard
  NAXIS   =                    0 / number of data axes
  EXTEND  =                    T / FITS dataset may contain extensions
  COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg
  KEYNEW1 = '2020-01-01T00:00:00' / record 5 changed to a DateTime type
  END

  Any[]
```
"""
function fits_edit_key(filename::String, hduindex::Int, key::String, val::Any, com::String)

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    key = _format_key(key)
    res = ["SIMPLE","BITPIX","NAXIS","NAXIS1","NAXIS2","NAXIS3","BZERO","END"]
    key ∈ res && return println("FitsError: '$key': cannot be edited (key protected under FITS standard)")

    h = FITS_headers[hduindex]
    i = Base.get(h.maps, key, 0)
    i == 0 && return println("FitsError: '$key': key not found")

    nold = length(h.keys)
    nobs = length(_fits_obsolete_records(h,i))
    newrecords = _fits_new_records(key, val, com)
    oldrecords = h.records[i+nobs:end]
   [Base.pop!(h.records) for j=i:nold]
   [Base.push!(h.records, newrecords[i]) for i ∈ eachindex(newrecords)]
   [Base.push!(h.records, oldrecords[i]) for i ∈ eachindex(oldrecords)]

    FITS_headers[hduindex] = _cast_header(h.records, hduindex)

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return _fits_save(FITS)

end
# test ...
function fits_edit_key()

    strExample="minimal.fits"
    fits_create(strExample;protect=false)
    fits_add_key(strExample, 1, "KEYNEW1", true, "FITS dataset may contain extension")
    fits_edit_key(strExample, 1, "KEYNEW1", false, "comment has changed")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW1",0)
    r = f[1].header.records;

    test = r[i] == "KEYNEW1 =                    F / comment has changed                            "

    rm(strExample)

    return test

end

"""
    fits_delete_key(filename, hduindex, key)

Delete a header record of given `key`, `value` and `comment` to `FITS_HDU[hduindex]` of file with name  'filename'
#### Examples:
```
strExample="minimal.fits"
fits_create(strExample;protect=false)
fits_add_key(strExample, 1, "KEYNEW1", true, "this is record 5")

f = fits_read(strExample)
get(f[1].header.maps,"KEYNEW1",0)
  5

fits_delete_key(strExample, 1, "KEYNEW1")

f = fits_read(strExample)
get(f[1].header.maps,"KEYNEW1",0)
  0

fits_delete_key(filnam, 1, "NAXIS")
 'NAXIS': cannot be deleted (key protected under FITS standard)
```
"""
function fits_delete_key(filename::String, hduindex::Int, key::String)

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    key = _format_key(key)
    res = ["SIMPLE","BITPIX","NAXIS","NAXIS1","NAXIS2","NAXIS3","BZERO","END"]
    key ∈ res && return println("FitsError: '$key': cannot be edited (key protected under FITS standard)")

    h = FITS_headers[hduindex]
    i = Base.get(h.maps, key, 0)
    i == 0 && return println("FitsError: '$key': key not found")

    nold = length(h.keys)
    nobs = length(_fits_obsolete_records(h,i))
    oldrecords = h.records[i+nobs:end]
   [Base.pop!(h.records) for j=i:nold]
   [Base.push!(h.records, oldrecords[i]) for i ∈ eachindex(oldrecords)]

    FITS_headers[hduindex] = _cast_header(h.records, hduindex)

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return _fits_save(FITS)

end
# test ...
function fits_delete_key()

    strExample="minimal.fits"
    fits_create(strExample;protect=false)
    fits_add_key(strExample, 1, "KEYNEW1", true, "FITS dataset may contain extension")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW1",0)

    test1 = i == 5

    fits_delete_key(strExample, 1, "KEYNEW1")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW1",0)

    test2 = i == 0

    test = .![test1, test2]

    rm(strExample)

    return !convert(Bool,sum(test))

end

"""
    fits_rename_key(filename, hduindex, keyold, kewnew)

Rename the key of a header record of file with name 'filename'
#### Example:
```
strExample="minimal.fits"
fits_create(strExample;protect=false)
fits_add_key(strExample, 1, "KEYNEW1", true, "this is record 5")
fits_rename_key(strExample, 1, "KEYNEW1",  "KEYNEW2")

f = fits_read(strExample)
fits_info(f[1])

  File: minimal.fits
  hdu: 1
  hdutype: PRIMARY
  DataType: Any
  Datasize: (0,)

  Metainformation:
  SIMPLE  =                    T / file does conform to FITS standard
  NAXIS   =                    0 / number of data axes
  EXTEND  =                    T / FITS dataset may contain extensions
  COMMENT    Primary FITS HDU    / http://fits.gsfc.nasa.gov/iaufwg
  KEYNEW2 =                    T / this is record 5
  END

  Any[]
```
"""
function fits_rename_key(filename::String, hduindex::Int, keyold::String, keynew::String)

    o = _fits_read_IO(filename)

    nhdu = _hdu_count(o)

    FITS_headers = [_read_header(o,i) for i=1:nhdu]
       FITS_data = [_read_data(o,i) for i=1:nhdu]

    keyold = _format_key(keyold)
       res = ["SIMPLE","BITPIX","NAXIS","NAXIS1","NAXIS2","NAXIS3","BZERO","END"]
    keyold ∈ res && return println("FitsWarning: '$keyold': cannot be renamed (key protected under FITS standard)")

    h = FITS_headers[hduindex]
    i = Base.get(h.maps, keyold, 0)
    i == 0 && return println("FitsError: '$keyold': key not found")
    Base.get(h.maps, keynew, 0) > 0 && return println("FitsWarning: '$keynew': key in use (use different name or edit key)")

    h.records[i] = rpad(keynew,8) * h.records[i][9:80]

    FITS_headers[hduindex] = _cast_header(h.records, hduindex)

    FITS = [FITS_HDU(filename, i, FITS_headers[i], FITS_data[i]) for i=1:nhdu]

    return _fits_save(FITS)

end
# test ...
function fits_rename_key()

    strExample="minimal.fits"
    fits_create(strExample;protect=false)
    fits_add_key(strExample, 1, "KEYNEW1", true, "this is record 5")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW1",0)

    test1 = i == 5

    fits_rename_key(strExample, 1,"KEYNEW1", "KEYNEW2")

    f = fits_read(strExample)
    i = get(f[1].header.maps,"KEYNEW2",0)

    test2 = i == 5

    test = .![test1, test2]

    rm(strExample)

    return !convert(Bool,sum(test))

end

# ....................... parse FITS_TABLE into a Vector of its columns .........................................

"""
    parse_FITS_TABLE(hdu)

Parse `FITS_TABLE` (ASCII table) into a Vector of its columns for further processing by the user.
Default formatting in ISO 2004 FORTRAN data format specified by keys "TFORMS1" - "TFORMSn")
Display formatting in ISO 2004 FORTRAN data format ("TDISP1" - "TDISPn") prepared for user editing.
#### Example:
```
strExample = "example.fits"
data = [10, 20, 30]
fits_create(strExample, data; protect=false)

t1 = Float16[1.01E-6,2.0E-6,3.0E-6,4.0E-6,5.0E-6]
t2 = [0x0000043e, 0x0000040c, 0x0000041f, 0x0000042e, 0x0000042f]
t3 = [1.23,2.12,3.,4.,5.]
t4 = ['a','b','c','d','e']
t5 = ["a","bb","ccc","dddd","ABCeeaeeEEEEEEEEEEEE"]
data = [t1,t2,t3,t4,t5]
fits_extend(strExample, data, "TABLE")

f = fits_read(strExample)
d = f[2].header.dict
d = [get(d,"TFORM\$i",0) for i=1:5]; println(strip.(d))
  SubString{String}["'E6.1    '", "'I4      '", "'F4.2    '", "'A1      '", "'A20     '"]

f[2].dataobject.data                            # this is the table hdu
  5-element Vector{String}:
   "1.0e-6 1086 1.23 a a                    "
   "2.0e-6 1036 2.12 b bb                   "
   "3.0e-6 1055 3.0  c ccc                  "
   "4.0e-6 1070 4.0  d dddd                 "
   "5.0e-6 1071 5.0  e ABCeeaeeEEEEEEEEEEEE "

parse_FITS_TABLE(f[2])
  5-element Vector{Vector{T} where T}:
   [1.0e-6, 2.0e-6, 3.0e-6, 4.0e-6, 5.0e-6]
   [1086, 1036, 1055, 1070, 1071]
   [1.23, 2.12, 3.0, 4.0, 5.0]
   ["a", "b", "c", "d", "e"]
   ["a                   ", "bb                  ", "ccc                 ", "dddd                ", "ABCeeaeeEEEEEEEEEEEE"]
```
"""
function parse_FITS_TABLE(hdu::FITS_HDU)

    dict = hdu.header.dict
    thdu = Base.strip(Base.get(dict,"XTENSION", "UNKNOWN") ,['\'',' '])

    thdu == "TABLE" || return error("Error: $thdu is not an ASCII TABLE HDU")

    ncols = Base.get(dict,"TFIELDS", 0)
    nrows = Base.get(dict,"NAXIS2", 0)
    tbcol = [Base.get(dict,"TBCOL$n", 0) for n=1:ncols]
    tform = [Base.get(dict,"TFORM$n", 0) for n=1:ncols]
    ttype = [cast_FORTRAN_format(tform[n]).Type for n=1:ncols]
    tchar = [cast_FORTRAN_format(tform[n]).TypeChar for n=1:ncols]
    width = [cast_FORTRAN_format(tform[n]).width for n=1:ncols]
      itr = [(tbcol[k]:tbcol[k]+width[k]-1) for k=1:ncols]

     data = hdu.dataobject.data
     data = [[data[i][itr[k]] for i=1:nrows] for k=1:ncols]
     data = [tchar[k] == 'D' ? Base.join.(Base.replace!.(Base.collect.(data[k]), 'D'=>'E')) : data[k] for k=1:ncols]
     Type = [ttype[k] == "Aw" ? (width[k] == 1 ? Char : String) : ttype[k] == "Iw" ? Int : Float64 for k=1:ncols]
     data = [ttype[k] == "Aw" ? data[k] : parse.(Type[k],(data[k])) for k=1:ncols]

    return data

end


text:
using QUIP
using Test, LinearAlgebra, BenchmarkTools
using JuLIP, JuLIP.Testing
using QUIP: soap

@testset "QUIP.jl" begin
    include("test.jl")
end

ext: jl


text:
# 32 bit unicode character
a = 'a';
println(typeof(a));
println(sizeof(a));

# Char <--> int conversions
println(Int('a'));
println(Char(97));

# Character comparisons

println('A' <= 'a');
println('a' > 'b');
println('a' - 'd');
println('A' + 2);
println(2 + 'A');

# Strings

str = "hello world, this is a string\n";
println(typeof(str));

# Extracting characters via index from 1 to n
println(str[begin]);
println(str[1]);
println(str[5]);
println(str[end]);

# Get first and last indexes of string
println(firstindex(str));
println(lastindex(str));

# Length of string
println(length(str));

# length and lastindex might not always match
# - because some chars might take up 2 code units

# arithmetic indexing
str[1 + 4];

# range indexing with inclusive bounds
str[4:9];

# single indexing vs range indexing

str[1:1] # returns string
str[1] # returns character

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
# Vortex detection and imprinting

"""
    P, Q = poles(u)

Return the Wirtinger derivatives of u, assuming a unit grid step
"""
function poles(u)
    N = size(u, 1)
    @assert N == size(u, 2)
    u = complex.(u)
    rs = (-1:1) .+ 1im * (-1:1)'
    rs /= sum(abs2.(rs))
    conv(u, rs) = [rs .* u[j:j+2, k:k+2] |> sum for j = 1:N-2, k = 1:N-2]
    P = zero(u)
    P[2:end-1, 2:end-1] .= conv(u, conj(rs))
    Q = zero(u)
    Q[2:end-1, 2:end-1] .= conv(u, rs)
    P, Q
end

"""
    unroll(θs)

Return θs adjusted mod 2π to be as continuous as possible
"""
function unroll(θ)
    Θ = similar(θ)
    Θ[1] = θ[1]
    poff = 0.0
    for j = 2:length(θ)
        jump = θ[j] - θ[j-1]
        if jump > π
            poff -= 2π
        elseif jump < -π
            poff += 2π
        end
        Θ[j] = θ[j] + poff
    end
    Θ
end

"""
    find_vortex(d, u)

Return the centre of the sole vortex in u
"""
function find_vortex(d::Discretisation, u, ixs = nothing)
    P, Q = poles(u)
    if isnothing(ixs)
        ixs = abs.(P) .> 0.5maximum(abs.(P))
    end
    regress_core(d, u, ixs)
end

find_vortex(u) = find_vortex(default(:discretisation), u)

function find_vortices(d::Discretisation, u)
    P, Q = poles(u)
    ixs = abs.(P) .> 0.5maximum(abs.(P))
    clusters = cluster_adjacent(adjacent_index, keys(ixs)[ixs])
    [find_vortex(d, u, [j ∈ C for j in keys(u)]) for C in clusters]
end

function find_moat(d::Discretisation, u)
    P, Q = poles(u)
    v = @. (R-w/2 < r < R+w/2)*abs(P+conj(Q))/abs(u)
    ixs = v .> 0.5maximum(v)
    regress_core(d, u, ixs)
end

function regress_core(d, u, ixs)
    z = argand(d)
    ixs = expand_box(box(ixs))[:]
    a, b, c = [z[ixs] conj(z[ixs]) ones(size(z[ixs]))] \ u[ixs]
    (b * conj(c) - conj(a) * c) / (abs2(a) - abs2(b))
end

function expand_box(bb::CartesianIndices{2})
    e1, e2 = ((1, 0), (0, 1)) .|> CartesianIndex
    # TODO raise issue on CartesianIndex half-pregnancy
    if length(bb) == 1
        bb = bb[]
        bb = [
            bb-e1-e2 bb-e1 bb-e1+e2
            bb-e2 bb bb+e2
            bb+e1-e2 bb+e1 bb+e1+e2
        ]
    elseif size(bb, 1) == 1
        bb = [bb .- e1; bb; bb .+ e1]
    elseif size(bb, 2) == 1
        bb = [bb .- e2 bb bb .+ e2]
    end
    bb
end

"Bounding box of cartesian indices"
function box(cixs::AbstractVector{CartesianIndex{2}})
    j1 = k1 = typemax(Int)
    j2 = k2 = typemin(Int)
    for c in cixs
        j1 = min(j1, c[1])
        j2 = max(j2, c[1])
        k1 = min(k1, c[2])
        k2 = max(k2, c[2])
    end
    CartesianIndices((j1:j2, k1:k2))
end

box(ixs::AbstractMatrix{Bool}) = box(keys(ixs)[ixs])

function cluster_adjacent(f, ixs)
    # function f determines adjacency
    clusters = Set()
    for i in ixs
        out = Set()
        ins = Set()
        for C in clusters
            if any(j -> f(i, j), C)
                push!(ins, C)
            else
                push!(out, C)
            end
        end
        push!(out, union(Set([i]), ins...))
        clusters = out
    end
    clusters
end

adjacent_index(j, k) = -1 ≤ j[1] - k[1] ≤ 1 && -1 ≤ j[2] - k[2] ≤ 1

"""
    PinnedVortices([d], rvs, as) <: Optim.Manifold

Constrain a field to have vortices centred at the rvs

TODO document the new point-estimates and Gaussians approach
"""
struct PinnedVortices <: Manifold
    # V'*ψ[:] is a vector of interpolated ψ.(rvs)
    # (1-U*V')*ψ projects the singularities back to the rvs
    # TODO make these Tensars
    # TODO solve a linear system in the u(rvs) to make the projection exact
    U::Matrix{Complex{Float64}}
    V::Matrix{Complex{Float64}}

    function PinnedVortices(
        d::Sampling{2},
        rvs::Vector{Complex{Float64}},
        as::Vector{Float64}
    )
        U = Array{Complex{Float64}}(undef, d.n^2, length(rvs))
        V = similar(U)
        for j in eachindex(rvs)
            U[:, j] = finterp(d, rvs[j], as[j])[:]
            V[:, j] = interpolant(d, rvs[j])[:]
        end
        new(U, V)
    end
end

"""
    finterp(d, r, a)

Return interpolant(d, r) convoluted with a Gaussian of width a.

Subtracting this is a stable way to shift a vortex by further than one pixel.
"""
function finterp end

# use coordinates instead of this

# PinnedVortices(d::Discretisation, rvs::Vararg{Number}; kwargs...) =
#     PinnedVortices(d, [convert(Complex{Float64}, r) for r in rvs]...; kwargs...)
# PinnedVortices(rvs::Vararg{Number}; kwargs...) =
#     PinnedVortices(default(:discretisation), rvs...; kwargs...)

function prjct!(M, q)
    w = q[:]
    q .-= reshape(M.U * (M.V' * w), size(q))
end

# The "vortex at R" space is invariant under normalisation
Optim.retract!(M::PinnedVortices, q) = Optim.retract!(Sphere(), prjct!(M, q))
Optim.project_tangent!(M::PinnedVortices, dq, q) =
    Optim.project_tangent!(Sphere(), prjct!(M, dq), q)


text:
function maxadjacent(M, pattern)
	H, W = size(M)
	h, w = size(pattern)

	bestprod = 0

	for i = 1:H-h+1
		for j = 1:W-w+1
			m = M[i:i+h-1,j:j+w-1] .* pattern
			println(m)
			p = prod(m[pattern.==1])
			if p > bestprod
				bestprod = p
			end
		end
	end

	return bestprod
end

function main()
	M = [08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08;
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00;
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65;
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91;
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80;
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50;
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70;
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21;
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72;
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95;
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92;
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57;
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58;
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40;
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66;
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69;
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36;
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16;
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54;
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48]

	N = 4

	patterns = {}
	push!(patterns, ones(eltype(M), N, 1))
	push!(patterns, ones(eltype(M), 1, N))
	push!(patterns, eye(eltype(M), N))
	push!(patterns, rotr90(eye(eltype(M), N)))

	bestprod = 0
	for pattern in patterns
		n = maxadjacent(M, pattern)
		if n > bestprod
			bestprod = n
		end
	end

	println(bestprod)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
# multisession support for Caesar.jl

function multisessionquery(conn, session::T, multisessions::Vector{T}) where {T <: AbstractString}
  len = length(multisessions)
  loadtx = transaction(conn)
  # construct the query
  query = "match (m:$(session):LANDMARK) "*
          "with m.label as mlb "*
          "order by mlb asc "*
          "match (n:LANDMARK) "*
          "where n.label = mlb "*
          "and (";
  for i in 1:len
    ots = multisessions[i]
    query = query*"n:$(ots) "
    query = i < len ? query*"or " : query
  end
  query = query*") "*
                "with distinct n.label as lbl, collect(id(n)) as others, collect(labels(n)) as nlb "*
                "return lbl, others, nlb";

  loadtx(query, submit=true)
end

# multisession must have session label starting with SESS* in each query response line
function parsemultisessionqueryresult!(lm2others::Dict{Symbol, Dict{Symbol, Int}}, cph)
  i = 0
  # @show cph.results[1]["data"]
  for i in 1:length(cph.results[1]["data"])
    res = cph.results[1]["data"][i]["row"]
    for j in 1:length(res[2])
      sess = ""
      for lb in res[3][j]
        if lb[1:4] == "SESS"
          sess = Symbol(lb)
        end
      end
      if !haskey(lm2others, sess)  lm2others[sess] = Dict{Symbol, Int}() end
      lm2others[sess][Symbol(res[1])] = res[2][j]
    end
  end
  nothing
end

"""
    getLandmOtherSessNeoIDs{T <: AbstractString}(::CloudGraph, session::T="", multisessions=Vector{T}())

Return dict of dict of Neo4j vertex IDs by session and landmark symbols.
"""
function getLandmOtherSessNeoIDs(cg::CloudGraph;
      session::T="",
      multisessions::Vector{T}=String[]  ) where {T <: AbstractString}
  #
  lm2others = Dict{Symbol, Dict{Symbol, Int}}()

  len = length(multisessions)
  len > 0 ? nothing : (return lm2others)

  if length(multisessions)==0
    cph = multisessionquery(cg.neo4j.connection, session, multisessions)
    parsemultisessionqueryresult!(lm2others, cph)
  else
    info("Ignoring multisession")
  end

  return lm2others
end

"""
    getprp2kde(::CloudGraph, neoids::Vector{Int}; N::Int=100)

Return PriorPoint2DensityNH with N points based on beliefs of neoids, and equal
share null hypothesis between length(neoids)+1 beliefs.
"""
function getprpt2kde(cloudGraph::CloudGraph, neoids::Vector{Int}; N::Int=100 )
  #
  # P = BallTreeDensity[]
  PTS = Array{Float64,2}()
  for fielem in neoids # lm2others[lm]
    cv = CloudGraphs.get_vertex(cloudGraph, fielem)
    vert = cloudVertex2ExVertex(cv)
    pts = getVal(vert)
    PTS = size(PTS,2) == 0 ? pts' : [PTS; pts']
    # push!(P, kde!(pts))
  end
  PTS = PTS';

  l = collect(1:size(PTS,2))
  shuffle!(l)

  # TODO -- not up sampling if N is bigger than size(PTS)
  ptsN = PTS[:,l[1:N]];
  @show size(ptsN)

  nh = 1.0/(length(neoids)+1.0)
  ph = Float64[nh;1.0-nh]
  return PriorPoint2DensityNH(kde!(ptsN), ph)
end

"""
    getAllLandmarkNeoIDs(::Dict{Symbol, Dict{Symbol, Int}}, ::Symbol)

Return Vector{Int} of Neo4j vertex IDs relating to symbol, as listed in lm2others.
"""
function getAllLandmarkNeoIDs(lm2others::Dict{Symbol, Dict{Symbol,Int}}, slm::Symbol)
  mslmneoids = Int[]
  for (sess, vals) in lm2others
    if haskey(vals, slm)
      push!(mslmneoids, vals[slm])
    end
  end
  return mslmneoids
end

"""
    getLocalSubGraphMultisession{T <: AbstractString}(cg::CloudGraph, lm2others; session::T="", numneighbors::Int=0)

Return subgraph copy of type FactorGraph contaning values from session in lm2others, and Vector{Symbol} of primary
key symbols used for graph exstraction.
"""
function getLocalSubGraphMultisession(cg::CloudGraph,
            lm2others;
            session::T="",
            numneighbors::Int=0  ) where {T <: AbstractString}
  #
  res = Dict{Symbol, Int}()
  sfg = Caesar.initfg(sessionname=session, cloudgraph=cg)
  if length(lm2others) > 0
    for (sess,ms) in lm2others
      for (sym, neoid) in ms
        res[sym] = 0
      end
    end
    getVertNeoIDs!(cg, res, session=session)
    fullcurrneolist = collect(values(res))
    fetchsubgraph!(sfg, fullcurrneolist, numneighbors=numneighbors) # can set numneighbors=0
  end
  return sfg, collect(keys(res))
end

"""
    findExistingMSConstraints(fgl::FactorGraph)

Return Dict{Symbol, Int} of vertex symbol to Neo4j node ID of MULTISESSION constraints in this `fgl.sessionname`.
"""
function findExistingMSConstraints(fgl::FactorGraph)
  loadtx = transaction(fgl.cg.neo4j.connection)
  query =  "match (n:$(fgl.sessionname):MULTISESSION)
            return id(n), n.exVertexId, n.label"
  cph = loadtx(query, submit=true)

  # parse the response into a dictionary
  existingms = Dict{Symbol, Int}() # symlbl => neoid
  for res in cph.results[1]["data"]
    neoid, exvid, lbl = res["row"][1], res["row"][2], res["row"][3]
    existingms[Symbol(lbl)] = neoid
  end

  return existingms
end

"""
    removeReinsertMultisessionPrior!{T <: FunctorSingleton}(fgl::FactorGraph, exims::Dict{Symbol, Int}, prp2::T, sym::Symbol, uid::Int)

Return new multisession factor, symbol sym, inserted using prp2::PriorPoint2DensityNH using user ExVertex id (uid),
after checking existing multisession dict (exims::(exvsym=>neoid)) if any nodes should be removed before inserting new ones
on that graph node. Assuming just one multisession prior per node.
"""
function removeReinsertMultisessionPrior!{T <: FunctorSingleton}(fgl::FactorGraph,
        exims::Dict{Symbol, Int},
        prp2::T,
        sym::Symbol,
        uid::Int  )
  #
  # remove previous Multi session constraint on this sym vertex
  if haskey(exims, sym)
    info("Removing previous multisession neoid=$(exims[sym]) on sym=$(sym).")
    removeNeo4jID(fgl.cg, neoid=exims[sym])
    delete!(exims, sym)
  end
  # add new multisession constraint to the graph
  addFactor!(fgl, [sym], prp2, labels=["MULTISESSION"], uid=uid)
end

"""
    rmInstMultisessionPriors!(::CloudGraph; session<:AbstractString=, multisessions::Vector{<:AbstractString}= )
"""
function rmInstMultisessionPriors!(cloudGraph::CloudGraph;
      session::T="NA",
      multisessions::Vector{T}=String[]  ) where {T <: AbstractString}
  #
  session!="NA" ? nothing : error("Please specify a valid session, currently = $(session)")

  multisessionsl = Vector{String}(setdiff(multisessions, [session]))
  length(multisessionsl) > 0 ? nothing : (return nothing)
  # get the landmarks of interest from neighboring sessions
  lm2others = getLandmOtherSessNeoIDs(cloudGraph,
                  session=session,multisessions=multisessionsl)
  #

  # grab local subgraph using NeoIDs
  sfg, lms = getLocalSubGraphMultisession(cloudGraph, lm2others,
                  session=session, numneighbors=1)
  #

  # get dict(sym => neoid) of exiting multisession constraints which may be removed during this update
  exims = findExistingMSConstraints(sfg)
  # get highest factor exvid
  mfn4jid = getmaxfactorid(cloudGraph.neo4j.connection, session)

  println("Multisession constraints in $(session), from $(multisessionsl), on $(lms)")
  for sym in lms
    neoids = getAllLandmarkNeoIDs(lm2others, sym)
    prp2 = getprpt2kde(cloudGraph, neoids, N=100 )
    mfn4jid += 1
    factorms = removeReinsertMultisessionPrior!(sfg, exims, prp2, sym, mfn4jid)
  end
  nothing
end

function removeMultisessions!(cloudGraph::CloudGraph; session::AbstractString="NA")
  session!="NA" ? nothing : error("Please specify a valid session, currently = $(session)")

  loadtx = transaction(cloudGraph.neo4j.connection)
  query =  "match (n:$(session):MULTISESSION)
            detach delete n
            return count(n)"
  cph = loadtx(query, submit=true)
  commit(loadtx)
  return cph
end

#


text:
using DashBootstrapComponents, DashHtmlComponents

toast = html_div([
    dbc_button(
        "Open toast",
        id = "auto-toast-toggle",
        color = "primary",
        className = "mb-3",
        n_clicks = 0,
    ),
    dbc_toast(
        [html_p("This is the content of the toast", className = "mb-0")],
        id = "auto-toast",
        header = "This is the header",
        icon = "primary",
        duration = 4000,
        is_open = false,
    ),
]);

callback!(app, Output("auto-toast", "is_open"), Input("auto-toast-toggle", "n_clicks")) do n
    if n > 0
        return true
    end
    return no_update()
end;

max_stars_repo_name: DrGFreeman/dash-bootstrap-components

max_issues_repo_name: DrGFreeman/dash-bootstrap-components

max_forks_repo_name: DrGFreeman/dash-bootstrap-components


text:
"""
    +   +---+---+
    | 1   2   3 |
    +---+   +   +
    | 4   5 | 6
    +---+---+---+

    julia> const graph = [
            0 1 0 0 0 0;
            1 0 1 0 1 0;
            0 1 0 0 0 1;
            0 0 0 0 1 0;
            0 1 0 1 0 0;
            0 0 1 0 0 0]

    julia> dist, path = dijkstra(graph, 1)
    (Dict(4=>3,2=>1,3=>2,5=>2,6=>3,1=>0), Dict(4=>5,2=>1,3=>2,5=>2,6=>3,1=>0))

    julia> printpath(path, 6) # Display solution of the maze
    1 -> 2 -> 3 -> 6

"""
function dijkstra(graph, source::Int=1)
    # ensure that the adjacency matrix is squared
    @assert size(graph, 1) == size(graph, 2)
    inf = typemax(Int64)
    n   = size(graph, 1)

    Q    = IntSet(1:n)                  # Set of unvisited nodes
    dist = Dict(n => inf for n in Q)    # Unknown distance function from source to v
    prev = Dict(n => 0   for n in Q)    # Previous node in optimal path from source
    dist[source] = 0                    # Distance from source to source

    function _minimumdist(nodes) # Find the less distant node among nodes
        kmin, vmin = nothing, inf
        for (k, v) in dist
            if k ∈ nodes && v ≤ vmin
                kmin, vmin = k, v
            end
        end
        return kmin
    end
    # Until all nodes are visited...
    while !isempty(Q)
        u = _minimumdist(Q)         # Vertex in Q with smallest dist[]
        pop!(Q, u)
        if dist[u] == inf break end # All remaining vertices are inaccessible from source
        for v in 1:n                # Each neighbor v of u
            if graph[u, v] != 0 && v ∈ Q # where v has not yet been visited
                alt = dist[u] + graph[u, v]
                if alt < dist[v]    # Relax (u, v, a)
                    dist[v] = alt
                    prev[v] = u
                end
            end
        end
    end

    return dist, prev
end

function printpath(prev::Dict, target::Int)
    path = "$target"
    while prev[target] != 0
        target = prev[target]
        path = "$target -> " * path
    end
    println(path)
end

const graph = [
    0 1 0 0 0 0;
    1 0 1 0 1 0;
    0 1 0 0 0 1;
    0 0 0 0 1 0;
    0 1 0 1 0 0;
    0 0 1 0 0 0]

dist, path = dijkstra(graph)
printpath(path, 6)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

for T in [Float64, Float32]

  #################
  # ProbabilityWeight semiring
  #################
  a = ProbabilityWeight(rand(T))
  b = ProbabilityWeight(rand(T))
  c = ProbabilityWeight(rand(T))
  d = ProbabilityWeight{T}(a) # check ProbabilityWeight(x::ProbabilityWeight)
  e = ProbabilityWeight{T}( ProbabilityWeight{T == Float64 ? Float32 : Float64}(get(a)) )
  f = ProbabilityWeight{T}( rand(T == Float64 ? Float32 : Float64) )
  @test typeof(e) == ProbabilityWeight{T}
  @test typeof(f) == ProbabilityWeight{T}

  @test (a*b).x == a.x*b.x
  @test (a + b).x == a.x + b.x
  @test (a/b).x == a.x/b.x
  @test reverse(a) == a
  @test FiniteStateTransducers.reverseback(a) == a
  o = one(ProbabilityWeight{T})
  @test o.x  == one(T)
  z = zero(ProbabilityWeight{T})
  @test z.x == zero(T)

  # + commutative monoid
  @test (a+b)+c ≈ a+(b+c)
  @test a+z ≈ z+a ≈ a
  @test a+b ≈ b+a
  # * monoid
  @test (a*b)*c ≈ a*(b*c)
  @test a*o == o*a == a # identity
  @test a*z == z*a == z # annihilator
  # distribution over addition
  @test a*(b+c) ≈ (a*b)+(a*c)
  @test (a+b)*c ≈ (a*c)+(b*c)
  # division
  c = a*b
  @test a*(c/a) ≈ c
  @test (c/b)*b ≈ c

  # properties
  W = ProbabilityWeight{T}
  @test FiniteStateTransducers.isidempotent(W) == false
  @test FiniteStateTransducers.iscommulative(W)
  @test FiniteStateTransducers.isleft(W)
  @test FiniteStateTransducers.isright(W)
  @test FiniteStateTransducers.isweaklydivisible(W)
  @test FiniteStateTransducers.iscomplete(W)
  @test FiniteStateTransducers.ispath(W) == false
  @test FiniteStateTransducers.reversetype(W) == W

  a = parse(W,"1.0")
  @test a.x == T(1.0)

  #################
  # LogWeight semiring
  #################
  a = LogWeight(1)
  b = LogWeight(2)
  @test (a + b).x ≈ log(exp(a.x)+exp(b.x))

  a = LogWeight(2)
  b = LogWeight(1)
  @test (a + b).x ≈ log(exp(a.x)+exp(b.x))

  a = LogWeight(rand(T))
  b = LogWeight(rand(T))
  c = LogWeight(rand(T))
  d = LogWeight{T}(a) # check ProbabilityWeight(x::ProbabilityWeight)
  e = LogWeight{T}( LogWeight{T == Float64 ? Float32 : Float64}(get(a)) )
  @test typeof(e) == LogWeight{T}

  @test (a*b).x ≈ a.x + b.x
  @test (a + b).x ≈ log(exp(a.x)+exp(b.x))
  @test (a / b).x ≈ log(exp(a.x)/exp(b.x))
  @test reverse(a) == a
  @test FiniteStateTransducers.reverseback(a) == a
  o = one(LogWeight{T})
  @test o.x  == 0
  z = zero(LogWeight{T})
  @test z.x == T(log(0))

  # + commutative monoid
  @test (a+b)+c ≈ a+(b+c)
  @test a+z ≈ z+a ≈ a
  @test a+b ≈ b+a
  # * monoid
  @test (a*b)*c ≈ a*(b*c)
  @test a*o == o*a == a # identity
  @test a*z == z*a == z # annihilator
  # distribution over addition
  @test a*(b+c) ≈ (a*b)+(a*c)
  @test (a+b)*c ≈ (a*c)+(b*c)
  # division
  c = a*b
  @test a*(c/a) ≈ c
  @test (c/b)*b ≈ c

  # properties
  W = LogWeight{T}
  @test FiniteStateTransducers.isidempotent(W) == false
  @test FiniteStateTransducers.iscommulative(W)
  @test FiniteStateTransducers.isleft(W)
  @test FiniteStateTransducers.isright(W)
  @test FiniteStateTransducers.isweaklydivisible(W)
  @test FiniteStateTransducers.iscomplete(W)
  @test FiniteStateTransducers.ispath(W) == false
  @test FiniteStateTransducers.reversetype(W) == W

  # parse
  a = parse(W,"1.0")
  @test a.x == T(1.0)

  #################
  # NLogWeight semiring
  #################
  a = NLogWeight(1)
  b = NLogWeight(2)
  @test (a + b).x ≈ -log(exp(-a.x)+exp(-b.x))

  a = NLogWeight(2)
  b = NLogWeight(1)
  @test (a + b).x ≈ -log(exp(-a.x)+exp(-b.x))

  a = NLogWeight(rand(T))
  b = NLogWeight(rand(T))
  c = NLogWeight(rand(T))
  d = NLogWeight{T}(a) # check ProbabilityWeight(x::ProbabilityWeight)
  e = NLogWeight{T}( NLogWeight{T == Float64 ? Float32 : Float64}(get(a)) )
  @test typeof(e) == NLogWeight{T}

  @test (a*b).x ≈ a.x + b.x
  @test (a + b).x ≈ -log(exp(-a.x)+exp(-b.x))
  @test (a / b).x ≈ log(exp(a.x)/exp(b.x))
  @test reverse(a) == a
  @test FiniteStateTransducers.reverseback(a) == a
  o = one(NLogWeight{T})
  @test o.x  == 0
  z = zero(NLogWeight{T})
  @test z.x == T(-log(0))

  # + commutative monoid
  @test (a+b)+c ≈ a+(b+c)
  @test a+z ≈ z+a ≈ a
  @test a+b ≈ b+a
  # * monoid
  @test (a*b)*c ≈ a*(b*c)
  @test a*o == o*a == a # identity
  @test a*z == z*a == z # annihilator
  # distribution over addition
  @test a*(b+c) ≈ (a*b)+(a*c)
  @test (a+b)*c ≈ (a*c)+(b*c)
  # division
  c = a*b
  @test a*(c/a) ≈ c
  @test (c/b)*b ≈ c

  # properties
  W = NLogWeight{T}
  @test FiniteStateTransducers.isidempotent(W) == false
  @test FiniteStateTransducers.iscommulative(W)
  @test FiniteStateTransducers.isleft(W)
  @test FiniteStateTransducers.isright(W)
  @test FiniteStateTransducers.isweaklydivisible(W)
  @test FiniteStateTransducers.iscomplete(W)
  @test FiniteStateTransducers.ispath(W) == false
  @test FiniteStateTransducers.reversetype(W) == W

  # parse
  a = parse(W,"1.0")
  @test a.x == T(1.0)

  #################
  # TropicalWeight semiring
  #################
  a = TropicalWeight(rand(T))
  b = TropicalWeight(rand(T))
  c = TropicalWeight(rand(T))
  d = TropicalWeight{T}(a) # check ProbabilityWeight(x::ProbabilityWeight)
  e = TropicalWeight{T}( TropicalWeight{T == Float64 ? Float32 : Float64}(get(a)) )
  @test typeof(e) == TropicalWeight{T}

  @test (a*b).x == a.x + b.x
  @test (a + b).x == min(a.x,b.x)
  @test (a / b).x ≈ log(exp(a.x)/exp(b.x))
  o = one(TropicalWeight{T})
  @test o.x  == 0
  z = zero(TropicalWeight{T})
  @test z.x == T(Inf)

  # + commutative monoid
  @test (a+b)+c ≈ a+(b+c)
  @test a+z ≈ z+a ≈ a
  @test a+b ≈ b+a
  # * monoid
  @test (a*b)*c ≈ a*(b*c)
  @test a*o == o*a == a # identity
  @test a*z == z*a == z # annihilator
  # distribution over addition
  @test a*(b+c) ≈ (a*b)+(a*c)
  @test (a+b)*c ≈ (a*c)+(b*c)
  # division
  c = a*b
  @test a*(c/a) ≈ c
  @test (c/b)*b ≈ c

  W = TropicalWeight{T}
  @test FiniteStateTransducers.isidempotent(W)
  @test FiniteStateTransducers.iscommulative(W)
  @test FiniteStateTransducers.isleft(W)
  @test FiniteStateTransducers.isright(W)
  @test FiniteStateTransducers.isweaklydivisible(W)
  @test FiniteStateTransducers.iscomplete(W)
  @test FiniteStateTransducers.ispath(W)
  @test FiniteStateTransducers.reversetype(W) == W

  # parse
  a = parse(W,"1.0")
  @test a.x == T(1.0)
end

a = BoolWeight(rand([true;false]))
b = BoolWeight(rand([true;false]))
c = BoolWeight(rand([true;false]))
d = BoolWeight(a) # check ProbabilityWeight(x::ProbabilityWeight)

@test (a*b).x == ( a.x && b.x)
@test (a + b).x == ( a.x || b.x )
o = one(BoolWeight)
@test o.x  == true
z = zero(BoolWeight)
@test z.x == false
@test_throws MethodError a/b

# + commutative monoid
@test (a+b)+c ≈ a+(b+c)
@test a+z ≈ z+a ≈ a
@test a+b ≈ b+a
# * monoid
@test (a*b)*c ≈ a*(b*c)
@test a*o == o*a == a # identity
@test a*z == z*a == z # annihilator
# distribution over addition
@test a*(b+c) ≈ (a*b)+(a*c)
@test (a+b)*c ≈ (a*c)+(b*c)

W = BoolWeight
@test FiniteStateTransducers.isidempotent(W)
@test FiniteStateTransducers.iscommulative(W)
@test FiniteStateTransducers.isleft(W)
@test FiniteStateTransducers.isright(W)
@test FiniteStateTransducers.isweaklydivisible(W) == false
@test FiniteStateTransducers.ispath(W)
@test FiniteStateTransducers.reversetype(W) == W

a = parse(W,"true")
@test a.x == true

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

a = LeftStringWeight("ciao")
b = LeftStringWeight("caro")
c = LeftStringWeight("mona")
d = LeftStringWeight(a) # check ProbabilityWeight(x::ProbabilityWeight)

@test (a*b).x == "ciaocaro"
@test (a*c).x == "ciaomona"
@test (a + b).x == "c"
@test (a + c).x == ""

o = one(LeftStringWeight)
@test o.x  == ""
@test o.iszero  == false
z = zero(LeftStringWeight)
@test z.x == ""
@test z.iszero

# + commutative monoid
@test (a+b)+c == a+(b+c)
@test a+z == z+a == a
@test a+b == b+a
# * monoid
@test (a*b)*c == a*(b*c)
@test a*o == o*a == a # identity
@test a*z == z*a == z # annihilator
# distribution over addition
@test a*(b+c) == (a*b)+(a*c)
# division
c = a*b
@test a*(c/a) == c
@test (c/b)*b != c

W = LeftStringWeight
@test FiniteStateTransducers.isidempotent(W)
@test FiniteStateTransducers.iscommulative(W) == false
@test FiniteStateTransducers.isleft(W)
@test FiniteStateTransducers.isright(W) == false
@test FiniteStateTransducers.isweaklydivisible(W)
@test FiniteStateTransducers.ispath(W) == false
@test FiniteStateTransducers.reversetype(W) == RightStringWeight

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

a = RightStringWeight("ciao")
b = RightStringWeight("caro")
c = RightStringWeight("mona")
d = RightStringWeight(a) # check ProbabilityWeight(x::ProbabilityWeight)

@test (a*b).x == "ciaocaro"
@test (a*c).x == "ciaomona"
@test (a + b).x == "o"
@test (a + c).x == ""

o = one(RightStringWeight)
@test o.x  == ""
@test o.iszero == false
z = zero(RightStringWeight)
@test z.x == ""
@test z.iszero

# + commutative monoid
@test (a+b)+c == a+(b+c)
@test a+z == z+a == a
@test a+b == b+a
# * monoid
@test (a*b)*c == a*(b*c)
@test a*o == o*a == a # identity
@test a*z == z*a == z # annihilator
# distribution over addition
@test (a+b)*c == (a*c)+(b*c)
# division
c = a*b
@test a*(c/a) != c
@test (c/b)*b == c

W = RightStringWeight
@test FiniteStateTransducers.isidempotent(W)
@test FiniteStateTransducers.iscommulative(W) == false
@test FiniteStateTransducers.isleft(W) == false
@test FiniteStateTransducers.isright(W)
@test FiniteStateTransducers.isweaklydivisible(W)
@test FiniteStateTransducers.ispath(W) == false
@test FiniteStateTransducers.reversetype(W) == LeftStringWeight

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

W1, W2 = LeftStringWeight, TropicalWeight{Float64}
a = ProductWeight(W1("ciao"),W2(1.0))
b = ProductWeight(W1("caro"),W2(2.0))
c = ProductWeight(W1("mona"),W2(3.0))
d = ProductWeight(a) # check ProbabilityWeight(x::ProbabilityWeight)
dd = typeof(d)(a)

@test (a*b) == ProductWeight(W1("ciaocaro"),W2(3.0))
@test (a*c) == ProductWeight(W1("ciaomona"),W2(4.0))
@test (a+b) == ProductWeight(W1("c"),W2(1.0))
@test (a+c) == ProductWeight(W1(""),W2(1.0))

W = typeof(a)
o = one(W)
@test get(o)[1]  == one(W1)
@test get(o)[2]  == one(W2)
z = zero(typeof(a))
@test get(z)[1] == zero(W1)
@test get(z)[2] == zero(W2)

# + commutative monoid
@test (a+b)+c == a+(b+c)
@test a+z == z+a == a
@test a+b == b+a
# * monoid
@test (a*b)*c == a*(b*c)
@test a*o == o*a == a # identity
@test a*z == z*a == z # annihilator
# distribution over addition
@test a*(b+c) == (a*b)+(a*c)
# division
c = a*b
@test a*(c/a) == c

@test FiniteStateTransducers.isidempotent(W)
@test FiniteStateTransducers.iscommulative(W) == false
@test FiniteStateTransducers.isleft(W)
@test FiniteStateTransducers.isright(W) == false
@test FiniteStateTransducers.isweaklydivisible(W)
@test FiniteStateTransducers.ispath(W) == false
@test FiniteStateTransducers.reversetype(W) == ProductWeight{2,Tuple{RightStringWeight,TropicalWeight{Float64}}}


text:
#==========================================================================================#

# TYPE

mutable struct Probit <: MLE

    sample::Microdata
    β::Vector{Float64}
    V::AbstractMatrix{Float64}

    Probit() = new()
end

#==========================================================================================#

# CONSTRUCTOR

function Probit(MD::Microdata)
    obj        = Probit()
    obj.sample = MD
    return obj
end

#==========================================================================================#

# ESTIMATION

function _fit!(obj::Probit, ::UnitWeights)

    y  = iszero.(getvector(obj, :response))
    x  = getmatrix(obj, :control)
    μ  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    p  = mean(y)
    p  = 1.0 / normpdf(norminvcdf(p))
    β₀ = lmul!(p, x \ y)

    function L(β::Vector)

        mul!(μ, x, β)
        ll = 0.0

        @inbounds for (yi, μi) in zip(y, μ)
            ll -= (yi ? normlogccdf(μi) : normlogcdf(μi))
        end

        return ll
    end

    function G!(g::Vector, β::Vector)

        mul!(μ, x, β)

        @inbounds for (i, (yi, μi)) in enumerate(zip(y, μ))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            μ[i] = - normpdf(μi) / ηi
        end

        mul!(g, transpose(x), μ)
    end

    function LG!(g::Vector, β::Vector)

        mul!(μ, x, β)
        ll = 0.0

        @inbounds for (i, (yi, μi)) in enumerate(zip(y, μ))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            μ[i] = - normpdf(μi) / ηi
            ll  -= (yi ? log(- ηi) : log(ηi))
        end

        mul!(g, transpose(x), μ)

        return ll
    end

    function H!(h::AbstractMatrix, β::Vector)

        mul!(μ, x, β)

        @inbounds for (i, (yi, μi)) in enumerate(zip(y, μ))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            ηi   = normpdf(μi) / ηi
            μ[i] = muladd(μi, ηi, abs2(ηi))
        end

        xx .= x .* μ
        mul!(h, transpose(x), xx)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("likelihood maximization did not converge")
    end
end

function _fit!(obj::Probit, w::AbstractWeights)

    y  = getvector(obj, :response)
    x  = getmatrix(obj, :control)
    μ  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    p  = mean(y)
    p  = 1.0 / normpdf(norminvcdf(p))
    β₀ = lmul!(p, x \ y)

    function L(β::Vector)

        mul!(μ, x, β)
        ll = 0.0

        @inbounds for (yi, μi, wi) in zip(y, μ, w)
            ll -= wi * (yi ? normlogccdf(μi) : normlogcdf(μi))
        end

        return ll
    end

    function G!(g::Vector, β::Vector)

        mul!(μ, x, β)

        @inbounds for (i, (yi, μi, wi)) in enumerate(zip(y, μ, w))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            μ[i] = - wi * normpdf(μi) / ηi
        end

        mul!(g, transpose(x), μ)
    end

    function LG!(g::Vector, β::Vector)

        mul!(μ, x, β)
        ll = 0.0

        @inbounds for (i, (yi, μi, wi)) in enumerate(zip(y, μ, w))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            μ[i] = - wi * normpdf(μi) / ηi
            ll  -= wi * (yi ? log(- ηi) : log(ηi))
        end

        mul!(g, transpose(x), μ)

        return ll
    end

    function H!(h::AbstractMatrix, β::Vector)

        mul!(μ, x, β)

        @inbounds for (i, (yi, μi, wi)) in enumerate(zip(y, μ, w))
            ηi   = (yi ? (normcdf(μi) - 1.0) : normcdf(μi))
            ηi   = normpdf(μi) / ηi
            μ[i] = wi * muladd(μi, ηi, abs2(ηi))
        end

        xx .= x .* μ
        mul!(h, transpose(x), xx)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("likelihood maximization did not converge")
    end
end

#==========================================================================================#

# SCORE (DERIVATIVE OF THE LIKELIHOOD FUNCTION)

function score(obj::Probit)

    y = getvector(obj, :response)
    x = getmatrix(obj, :control)
    v = x * obj.β

    @inbounds for (i, (yi, vi)) in enumerate(zip(y, v))
        ηi   = (iszero(yi) ? (normcdf(vi) - 1.0) : normcdf(vi))
        v[i] = normpdf(vi) / ηi
    end

    return Diagonal(v) * x
end

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

function jacobian(obj::Probit, ::UnitWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :control)
    v = x * obj.β

    @inbounds for (i, (yi, vi)) in enumerate(zip(y, v))
        ηi   = (iszero(yi) ? (normcdf(vi) - 1.0) : normcdf(vi))
        ηi   = normpdf(vi) / ηi
        v[i] = muladd(vi, ηi, abs2(ηi))
    end

    return - crossprod(x, v)
end

function jacobian(obj::Probit, w::AbstractWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :control)
    v = x * obj.β

    @inbounds for (i, (yi, vi, wi)) in enumerate(zip(y, v, w))
        ηi   = (iszero(yi) ? (normcdf(vi) - 1.0) : normcdf(vi))
        ηi   = normpdf(vi) / ηi
        v[i] = wi * muladd(vi, ηi, abs2(ηi))
    end

    return - crossprod(x, v)
end

#==========================================================================================#

# LINEAR PREDICTOR

function linear_predictor(obj::Probit, MD::Microdata)
    (getnames(obj, :control) != getnames(MD, :control)) && throw("missing variables")
    getmatrix(MD, :control) * obj.β
end

# FITTED VALUES

predict(obj::Probit, MD::Microdata) = normcdf.(linear_predictor(obj, MD))

# DERIVATIVE OF FITTED VALUES

function jacobexp(obj::Probit)
    x  = getmatrix(obj, :control)
    v  = x * obj.β
    v .= normpdf.(v)
    return Diagonal(v) * x
end

#==========================================================================================#

# UTILITIES

coefnames(obj::Probit) = getnames(obj, :control)

# LIKELIHOOD FUNCTION

function _loglikelihood(obj::Probit, ::UnitWeights)

    y  = getvector(obj, :response)
    μ  = linear_predictor(obj)
    ll = 0.0

    @inbounds for (yi, μi) in zip(y, μ)
        ll += (iszero(yi) ? normlogccdf(μi) : normlogcdf(μi))
    end

    return ll
end

function _loglikelihood(obj::Probit, w::AbstractWeights)

    y  = getvector(obj, :response)
    μ  = linear_predictor(obj)
    ll = 0.0

    @inbounds for (yi, μi, wi) in zip(y, μ, w)
        ll += wi * (iszero(yi) ? normlogccdf(μi) : normlogcdf(μi))
    end

    return ll
end

# LIKELIHOOD FUNCTION UNDER NULL MODEL

function _nullloglikelihood(obj::Probit, w::AbstractWeights)
    y = getvector(obj, :response)
    μ = mean(y, w)
    return nobs(obj) * (μ * log(μ) + (1.0 - μ) * log(1.0 - μ))
end

# DEVIANCE

deviance(obj::Probit) = - 2.0 * _loglikelihood(obj, getweights(obj))

# DEVIANCE UNDER NULL MODEL

nulldeviance(obj::Probit) = - 2.0 * _nullloglikelihood(obj, getweights(obj))

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
using Graphs
using DataStructures: Queue, enqueue!, dequeue!, iterate
include("Util.jl")

@enum MESSAGE begin
    explorer
    confirmation
    echo
    request
    token
end

abstract type MessageContent end
abstract type StopContent <: MessageContent end

struct Letter{T<:Integer,S<:MessageContent}
    e::Edge{T}
    m::MESSAGE
    content::S
end

struct Network{T<:Integer}
    graph::AbstractGraph{T}
    messages::Queue{Letter}
    Network{T}(graph::AbstractGraph{T}) where {T<:Integer} = new(graph, Queue{Letter}())
end

add_letter!(nw::Network, l::Letter) = enqueue!(nw.messages, l)
pop_letter!(nw::Network) = dequeue!(nw.messages)
next_letter(nw::Network) = first(nw.messages)
all_letters(nw::Network) = nw.messages
has_letters(nw::Network) = !isempty(nw.messages)

send!(nw::Network{T}, src::T, dst::T, m::MESSAGE, c::S) where {T<:Integer,S<:MessageContent} = add_letter!(nw, Letter(Edge(src, dst), m, c))
send!(nw::Network{T}, src::T, dsts::AbstractArray{T}, m::MESSAGE, c::S) where {T<:Integer,S<:MessageContent} = foreach(dst -> send!(nw, src, dst, m, c), dsts)

max_stars_repo_name: tomtuamnuq/dist-algorithms

max_issues_repo_name: tomtuamnuq/dist-algorithms

max_forks_repo_name: tomtuamnuq/dist-algorithms


text:
############################################################################ hyperparameter updates
function getParameters(data::StateSpace{N,T}; maxIter::Int=800, tol=1e-12, kwargs...) where {N,T<:Number}
	newtonOptions = NewtonPar(verbose=false,maxIter=maxIter,tol=T(tol))

	# support for StaticArrays github.com/JuliaArrays/StaticArrays.jl/issues/73
	newtonOptions = @set newtonOptions.linsolver.useFactorization = false
	return ContinuationPar(

        pMin=minimum(data.parameter), pMax=maximum(data.parameter), maxSteps=10*length(data.parameter),
        ds=step(data.parameter), dsmax=step(data.parameter), dsmin=step(data.parameter),

		newtonOptions=newtonOptions, detectFold=false, detectBifurcation=true,
		saveEigenvectors=false, nev=N )
end

function updateParameters!(parameters::ContinuationPar{T, S, E}, steady_states::Vector{Branch{V,T}};
    resolution=400 ) where {T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver}

    # estimate scale from steady state curves
    branch_points = map(length,steady_states)
    ds = maximum(branch_points)*parameters.ds/resolution
    parameters = setproperties(parameters;ds=ds,dsmin=ds,dsmax=ds)
end

function mean(f::Function,x::AbstractArray{T}; type::Symbol=:geometric, kwargs...) where T
	return type == :geometric ? mapreduce(y->f(y)^(1/length(x)),*,x;init=1,kwargs...) : mapreduce(y->f(y)/length(x),+,x;init=zero(f(zero(T))),kwargs...)
end

############################################################################# training loop
function train!( F::Function, parameters::NamedTuple, data::StateSpace;
				 iter::Int=200, optimiser=Momentum(0.001), plot_solution = false,
				 convergence = 0.02, kwargs...)

	hyperparameters = getParameters(data;kwargs...)
	Loss = steady_states = nothing

	∇Loss = similar(parameters.θ)
	trajectory = typeof(parameters.θ)[]

	for i=1:iter
		try
			steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...)
			Loss,∇Loss = ∇loss(F,steady_states,parameters.θ,data;kwargs...)

		catch error
			printstyled(color=:red,   "Iteration $i\tError = $error\n") end
			printstyled(color=:yellow,"Iteration $i\tLoss = $Loss\n")

		printstyled(color=:blue,"$steady_states\n")
		println("Parameters\t$(parameters.θ)")
		println("Gradients\t$(∇Loss)")

		if isinf(Loss) throw("infinite loss") end
		push!(trajectory,copy(parameters.θ))

		predictions = unique([ s.z for branch ∈ steady_states for s ∈ branch if s.bif ]; atol=2*step(data.parameter))
		if (Loss < convergence) & (length(predictions) == length(data.targets))
			break
		end

		if plot_solution>0 if i%plot_solution==0 plot(steady_states,data) end end
		update!(optimiser, parameters.θ, ∇Loss )
	end

	return trajectory
end

############################################################################## loss evaluation helper
function loss(F::Function, θ::AbstractVector, data::StateSpace, hyperparameters::ContinuationPar; kwargs...)

	try
		parameters = (θ=θ,p=minimum(data.parameter))
		steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...)
		return loss(F,steady_states,θ,data;kwargs...)

	catch error
		printstyled(color=:red,"$error\n")
		return NaN
	end
end

function loss(F::Function, θ::AbstractVector, data::StateSpace; kwargs...)
	return loss(F,θ,data,getParameters(data); kwargs...)
end

function ∇loss(F::Function, θ::AbstractVector, data::StateSpace, hyperparameters::ContinuationPar; kwargs...)
	try
		parameters = (θ=θ,p=minimum(data.parameter))
		steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...)
		return ∇loss(F,steady_states,θ,data;kwargs...)

	catch error
		printstyled(color=:red,"$error\n")
		return NaN,fill(NaN,length(θ))
	end
end

function ∇loss(F::Function, θ::AbstractVector, data::StateSpace; kwargs...)
	return ∇loss(F,θ,data,getParameters(data); kwargs...)
end


text:
struct MultiFieldCellField{DS<:DomainStyle} <: CellField
  single_fields::Vector{<:CellField}
  trian::Triangulation
  domain_style::DS

  function MultiFieldCellField(single_fields::Vector{<:CellField})
    @assert length(single_fields) > 0
    f1 = first(single_fields)

    # Find a suitable domain style
    if any( map(i->DomainStyle(i)==ReferenceDomain(),single_fields) )
      domain_style = ReferenceDomain()
    else
      domain_style = PhysicalDomain()
    end

    trian = get_triangulation(f1)
    #@check all(map(i-> get_triangulation(i) === trian,single_fields))

    new{typeof(domain_style)}(single_fields,trian,domain_style)
  end
end

function CellData.get_data(f::MultiFieldCellField)
  s = """
  Function get_data is not implemented for MultiFieldCellField at this moment.
  You need to extract the individual fields and then evaluate them separatelly.

  If ever implement this, evaluating a `MultiFieldCellField` directly would provide,
  at each evaluation point, a tuple with the value of the different fields.
  """
  @notimplemented s
end

CellData.get_triangulation(f::MultiFieldCellField) = f.trian
CellData.DomainStyle(::Type{MultiFieldCellField{DS}}) where DS = DS()
num_fields(a::MultiFieldCellField) = length(a.single_fields)
Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i]
Base.iterate(a::MultiFieldCellField)  = iterate(a.single_fields)
Base.iterate(a::MultiFieldCellField,state)  = iterate(a.single_fields,state)


text:
"""
    selinv_from_ldlt(Fp,Fi,Fv; conj = conj) -> Bv

Compute the entries `Bv` of the inverse of `A = F.L*F.D*F.Lt` contained in
the sparsity pattern of the (incomplete) LDLt factorization `F`.
"""
function selinv_from_ldlt(Fp,Fi,Fv; conj = Base.conj)
    Tv = eltype(Fv)
    n = length(Fp)-1

    @assert length(Fp) >= 1
    @assert issorted(Fp)
    @assert length(Fi) >= Fp[end]-1
    @assert length(Fv) >= Fp[end]-1
    @assert all(1 .<= Fi .<= n)

    @inbounds begin
        # Return variables
        Bv = Vector{Tv}(undef,length(Fi))

        # Workspace for a single column
        Fv_aj = Vector{Tv}(undef,n)
        Bv_aj = Vector{Tv}(undef,n)

        # Main algorithm
        for j in reverse(1:n)
            # Initialise column
            for p in Fp[j]+1:Fp[j+1]-1
                i = Fi[p]
                Fv_aj[i] = Fv[p]
                Bv_aj[i] = zero(Tv)
            end

            # Pull updates into B[j+1:n,j]
            for p in Fp[j]+1:Fp[j+1]-1
                k = Fi[p]
                Fv_kj = Fv_aj[k]
                Bv_kj = Bv_aj[k] - Bv[Fp[k]]*Fv_kj
                for p in Fp[k]+1:Fp[k+1]-1
                    i = Fi[p]
                    Fv_ij = Fv_aj[i]
                    Bv_ik = Bv[p]
                    Bv_aj[i] -=      Bv_ik *Fv_kj
                    Bv_kj    -= conj(Bv_ik)*Fv_ij
                end
                Bv_aj[k] = Bv_kj
            end

            # Copy temporary column into B
            for p in Fp[j]+1:Fp[j+1]-1
                i = Fi[p]
                Fv_aj[i] = zero(Tv)
                Bv[p] = Bv_aj[i]
            end

            # Deal with diagonal
            d = inv(Fv[Fp[j]])
            for p in Fp[j]+1:Fp[j+1]-1
                d -= conj(Bv[p])*Fv[p]
            end
            Bv[Fp[j]] = d
        end
    end
    return Bv
end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text: A = Mmap.mmap(open("/devel/myfile.txt"), Array{UInt8,1})


text:
module Lines

using ..Ahorn, Maple

displayName = "Lines"
group = "Brushes"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing
material = '0'

line = nothing
lineBrush = nothing

function drawLine(line::Ahorn.Line, layer::Ahorn.Layer)
    points = Ahorn.pointsOnLine(line)
    pixels, ox, oy = Ahorn.nodesToBrushPixels(points)
    global lineBrush = Ahorn.Brush("Line", pixels, (ox, oy))

    Ahorn.drawBrush(lineBrush, layer, 1, 1)
end

function drawLines(layer::Ahorn.Layer, room::Ahorn.DrawableRoom, camera::Ahorn.Camera)
    if line !== nothing
        drawLine(line, toolsLayer)
    end
end

function applyLine!(line::Ahorn.Line, material::Char)
    Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Line $(material)", Ahorn.loadedState.room))

    Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    Ahorn.applyBrush!(lineBrush, tiles, material, 1, 1)
end

function cleanup()
    global line = nothing

    Ahorn.redrawLayer!(toolsLayer)
end

function setMaterials!(layer::Ahorn.Layer)
    tileNames = Ahorn.tileNames(layer)
    displayNames = Ahorn.getBrushMaterialsNames(layer, get(Ahorn.config, "sort_brush_tiles_by_name", true))

    Ahorn.setMaterialList!(displayNames, row -> row[1] == get(tileNames, material, nothing))
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    layerName = Ahorn.layerName(targetLayer)
    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(layerName)", tileNames["Air"])[1]

    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")
    Ahorn.updateLayerList!(["fgTiles", "bgTiles"], row -> row[1] == wantedLayer)

    Ahorn.redrawingFuncs["tools"] = drawLines
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)
    Ahorn.persistence["brushes_layer"] = selected

    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(selected)", tileNames["Air"])[1]
    setMaterials!(targetLayer)
end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    tileNames = Ahorn.tileNames(targetLayer)
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = tileNames[selected]
    global material = tileNames[selected]
end

function materialFiltered(list::Ahorn.ListContainer)
    tileNames = Ahorn.tileNames(targetLayer)
    Ahorn.selectRow!(list, row -> row[1] == get(tileNames, material, '0'))
end

function getFavorites()
    if targetLayer !== nothing
        key = "favorites_brushes_$(targetLayer.name)"

        return Ahorn.getFavorites(Ahorn.persistence, key)
    end

    return []
end

function materialDoubleClicked(material::String)
    if targetLayer !== nothing
        key = "favorites_brushes_$(targetLayer.name)"
        tileNames = Ahorn.tileNames(targetLayer)

        Ahorn.toggleFavorite(Ahorn.persistence, key, material)
        setMaterials!(targetLayer)
        Ahorn.updateMaterialFilter!(targetLayer.name)
    end
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "fgTiles")
end

function leftClick(x::Number, y::Number)
    if line !== nothing
        applyLine!(line, material)

        global line = nothing

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

function middleClick(x::Number, y::Number)
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    tileNames = Ahorn.tileNames(targetLayer)
    target = get(tiles.data, (y, x), '0')

    global material = target
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = material
    Ahorn.selectMaterialList!(tileNames[target])
end

function mouseMotion(x::Number, y::Number)
    # Only make a preview square if we aren't dragging
    if !Ahorn.mouseButtonHeld(0x1)
        newLine = Ahorn.Line(x, y, x, y)

        if newLine != line
            global line = newLine

            Ahorn.redrawLayer!(toolsLayer)
        end
    end
end

function selectionMotion(x1::Number, y1::Number, x2::Number, y2::Number)
    newLine = Ahorn.Line(x1, y1, x2, y2)

    if newLine != line
        global line = newLine

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function selectionFinish(x1::Number, y1::Number, x2::Number, y2::Number)
    if line !== nothing
        applyLine!(line, material)

        global line = nothing

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
module TestKinematics

using Test
using LinearAlgebra
using BenchmarkTools

using Flight.Kinematics
using Flight.Attitude
using Flight.Geodesy

export test_kinematics

function test_kinematics()
    @testset verbose = true "Geodesy Kinematics" begin
        @testset verbose = true "Initialization" begin test_init() end
        @testset verbose = true "Position Update" begin test_init() end
    end
end

function test_init()

    init = KinInit(
        q_nb = RQuat([1, 2, 3, -2]),
        Ob = NVectorAlt(ϕ = π/3, λ = -π/6, h = 1500),
        ω_lb_b = [0.1, 0.01, -0.4],
        v_eOb_b = [100, 5, -10])

    x = init_x0(init)
    #now rebuild the initializer from the kinematic state vector and check it
    #against the original
    q_lb = RQuat(x.pos.q_lb, normalization = false)
    q_el = RQuat(x.pos.q_el, normalization = false)
    h = x.pos.h[1]
    ω_eb_b = x.vel.ω_eb_b
    v_eOb_b = x.vel.v_eOb_b

    Ob = NVectorAlt(NVector(q_el), h)
    (R_N, R_E) = radii(Ob)
    q_nb = q_lb
    v_eOb_n = q_nb * v_eOb_b
    ω_el_n = [
        v_eOb_n[2] / (R_E + h),
        -v_eOb_n[1] / (R_N + h),
        0.0]
    ω_el_b = q_nb' * ω_el_n
    ω_lb_b = ω_eb_b - ω_el_b

    init_test = KinInit(q_nb = q_lb, Ob = Ob, ω_lb_b = ω_lb_b, v_eOb_b = v_eOb_b )

    @test init.q_nb ≈ init_test.q_nb
    @test init.Ob ≈ init_test.Ob
    @test init.ω_lb_b ≈ init_test.ω_lb_b
    @test init.v_eOb_b ≈ init_test.v_eOb_b

end

function test_fpos()

    init = KinInit(
        q_nb = RQuat([1, 2, 3, -2]),
        Ob = NVectorAlt(ϕ = π/3, λ = -π/6, h = 1500),
        ω_lb_b = [0.1, 0.01, -0.4],
        v_eOb_b = [100, 5, -10])

    x_kin = init_x0(init)

    kin_data = Y(Kin())
    ẋ_pos = init_x0(Pos())

    @btime f_kin!($kin_data, $ẋ_pos, $x_kin)

end

end #module


text:
acosh(X::Array{Float32}) = Vectorize.Accelerate.acosh(X)

sub!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.sub!(X, Y, Z)

atanh(X::Array{Float32}) = Vectorize.Accelerate.atanh(X)

round(X::Array{Float32}) = Vectorize.Accelerate.round(X)

sinpi(X::Array{Float32}) = Vectorize.Accelerate.sinpi(X)

acos!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.acos!(X, Y)

log(X::Array{Float32}) = Vectorize.Accelerate.log(X)

log!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.log!(X, Y)

log2!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.log2!(X, Y)

sub(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.sub(X, Y)

expm1(X::Array{Float64}) = Vectorize.Accelerate.expm1(X)

mul(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.mul(X, Y)

atan!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.atan!(X, Y)

sqrt!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.sqrt!(X, Y)

exp!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.exp!(X, Y)

max!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.max!(X, Y, Z)

log2(X::Array{Float64}) = Vectorize.Accelerate.log2(X)

abs!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.abs!(X, Y)

sinpi!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.sinpi!(X, Y)

asin!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.asin!(X, Y)

rec(X::Array{Float64}) = Vectorize.Accelerate.rec(X)

sqrt(X::Array{Float64}) = Vectorize.Accelerate.sqrt(X)

div(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.div(X, Y)

log2(X::Array{Float32}) = Vectorize.Accelerate.log2(X)

asinh(X::Array{Float32}) = Vectorize.Accelerate.asinh(X)

tanpi!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.tanpi!(X, Y)

floor(X::Array{Float32}) = Vectorize.Accelerate.floor(X)

exp2(X::Array{Float64}) = Vectorize.Accelerate.exp2(X)

expm1(X::Array{Float32}) = Vectorize.Accelerate.expm1(X)

cosh(X::Array{Float32}) = Vectorize.Accelerate.cosh(X)

sub(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.sub(X, Y)

asin(X::Array{Float64}) = Vectorize.Accelerate.asin(X)

rec!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.rec!(X, Y)

log1p(X::Array{Float32}) = Vectorize.Accelerate.log1p(X)

floor(X::Array{Float64}) = Vectorize.Accelerate.floor(X)

invsqrt!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.invsqrt!(X, Y)

min!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.min!(X, Y, Z)

min(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.min(X, Y)

div!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.div!(X, Y, Z)

cos(X::Array{Float64}) = Vectorize.Accelerate.cos(X)

cospi(X::Array{Float32}) = Vectorize.Accelerate.cospi(X)

round!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.round!(X, Y)

mean(X::Array{Float64}) = Vectorize.Accelerate.mean(X)

rec!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.rec!(X, Y)

sqrt!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.sqrt!(X, Y)

summag(X::Array{Float32}) = Vectorize.Accelerate.summag(X)

add(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.add(X, Y)

sum(X::Array{Float64}) = Vectorize.Accelerate.sum(X)

exp!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.exp!(X, Y)

ceil(X::Array{Float64}) = Vectorize.Accelerate.ceil(X)

ceil!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.ceil!(X, Y)

trunc!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.trunc!(X, Y)

ceil(X::Array{Float32}) = Vectorize.Accelerate.ceil(X)

asinh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.asinh!(X, Y)

exponent(X::Array{Float32}) = Vectorize.Accelerate.exponent(X)

cosh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.cosh!(X, Y)

asinh(X::Array{Float64}) = Vectorize.Accelerate.asinh(X)

abs!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.abs!(X, Y)

exponent!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.exponent!(X, Y)

log1p!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.log1p!(X, Y)

acos!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.acos!(X, Y)

cospi!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.cospi!(X, Y)

cospi(X::Array{Float64}) = Vectorize.Accelerate.cospi(X)

add!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.add!(X, Y, Z)

exp2(X::Array{Float32}) = Vectorize.Accelerate.exp2(X)

acos(X::Array{Float64}) = Vectorize.Accelerate.acos(X)

asin!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.asin!(X, Y)

log10(X::Array{Float32}) = Vectorize.Accelerate.log10(X)

sum(X::Array{Float32}) = Vectorize.Accelerate.sum(X)

minimum(X::Array{Float32}) = Vectorize.Accelerate.minimum(X)

log10!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.log10!(X, Y)

trunc(X::Array{Float32}) = Vectorize.Accelerate.trunc(X)

log!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.log!(X, Y)

invsqrt(X::Array{Float32}) = Vectorize.Accelerate.invsqrt(X)

atan(X::Array{Float64}) = Vectorize.Accelerate.atan(X)

sinh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.sinh!(X, Y)

sqrt(X::Array{Float32}) = Vectorize.Accelerate.sqrt(X)

tanpi!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.tanpi!(X, Y)

abs(X::Array{Float64}) = Vectorize.Accelerate.abs(X)

atan!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.atan!(X, Y)

invsqrt!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.invsqrt!(X, Y)

log10!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.log10!(X, Y)

atanh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.atanh!(X, Y)

atanh(X::Array{Float64}) = Vectorize.Accelerate.atanh(X)

sinpi(X::Array{Float64}) = Vectorize.Accelerate.sinpi(X)

cospi!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.cospi!(X, Y)

tan!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Yeppp.tan!(X, Y)

rec(X::Array{Float32}) = Vectorize.Accelerate.rec(X)

log(X::Array{Float64}) = Vectorize.Accelerate.log(X)

sin(X::Array{Float32}) = Vectorize.Accelerate.sin(X)

asin(X::Array{Float32}) = Vectorize.Accelerate.asin(X)

div!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.div!(X, Y, Z)

exp2!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.exp2!(X, Y)

round!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.round!(X, Y)

tanh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.tanh!(X, Y)

max(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.max(X, Y)

mul!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.mul!(X, Y, Z)

sub!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.sub!(X, Y, Z)

tan!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.tan!(X, Y)

sin!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.sin!(X, Y)

sumsqr(X::Array{Float32}) = Vectorize.Accelerate.sumsqr(X)

asinh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.asinh!(X, Y)

atanh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.atanh!(X, Y)

ceil!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.ceil!(X, Y)

cos!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.cos!(X, Y)

mul(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.mul(X, Y)

expm1!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.expm1!(X, Y)

add(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Yeppp.add(X, Y)

summag(X::Array{Float64}) = Vectorize.Accelerate.summag(X)

sinpi!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.sinpi!(X, Y)

abs(X::Array{Float32}) = Vectorize.Accelerate.abs(X)

log1p!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.log1p!(X, Y)

cos(X::Array{Float32}) = Vectorize.Accelerate.cos(X)

mean(X::Array{Float32}) = Vectorize.Accelerate.mean(X)

cosh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.cosh!(X, Y)

sin(X::Array{Float64}) = Vectorize.Accelerate.sin(X)

tan(X::Array{Float64}) = Vectorize.Yeppp.tan(X)

max!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.max!(X, Y, Z)

trunc(X::Array{Float64}) = Vectorize.Accelerate.trunc(X)

tanh(X::Array{Float32}) = Vectorize.Accelerate.tanh(X)

tanpi(X::Array{Float32}) = Vectorize.Accelerate.tanpi(X)

div(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.div(X, Y)

sumsqr(X::Array{Float64}) = Vectorize.Accelerate.sumsqr(X)

minimum(X::Array{Float64}) = Vectorize.Accelerate.minimum(X)

cosh(X::Array{Float64}) = Vectorize.Accelerate.cosh(X)

max(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.max(X, Y)

sinh(X::Array{Float32}) = Vectorize.Accelerate.sinh(X)

maximum(X::Array{Float64}) = Vectorize.Accelerate.maximum(X)

exponent!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.exponent!(X, Y)

exp2!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.exp2!(X, Y)

invsqrt(X::Array{Float64}) = Vectorize.Accelerate.invsqrt(X)

acos(X::Array{Float32}) = Vectorize.Accelerate.acos(X)

mul!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.mul!(X, Y, Z)

exp(X::Array{Float32}) = Vectorize.Accelerate.exp(X)

add!(X::Array{Float64}, Y::Array{Float64}, Z::Array{Float64}) = Vectorize.Accelerate.add!(X, Y, Z)

log1p(X::Array{Float64}) = Vectorize.Accelerate.log1p(X)

acosh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.acosh!(X, Y)

log2!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.log2!(X, Y)

sinh(X::Array{Float64}) = Vectorize.Accelerate.sinh(X)

acosh(X::Array{Float64}) = Vectorize.Accelerate.acosh(X)

sin!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.sin!(X, Y)

exp(X::Array{Float64}) = Vectorize.Accelerate.exp(X)

exponent(X::Array{Float64}) = Vectorize.Accelerate.exponent(X)

maximum(X::Array{Float32}) = Vectorize.Accelerate.maximum(X)

min!(X::Array{Float32}, Y::Array{Float32}, Z::Array{Float32}) = Vectorize.Accelerate.min!(X, Y, Z)

tanh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.tanh!(X, Y)

log10(X::Array{Float64}) = Vectorize.Accelerate.log10(X)

floor!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.floor!(X, Y)

expm1!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.expm1!(X, Y)

tanh(X::Array{Float64}) = Vectorize.Accelerate.tanh(X)

acosh!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.acosh!(X, Y)

round(X::Array{Float64}) = Vectorize.Accelerate.round(X)

cos!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Yeppp.cos!(X, Y)

min(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.min(X, Y)

tanpi(X::Array{Float64}) = Vectorize.Accelerate.tanpi(X)

atan(X::Array{Float32}) = Vectorize.Accelerate.atan(X)

floor!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.floor!(X, Y)

tan(X::Array{Float32}) = Vectorize.Accelerate.tan(X)

sinh!(X::Array{Float64}, Y::Array{Float64}) = Vectorize.Accelerate.sinh!(X, Y)

trunc!(X::Array{Float32}, Y::Array{Float32}) = Vectorize.Accelerate.trunc!(X, Y)


text:
## ----------------------------------------------------------------------------
using ProjAssistant
@quickactivate

@time begin
    import DataFrames: DataFrame
    import CSV

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
end

## ----------------------------------------------------------------------------
# The human glycoprotein
# A1AT has a weight of 52 kDa. It consists of a single polypeptide chain with 394
# amino acids, and it has about 14 % carbohydrates by weight.

# from https=>//www.drugbank.ca/polypeptides/P01009
a1at_sec = string(
    "MPSSVSWGILLLAGLCCLVPVSLAEDPQGDAAQKTDTSHHDQDHPTFNKITPNLAEFAFS",
    "LYRQLAHQSNSTNIFFSPVSIATAFAMLSLGTKADTHDEILEGLNFNLTEIPEAQIHEGF",
    "QELLRTLNQPDSQLQLTTGNGLFLSEGLKLVDKFLEDVKKLYHSEAFTVNFGDTEEAKKQ",
    "INDYVEKGTQGKIVDLVKELDRDTVFALVNYIFFKGKWERPFEVKDTEEEDFHVDQVTTV",
    "KVPMMKRLGMFNIQHCKKLSSWVLLMKYLGNATAIFFLPDEGKLQHLENELTHDIITKFL",
    "ENEDRRSASLHLPKLSITGTYDLKSVLGQLGITKVFSNGADLSGVTEEAPLKLSKAVHKA",
    "VLTIDEKGTEAAGAMFLEAIPMSIPPEVKFNKPFVFLMIEQNTKSPLFMGKVVNPTQK"
)

aa_map = Dict("CYS"=> "C", "ASP"=> "D", "SER"=> "S", "GLN"=> "Q", "LYS"=> "K",
     "ILE"=> "I", "PRO"=> "P", "THR"=> "T", "PHE"=> "F", "ASN"=> "N",
     "GLY"=> "G", "HIS"=> "H", "LEU"=> "L", "ARG"=> "R", "TRP"=> "W",
     "ALA"=> "A", "VAL"=>"V", "GLU"=> "E", "TYR"=> "Y", "MET"=> "M")
for (k, v) in aa_map
    aa_map[v] = k
end

# relative abundance
a1at_aa_rel_ab = Dict()
aa_count = length(a1at_sec)
for i in a1at_sec
    aa = aa_map[string(i)]
    if haskey(a1at_aa_rel_ab, aa)
        a1at_aa_rel_ab[aa] += 1 / aa_count
    else
        a1at_aa_rel_ab[aa] = 1 / aa_count
    end
end

## ----------------------------------------------------------------------------s
df = DataFrame(id = collect(keys(a1at_aa_rel_ab)), rel_ab = collect(values(a1at_aa_rel_ab)))
CSV.write(Rd._a1at_aa_rel_abundance_proc_file(), df)
println("created $(relpath(Rd._a1at_aa_rel_abundance_proc_file()))")


text:
"""
# module Inputs

# Examples

```jldoctest
julia>
```
"""
module Inputs

using AbInitioSoftwareBase.Inputs: Namelist, groupname
using QuantumESPRESSOBase.Inputs: QuantumESPRESSOInput
using PyFortran90Namelists: Parser

struct InvalidInput
    msg::String
end

function Base.tryparse(::Type{T}, str::AbstractString) where {T<:Namelist}
    d::Dict{String,Any} = Parser().reads(str)
    return if haskey(d, lowercase(groupname(T)))
        dict = Dict(Symbol(k) => v for (k, v) in d[lowercase(groupname(T))])
        T(; dict...)
    end
end # function Base.tryparse

function Base.parse(::Type{T}, str::AbstractString) where {T<:Namelist}
    x = tryparse(T, str)
    if x === nothing
        throw(Meta.ParseError("cannot find namelist `$(groupname(T))`!"))
    else
        return x
    end
end # function Base.parse

include("PWscf.jl")
include("PHonon.jl")

# Idea from https://github.com/JuliaData/CSV.jl/blob/c3af297/src/CSV.jl#L64-L69
function Base.read(io::IO, ::Type{T}) where {T<:QuantumESPRESSOInput}
    str = read(io, String)
    return parse(T, str)
end
function Base.read(filename::AbstractString, ::Type{T}) where {T<:QuantumESPRESSOInput}
    str = read(filename, String)
    return parse(T, str)
end

end


text:
function dubins_LSL(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef,3)
    tmp0 = intermediate_results.d + intermediate_results.sa - intermediate_results.sb
    p_sq = 2 + intermediate_results.d_sq - (2*intermediate_results.c_ab) +
            (2 * intermediate_results.d * (intermediate_results.sa - intermediate_results.sb))
    tmp0 = round(tmp0; digits=10)
    p_sq = round(p_sq; digits=10)

    if p_sq >= 0
        tmp1 = atan((intermediate_results.cb - intermediate_results.ca), tmp0)
        tmp1 = round(tmp1; digits=10)
        out[1] = mod2pi(tmp1 - intermediate_results.α)
        out[2] = sqrt(p_sq)
        out[3] = mod2pi(intermediate_results.β - tmp1)
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end

function dubins_RSR(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef, 3)
    tmp0 = intermediate_results.d - intermediate_results.sa + intermediate_results.sb
    p_sq = 2 + intermediate_results.d_sq - (2 * intermediate_results.c_ab) +
            (2 * intermediate_results.d * (intermediate_results.sb - intermediate_results.sa))
    tmp0 = round(tmp0; digits=10)
    p_sq = round(p_sq; digits=10)

    if p_sq >= 0
        tmp1 = atan((intermediate_results.ca - intermediate_results.cb), tmp0)
        tmp1 = round(tmp1; digits=10)
        out[1] = mod2pi(intermediate_results.α - tmp1)
        out[2] = sqrt(p_sq)
        out[3] = mod2pi(tmp1 -intermediate_results.β)
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end

function dubins_LSR(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef,3)
    p_sq = -2 + (intermediate_results.d_sq) + (2 * intermediate_results.c_ab) +
                    (2 * intermediate_results.d * (intermediate_results.sa + intermediate_results.sb))
    p_sq = round(p_sq; digits=10)

    if p_sq >= 0
        p = sqrt(p_sq)
        tmp0 = atan((-intermediate_results.ca - intermediate_results.cb),
                    (intermediate_results.d + intermediate_results.sa + intermediate_results.sb)) - atan(-2.0, p)
        tmp0 = round(tmp0; digits=10)
        out[1] = mod2pi(tmp0 - intermediate_results.α)
        out[2] = p
        out[3] = mod2pi(tmp0 - mod2pi(intermediate_results.β))
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end

function dubins_RSL(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef,3)
    p_sq = -2 + intermediate_results.d_sq + (2 * intermediate_results.c_ab) -
            (2 * intermediate_results.d * (intermediate_results.sa + intermediate_results.sb))
    p_sq = round(p_sq; digits=10)

    if p_sq >= 0
        p = sqrt(p_sq)
        tmp0 = atan((intermediate_results.ca + intermediate_results.cb),
                        (intermediate_results.d - intermediate_results.sa - intermediate_results.sb)) - atan(2.0, p)
        tmp0 = round(tmp0; digits=10)
        out[1] = mod2pi(intermediate_results.α - tmp0)
        out[2] = p
        out[3] = mod2pi(intermediate_results.β - tmp0)
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end

function dubins_RLR(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef,3)
    tmp0 = (6.0 - intermediate_results.d_sq + 2*intermediate_results.c_ab +
            2*intermediate_results.d*(intermediate_results.sa - intermediate_results.sb)) / 8.0
    phi  = atan(intermediate_results.ca - intermediate_results.cb, intermediate_results.d - intermediate_results.sa + intermediate_results.sb)
    tmp0 = round(tmp0; digits=10)
    phi = round(phi; digits=10)

    if abs(tmp0) <= 1
        p = mod2pi((2*π) - acos(tmp0))
        t = mod2pi(intermediate_results.α - phi + mod2pi(p/2.0))
        p = round(p; digits=10)
        t = round(t; digits=10)
        out[1] = t
        out[2] = p
        out[3] = mod2pi(intermediate_results.α - intermediate_results.β - t + mod2pi(p))
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end

function dubins_LRL(intermediate_results::DubinsIntermediateResults)

    out = Vector{Float64}(undef,3)
    tmp0 = (6.0 - intermediate_results.d_sq + 2*intermediate_results.c_ab +
            2*intermediate_results.d*(intermediate_results.sb - intermediate_results.sa)) / 8.0
    phi = atan(intermediate_results.ca - intermediate_results.cb, intermediate_results.d + intermediate_results.sa - intermediate_results.sb)
    tmp0 = round(tmp0; digits=10)
    phi = round(phi; digits=10)

    if abs(tmp0) <= 1
        p = mod2pi(2*π - acos(tmp0))
        t = mod2pi(-intermediate_results.α - phi + p/2.0)
        p = round(p; digits=10)
        t = round(t; digits=10)
        out[1] = t
        out[2] = p
        out[3] = mod2pi(mod2pi(intermediate_results.β) - intermediate_results.α - t + mod2pi(p))
        return EDUBOK, out
    end

    return EDUBNOPATH, nothing
end


text:
function read_graphs(n, file)
    G1 = SimpleGraph(n)
    G2 = SimpleGraph(n)
    const_G2 = 0
    for line in eachline(file)
        #println(line)
        if line == "graph2"
            const_G2 += 1
        elseif line == "matrice log likelihood ratio"
            const_G2 += 1
        end
        if const_G2 == 0
            u = split(line, " : ")
            if length(u) >1
                add_edge!(G1, parse(Int, u[1]),parse(Int, u[2]))
                #println(parse(Int, u[1]), "|",  parse(Int, u[2]))
            end
        elseif const_G2 == 1
            u = split(line, " : ")
            if length(u) >1
                add_edge!(G2, parse(Int, u[1]),parse(Int, u[2]))
                #println(parse(Int, u[1]), "|",  parse(Int, u[2]))
            end
        end
    end
    return G1, G2
end

function Pair_IO(G1,G2,n,λ,s)
    N1 = Network(G1)
    N2 = Network(G2)
    #neig1 = create_neig_r(G1)
    #neig2 = create_neig_r(G2)
    return Pair_ER(N1,N2,SimpleGraph(n),n,λ,s,1:n,1:n,zeros(Int64,n,2))
end

function create_matrix_check(data)
    n = length(data)
    M_data = zeros(n,n)
    for (i,l) in enumerate(data)
        for (j,e) in enumerate(l)
            M_data[i,j] = e
        end
    end
    M_data
end

function check_data(PG,M_data,n_iter)
    n,_ = size(M_data)
    mp_fin = run_bp(PG,n_iter)
    M_lr = create_matrix_lr(PG,mp_fin)
    M_loglr = log.(M_lr)
    check = true
    for i in 1:n
        for j in 1:n
            if (M_data[i,j] - M_loglr[i,j]) > 0.00001
                #println(i,"|",j,"|", M_data[i,j],"|",M_loglr[i,j])
                check = false
            end
        end
    end
    check
end