text:
# TRANSMISSION OF WRATTEN FILTERS

# Compiled by Allie C. Peed, Jr. for The Eastman Kodak Company

# Data condensed from Kodak Wratten Filters for Scientific and Technical Use published by the Eastman Kodak Company, manufacturers of the filters.

# The following pages give (1) percentage luminous transmittance at wave lengths from 400 to 700μ at intervals of 10μ for the standard illuminant “C” adopted by the International Commission of Illumination, (2) dominant wavelength in millimicrons, and (3) percentage of excitation purity. Values of wave length followed by “c” indicate the complementary wave lengths of purple filters which do not have a dominant wave length.

# All colorimetric specifications are based on the 1931 standard ICI colorimetric and luminosity data.

# The transmittance data are given as representing standard samples of the filters. They are intended only for the information of users in choosing filters which will meet their requirements. Values taken from the tables of data should not be used by research workers as representing precisely the absorption characteristics of a particular filter. If such precise data are needed, they should be determined for the particular filter being used.

# Where the spectra extend into the ultraviolet this fact is indicated by an asterisk (*) in the transmission tables immediately beneath the filter number, and quantitative data are not given. The manufacturer should be consulted for this information. Transmission in the ultraviolet of wave lengths less than 330μ will be eliminated in the case of cemented filters, as glass absorbs ultraviolet radiation of wave lengths shorter than about 330μ.

# Stability ratings are given as three letter combinations following the filter description in the table below. In establishing the stability classifications each filter is exposed to a selected light source for a specific time interval. The following grading system is used to describe the result:

#       Class A — stable
#       Class B — relatively stable
#       Class C — somewhat unstable
#       Class D — unstable

# The classification letters, for example, AAA, describe the stability to the following three exposure tests in this order:

# 1.    Two weeks’ exposure to daylight in a south window
# 2.    Twenty-four hours’ exposure to a “Fade-Ometer”
# 3.    Two weeks’ exposure at two feet from a 1000-watt tungsten lamp.

# Filters are supplied in two forms: as lacquered gelatin film, or as a gelatin film cemented between pieces of optical glass. Filters in glass are cemented between sheets of plane-parallel glass, which is surfaced in quantities and is of sufficient accuracy for general photographic work, and for most scientific purposes.

# Most Wratten Gelatin Filters are stocked in 2- or 3-inch squares. Stocks of 2- or 3-inch square filters cemented in glass are maintained only in filters usually used for general photographic work.

# The booklet “Kodak Filters and Lens Attachments” gives more valuable information on this subject.

# No.	Description, use, and stability

#       Colorless

# 0     For compensating thickness of other gelatin filters in optical systems, AAA.
# 1     Absorbs ultraviolet below 360 mμ, DDD.
# 1A    Kodak Skylight Filter — Reduces excess bluishness in outdoor color photographs in open shade under a clear, blue sky, ACA.

#       Yellow

# 2B    Absorbs ultraviolet below 410 mμ, ACA.
# 3     Light yellow, CCD.
# 3N5   No. 3 plus 0.5 neutral density, AAA.
# 4     Light yellow — Approximate correction on panchromatic materials for outdoor scenes, including sky, CCC.
# 6     K1 — Light yellow — Partial correction outdoors, BBA.
# 8     K2 — Yellow — Full correction outdoors on Type B panchromatic materials.
#                       Widely used for proper sky, cloud, and foliage rendering. Green separation for Fluorescence Process, AAA.
# 8N5   No. 8 plus 0.5 neutral density, AAA.
# 9     K3 — Deep yellow. Moderate contrast in outdoor photography (with black-and-white films),AAA.
# 11    X1 — Greenish yellow. Correction for tungsten light on Type B panchromatic materials;
#                       also for daylight correction with Type C panchromatic materials in making # outdoor portraits, darkening skies, or lightening foliage, AAA.
# 12 	Minus blue. Haze cutting in aerial photography, AAA.
# 13 	X2 — Yellow green. Correction for Type C panchromatic materials in tungsten light, ABA.
# 15 	G — Deep yellow. Overcorrection in landscape photography. Contrast control in copying and in aerial infrared photography, AAA.
# 16 	Blue absorption, AAB.
# 18A	Transmits ultraviolet and infrared only (glass), AAA.

# 	    Oranges and Reds

# 21    Blue and blue-green absorption, CBB.
# 22    Yellow-orange. For increasing contrast in blue preparations in microscopy. Mercury yellow, BAC.
# 23A   Light red. Two-color projection — contrast effects, BAB.
# 24    Red for two-color photography (daylight or tungsten). White-flame-arc tricolor projection, AAB.
# 25 	A — Tricolor red for direct color separation. Contrast effects in commercial photography and in outdoor scenes. Two-color general viewing. Aerial infrared photography and haze cutting, AAA.
# 26 	Stereo red, AAA.
# 29 	Red color separation from transparencies and for the Kodak Fluorescence Process. Strong contrast effects. Copying blueprints. Tungsten tricolor projection, AAA.

# Magentas and Violets

# 30 	Green absorption, BBC.
# 31 	Green absorption, CCA.
# 32 	Minus green, CCD.
# 33 	Strong	green absorption, CCB.
# 34 	Violet, CDD.
# 34A 	Blue separation — Kodak Fluorescence Process, DCC.
# 35 	Contrast in microscopy, CDD.
# 36 	Dark violet, CCC.

# Blues and Blue-greens

# 38	Red absorption, BCA.
# 38A	Red absorption. Increasing contrast in visual microscopy, BBB
# 39 	Contrast control in printing motion-picture duplicates (glass) AAA.
# 40 	Green for two-olor photography (tungsten), CBC.
# 44	Minus red — Two-color general viewing, DDD.
# 44A	Minus red, DDD.
# 45	Contrast in microscopy, DDD.
# 45A 	Highest resolving power in visual microscopy, CDC.
# 46	Blue projection (experimental), DDD.
# 47 	Tricolor blue for direct color separation and from Kodak Ektacolor Film for Dye Transfer.
#                       Contrast effects in commercial photography. Tungsten and white- flame-arc tricolor projection, BBC.
# 47B	Tricolor blue for color separation from transparencies and from Kodak Ektacolor Film for Graphic Arts, BBB.
# 48	Green and red absorption, CBC.
# 48A	Green and red absorption, AAB.
# 49 	Dark blue, BCB.
# 49B	Very dark blue, BBB.
# 50 	Very dark blue. Mercury violet, CCC.

# Greens

# 52    Light green, AAB.
# 53    Medium green, CCB
# 54    Very dark green, AAA.
# 55    Stereo green, BBC.
# 56    Very light green, CBC.
# 57    Green for two-color photography (daylight), CBC.
# 57A   Light green, BBC.
# 58    Tricolor green for direct color separation. Contrast effects in commercial photography and microscopy, BBC.
# 59    Green for tricolor projection (white-flame-arc), BBB.
# 59A   Very light green, BBB.
# 60    Green for two-color photography (tungsten), BDC.
# 61    Green color separation from transparencies and Kodak Ektacolor Film. Tricolor projection (tungsten), ABC.
# 64    Red absorption (light), CDB.
# 65    Red absorption, ADB.
# 65A   Red absorption, CCD.
# 66    Contrast effects in microscopy and medical photography, DDC.
# 67A   Red absorption (light). Two-color projection, CDC.

# Narrow-band

# 70    Dark red. Infrared photography. Color separation for Kodak Ektacolor Film (with tungsten), ABC.
# 72B   Dark orange-yellow, CCC.
# 73    Dark yellow-green, ABB.
# 74    Dark green. Mercury green,	BBC.
# 75    Dark blue-green, ACC.
# 76    Dark violet (compound filter),	DDD.
# 77    Transmits 546 mμ mercury line (glass plus gelatin), AAA.
# 77A   Transmits 546 mμ mercury line (glass plus gelatin), AAA.

# Photometries

# 78    Bluish. Photometric filter (visual), BAB.
# 78AA  Bluish. Photometric filter (visual), BAA.
# 78A   Bluish. Photometric filter (visual), AAA.
# 78B   Bluish. Photometric filter (visual), AAA.
# 78C   Bluish. Photometric filter (visual), BAA.
# 86    Yellowish. Photometric filter (visual), BBA.
# 86A   Yellowish. Photometric filter (visual), AAA.
# 6B    Yellowish. Photometric filter (visual), BCA.
# 86C   Yellowish. Photometric filter (visual), AAA.

# Light Balancing

# 80A   For Kodachrome Film, Daylight Type, and photographic flood lamps, ABA.
# 81    Yellowish. For warmer color rendering.
# 81A   Yellowish. For Kodak Ektachrome Film, Type B, with photographic flood lamps.
# 81B   Yellowish. For warmer color rendering.
# 81C   Yellowish. For Kodachrome Film, Type A, with flash lamps.
# 81D   Yellowish. For Kodachrome Film, Type A, with flash lamps.
# 81EF  Yellowish. For Kodak Ektachrome Film, Type B, with flash lamps.
# 82    Bluish. For cooler color rendering.
# 82A   Bluish. For Kodachrome Film, Type A, with 3200 K lamps.
# 82B   Bluish. For cooler color tendering.
# 82C   Bluish. For cooler color rendering.
# 83    Yellowish. For 16 mm Commercial Kodachrome Film and daylight exposure, BBB.
# 85    Orange. For Type A Kodak color films and daylight exposure, BAA.
# 85B   Orange. For Kodak Ektachrome Film, Type B, and daylight exposure, BAB.

# Miscellaneous

# 79    Photographic sensitometry. Corrects 2360 K to 5500 K, AAA.
# 87    For infrared photography. Absorbs	visual.
# 87C   Absorbs visual, transmits infrared.
# 88A   For infrared photography. Absorbs	visual.
# 89B   For infrared photography, AAA.
# 90    Narrow-band viewing filter for judging brightness scale of scenes, CCD.
# 96    Neutral filters for controlling luminance, AAB.
# 97    Dichroic absorption, AAA.
# 102   Correction filter for Barrier-layer photocell, ABA.
# 106   Correction filter for S-4 type photocell, AAA.

"""
    wratten(f::Int)

    Returns `TSpec` of one of the Wratten filters of a given index `f` in the range `0 < f < 128`.

# Examples
```jldoctest
julia> wratten(77)
TSpec(Real[400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, 480.0, 490.0  …  610.0,
 620.0, 630.0, 640.0, 650.0, 660.0, 670.0, 680.0, 690.0, 700.0], Real[0.6409999999999999,
  0.665, 0.684, 0.695, 0.7, 0.6940000000000001, 0.675, 0.654, 0.629, 0.598  …  0.4, 0.376,
   0.355, 0.342, 0.336, 0.34, 0.341, 0.345, 0.348, 0.36], 1.0)
```
"""
function wratten(f::Int)
    if 0 < f < 128
        TSpec(WrattenVis[:,1],WrattenVis[:,f+1]/100.0,1.0)
    elseif 127 < f < 132
        TSpec(WrattenIR[:,1],WrattenIR[:,f-126]/100.0, 1.0)
    else
        error("Index $f does not exist. The allowed index range is 1 to 131.")
    end
end

"""
    wratten(f::Symbol)

Returns `TSpec` of one of the Wratten filters of a given name `f`.

# Examples

```jldoctest
julia> wratten(:n44)
TSpec(Real[400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, 480.0, 490.0  …  610.0,
620.0, 630.0, 640.0, 650.0, 660.0, 670.0, 680.0, 690.0, 700.0], Real[0.0044, 0.0036, 0.0063,
 0.0363, 0.131, 0.254, 0.365, 0.465, 0.536, 0.568…  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  0.0018, 0.016], 1.0)
```
"""
function wratten(f::Symbol)
    @inbounds for i in 1:length(Wrattenfilters)
        Wrattenfilters[i] == f ? (return wratten(i)) : nothing
    end
    error("Filter $f does not exist")
end

"""
    wratten_filter(f::Int, x::Real)

Returns `TSpec` of one of the Wratten filters of a given name adapted to `SPECENV` settings.
`x` is relative to unit thickness.

```jldoctest

```
"""
function wratten_filter(f::Int,x::Real, env = SPECENV)
    wr = wratten(f)
    extr = extrap(hermitespline(wr.λ, wr.t.^x, :kruger), collect(env.λmin:env.Δλ:env.λmax),:smooth)
    intr = interp(extr, env.λmin:env.Δλ:env.λmax)
    intr[2][intr[2] .< 0.0] .= 0.0
	intr[2][intr[2] .> 1.0] .= 1.0
    transmittance_spec(intr[1], intr[2], x)
end

"""
`wratten_filter(f::Int)`

returns `TSpec` of one of the Wratten filters at unit thickness.
"""
wratten_filter(f::Int) = wratten_filter(f, 1.0)

"""
`wratten_filter(f::Symbol)`

returns `TSpec` of one of the Wratten filters of a given name adapted to `SPECENV` settings.
"""
function wratten_filter(f::Symbol, x::Real)
    @inbounds for i in 1:length(Wrattenfilters)
        Wrattenfilters[i] == f ? (return wratten_filter(i, x)) : nothing
    end
    error("Filter $f does not exist")
end

"""
`wratten_filter(f::Symbol)`

returns `TSpec` of one of the Wratten filters of a given name at unit thickness `x` = 1.0

"""
wratten_filter(f::Symbol) = wratten_filter(f, 1.0)

const Wrattenfilters= [:n0, :n1, :n1A, :n2B, :n3, :n3N5, :n4, :n6, :n8, :n8N5, :n9, :n11, :n12, :n13, :n15, :n16, :n18A, :n21, :n22, :n23A, :n24, :n25, :n26, :n29, :n30, :n31, :n32, :n33, :n34, :n34A, :n35, :n36, :n38, :n38A, :n39, :n40, :n44, :n44A, :n45, :n45A, :n46, :n47, :n47B, :n48A, :n49, :n49B, :n50, :n52, :n53, :n54, :n55, :n56, :n57, :n57A, :n58, :n59, :n59A, :n60, :n61, :n64, :n65, :n65A, :n66, :n67A, :n70, :n72B, :n73, :n74, :n75, :n76, :n77, :n77A, :n78, :n78AA, :n78A, :n78B, :n78C, :n79, :n80A, :n81, :n81A, :n81B, :n81C, :n81D, :n81EF, :n82, :n82A, :n82B, :n82C, :n83, :n85, :n85B, :n86, :n86A, :n86B, :n86C, :n89B, :n90, :n96, :n97, :n102, :n106, :CC05R, :CC10R, :CC20R, :CC30R, :CC40R, :CC50R, :CC05B, :CC10B, :CC20B, :CC30B, :CC40B, :CC50B, :CC05G, :CC10G, :CC20G, :CC30G, :CC40G, :CC50G, :CC05Y, :CC10Y, :CC20Y, :CC30Y, :CC40Y, :CC50Y, :CC05M, :CC10M, :CC20M, :CC30M, :CC40M, :CC50M, :CC05C, :CC10C, :CC20C, :CC30C, :CC40C, :CC50C, :n87, :n87C, :n88A, :n89B]

# All Kodak Wratten filter transmittances in this array are given in %
# checked for correctness and errors/typos in the original table corrected.

"""
`const WrattenVis`

Array containing all 127 Wratten filters for the visible spectrum
"""
const WrattenVis=
   [400 88.0 85.0 59.0 19.0 0.00 0.00 0.00 7.40 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 48.6 13.8 38.0 0.85 64.0 0.00 48.0 36.5 60.5 33.4 85.2 0.00 0.44 2.52 0.00 0.00 1.20 7.80 16.0 0.96 5.65 3.30 1.70 0.45 2.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.00 0.00 0.00 12.3 1.10 0.00 0.00 0.00 0.00 0.00 0.22 0.00 0.00 37.2 43.0 56.0 64.1 74.9 24.0 67.6 77.7 65.1 55.1 46.1 38.2 30.7 83.0 80.1 76.7 73.4 13.5 6.00 1.59 0.50 8.00 20.0 44.0 0.00 0.00 4.28 0.00 1.12 0.00 81.0 73.0 61.5 51.6 42.5 36.4 87.0 85.5 82.2 80.2 77.0 74.1 87.6 86.6 85.6 84.2 82.3 80.9 87.3 86.0 83.9 82.3 80.4 78.8;
    410 88.5 85.5 76.0 48.0 0.00 0.00 0.00 8.32 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 47.4 14.5 37.9 0.71 70.1 0.10 57.0 45.5 66.5 41.2 78.2 0.00 0.36 3.39 0.00 0.00 0.60 17.4 29.5 3.16 10.0 4.28 2.00 0.39 1.51 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 9.20 0.00 0.00 13.0 0.93 0.00 0.00 0.00 0.00 0.00 0.18 0.00 0.00 41.7 46.0 58.6 66.5 76.6 26.0 73.1 78.1 65.9 55.8 46.6 38.4 31.5 83.7 80.8 78.0 75.0 13.1 18.0 9.32 0.81 12.2 26.1 55.0 0.00 0.00 4.91 0.00 0.96 0.00 81.0 72.4 60.0 50.0 40.0 33.9 87.5 86.4 84.0 82.5 80.3 78.4 88.2 87.7 86.6 85.7 84.6 83.6 88.2 87.5 85.2 84.5 83.4 82.7;
    420 88.9 86.0 82.0 67.0 0.00 0.00 0.00 10.4 0.00 0.00 0.00 0.16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 48.5 16.4 40.0 1.17 72.0 40.0 57.6 45.5 72.5 53.0 70.5 0.00 0.63 6.30 0.00 0.00 0.80 34.0 43.6 8.25 16.0 6.93 3.55 0.59 0.80 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 8.75 0.23 0.00 15.0 1.28 0.00 0.00 0.00 0.00 0.00 0.29 0.00 0.00 44.2 48.7 61.0 68.4 77.9 29.0 76.8 79.0 67.6 57.7 49.0 41.0 34.3 84.6 81.6 79.2 76.4 13.5 28.4 15.5 1.55 16.7 31.6 62.0 0.00 0.00 5.50 0.00 0.89 0.10 81.0 72.0 58.6 48.2 38.2 31.9 87.7 87.2 85.0 84.0 82.2 80.7 88.6 88.0 87.0 85.9 85.2 84.4 88.7 88.1 86.5 86.0 85.3 84.8;
    430 89.3 86.5 84.6 75.3 0.36 0.00 0.00 13.5 0.00 0.00 0.00 0.29 0.00 0.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 50.1 25.5 43.0 1.69 68.4 69.7 47.5 32.7 75.3 58.0 63.3 0.00 3.63 17.4 0.00 1.00 5.98 47.0 50.0 15.0 21.0 11.2 7.00 2.63 0.44 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.16 0.00 0.00 9.20 0.61 0.16 18.4 3.16 0.00 0.00 0.00 0.00 0.00 1.38 0.00 0.00 44.6 49.8 61.8 69.5 78.9 31.0 77.7 80.5 70.2 61.0 52.5 45.5 38.6 85.1 82.2 79.7 77.2 14.1 33.4 19.0 2.88 21.5 37.5 66.6 0.00 0.00 6.17 0.00 0.96 0.20 81.1 71.6 57.7 47.0 36.8 30.5 88.0 87.5 85.3 84.3 82.5 81.1 88.7 88.0 86.9 85.6 84.4 83.6 89.0 88.6 87.5 87.0 86.3 85.9;
    440 89.6 87.0 86.0 80.0 1.78 0.00 0.00 18.9 0.00 0.00 0.00 0.56 0.00 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 49.4 42.7 55.5 5.36 58.2 68.7 29.5 15.2 76.2 58.8 53.6 0.00 13.1 32.7 5.00 8.81 19.0 50.3 47.2 22.6 25.0 18.9 13.0 8.90 0.41 0.00 0.00 0.00 0.00 0.00 0.19 0.00 0.00 0.37 0.00 0.00 11.3 1.58 1.32 23.2 6.40 0.00 0.00 0.00 0.00 0.00 3.50 0.00 0.00 44.2 49.7 61.8 70.0 79.4 32.2 76.5 81.9 72.8 64.5 57.2 50.0 43.2 85.4 82.4 79.7 77.2 15.6 36.2 20.8 5.50 27.8 44.0 70.8 0.00 0.00 6.92 0.00 1.23 0.35 81.2 71.5 57.2 46.4 36.0 29.7 88.1 87.5 85.0 83.5 81.3 79.8 88.7 87.9 86.0 84.7 82.5 81.4 89.3 89.0 87.7 87.3 86.6 86.1;
    450 89.8 87.4 86.8 83.0 11.5 1.59 0.00 27.6 0.00 0.00 0.00 1.32 0.00 1.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 43.0 50.2 66.0 14.3 42.3 56.2 12.3 3.70 75.9 57.6 42.5 0.00 25.4 41.8 19.0 17.4 30.1 48.3 36.0 30.3 26.2 25.6 17.4 14.0 0.69 0.00 0.00 0.00 0.00 0.00 0.87 0.00 0.40 1.26 0.19 0.00 15.5 4.10 5.50 31.2 10.5 0.00 0.00 0.00 0.00 0.00 3.50 0.00 0.00 41.7 48.0 61.0 69.4 79.5 32.7 73.0 83.0 74.8 67.2 60.5 53.9 47.4 85.4 82.4 79.2 76.6 17.8 38.1 22.1 9.10 34.2 50.1 74.3 0.00 0.00 7.50 0.00 1.86 0.58 81.4 71.6 57.2 46.4 36.1 29.6 88.1 87.2 83.9 81.9 78.7 76.6 88.6 87.5 84.9 82.8 80.0 78.1 89.5 89.1 87.8 87.5 86.6 86.0;
    460 89.9 87.8 87.2 85.2 38.0 9.40 6.90 39.0 0.25 0.16 0.00 4.00 0.00 4.08 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 26.5 40.4 66.0 12.4 25.2 40.5 3.50 0.35 74.8 55.2 28.5 3.16 36.5 48.1 29.5 20.9 33.8 43.4 25.0 33.2 22.9 24.0 14.8 12.3 1.45 0.00 0.00 0.20 0.16 0.44 2.56 0.00 1.90 4.57 1.38 0.00 23.3 9.00 13.0 42.2 17.7 0.00 0.00 0.00 0.00 1.97 1.92 0.00 0.00 38.0 44.9 58.7 67.5 79.3 31.4 69.0 83.7 76.0 69.1 63.0 56.5 50.2 85.0 81.7 78.0 75.2 21.0 40.4 24.3 13.5 40.4 55.4 76.8 0.00 0.00 7.81 0.00 3.23 0.98 81.7 72.4 58.5 47.5 37.5 31.0 87.9 86.4 82.5 79.5 75.9 72.9 88.4 86.5 83.1 80.0 76.1 73.7 89.6 89.1 87.7 87.3 86.4 85.7;
    470 90.1 88.2 87.5 86.7 68.0 18.5 42.0 52.3 5.50 2.00 1.78 12.0 0.00 11.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 13.8 22.6 57.0 5.00 12.1 23.8 0.25 0.00 73.4 51.9 17.3 21.6 46.5 51.7 34.4 21.6 32.1 36.2 13.2 29.6 16.5 15.7 7.60 5.36 2.70 0.10 0.00 2.90 3.12 3.10 7.80 0.23 7.70 13.2 5.38 0.00 34.4 16.8 24.9 55.5 28.5 0.00 0.00 0.00 0.00 10.0 0.51 0.00 0.00 33.8 40.3 55.0 65.4 78.6 28.8 63.6 84.3 77.1 70.6 64.2 58.1 52.0 84.6 80.7 76.3 73.2 25.5 43.0 27.5 17.8 45.0 59.5 78.7 0.00 0.00 8.15 0.22 6.45 1.50 82.3 73.7 60.6 49.9 40.0 33.6 87.5 85.3 80.3 76.2 72.0 67.9 87.8 85.2 80.8 76.4 71.3 68.0 89.7 89.0 87.5 87.0 85.8 85.2;
    480 90.3 88.5 87.3 88.1 80.8 23.5 74.0 65.8 19.0 6.30 8.31 26.0 0.00 23.5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 8.20 40.0 0.50 2.70 9.20 0.00 0.00 71.6 48.5 10.2 44.7 53.6 52.9 35.7 20.5 27.0 28.5 4.50 22.4 9.55 6.93 2.76 1.55 4.90 0.70 0.00 13.1 13.0 13.1 21.6 1.38 21.0 30.0 15.0 0.33 46.8 24.9 36.6 68.4 41.4 0.00 0.00 0.00 0.00 17.4 0.00 0.00 0.00 27.5 35.6 51.0 62.9 77.8 25.6 57.6 84.6 77.8 71.3 65.0 59.0 53.0 84.0 79.3 74.4 70.7 30.2 45.3 30.9 21.3 48.7 62.5 80.2 0.00 0.00 8.47 0.43 14.0 2.30 82.8 74.9 62.0 52.0 42.5 35.9 87.0 84.0 77.8 72.5 67.5 62.7 87.0 83.6 77.9 72.1 65.8 61.7 89.7 89.0 87.2 86.5 85.3 84.3;
    490 90.4 88.7 86.8 88.8 85.2 25.5 84.7 76.8 41.0 13.2 20.7 43.7 0.00 39.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.63 1.85 21.0 0.00 0.20 2.30 0.00 0.00 69.5 44.6 4.00 61.4 56.8 52.2 34.5 18.0 20.2 19.6 1.30 14.1 4.27 2.14 0.40 0.10 8.50 2.14 0.00 34.2 34.5 31.9 41.7 4.90 41.5 50.8 32.0 4.00 56.6 31.3 45.1 77.6 52.1 0.00 0.00 0.00 0.00 18.0 0.00 0.00 0.00 23.5 30.9 47.1 59.8 76.7 22.2 51.3 84.9 78.3 71.8 65.7 60.0 54.0 83.3 78.0 72.1 68.1 35.8 47.2 34.3 24.5 51.2 64.6 81.2 0.00 0.00 8.60 0.39 21.6 3.50 83.2 75.8 63.5 53.9 44.8 37.9 86.2 82.4 74.2 68.3 62.3 56.6 86.0 81.8 74.4 67.0 60.0 55.0 89.7 89.0 87.0 86.0 84.4 83.4;
    500 90.5 88.9 86.3 89.5 86.9 26.3 87.5 83.5 63.5 20.0 34.5 55.0 1.50 50.8 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.12 9.56 0.00 0.00 0.33 0.00 0.00 66.7 40.2 1.33 70.2 55.8 49.8 29.7 14.4 11.1 0.36 0.17 7.30 1.58 0.46 0.00 0.00 13.3 4.47 0.00 53.4 59.0 50.5 58.8 17.7 59.0 66.0 48.4 16.6 62.1 33.7 45.8 82.7 57.9 0.00 0.00 0.00 0.00 13.0 0.00 0.00 0.00 19.5 26.5 43.5 57.0 75.5 19.3 45.2 85.3 78.6 72.6 66.4 60.8 55.4 82.6 76.6 70.2 65.7 43.5 48.9 38.3 26.8 52.8 66.0 81.7 0.00 0.00 8.73 0.15 30.7 5.20 83.3 76.6 64.6 55.2 46.1 39.4 85.2 80.5 71.2 63.8 56.7 50.1 85.0 79.7 70.5 61.7 53.7 48.1 89.6 89.0 86.5 85.2 83.5 82.3;
    510 90.6 89.1 85.5 89.9 87.8 26.7 88.5 87.0 78.0 24.3 48.8 60.0 17.3 55.2 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.51 0.00 0.00 0.00 0.00 0.00 63.9 35.8 0.35 72.4 50.9 44.8 21.5 10.1 4.39 0.00 0.00 2.64 0.48 0.00 0.00 0.00 18.2 7.24 0.10 67.0 73.0 60.6 67.9 38.8 67.7 73.0 57.2 32.3 62.9 32.4 39.7 84.6 58.8 0.00 0.00 0.00 0.96 7.35 0.00 0.30 0.10 15.8 23.4 40.0 54.2 74.2 16.8 39.4 85.4 79.0 72.9 66.5 61.1 56.2 82.0 75.3 68.3 63.5 46.3 49.2 40.7 27.9 53.4 66.4 81.9 0.00 0.00 8.85 0.00 41.4 7.70 83.0 76.1 64.0 54.5 46.0 38.5 84.4 78.6 67.7 58.7 51.0 44.5 83.8 77.5 66.7 56.5 47.7 41.6 89.6 88.7 86.0 84.4 82.4 80.8;
    520 90.7 89.3 84.8 90.3 88.4 27.0 89.1 88.4 84.1 26.7 62.0 60.2 55.0 56.5 16.0 3.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.13 0.00 0.00 0.00 0.00 0.00 60.8 31.7 0.00 70.5 42.1 36.8 11.5 5.60 1.66 0.00 0.00 0.50 0.00 0.00 0.00 0.00 23.7 10.7 0.31 69.3 79.0 63.3 70.1 52.2 69.8 75.1 59.2 40.0 59.1 27.5 29.7 84.0 55.4 0.00 0.00 0.00 7.95 3.20 0.00 9.10 5.35 13.8 20.3 36.9 51.4 73.0 14.2 34.2 85.5 79.5 73.2 67.0 61.6 57.0 81.4 74.0 66.5 61.5 47.2 48.2 40.6 28.6 53.7 66.6 82.0 0.00 0.00 8.90 0.00 51.3 10.7 82.4 74.9 61.5 51.6 42.5 35.0 83.5 77.0 64.0 54.4 46.0 38.6 82.7 75.3 63.4 52.0 42.8 36.3 89.5 88.5 85.2 83.5 81.1 79.2;
    530 90.7 89.5 84.3 90.5 89.0 27.2 89.4 89.0 86.5 28.0 76.0 57.8 77.8 55.0 52.1 22.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 57.0 27.2 0.00 64.8 30.5 26.8 3.80 2.52 0.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 28.5 14.0 0.64 65.1 79.9 61.0 67.6 53.6 67.2 73.2 55.5 39.6 51.6 20.7 17.8 82.6 47.5 0.00 0.00 0.00 14.6 0.83 0.00 13.5 1.90 11.8 17.8 34.4 49.3 72.1 12.7 30.0 86.0 80.4 74.5 68.8 62.5 59.5 81.0 73.1 65.5 59.9 48.3 48.3 40.7 30.4 55.0 67.6 82.4 0.00 0.00 9.01 0.00 59.4 15.1 81.6 73.5 59.4 48.5 38.5 31.5 82.6 75.2 61.5 50.7 41.6 34.1 81.8 73.7 60.5 48.6 39.0 31.9 89.4 88.0 84.3 82.4 79.6 77.3;
    540 90.8 89.7 84.0 90.6 89.5 27.5 89.6 89.4 87.7 28.6 83.8 54.2 86.0 51.0 70.7 48.0 0.00 2.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 52.6 22.3 0.00 55.5 18.6 16.8 0.85 0.40 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 32.1 16.6 0.89 56.7 77.5 55.0 61.8 47.6 61.5 68.5 47.5 34.5 41.3 13.7 7.90 79.1 36.6 0.00 0.00 0.00 12.9 0.14 0.00 46.1 35.0 10.5 16.6 32.7 48.1 71.5 11.8 27.1 86.5 81.5 76.0 71.0 66.1 62.7 80.8 72.7 65.0 59.1 49.6 49.2 41.6 32.5 56.5 69.0 83.0 0.00 0.00 9.07 0.00 64.2 20.2 81.2 72.5 57.8 46.5 36.6 29.4 82.1 73.9 59.5 48.3 39.0 31.3 81.3 72.5 58.6 46.6 36.7 29.8 89.2 87.5 83.4 81.0 77.7 75.0;
    550 90.8 89.9 83.9 90.7 89.8 27.8 89.8 89.7 88.4 29.0 87.0 50.0 88.4 46.0 84.3 69.5 0.00 29.0 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 48.0 17.6 0.00 44.2 8.99 8.20 0.00 0.10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 33.1 17.3 0.93 45.0 72.6 47.1 53.5 38.4 54.0 62.0 36.8 26.3 28.0 6.50 2.40 73.7 25.0 0.00 0.00 0.00 7.60 0.00 0.00 78.0 71.8 9.56 14.9 31.2 46.7 70.7 11.0 24.8 86.8 82.3 77.0 72.0 67.3 64.5 80.6 72.4 64.5 58.3 51.8 51.0 43.2 35.0 58.5 70.2 83.5 0.00 0.00 9.20 0.00 66.7 25.7 81.1 72.4 57.1 45.6 35.6 28.7 81.5 73.0 58.0 46.6 36.9 29.5 81.2 72.2 58.0 46.0 36.0 29.1 88.9 87.0 82.3 79.0 75.3 72.2;
    560 90.9 90.1 84.1 90.8 90.1 27.9 90.0 89.9 88.8 29.3 88.3 44.8 89.4 39.2 87.5 79.5 0.00 65.0 19.0 0.00 0.00 0.00 0.00 0.00 0.10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 42.8 12.9 0.00 32.5 3.59 2.95 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 31.0 15.4 0.62 33.1 66.1 37.3 43.3 27.8 45.0 54.4 25.2 17.3 16.2 1.66 0.32 67.1 14.2 0.00 0.00 2.24 3.06 0.00 0.00 75.8 63.1 8.53 13.2 29.4 45.0 69.8 9.76 23.5 87.0 82.6 77.6 72.5 68.0 65.3 80.4 71.8 63.8 57.2 56.5 55.8 47.1 41.2 63.0 73.0 84.6 0.00 9.00 9.30 0.00 66.3 31.0 81.4 72.8 58.0 46.9 36.4 29.7 81.4 72.7 57.5 45.9 35.9 28.6 81.5 72.8 58.3 46.5 36.7 29.7 88.5 86.1 80.5 76.7 72.7 69.0;
    570 90.9 90.2 84.8 90.9 90.4 28.0 90.2 90.1 89.2 29.5 88.8 38.9 89.7 32.0 88.7 84.0 0.00 80.6 60.0 11.0 0.00 0.00 0.00 0.00 10.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 37.0 8.78 0.00 20.3 0.80 0.91 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 25.6 11.4 0.21 20.7 58.0 26.5 31.6 17.4 35.0 44.5 14.4 9.70 7.95 0.40 0.00 58.8 5.50 0.00 0.00 5.97 0.83 0.00 0.00 8.00 0.00 7.77 12.1 28.0 43.6 69.0 8.81 22.6 87.1 82.7 77.8 72.7 68.3 65.8 80.2 71.5 63.2 56.2 65.0 64.5 56.0 53.0 70.9 78.1 86.8 0.00 30.5 9.20 0.00 63.0 35.6 82.5 74.5 60.6 49.2 39.2 32.7 81.4 73.0 57.9 46.3 36.1 28.7 82.5 74.6 60.5 49.8 40.2 32.3 88.0 85.0 78.5 74.0 69.3 65.0;
    580 90.9 90.3 86.0 90.9 90.6 28.4 90.4 90.3 89.5 29.6 89.1 33.1 90.1 25.1 89.3 86.3 0.00 85.4 81.0 47.0 4.55 0.00 0.00 0.00 45.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 30.6 5.65 0.00 9.56 0.00 0.10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 19.1 6.90 0.00 9.00 46.1 16.6 19.4 9.00 24.0 33.0 6.30 4.40 3.10 0.00 0.00 47.2 1.40 0.00 0.00 4.56 0.12 0.00 0.00 1.00 0.00 7.41 11.6 27.5 43.1 68.8 8.50 22.6 87.1 82.8 78.0 73.0 68.5 66.0 80.2 71.5 63.2 56.1 75.5 75.0 68.1 67.5 79.0 84.0 88.9 0.00 34.3 9.19 0.00 58.0 43.2 83.9 77.3 65.0 54.6 45.0 38.6 81.9 73.9 59.3 47.9 37.8 30.5 84.0 77.3 64.9 55.6 46.2 39.0 87.5 83.8 76.1 70.9 65.4 60.5;
    590 91.0 90.4 87.4 91.0 90.7 29.0 90.6 90.5 89.8 29.8 89.3 27.6 90.3 18.2 89.7 87.8 0.00 87.3 87.0 69.6 37.3 12.6 2.90 0.00 76.0 0.63 0.00 0.00 0.00 0.00 0.00 0.00 25.5 3.48 0.00 3.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 12.6 3.60 0.00 2.70 33.8 8.69 9.70 3.50 14.0 22.0 1.82 1.66 0.80 0.00 0.00 34.5 0.28 0.00 1.26 2.00 0.00 0.00 0.00 0.32 0.00 6.93 11.1 27.0 42.9 68.6 8.29 23.2 87.4 83.1 78.2 74.0 69.5 66.5 80.3 71.7 63.4 56.0 83.0 83.0 78.1 76.5 85.2 87.5 89.9 0.00 25.2 9.54 0.00 51.9 53.8 85.7 80.7 71.0 61.8 53.5 47.6 82.7 75.1 61.6 50.3 40.8 33.4 85.8 80.8 70.6 63.3 54.9 48.7 87.0 82.5 73.9 67.5 61.6 55.8;
    600 91.0 90.5 88.5 91.1 90.8 29.5 90.8 90.6 90.1 29.9 89.5 22.7 90.4 13.5 90.0 89.0 0.00 88.1 88.5 82.7 72.3 50.0 30.0 0.00 87.4 26.0 6.04 0.00 0.00 0.00 0.00 0.00 20.9 2.09 0.00 1.10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 7.78 1.41 0.00 0.40 24.0 3.70 4.50 1.50 7.95 14.6 0.48 0.38 0.00 0.00 0.00 24.4 0.00 0.00 5.89 0.56 0.00 0.00 0.00 16.2 1.60 6.45 10.4 26.0 41.8 68.0 7.56 23.7 87.6 84.0 79.1 75.6 72.0 68.1 80.2 71.5 63.0 55.0 87.3 87.2 85.0 85.0 88.1 89.3 90.6 0.00 16.1 9.64 0.00 45.2 65.6 87.6 84.0 77.0 70.5 64.0 58.7 83.4 76.3 63.5 53.0 43.5 36.0 88.0 84.5 77.1 71.6 64.9 59.9 86.4 81.0 71.2 64.1 57.6 51.3;
    610 91.0 90.5 89.5 91.2 90.9 29.5 90.9 90.7 90.3 29.6 89.7 19.0 90.5 9.60 90.1 89.6 0.00 88.7 89.0 85.8 82.9 75.0 63.2 10.0 89.5 67.2 41.0 0.80 0.00 0.13 0.00 0.00 16.8 1.15 0.00 0.32 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.17 0.40 0.00 0.00 18.7 1.60 2.00 0.41 4.90 10.5 0.10 0.00 0.00 0.00 0.00 18.5 0.00 0.00 5.25 0.10 0.00 0.00 0.00 52.1 32.1 5.50 9.20 24.1 40.0 66.7 6.45 23.2 88.1 85.0 81.0 78.5 75.0 71.6 79.3 70.3 61.5 53.0 89.3 88.9 88.0 88.1 89.8 90.3 91.0 0.00 11.3 9.73 0.00 37.8 77.0 89.0 86.5 82.0 77.8 73.0 69.8 83.6 76.7 64.5 54.6 44.7 37.8 89.3 87.0 82.2 79.2 74.9 70.7 85.5 79.5 68.5 60.4 53.4 46.2;
    620 91.0 90.6 90.2 91.3 91.0 29.3 91.0 90.8 90.5 29.4 89.8 14.9 90.7 6.40 90.2 90.0 0.00 89.0 89.5 87.2 86.4 82.6 78.9 45.3 90.2 84.0 75.0 24.9 0.00 1.00 0.00 0.00 12.9 0.59 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.34 0.15 0.00 0.00 13.2 0.49 0.87 0.00 2.70 6.92 0.00 0.00 0.00 0.00 0.00 13.7 0.00 0.00 2.88 0.00 0.00 0.00 0.00 83.0 78.0 4.80 7.70 21.8 37.6 65.0 5.13 21.0 88.8 86.1 83.1 80.8 78.0 74.7 78.4 68.5 59.0 50.2 90.4 90.0 89.6 89.6 90.5 90.7 91.1 0.00 7.40 9.56 0.00 30.5 82.8 90.0 88.9 86.2 83.5 80.8 78.3 83.5 76.5 64.3 54.4 44.3 37.5 90.2 88.9 86.1 84.1 81.4 79.2 84.5 77.9 65.5 56.7 49.2 42.3;
    630 91.0 90.6 90.6 91.3 91.0 29.1 91.1 90.9 90.7 29.1 89.9 11.4 90.8 3.66 90.3 90.2 0.00 89.5 89.8 87.9 87.8 85.5 84.0 71.4 90.5 88.1 86.1 60.8 0.00 6.30 0.00 0.00 10.0 0.28 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.38 0.00 0.00 0.00 7.22 0.00 0.22 0.00 1.00 3.16 0.00 0.00 0.00 0.00 0.00 7.70 0.00 0.00 1.26 0.00 0.00 0.00 0.00 84.9 79.5 3.94 6.50 19.7 35.5 63.8 4.17 18.2 89.2 87.0 84.2 82.1 79.8 77.0 77.5 66.9 56.9 47.4 90.8 90.5 90.3 90.4 90.8 90.9 91.2 0.00 2.91 9.27 0.00 25.0 86.0 90.6 89.9 88.1 87.2 85.1 84.2 83.2 76.6 63.1 53.2 42.5 35.6 90.6 90.0 88.7 87.4 86.0 84.5 83.8 76.3 62.7 53.1 45.0 38.0;
    640 91.1 90.7 90.8 91.4 91.1 29.0 91.2 91.0 90.9 28.8 90.0 9.10 90.9 2.20 90.4 90.3 0.00 89.9 90.0 88.5 88.5 86.7 86.1 82.7 90.7 89.8 89.0 78.0 0.40 22.0 0.00 0.00 7.79 0.13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.80 0.00 0.00 0.00 3.02 0.00 0.00 0.00 0.17 1.07 0.00 0.00 0.00 0.00 0.00 3.00 0.00 0.00 0.48 0.00 0.00 0.00 0.00 88.1 86.5 3.46 5.60 18.6 34.2 63.0 3.47 15.8 89.4 87.4 85.1 83.0 80.8 78.4 76.8 65.5 55.0 45.2 91.0 90.7 90.7 90.7 91.1 91.1 91.3 0.00 0.76 9.10 0.00 20.6 87.6 91.1 90.5 89.8 89.2 88.0 87.5 82.8 74.5 61.6 51.5 40.2 33.7 90.8 90.5 90.0 89.3 88.7 87.6 83.3 75.1 60.8 50.4 42.0 34.9;
    650 91.1 90.7 91.0 91.4 91.2 29.4 91.3 91.1 91.0 28.9 90.1 8.05 91.0 1.58 90.5 90.4 0.00 90.2 90.1 89.0 89.0 87.6 87.2 86.6 90.8 90.2 90.0 85.0 4.00 45.0 0.10 0.00 6.68 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.54 0.00 0.00 0.00 1.48 0.00 0.00 0.00 0.00 0.50 0.00 0.00 0.00 0.00 0.00 1.46 0.00 0.63 0.14 0.00 0.00 0.00 0.00 89.8 89.2 3.24 5.50 18.4 33.6 62.7 3.16 14.5 89.5 87.7 85.6 83.5 81.5 79.2 76.5 64.8 54.1 44.1 91.1 90.9 90.9 91.0 91.2 91.2 91.4 0.00 0.29 9.07 0.00 17.5 88.7 91.2 90.8 90.5 90.3 89.5 89.2 82.5 74.0 60.6 50.3 39.0 32.4 91.0 90.7 90.5 90.2 90.0 89.7 82.8 74.4 59.5 48.8 40.2 32.9;
    660 91.1 90.8 91.1 91.5 91.3 29.6 91.4 91.2 91.1 29.2 90.1 7.50 91.1 1.74 90.6 90.5 0.00 90.4 90.2 89.4 89.3 88.2 88.1 88.4 90.9 90.4 90.6 87.5 20.7 65.0 3.00 0.21 6.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.36 0.00 0.00 0.00 1.91 0.00 0.00 0.00 0.00 0.91 0.00 0.00 0.00 0.00 0.00 1.91 0.00 10.5 0.00 0.00 0.00 0.00 0.00 89.8 89.0 3.16 5.60 18.5 34.0 63.0 3.09 13.8 89.8 88.0 86.0 84.1 82.1 80.1 76.2 64.6 53.7 43.6 91.3 91.0 91.0 91.1 91.3 91.3 91.5 0.00 0.41 9.00 0.00 15.2 89.5 91.3 91.1 90.8 90.6 90.4 90.1 82.5 73.8 60.1 49.5 38.4 31.8 91.1 91.0 90.8 90.8 90.4 90.4 82.5 74.0 58.5 48.0 39.4 32.0;
    670 91.1 90.8 91.1 91.5 91.4 29.8 91.5 91.2 91.2 29.4 90.2 7.05 91.2 2.62 90.6 90.6 0.00 90.5 90.3 89.6 89.7 88.5 88.5 89.4 91.0 90.5 90.7 88.7 45.2 77.3 19.0 7.50 5.91 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.27 0.00 0.00 0.00 7.95 0.00 0.00 0.00 0.63 3.00 0.00 0.00 0.00 0.00 0.00 6.17 0.00 35.0 0.00 0.00 0.00 0.00 0.00 85.5 79.5 3.39 5.80 18.7 34.1 63.3 3.16 13.4 90.0 88.2 86.5 84.8 83.0 80.9 76.1 64.5 53.7 43.5 91.5 91.0 91.2 91.2 91.4 91.4 91.6 0.00 2.30 9.13 0.00 13.7 90.0 91.4 91.3 91.0 90.9 90.8 90.7 82.3 73.3 59.6 49.0 37.7 31.0 91.2 91.2 91.0 91.0 90.7 90.7 82.4 73.6 57.9 47.2 38.6 31.0;
    680 91.1 90.9 91.1 91.6 91.5 30.0 91.5 91.3 91.3 29.5 90.2 6.50 91.2 3.55 90.7 90.7 0.00 90.5 90.4 89.8 89.9 89.0 88.9 90.0 91.1 90.7 90.8 89.4 66.5 85.0 43.5 29.0 5.41 0.00 0.50 0.80 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.23 0.00 0.00 0.66 23.0 0.00 0.16 0.00 4.00 10.0 0.00 0.00 0.00 0.00 0.00 19.9 0.00 55.2 0.00 0.00 0.00 0.00 0.00 76.1 62.5 3.45 6.10 19.0 34.5 63.4 3.16 12.7 90.1 88.5 87.0 85.5 83.7 81.8 76.1 64.4 53.5 43.1 91.5 91.0 91.3 91.3 91.4 91.5 91.6 0.10 9.52 9.08 0.44 12.8 90.5 91.6 91.5 91.2 91.1 91.0 91.1 82.0 72.8 58.6 48.2 36.5 30.0 91.3 91.3 91.2 91.1 91.0 91.0 82.0 73.0 57.5 46.0 37.5 29.9;
    690 91.1 90.9 91.1 91.7 91.6 30.2 91.6 91.4 91.4 29.7 90.3 6.10 91.2 4.48 90.7 90.8 0.25 90.6 90.5 90.0 90.2 89.3 89.2 90.3 91.1 90.8 90.9 89.8 78.8 88.2 66.0 55.0 4.90 0.00 4.06 6.99 0.18 0.00 0.00 0.20 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.19 0.00 0.00 6.90 44.1 0.00 1.15 0.00 12.0 20.0 2.10 0.00 0.10 0.00 0.20 42.6 0.00 70.0 0.00 0.00 0.00 0.00 0.13 75.0 62.4 3.51 6.10 19.3 34.8 63.6 3.16 11.7 90.3 89.0 87.5 86.1 84.6 82.9 76.2 64.2 53.4 42.8 91.5 91.0 91.3 91.3 91.5 91.6 91.6 1.58 28.5 9.21 5.02 12.1 90.8 91.7 91.7 91.4 91.4 91.3 91.1 81.9 72.5 58.1 47.2 35.4 29.0 91.4 91.4 91.4 91.3 91.3 91.3 82.0 72.8 57.4 45.5 36.7 29.1;
    700 91.1 91.0 91.1 91.8 91.7 31.0 91.6 91.3 91.5 30.2 90.3 6.20 91.3 5.25 90.8 90.8 1.20 90.6 90.6 90.2 90.3 89.5 89.5 90.4 91.1 91.0 91.0 90.0 85.0 89.8 77.7 71.3 5.00 0.00 17.8 23.5 1.60 0.00 1.00 2.24 0.85 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.17 0.00 0.00 27.8 64.8 0.00 3.17 0.53 22.6 30.0 8.70 0.00 4.50 0.00 2.18 63.1 0.40 79.0 0.00 0.00 0.00 0.14 1.24 86.5 83.0 3.90 6.50 20.2 36.0 65.0 3.31 11.5 90.5 89.2 88.0 86.8 85.5 84.0 77.1 64.6 54.1 43.5 91.5 91.0 91.3 91.3 91.5 91.6 91.6 11.2 51.9 9.52 18.7 12.0 91.0 91.9 91.9 91.5 91.5 91.4 91.2 82.2 73.0 58.5 47.5 35.6 29.0 91.5 91.5 91.5 91.5 91.5 91.5 82.5 74.0 58.5 46.4 37.3 29.8]

# Wratten IR filters
"""
`const WrattenIR`

Array containing the 4 Wratten filters for the IR spectrum
"""
const WrattenIR =
   [700 0.00 0.00 0.00 11.2;
    710 0.00 0.00 0.00 32.4;
    720 0.00 0.00 0.00 57.6;
    730 0.00 0.00 7.40 69.1;
    740 0.10 0.00 32.8 77.6;
    750 2.19 0.00 56.3 83.1;
    760 7.95 0.00 69.2 85.0;
    770 17.4 0.00 74.2 86.1;
    780 31.6 0.00 77.6 87.0;
    790 43.7 0.00 79.7 87.7;
    800 53.8 0.32 81.4 88.1;
    810 61.7 3.20 82.6 88.4;
    820 69.2 8.90 83.7 88.6;
    830 74.1 17.8 84.7 88.8;
    840 77.7 28.2 85.5 89.0;
    850 81.4 41.0 86.1 89.2;
    860 84.0 53.8 86.6 89.4;
    870 85.4 61.6 87.2 89.6;
    880 86.8 69.2 87.5 89.8;
    890 87.8 74.1 87.8 89.9;
    900 88.4 78.5 88.0 90.0;
    910 88.8 81.5 88.2 90.1;
    920 89.1 83.6 88.4 90.2;
    930 89.1 85.1 88.6 90.3;
    940 89.1 86.0 88.8 90.4;
    950 89.1 87.0 89.0 90.5]


text:
export CPU, simulate!
export simcallback

mutable struct CPU
   accumulator::Int
   pc::Int
end

"""
    simcallback(cpu::CPU, instruction::Integer)
A callback function for `simulate!` to visualize
instructions executed.
"""
function simcallback(cpu::CPU, instruction::Integer)
    println(cpu.pc, ":", instruction,
            "\tAccu: ", rpad(cpu.accumulator, 3),
            "// ", disassemble(instruction))
end

function donothing(cpu::CPU, instruction::Integer) end

"""
    simulate!(mem = [0], inputs = Int[], callback = donothing)
Executes a program stored in `mem`, where each entry in an integer array is a numeric code
for the little man computer. The callback is a function with the following signature:

    callback(cpu::CPU, instruction::Integer)

It tells what the current state of the virtual LMC CPU is before the `instruction`
is executed.
"""
function simulate!(mem::Vector{Int}=[0], inputs::Vector{Int}=Int[]; callback=donothing)
    cpu = CPU(0, 0)
    outputs = Int[]
    j = 1 # keep track of what input data should be read

    address(IR::Int) = rem(IR, 100)
    data(IR::Int) = mem[begin + address(IR)]

    # limit to 1000 instructions to avoid getting trapped
    for i in 1:50
        IR = mem[begin + cpu.pc]
        callback(cpu, IR)
        cpu.pc += 1

        opcode = div(IR, 100)
        if     opcode == 1
            cpu.accumulator += data(IR)
        elseif opcode == 2
            cpu.accumulator -= data(IR)
        elseif opcode == 3
            mem[begin + address(IR)] = cpu.accumulator
        elseif opcode == 5
            cpu.accumulator = data(IR)
        elseif opcode == 6
            cpu.pc = address(IR)
        elseif opcode == 7
            if cpu.accumulator == 0
                cpu.pc = address(IR)
            end
        elseif opcode == 8
            if cpu.accumulator >= 0
                cpu.pc = address(IR)
            end
        elseif opcode == 0
            break
        elseif IR == 901
            if j > length(inputs)
                break
            end
            cpu.accumulator = inputs[j]
            j += 1
        elseif IR == 902
            push!(outputs, cpu.accumulator)
        else
            error("Unknown instruction found on address $(cpu.pc)")
        end
    end
    outputs
end

ext: jl


text:
function DSSTThirdBody(arg0::CelestialBody, arg1::jdouble)
    return DSSTThirdBody((CelestialBody, jdouble), arg0, arg1)
end

function get_body(obj::DSSTThirdBody)
    return jcall(obj, "getBody", CelestialBody, ())
end

function get_events_detectors(obj::DSSTThirdBody)
    return jcall(obj, "getEventsDetectors", Vector{EventDetector}, ())
end

function get_field_events_detectors(obj::DSSTThirdBody, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Vector{FieldEventDetector}, (Field,), arg0)
end

function get_mean_element_rate(obj::DSSTThirdBody, arg0::FieldSpacecraftState, arg1::FieldAuxiliaryElements, arg2::Vector{RealFieldElement})
    return jcall(obj, "getMeanElementRate", Vector{RealFieldElement}, (FieldSpacecraftState, FieldAuxiliaryElements, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function get_mean_element_rate(obj::DSSTThirdBody, arg0::SpacecraftState, arg1::AuxiliaryElements, arg2::Vector{jdouble})
    return jcall(obj, "getMeanElementRate", Vector{jdouble}, (SpacecraftState, AuxiliaryElements, Vector{jdouble}), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::DSSTThirdBody)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function initialize(obj::DSSTThirdBody, arg0::AuxiliaryElements, arg1::PropagationType, arg2::Vector{jdouble})
    return jcall(obj, "initialize", List, (AuxiliaryElements, PropagationType, Vector{jdouble}), arg0, arg1, arg2)
end

function initialize(obj::DSSTThirdBody, arg0::FieldAuxiliaryElements, arg1::PropagationType, arg2::Vector{RealFieldElement})
    return jcall(obj, "initialize", List, (FieldAuxiliaryElements, PropagationType, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function register_attitude_provider(obj::DSSTThirdBody, arg0::AttitudeProvider)
    return jcall(obj, "registerAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function update_short_period_terms(obj::DSSTThirdBody, arg0::Vector{RealFieldElement}, arg1::Vector{FieldSpacecraftState})
    return jcall(obj, "updateShortPeriodTerms", void, (Vector{RealFieldElement}, Vector{FieldSpacecraftState}), arg0, arg1)
end

function update_short_period_terms(obj::DSSTThirdBody, arg0::Vector{jdouble}, arg1::Vector{SpacecraftState})
    return jcall(obj, "updateShortPeriodTerms", void, (Vector{jdouble}, Vector{SpacecraftState}), arg0, arg1)
end


text:
"""
Description

      The kplot() method takes a TimeSeries.TimeArray object and plots a kernel density plot using the Winston plotting package.

Usage

      using MarketData
      kplot(cl)
      kplot(cl, caption="Apple Closing prices")
      kplot(ohlc, 4) # the fourth column ("Close") will be transformed and plotted

Method Signature(s)

      kplot(t::TimeSeries.TimeArray;
            n::Int=1,
            LHS=colorant"goldenrod",
            RHS=colorant"khaki",
            center::Float64=0.0,
            title::ASCIIString="")

Details

      The kplot() method plots a kernel density using the Winston plotting package. The plot is colored around a defined
      center, which defaults to 0.0. Left- (LHS) and right-side (RHS) colors default to "red" and "green" respectively. The first
      argument is a TimeSeries.TimeArray, while the second argument defines the column to be plotted and defaults to n=1. When the meta
      field is occupied with a string, it, along with the column name, becomes the plot's caption. Since the caption is a keyword
      argument, it can be selected be any string with the function call.

References

      https://github.com/nolta/Winston.jl

See Also

      uhist(x) which plots time series in REPL
"""

function kplot(t::TimeSeries.TimeArray; n::Int=1,
               LHS=colorant"goldenrod",
               RHS=colorant"khaki",
               center::Float64=0.0, title::ASCIIString="")

    # initialize Winston
    # KernelDensity.Winston_init()

    v = t.values[:,n]
    c = t.colnames[n]
    isa(t.meta, AbstractString) ? m = t.meta :  m = ""

    # title auto-fill
    if title == ""
        title = string(m, " ", c)
    end

    k  = kde(v)

    x1 = k.x[k.x .<= center]
    x2 = k.x[k.x .> center]
    y1 = k.density[1:length(x1)]
    y2 = k.density[length(x1)+1:end]

    f1 = FillBelow(x1, y1, color=LHS)
    f2 = FillBelow(x2, y2, color=RHS)

    #clean up below
    x3 = linspace(x1[1], last(x2))
    y3 = zeros(x3)
    f3 = FillBelow(x3, y3)

    # plotting
    p = Winston.plot(title=title)
    add(p, f1, f2, f3)

    # add some statistics
    pos, neg = length(findall(t .>= 0)), length(findall(t .< 0))
    Winston.text(minimum(x1), min(maximum(y1), maximum(y2)), "Positive = $pos \nNegative = $neg")
end


text:
using REPL.TerminalMenus
!endswith(pwd(), "test") && cd("test")
include("generalize_calls.jl")

is_loaded(foo) = isdefined(Main, foo)

function pick_modules(;verb_vec = verb_vec)
    !isdefined(Main, :Spotify) && return "-> using Spotify"
    submods = [:Albums, :Artists, :Browse, :Episodes,
        :Follow, :Library, #:Markets,#
        :Personalization, :Player, :Playlists,
        :Search, :Shows, :Tracks,
        :UsersProfile, :Objects]

    foodi = Dict{Int, UnitRange}(
            1 => 1:2,
            2 => 3:6,
            3 => 7:13,
            4 => 14:14,
            5 => 15:21,
            6 => 22:33,
            7 => 34:35,
            8 => 36:37,
            9 => 38:40)
    menuitems = String[]
    for ke in sort(collect(keys(foodi)))
        va = foodi[ke]
        str = "$va \t $(submods[ke])"
        if is_loaded(verb_vec[first(va)])
            push!(menuitems, str )
        else
            # Red color
            push!(menuitems, "\e[31m" * str * "\e[39m")
        end
    end

    menu = MultiSelectMenu(menuitems; pagesize = length(menuitems),charset = :unicode)
    selset = request("Pick group (red: module not loaded)", menu)
    feedback = String[]
    foos = []
    for i in selset
        addfoos = collect(foodi[i])
        if !is_loaded(verb_vec[first(addfoos)])
            m = submods[i]
            push!(feedback, " Spotify.$m")
        end
        foos = vcat(foos, addfoos)
    end
    if feedback != String[]
        fb = join(feedback, ", ")
        @error "Some submodules are not in scope. Try:\n using $fb"
        return Int[]
    else
        sort(foos)
    end
end
function pick_function(;verb_vec = verb_vec)
    subset = pick_modules(;verb_vec)
    if subset == Int[]
        return
    end
    list = String[]
    for i in subset
        buf = IOBuffer()
        show_default_call(IOContext(buf, :color=>true)   , i)
        push!(list, lpad(string(i), 2) * " " * String(take!(buf)))
    end
    menu = RadioMenu(list; charset = :unicode, pagesize= length(list))
    listno =1
    while true
        listno = request("Pick function\n[press: enter = select, q=quit]", menu;cursor = listno)
        listno < 1 && break
        no = subset[listno]
        println()
        print("Now calling:\n\t  ")
        show_default_call(no)
        t0 = now()
        r = make_default_call(no)
        t1 = now()
        display(r)
        println()
        println("Response after $(t1-t0)")
        println()
        r == () && break
    end
end


text:
"""
    abstract type Grid{Dc,Dp} <: Triangulation{Dc,Dp}

Abstract type that represents conforming triangulations, whose cell-wise
nodal coordinates are defined with a vector of nodal coordinates, plus
a cell-wise vector of node ids.

The interface of `Grid` is defined by overloading the
methods in `Triangulation` plus the following ones:

- [`get_node_coordinates(trian::Grid)`](@ref)
- [`get_cell_node_ids(trian::Grid)`](@ref)

From these two methods a default implementation of [`get_cell_coordinates(trian::Triangulation)`](@ref)
is available.

The `Grid`  interface has the following traits

- [`OrientationStyle(::Type{<:Grid})`](@ref)
- [`RegularityStyle(::Type{<:Grid})`](@ref)

The interface of `Grid` is tested with

- [`test_grid`](@ref)

"""
abstract type Grid{Dc,Dp} <: Triangulation{Dc,Dp} end

# Traits

abstract type OrientationStyle end
struct Oriented <: OrientationStyle end
struct NonOriented <: OrientationStyle end

"""
    OrientationStyle(::Type{<:Grid})
    OrientationStyle(::Grid)

`Oriented()` if has oriented faces, `NonOriented()` otherwise (default).
"""
OrientationStyle(a::Grid) = OrientationStyle(typeof(a))
OrientationStyle(::Type{<:Grid}) = NonOriented()

abstract type RegularityStyle end
struct Regular <: RegularityStyle end
struct Irregular <: RegularityStyle end

"""
    RegularityStyle(::Type{<:Grid})
    RegularityStyle(::Grid)

`Regular()` if no hanging-nodes default), `Irregular()` otherwise.
"""
RegularityStyle(a::Grid) = RegularityStyle(typeof(a))
RegularityStyle(::Type{<:Grid}) = Regular()

# Interface

"""
    get_node_coordinates(trian::Grid) -> AbstractArray{<:Point{Dp}}
"""
function get_node_coordinates(trian::Grid)
  @abstractmethod
end

"""
    get_cell_node_ids(trian::Grid)
"""
function get_cell_node_ids(trian::Grid)
  @abstractmethod
end

"""
    test_grid(trian::Grid)
"""
function test_grid(trian::Grid)
  test_triangulation(trian)
  nodes_coords = get_node_coordinates(trian)
  @test isa(nodes_coords,AbstractArray{<:Point})
  cell_node_ids = get_cell_node_ids(trian)
  @test isa(cell_node_ids,AbstractArray{<:AbstractArray{<:Integer}})
  @test num_nodes(trian) == length(nodes_coords)
  @test isa(is_oriented(trian),Bool)
  @test isa(is_regular(trian),Bool)
  @test OrientationStyle(trian) in (Oriented(), NonOriented())
  @test RegularityStyle(trian) in (Regular(), Irregular())
  @test is_oriented(trian) == (OrientationStyle(trian) == Oriented())
  @test is_regular(trian) == (RegularityStyle(trian) == Regular())
end

# Methods from triangulation

function get_cell_coordinates(trian::Grid)
  node_to_coords = get_node_coordinates(trian)
  cell_to_nodes = get_cell_node_ids(trian)
  lazy_map(Broadcasting(Reindex(node_to_coords)),cell_to_nodes)
end

# Some API

"""
    is_oriented(::Type{<:Grid}) -> Bool
    is_oriented(a::Grid) -> Bool
"""
is_oriented(a::Grid) = is_oriented(typeof(a))
is_oriented(a::Type{T}) where T<:Grid = OrientationStyle(T) == Oriented()

"""
    is_regular(::Type{<:Grid}) -> Bool
    is_regular(a::Grid) -> Bool
"""
is_regular(a::Grid) = is_regular(typeof(a))
is_regular(a::Type{T}) where T<:Grid = RegularityStyle(T) == Regular()

"""
    Grid(reffe::LagrangianRefFE)
"""
function Grid(reffe::LagrangianRefFE)
  UnstructuredGrid(reffe)
end

"""
    compute_linear_grid(reffe::LagrangianRefFE)
"""
function compute_linear_grid(reffe::LagrangianRefFE)
  @notimplementedif ! is_n_cube(get_polytope(reffe)) "linear grid only implemented for n-cubes at the moment"
  if get_order(reffe) == 0
    D = num_cell_dims(reffe)
    partition = tfill(1,Val{D}())
  else
    partition = get_orders(reffe)
  end
  pmin, pmax = get_bounding_box(get_polytope(reffe))
  desc = CartesianDescriptor(pmin,pmax,partition)
  CartesianGrid(desc)
end

"""
    compute_reference_grid(p::LagrangianRefFE, nelems::Integer)
"""
function compute_reference_grid(reffe::LagrangianRefFE, nelems::Integer)
  p = get_polytope(reffe)
  r = LagrangianRefFE(Float64,p,nelems)
  compute_linear_grid(r)
end

"""
    Grid(::Type{ReferenceFE{d}},p::Polytope) where d
"""
function Grid(::Type{ReferenceFE{d}},p::Polytope) where d
  UnstructuredGrid(ReferenceFE{d},p)
end

"""
    simplexify(grid::Grid)
"""
function simplexify(grid::Grid)
  simplexify(UnstructuredGrid(grid))
end


text:
function Gaussian2DKernel(sx::Float64,sy::Float64,t::Float64)::Matrix{Float64}
	a=cos(t)^2.0/(2.0*sx^2.0)+sin(t)^2.0/(2.0*sy^2.0)
	b=sin(2.0*t)/(2.0*sx^2.0)-sin(2.0*t)/(2.0*sy^2.0)
	c=sin(t)^2.0/(2.0*sx^2.0)+cos(t)^2.0/(2.0*sy^2.0)
	A=1.0/(2.0*pi*sx*sy)
	L=4.0*max(sx,sy)
	xx=collect(-L:1:L)
	yy=xx
	return [A*exp(-a*x^2.0-b*x*y-c*y^2.0) for x in xx, y in yy]
end

function load_cube(cube_fits,name;z=0.0,rms=0.0,centerx="median",centery="median",beam="header",zaxis="freq",line_rest=230.538e9,vel_rest=0.0,x0=-5.0,x1=5.0,y0=-5.0,y1=5.0,v0=-500.0,v1=500.0)
    f = FITS(cube_fits)
    header = read_header(f[1])
    cosmo = cosmology()
    D = luminosity_dist(cosmo, z)

    data = read(f[1])
    if size(data)[end]==1
        data=data[:, :, :, 1]
    end
    rad_Mpc = angular_diameter_dist(cosmo, z) / u"rad"
    deg2kpc = uconvert(u"kpc" / u"°", rad_Mpc)

    Nx = header["NAXIS1"]
    dx = abs(header["CDELT1"])
    b_x =header["CRVAL1"]-header["CRPIX1"]*dx
    x=LinRange(b_x+dx,b_x+dx*Nx,Nx)u"°"
    if isnothing(centerx)
        x_center=header["RA"]u"°"
    elseif centerx=="median"
        x_center=median(x)
    else
        x_center=centerx
    end
    x=x.-x_center
    x_kpc=x .*deg2kpc
    xx = ustrip(x_kpc)

    Ny = header["NAXIS2"]
    dy = abs(header["CDELT2"])
    b_y =header["CRVAL2"]-header["CRPIX2"]*dy
    y=LinRange(b_y+dy,b_y+dy*Ny,Ny)u"°"
    if isnothing(centerx)
        y_center=header["DEC"]u"°"
    elseif centery=="median"
        y_center=median(y)
    else
        y_center=centery
    end
    y=y.-y_center
    y_kpc=y .*deg2kpc
    yy = ustrip(y_kpc)
    if beam == "header"
        bmaj=(header["BMAJ"]u"°")*deg2kpc
        bmin=(header["BMIN"]u"°")*deg2kpc
        bpa=header["BPA"]
        beamr=[ustrip(bmaj),ustrip(bmin),bpa]
    else
        beamr=beam #in kpc
    end
    dx=diff(xx)[1]
    b1=round(beamr[1]/dx /2.35)
	b2=round(beamr[2]/dx /2.35)
	b_pa_rad=deg2rad(beamr[3]+90.0) #correction for north
	psf= Gaussian2DKernel(b1,b2,b_pa_rad);

    if zaxis == "freq"
        dnu=header["CDELT3"]u"Hz"
        b_nu=header["CRVAL3"]u"Hz"-header["CRPIX3"]*dnu
        Nnu=header["NAXIS3"]
        nu=LinRange(b_nu+dnu,b_nu+dnu*Nnu,Nnu)
        line_rest = (line_rest)u"Hz" / (1.0 + z)
        v =((line_rest^2.0 .- nu .^ 2.0) ./ (line_rest^2.0 .+ nu .^ 2.0) .*c_0)#u"km" / u"s"
        v=v .|> u"km"/u"s"
        vv=[ustrip(vi) for vi ∈ v]
    elseif (zaxis == "vel")
        Nv = header["NAXIS3"]
        dv = abs(header["CDELT3"])
        b_v =header["CRVAL3"]-header["CRPIX3"]*dv
        v=LinRange(b_v+dv,b_v+dv*Nv,Nv)
        v=v.-vel_rest
        #if !(header["CUNIT3"]=="km/s")
        vv=v ./ 1000.0
        #else
        #    vv=v
        #end
    else
        @error "Cannot understand the z coordinate"
    end
    #masking
    fv=findall(x -> (x .> v0) .& (x .< v1), vv)
    fx=findall(x -> (x .> x0) .& (x .< x1), xx)
    fy=findall(x -> (x .> y0) .& (x .< y1), yy)

    return Dict("name"=>name,"X"=>xx[fx],"Y"=>yy[fy],"V"=>vv[fv],"dx"=>diff(xx)[1],"dv"=>diff(vv)[1],"data"=>Float64.(data[fx,fy,fv]),"rms"=>rms,"beam"=>beamr,"BEAM"=>psf,"header"=>header)
end

F(θ,X,Y,var,;slit=0.1)= var[(Y .> -slit/(2.0*cos(θ-pi/2.0)) .+ X .*tan(θ-pi/2.0)) .& (Y .< slit/(2.0*cos(θ-pi/2.0)) .+ X .*tan(θ-pi/2.0))]

function pvd(ang,data;slit=0.1)
    slit_pix=Int(round(slit/data["dx"]))
	icx=Int(round(size(data["data"])[1]/2))
    A=zeros(size(data["data"]))
    for i in 1:size(data["data"])[3]
        A[:,:,i]=collect(imrotate(data["data"][:,:,i],deg2rad(ang),size(data["data"][:,:,1]),fill=0.0))[:,:,1]
    end
    return sum(A[icx-slit_pix:icx+slit_pix,:,:], dims=1)[1,:,:] ./(2.0*slit_pix*data["rms"])
end

function Clouds_Mass(I,data;n=2,v0=115.271,aco=4.2)
    beam_area=data["beam"][1]*data["beam"][2]*pi/(4.0*log(2.0)) #kpc^2
    pixel_area=data["dx"]^2.0
    #Ip=clouds.Fp .* clouds.Sp .* sqrt(2.0*pi) .*(pixel_area/beam_area) #jy km/s
    Mp=(3.25e7*aco*D^2.0/(n*v0)^2.0/(1.0+z)) .*I
    return Mp
end

function Outflow_Rate(M,V,R)
    Vkm_yr=V/3.168808781402895e-08
    R_km=R*3.0856775814671916e+16
    return Vkm_yr*M/R_km
end

function Outflow_Power(M,V,R)
    dotM=Outflow_Rate(M,V,R)
    return 0.5*V^2.0*dotM*6.300890659296177e+35
end
function Outflow_Momentum(M,V,R)
    dotM=Outflow_Rate(M,V,R)
    return V*dotM*6.300890659296178e+30
end


text:
#=
# Semidefinite Program Decomposition
This example illustrates how to decompose a positive semidefinite constraint with
`Chordal.jl`. The decomposed constraint can then be plugged into your favorite
semidefinite program (SDP) solver.
=#

#=
## Problem Setup
We first generate a random SDP using `generate_random_sdp(n)`, which
constructs problems of the form

$$
\begin{array}{ll}
\text{minimize} &c^Tx \\
\text{subject to} & \sum_{i=1}^m F_ix_i + G \succeq 0.
\end{array}
$$
=#

#-
using Chordal
using JuMP, SCS
using SparseArrays, LinearAlgebra, Random

rand_seed = 1234
n = 500
optimizer = SCS.Optimizer

## D is the dual variable and xstar is an optimal solution
c, F, G, xstar, D = Chordal.generate_random_sdp(n, rand_seed=rand_seed);

#=
## Standard Solve
First, we build the SDP using `JuMP` and solve it using `SCS` without chordal
decomposition.
=#
#-
## Without decomposition
m = Model(optimizer)
JuMP.set_silent(m)
@variable(m, x[1:n])
@constraint(m, sum(F[i]*x[i] for i in 1:n) + G ∈ PSDCone())
@objective(m, Min, dot(c, x))
time_non_chordal = @elapsed optimize!(m)

xv = value.(x)
pstar = dot(c,xv)
@info "termination status: $(termination_status(m))"
@info "solution status: $(primal_status(m))"
@info "Difference with true optimal: $(norm(xv - xstar))"
@info "Optimal value: $pstar"

#=
## Solving after Chordal Decomposition
Next, we use chordal decomposition tools from `Chordal.jl` to decompose the PSD
constraint into smaller blocks.
=#

## Chordal Decomposition Setup
m2 = Model(optimizer)
JuMP.set_silent(m2)
@variable(m2, y[1:n])

## build_A is a helper function to construct the SDP
## A = sum(F[i]*y[i] for i in 1:n) + G
## A + S ∈ PSDCone(), where S ⪰ 0
A = Chordal.build_A(y, F, G)
nnzA = nnz(A)
@info "There are $nnzA nonzeros; density = $(round(nnzA/n^2, digits=3))"
time_constraints = @elapsed build_constraints_lmi!(m2, A)
## NOTE: Can also call build_constraints_lmi!(m2, y, F, G)
@info "Finished adding constraints, time = $(round(time_constraints, digits=3))s"
@objective(m2, Min, dot(c,y))
@info "Finished setup"

#=
We can solve the decomposed model with any SDP solver (here we use SCS).
=#
#-
## Chordal Decomposition Solve
time_chordal = @elapsed optimize!(m2)
yv = value.(y)
pstar_chordal = dot(c,yv)

@info "termination status: $(termination_status(m2))"
@info "solution status: $(primal_status(m2))"
@info "difference with non-chordal: $(norm(yv - xv))"
@info "difference with true variable: $(norm(yv - xstar))"
@info "Time with chordal is $(round(time_chordal, digits=3))s vs $(round(time_non_chordal, digits=3))s"


text:
# define new aggregator algorithms here and specify their properties

# Direct type methods

"""
Gillespie's Direct method. `ConstantRateJump` rates and affects are stored in
tuples. Fastest for a small (total) number of `ConstantRateJump`s or
`MassActionJump`s (~10). For larger numbers of possible jumps use other
methods.

Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct Direct <: AbstractAggregatorAlgorithm end

"""
Gillespie's Direct method. `ConstantRateJump` rates and affects are stored via
`FunctionWrappers`, which is more performant than `Direct` for very large
numbers of `ConstantRateJump`s. However, for such large numbers of jump
different classes of aggregators are usually much more performant (i.e.
`SortingDirect`, `DirectCR`, `RSSA` or `RSSACR`).

Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct DirectFW <: AbstractAggregatorAlgorithm end

"""
The Composition-Rejection Direct method. Performs best relative to other methods
for systems with large numbers of jumps with special structure (for example a
linear chain of reactions, or jumps corresponding to particles hopping on a grid
or graph).

- A. Slepoy, A.P. Thompson and S.J. Plimpton, A constant-time kinetic Monte
  Carlo algorithm for simulation of large biochemical reaction networks, Journal
  of Chemical Physics, 128 (20), 205101 (2008). doi:10.1063/1.2919546

- S. Mauch and M. Stalzer, Efficient formulations for exact stochastic
  simulation of chemical systems, ACM Transactions on Computational Biology and
  Bioinformatics, 8 (1), 27-35 (2010). doi:10.1109/TCBB.2009.47
"""
struct DirectCR <: AbstractAggregatorAlgorithm end

"""
The Sorting Direct method. Often the fastest algorithm for smaller to moderate
sized systems (tens of jumps), or systems where a few jumps occur much more
frequently than others.

J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson and N. F. Samatova, The
  sorting direct method for stochastic simulation of biochemical systems with
  varying reaction execution behavior, Computational Biology and Chemistry, 30
  (1), 39049 (2006). doi:10.1016/j.compbiolchem.2005.10.007
"""
struct SortingDirect <: AbstractAggregatorAlgorithm end

"""
The Rejection SSA method. One of the best methods for systems with hundreds to
many thousands of jumps (along with `RSSACR`) and sparse dependency graphs.

- V. H. Thanh, C. Priami and R. Zunino, Efficient rejection-based simulation of
  biochemical reactions with stochastic noise and delays, Journal of Chemical
  Physics, 141 (13), 134116 (2014). doi:10.1063/1.4896985

- V. H. Thanh, R. Zunino and C. Priami, On the rejection-based algorithm for
  simulation and analysis of large-scale reaction networks, Journal of Chemical
  Physics, 142 (24), 244106 (2015). doi:10.1063/1.4922923
"""
struct RSSA <: AbstractAggregatorAlgorithm end

"""
The Rejection SSA Composition-Rejection method. Often the best performer for
systems with tens of thousands of jumps and sparse depedency graphs.

V. H. Thanh, R. Zunino, and C. Priami, Efficient Constant-Time Complexity
Algorithm for Stochastic Simulation of Large Reaction Networks, IEEE/ACM
Transactions on Computational Biology and Bioinformatics, Vol. 14, No. 3,
657-667 (2017).
"""
struct RSSACR <: AbstractAggregatorAlgorithm end

"""
A rejection-based direct method.
"""
struct RDirect <: AbstractAggregatorAlgorithm end

# NRM-based methods

"""
Gillespie's First Reaction Method. Should not be used for practical applications
due to slow performance relative to all other methods.

Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct FRM <: AbstractAggregatorAlgorithm end

"""
Gillespie's First Reaction Method with `FunctionWrappers` for handling
`ConstantRateJump`s. Should not be used for practical applications
due to slow performance relative to all other methods.

Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct FRMFW <: AbstractAggregatorAlgorithm end

"""
The Next Reaction Method. Can significantly outperform Direct for systems with
large numbers of jumps and sparse dependency graphs, but is usually slower than
one of `DirectCR`, `RSSA`, or `RSSACR` for such systems.

M. A. Gibson and J. Bruck, Efficient exact stochastic simulation of chemical
systems with many species and many channels, Journal of Physical Chemistry A,
104 (9), 1876-1889 (2000). doi:10.1021/jp993732q
"""
struct NRM <: AbstractAggregatorAlgorithm end

# spatial methods

"""
The Next Subvolume Method for spatial jump process simulations. Usually slower
than `DirectCRDirect`. Uses an indexed priority queue tree structure to
determine where on the grid/graph the next jump occurs, and then the `Direct`
method to determine which jump at the given location occurs.

Elf, Johan and Ehrenberg, M, Spontaneous separation of bi-stable biochemical
systems into spatial domains of opposite phases,Systems Biology, 2004 vol. 1(2)
pp. 230-236. doi:10.1049/sb:20045021
"""
struct NSM <: AbstractAggregatorAlgorithm end

"""
The Direct Composition-Rejection Direct method. Uses the `DirectCR` method to
determine where on the grid/graph a jump occurs, and the `Direct` method to
determine which jump occurs at the sampled location.

Constant-complexity stochastic simulation algorithm with optimal binning, Kevin
R. Sanft and Hans G. Othmer, Journal of Chemical Physics 143, 074108 (2015);
doi: 10.1063/1.4928635
"""
struct DirectCRDirect <: AbstractAggregatorAlgorithm end

const JUMP_AGGREGATORS = (Direct(),DirectFW(),DirectCR(),SortingDirect(),RSSA(),FRM(),FRMFW(),NRM(),RSSACR(),RDirect())

# For JumpProblem construction without an aggregator
struct NullAggregator <: AbstractAggregatorAlgorithm end

# true if aggregator requires a jump dependency graph
needs_depgraph(aggregator::AbstractAggregatorAlgorithm) = false
needs_depgraph(aggregator::DirectCR) = true
needs_depgraph(aggregator::SortingDirect) = true
needs_depgraph(aggregator::NRM) = true
needs_depgraph(aggregator::RDirect) = true

# true if aggregator requires a map from solution variable to dependent jumps.
# It is implicitly assumed these aggregators also require the reverse map, from
# jumps to variables they update.
needs_vartojumps_map(aggregator::AbstractAggregatorAlgorithm) = false
needs_vartojumps_map(aggregator::RSSA) = true
needs_vartojumps_map(aggregator::RSSACR) = true

is_spatial(aggregator::AbstractAggregatorAlgorithm) = false
is_spatial(aggregator::NSM) = true
is_spatial(aggregator::DirectCRDirect) = true


text:
scale_spec =
[
	(:name, AbstractString, "x"),
	(:_type, AbstractString, "linear"),
	(:domain, Union{Vector, VegaDataRef}, nothing),
	(:domainMin, Union{Number, VegaDataRef, VegaSignal}, nothing),
	(:domainMax, Union{Number, VegaDataRef, VegaSignal}, nothing),
	(:range, Union{Vector, AbstractString, VegaDataRef}, nothing), #TODO: What are Scale Range Literal docs trying to tell us?
	(:rangeMin, Union{Number, VegaDataRef, VegaSignal}, nothing),
	(:rangeMax, Union{Number, VegaDataRef, VegaSignal}, nothing),
	(:reverse, Union{Bool, VegaDataRef}, nothing),
	(:round, Bool, nothing),
	(:points, Bool, nothing),
	(:padding, Number, nothing),
	(:outerPadding, Number, nothing),
	(:bandWidth, Union{Number, VegaSignal}, nothing),
	(:sort, Union{Bool, Dict{Any, Any}}, nothing),  # TODO: docs say sort should be object, union with bool until we're sure
	(:clamp, Bool, nothing),
	(:nice, Union{Bool, AbstractString}, nothing), #TODO: Docs say string now, union until we know
	(:exponent, Number, nothing),
	(:zero, Bool, nothing)
]

primitivefactory(:VegaScale, scale_spec)

#isordinal(x::VegaScale) = x._type == :ordinal

max_stars_repo_name: jpthor/vega


text:
using Plots
plotlyjs();

include("../julia/fitgamma.jl");
include("../julia/utils.jl");

x_true = readcsv("../data/spikes.csv")[:];
isi = diff(x_true);

lambda, k = fitgamma(isi);

@printf "Estimated parameters:\n"
@printf "lambda: %.4f, k: %.4f\n" lambda k

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
module NamedVecs

export NamedVec

"""
A `NamedVec` is a kind of hybrid between a `Vector` and a `NamedTuple`. This can be useful,
for example, in mathematical modelling where collections of parameters are some times treated
as a single vector but at other times acessed individually.

`NamedVec` supports vector operations, such as `+` and `*`, that are useful in linear algebra.

At the same time, a `NamedVec` can be indexed by field names to obtain reshaped views into
parts of the vector.

Example:
```
v = NamedVec((a=[1], b=[2, 3]))
v' * v # returns 14
f((;a, b)) = 2a .+ b # This notation requires Julia 1.7
f(v) # returns [4, 5]
```
"""
struct NamedVec{T,Names,D<:AbstractVector{T},M} <: AbstractVector{T}
    data::D
    maps::NamedTuple{Names, M}
end

# Constructors

function NamedVec(data::AbstractVector, names::NTuple{N, Symbol}, maps::Tuple{Vararg{Any,N}}) where {N}
    NamedVec(data, NamedTuple{names}(maps))
end

"""
A `NamedVec` can be created from a `NamedTuple` of arrays. The arrays are copied (not aliased) into the new object.
"""
function NamedVec(::Type{T}, xs::NamedTuple{Names, <:Tuple{Any, Vararg{Any}}}) where {T<:AbstractVector, Names}
    ns = map(length, Tuple(xs)) # TODO: should look at number of elements actually used.
    lastix = cumsum(ns)
    firstix = lastix .- ns .+ 1
    data = similar(T isa UnionAll ? T{promote_type(eltype.(Tuple(xs))...)} : T, (sum(ns),))
    maps = ntuple(length(xs)) do i
        (xs[i] isa Real ? firstix[i] : firstix[i]:lastix[i]) => viewfun(xs[i])
    end
    v = NamedVec(data, Names, maps)
    for n in Names
        setproperty!(v, n, getproperty(xs, n))
    end
    v
end

"NamedVec, if the vector type is not specified, uses `Vector`."
function NamedVec(xs::NamedTuple{Names, <:Tuple{Any, Vararg{Any}}}) where {Names}
    E = promote_type(eltype.(Tuple(xs))...)
    NamedVec(Vector{E}, xs)
end

NamedVec(::Type{T}; kwargs...) where {T<:AbstractVector} = NamedVec(T, NamedTuple(kwargs))
NamedVec(; kwargs...) = NamedVec(NamedTuple(kwargs))

"""
`viewfun(x)` - Returns a function that create an object that is like `x`, from a vector
with the correct number of elements.
"""
function viewfun(x::AbstractArray)
    let sz = size(x)
        v -> reshape(v, sz)
    end
end

function viewfun(x::NamedVec)
    let m = maps(x)
        v -> NamedVec(v, m)
    end
end

viewfun(::AbstractVector{<:Real}) = identity

viewfun(::Real) = getindex

# Conversion

"`vec(v::NamedVec)` returns the underlying data vector of a NamedVec"
Base.vec(v::NamedVec) = getfield(v, :data)

"oftype(x::NamedVec, y) converts a vector y into a NamedVec with the names and indices of x"
function Base.oftype(v::NamedVec{Names, <:Tuple{Vararg{AbstractArray}}}, w::AbstractVector) where {Names}
    @boundscheck eachindex(v) == eachindex(w) || error("Size missmatch.")
    NamedVec(w, maps(v))
end

"""
`Tuple(v::NamedVec)` returns a tuple of the views into `v`.
For a tuple of the elements of the vector, use `Tuple(vec(v))`.
"""
function Base.Tuple(v::NamedVec)
    let data=vec(v)
        map(m->m[2](view(data, m[1])), maps(v))
    end
end

"`NamedTuple(v::NamedVec)` returns a `NamedTuple` of the views into `v`."
Base.NamedTuple(v::NamedVec{<:Any,Names}) where {Names} = NamedTuple{Names}(Tuple(v))

# Display

function Base.show(io::IO, ::MIME"text/plain", v::NamedVec)
    println(io, length(v), "-element NamedVec{", eltype(v), "}:")
    show(io, NamedTuple(v))
end

# Vector-like

Base.eltype(::NamedVec{T}) where {T} = T
Base.length(v::NamedVec) = length(vec(v))
Base.eachindex(v::NamedVec) = eachindex(vec(v))
Base.size(v::NamedVec) = size(vec(v))
Base.size(v::NamedVec, dim::Integer) = size(vec(v), dim)
Base.getindex(v::NamedVec, ix::Integer) = getindex(vec(v), ix)
Base.setindex!(v::NamedVec, y, ix::Integer) = setindex!(vec(v), y, ix)

Base.:*(a::Number, v::NamedVec) = NamedVec(a * vec(v), maps(v))
Base.:+(a::NamedVec{Names}, b::NamedVec{Names}) where {Names} = NamedVec(vec(a) + vec(b), maps(a))

Base.similar(a::NamedVec) = NamedVec(similar(vec(a)), maps(a))
Base.copy(a::NamedVec) = NamedVec(copy(vec(a)), maps(a))

# Vector-like broadcasting

struct NamedVecStyle <: Broadcast.AbstractArrayStyle{1}
end

Base.Broadcast.BroadcastStyle(::Type{<:NamedVec}) = NamedVecStyle()
Base.Broadcast.BroadcastStyle(a::NamedVecs.NamedVecStyle, ::Base.Broadcast.DefaultArrayStyle{0}) = a
Base.Broadcast.BroadcastStyle(a::NamedVecs.NamedVecStyle, ::Base.Broadcast.DefaultArrayStyle{1}) = a
Base.Broadcast.BroadcastStyle(::NamedVecs.NamedVecStyle, b::Base.Broadcast.DefaultArrayStyle) = b

function Base.similar(bc::Broadcast.Broadcasted{NamedVecStyle}, ::Type{T}) where {T}
    # TOOD: We should maybe check that there are not conflicting sets of names.
    v = getnv(bc)
    NamedVec(similar(vec(v), T), maps(v))
end

"getnv extracts the first NamedVec it can find in a Broadcasted object"
getnv(::Any) = nothing
getnv(x::NamedVec) = x
function getnv(bc::Broadcast.Broadcasted{NamedVecStyle})
    for a in bc.args
        y = getnv(a)
        y isa NamedVec && return y
    end
    nothing
end

# NamedTuple-like

Base.propertynames(::NamedVec{<:Any, Names}, ::Bool=false) where {Names} = Names

function Base.getproperty(v::NamedVec, s::Symbol)
    ix, f = getproperty(maps(v), s)
    f(view(vec(v), ix))
end

function Base.setproperty!(v::NamedVec, s::Symbol, y)
    ix, f = getproperty(maps(v), s)
    if f === getindex
        vec(v)[ix] = y
    else
        f(view(vec(v), ix)) .= y
    end
end

# Utility functions
# (Unexported and mainly for internal use)

maps(v::NamedVec) = getfield(v, :maps)

end # modulev


text:
function toMagma(f::IOStream, a::Array{T, 1}; name::String="R") where T
  print(f, name, " := [\n")
  for i=1:(length(a)-1)
    try
      toMagma(f, a[i]);
    catch b
      print(f, a[i])
    end
    print(f, ",\n")
  end
  try
    toMagma(f, a[end])
  catch b
    print(f, a[end])
  end
  print(f, "];\n")
end

function toMagma(f::IOStream, t::Tuple)
  print(f, "<")
  for i=1:length(t)
    try
      toMagma(f, t[i])
    catch e
      print(f, t[i])
    end
    if i<length(t)
      print(f, ", ")
    else
      print(f, ">\n")
    end
  end
end

function toMagma(s::String, a::Array{T, 1}; name::String="R", mode::String ="w") where T
  f = open(s, mode)
  toMagma(f, a, name = name)
  close(f)
end


text:
"""
    SGLD(n_iters::Int, step_size::Float64)

 Stochastic Gradient Langevin Dynamics sampler.

Usage:

```julia
SGLD(1000, 0.5)
```

Example:

```julia
@model example begin
  ...
end

sample(example, SGLD(1000, 0.5))
```
"""
mutable struct SGLD{T} <: Hamiltonian
  n_iters   ::  Int       # number of samples
  step_size ::  Float64   # constant scale factor of learning rate
  space     ::  Set{T}    # sampling space, emtpy means all
  gid       ::  Int
end
function SGLD(step_size::Float64, space...)
  SGLD(1, step_size, space..., 0)
end
function SGLD(n_iters, step_size)
  SGLD(n_iters, step_size, Set(), 0)
end
function SGLD(n_iters, step_size, space...)
  _space = isa(space, Symbol) ? Set([space]) : Set(space)
  SGLD(n_iters, step_size, _space, 0)
end
SGLD(alg::SGLD, new_gid::Int) = SGLD(alg.n_iters, alg.step_size, alg.space, new_gid)

function step(model, spl::Sampler{<:SGLD}, vi::VarInfo, is_first::Bool)
  if is_first
    if ~haskey(spl.info, :wum)
      if spl.alg.gid != 0 link!(vi, spl) end    # X -> R

      D = length(vi[spl])
      ve = VarEstimator{Float64}(0, zeros(D), zeros(D))
      wum = WarmUpManager(1, Dict(), ve)
      wum[:ϵ] = [spl.alg.step_size]
      wum[:stds] = ones(D)
      spl.info[:wum] = wum

      oldθ = realpart(vi[spl])
      vi[spl] = oldθ

      # Initialize iteration counter
      spl.info[:t] = 0

      if spl.alg.gid != 0 invlink!(vi, spl) end # R -> X
    end

    push!(spl.info[:accept_his], true)

    vi
  else
    # Update iteration counter
    t = deepcopy(spl.info[:t]) + 1
    spl.info[:t] = deepcopy(t)

    @debug "compute current step size..."
    γ = .35
    ϵ_t = spl.alg.step_size / t^γ # NOTE: Choose γ=.55 in paper
    push!(spl.info[:wum][:ϵ], ϵ_t)

    @debug "X-> R..."
    if spl.alg.gid != 0
      link!(vi, spl)
      runmodel(model, vi, spl)
    end

    @debug "recording old variables..."
    old_θ = realpart(vi[spl])
    θ = deepcopy(old_θ)
    _, grad = gradient_forward(old_θ, vi, model, spl)

    if verifygrad(grad)
      @debug "update latent variables..."
      v = zeros(Float64, size(old_θ))
      for k in 1:size(old_θ, 1)
        noise = rand(MvNormal(zeros(length(old_θ[k,:])), sqrt.(ϵ_t)*ones(length(old_θ[k,:]))))
        θ[k,:] = old_θ[k,:] - 0.5 * ϵ_t * grad[k,:] + noise
      end
    end

    @debug "always accept..."
    push!(spl.info[:accept_his], true)
    vi[spl] = θ

    @debug "R -> X..."
    if spl.alg.gid != 0 invlink!(vi, spl); cleandual!(vi) end

    vi
  end
end


text:
using VLKeggSDK
using DataFrames
using ProgressMeter
using BSON

# setup list of reactions -
rn_number_array = [

    # upper glycolysis -
    "rn:R00299" # 1
    "rn:R00771" # 2
    "rn:R00756" # 3
    "rn:R00762" # 4
    "rn:R01068" # 5
    "rn:R01015" # 6
    "rn:R01061" # 7
    "rn:R01512" # 8
    "rn:R01518" # 9
    "rn:R00658" # 10
    "rn:R00200" # 11

    # AcCoA formation -
    "rn:R00199" # 12
    "rn:R01699" # 13
    "rn:R02569" # 14
    "rn:R08550" # 15
    "rn:R00235" # 16
    "rn:R00229" # 17

    # TCA cycle -
    "rn:R00351" # 18
    "rn:R01324" # 19
    "rn:R00267" # 20
    "rn:R01700" # 21
    "rn:R02570" # 22
    "rn:R08550" # 23
    "rn:R00405" # 24
    "rn:R01082" # 25
    "rn:R00342" # 26

    # OxPhos -
    "rn:R02164" # 27
    "rn:R00081" # 28
    "rn:R00004" # 29
]

# download these reactions -
reaction_kegg_metabolite_markup_array = Array{KEGGReaction,1}()
@showprogress "Reactions: " for kegg_reaction_number in rn_number_array

    # get the reaction in KEGG metabolite format -
    kegg_reaction_object = get_reaction_for_rn_number(kegg_reaction_number) |> check
    if (isnothing(kegg_reaction_object) == false)
        push!(reaction_kegg_metabolite_markup_array, kegg_reaction_object)
    end
end

# ok: lets add some reactions by hand -
# ec_number::String
# reaction_number::String
# enzyme_name::String
# reaction_markup::String
# reaction_forward::Union{Missing,String}
# reaction_reverse::Union{Missing,String}
# stoichiometric_dictionary::Union{Missing,Dict{String,Number}}

# complex V
CV = KEGGReaction() # 30
CV.ec_number = "ec:7.1.2.1"
CV.reaction_number = "rn:74aa5235-32e2-403a-9043-5b7e59d56255"
CV.enzyme_name = "proton-translocating ATPase"
CV.reaction_markup = "C00008 + C00009 + 3 C00080 <=> C00002 + C00001 + 3 C00080"
CV.reaction_forward = "C00008 + C00009 + 3 C00080"
CV.reaction_reverse = "C00002 + C00001 + 3 C00080"
CV.stoichiometric_dictionary = Dict("C00008" => -1, "C00009" => -1, "C00002" => 1, "C00001" => 1)
push!(reaction_kegg_metabolite_markup_array, CV)

# complex III
CIII = KEGGReaction() # 31
CIII.ec_number = "ec:7.1.1.8"
CIII.reaction_number = "rn:a484c567-c328-42ea-9c39-3ebac4a587c8"
CIII.enzyme_name = "complex III"
CIII.reaction_markup = "C00530 +  2 C00125 <=> C00472 +  2 C00126 + 2 C00080"
CIII.reaction_forward = "C00530 +  2 C00125"
CIII.reaction_reverse = "C00472 +  2 C00126 + 2 C00080"
CIII.stoichiometric_dictionary = Dict("C00530" => -1, "C00125" => -2, "C00472" => 1, "C00126" => 2, "C00080" => 2)
push!(reaction_kegg_metabolite_markup_array, CIII)

# CI -
CI = KEGGReaction() # 32
CI.ec_number = "ec:7.1.1.2"
CI.reaction_number = "rn:R11945"
CI.enzyme_name = "complex I"
CI.reaction_markup = "C00399 + C00004 + 6 C00080 <=> C00390 + C00003 + 6 C00080"
CI.reaction_forward = "C00399 + C00004 + 6 C00080"
CI.reaction_reverse = "C00390 + C00003 + 6 C00080"
CI.stoichiometric_dictionary = Dict("C00399" => -1, "C00004" => -1, "C00390" => 1, "C00003" => 1)
push!(reaction_kegg_metabolite_markup_array, CI)

# build the metabolites table -
compound_record_array = Array{KEGGCompound,1}()
@showprogress "Metabolites: " for reaction_obj in reaction_kegg_metabolite_markup_array

    # get the reaction number -
    kegg_reaction_markup = reaction_obj.reaction_markup

    # get the compound symbols -
    tmp_compound_array = extract_metabolite_symbols(kegg_reaction_markup)
    compound_record = get_compound_records(tmp_compound_array) |> check
    if (isnothing(compound_record) == false)
        append!(compound_record_array, compound_record)
    end
end

# for some reason unique is not working - so ...
unique_compound_array = Array{KEGGCompound,1}()
for compound_object in compound_record_array
    if (in(compound_object, unique_compound_array) == false)
        push!(unique_compound_array, compound_object)
    end
end

# Dump model to disk -
model = Dict{Symbol,Any}()
model[:compounds] = unique!((unique_compound_array |> DataFrame))
model[:reactions] = unique!((reaction_kegg_metabolite_markup_array |> DataFrame))

# setup path -
_PATH_TO_MODEL_FILE = "./test/model/model.bson"
bson(_PATH_TO_MODEL_FILE, model)


text:
# only for the first u-domain, since we need to separate the u=0 point
function compute_bulkevolved_t_1st!(bulkevol_t::BulkEvolved,
                                    bulkconstrain::BulkConstrained,
                                    bulkevol::BulkEvolved, deriv::BulkDeriv,
                                    sys::System, evoleq::AffineNull)
    phiGF   = getphi(bulkevol)
    phiGF_t = getphi(bulkevol_t)

    phidGF = getphid(bulkconstrain)
    # SdGF   = getSd(bulkconstrain)
    AGF    = getA(bulkconstrain)

    Du_phi  = deriv.Du_phi
    Du_phid = deriv.Du_phid

    uu  = sys.ucoord
    # Du  = sys.Du
    # Dx  = sys.Dx
    # Dy  = sys.Dy

    Nu, Nx, Ny = size(sys)

    # u = 0
    @fastmath @inbounds for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            phi_u  = Du_phi[1,i,j]
            phid_u = Du_phid[1,i,j]

            phiGF_t[1,i,j] = 3//2 * (phi_u - phid_u) + phi_u / 2
        end
    end

    # remaining inner grid points
    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            @inbounds @simd for a in 2:Nu
                u      = uu[a]
                u3     = u * u * u
                u4     = u * u3

                phi    = phiGF[a,i,j]
                phid   = phidGF[a,i,j]
                A      = AGF[a,i,j]

                phi_u  = Du_phi[a,i,j]
                # phid_u = Du_phid[a,i,j]

                phiGF_t[a,i,j] =  3//2 * (phi - phid) / u + phi_u / 2 + 3//2 * u3 * A * phi +
                    A * u4 * phi_u / 2
            end
        end
    end

    nothing
end

# remaining u-domains
function compute_bulkevolved_t!(bulkevol_t::BulkEvolved,
                                bulkconstrain::BulkConstrained,
                                bulkevol::BulkEvolved, deriv::BulkDeriv,
                                sys::System, evoleq::AffineNull)
    phiGF   = getphi(bulkevol)
    phiGF_t = getphi(bulkevol_t)

    phidGF = getphid(bulkconstrain)
    # SdGF   = getSd(bulkconstrain)
    AGF    = getA(bulkconstrain)

    Du_phi  = deriv.Du_phi
    Du_phid = deriv.Du_phid

    uu  = sys.ucoord
    # Du  = sys.Du
    # Dx  = sys.Dx
    # Dy  = sys.Dy

    Nu, Nx, Ny = size(sys)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            @inbounds @simd for a in 1:Nu
                u      = uu[a]
                u3     = u * u * u
                u4     = u * u3

                phi    = phiGF[a,i,j]
                phid   = phidGF[a,i,j]
                A      = AGF[a,i,j]

                phi_u  = Du_phi[a,i,j]
                # phid_u = Du_phid[a,i,j]

                phiGF_t[a,i,j] =  3//2 * (phi - phid) / u + phi_u / 2 + 3//2 * u3 * A * phi +
                    A * u4 * phi_u / 2
            end
        end
    end

    nothing
end

function sync_bulkevolved!(bulkevols_t, bulkconstrains,
                           systems::SystemPartition, evoleq)
    Nsys = length(systems)

    @inbounds for i in 1:Nsys-1
        sync_bulkevolved!(bulkevols_t[i], bulkevols_t[i+1], bulkconstrains[i+1],
                          systems[i], systems[i+1], evoleq)
    end

    nothing
end

function sync_bulkevolved!(bulkevol1_t::BulkEvolved, bulkevol2_t::BulkEvolved,
                           bulkconstrain2::BulkConstrained,
                           sys1::System, sys2::System, evoleq::AffineNull)
    u0 = sys2.ucoord[1]

    _, Nx, Ny = size(sys2)

    phiGF1_t = getphi(bulkevol1_t)
    phiGF2_t = getphi(bulkevol2_t)
    AGF      = getA(bulkconstrain2)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            A = AGF[1,i,j]

            # characteristic speed
            c = u0 * u0 * A/2

            # if c > 0, mode is entering grid1 from grid2;
            # if c < 0, mode is entering grid2 from grid1.
            # we assume here that grids merely touch at the interface
            if c > 0
                phiGF1_t[end,i,j] = phiGF2_t[1,i,j]
            elseif c < 0
                phiGF2_t[1,i,j] = phiGF1_t[end,i,j]
            end
        end
    end

    nothing
end


text:
__precompile__()

module FastqSubsample

using Libz

export fastqsubsample

"""
    fastqsubsample(ifastq, ofastq, nSubsample [, seed = seed])

Subsample fixed number of reads from input FASTQ file.

- ifastq: Input FASTQ. File names ending with ".gz" will be treated as gzipped files.
- ofastq: Output FASTQ. File names ending with ".gz" will be treated as gzipped files.
- nSubsample: Number of reads after subsampling. Integer.
- seed: Seed for random number generation. The results of subsampling is the same when using the same seed. Integer.

"""

function fastqsubsample(ifastq::String, ofastq::String, nSubsample::Int; kwargs...)
    S = zeros(Int, nSubsample)
    n = countFastq(ifastq)

    println("Reads in $ifastq: $n")
    println("Reads will be subsampled to: $nSubsample")
    if n < nSubsample
        println("Cannot subsample more than the number of the original reads")
        exit()
    end

    isSeeded, seed = checkSeed(kwargs, :seed, 12345)
    rng = ifelse(isSeeded, MersenneTwister(seed), MersenneTwister())

    for i in 1:nSubsample
        S[i] = i
    end

    for i in (nSubsample+1):n
        j = rand(rng, 1:i)
        if j <= nSubsample
            S[j] = i
        end
    end

    # accepted = IntSet(shuffle(rng, 1:n)[1:nSubsample])
    accepted = IntSet(S)

    writeFastqSubsample(ifastq, ofastq, accepted)
    println("Subsampled FASTQ is written: $ofastq")
end

function checkSeed(kwargs::Array{Any,1}, s::Symbol, val::Int)
    for element in kwargs
        if element[1] == s
            return true, element[2]
        end
    end
    return false, val
end

function countFastq(ifastq::String)
    n = 0

    if ismatch(r".gz$", ifastq)
        f = ZlibInflateInputStream(open(ifastq))
    else
        f = open(ifastq)
    end

    try # For avoid EOFerror of GZip.open
        for line in eachline(f)
            n += 1
            readline(f)
            readline(f)
            readline(f)
        end
    end

    close(f)

    n
end

function writeFastqSubsample(ifastq::String, ofastq::String, accepted::IntSet)
    i = 0
    if ismatch(r".gz$", ifastq)
        f = ZlibInflateInputStream(open(ifastq))
    else
        f = open(ifastq)
    end

    if ismatch(r".gz$", ofastq)
        fw = ZlibDeflateOutputStream(open(ofastq, "w"))
    else
        fw = open(ofastq, "w")
    end

    try
        for line in eachline(f)
            i += 1
            if in(i, accepted)
                write(fw, line)
                write(fw, readline(f))
                write(fw, readline(f))
                write(fw, readline(f))
            else
                readline(f)
                readline(f)
                readline(f)
            end
        end
    end

    close(f)
    close(fw)
end

end # module

max_stars_repo_name: yuifu/FastqSubsample

max_issues_repo_name: yuifu/FastqSubsample

max_forks_repo_name: yuifu/FastqSubsample


text:
module Feather

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Feather" => Ahorn.EntityPlacement(
        Maple.Feather
    ),
    "Feather (Shielded)" => Ahorn.EntityPlacement(
        Maple.Feather,
        "point",
        Dict{String, Any}(
            "shielded" => true
        )
    )
)

function Ahorn.selection(entity::Maple.Feather)
    x, y = Ahorn.position(entity)
    shielded = get(entity.data, "shielded", false)

    return shielded ? Ahorn.Rectangle(x - 12, y - 12, 24, 24) : Ahorn.Rectangle(x - 8, y - 8, 16, 16)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Feather, room::Maple.Room)
    if get(entity.data, "shielded", false)
        Ahorn.Cairo.save(ctx)

        Ahorn.set_antialias(ctx, 1)
        Ahorn.set_line_width(ctx, 1);

        Ahorn.drawCircle(ctx, 0, 0, 12, (1.0, 1.0, 1.0, 1.0))

        Ahorn.Cairo.restore(ctx)
    end

    Ahorn.drawSprite(ctx, "objects/flyFeather/idle00.png", 0, 0)
end

end


text:
import InteractiveUtils

g_cached_subtypes = Dict{DataType, Vector{DataType}}()

"""
Returns an array of all concrete subtypes of T.
Note that this does not find parameterized types.
"""
function get_all_concrete_subtypes(::Type{T}) where {T}
    if haskey(g_cached_subtypes, T)
        return g_cached_subtypes[T]
    end

    sub_types = Vector{DataType}()
    _get_all_concrete_subtypes(T, sub_types)
    g_cached_subtypes[T] = sub_types
    return sub_types
end

"""Recursively builds a vector of subtypes."""
function _get_all_concrete_subtypes(::Type{T}, sub_types::Vector{DataType}) where {T}
    for sub_type in InteractiveUtils.subtypes(T)
        if isconcretetype(sub_type)
            push!(sub_types, sub_type)
        elseif isabstracttype(sub_type)
            _get_all_concrete_subtypes(sub_type, sub_types)
        end
    end

    return nothing
end

"""Returns an array of concrete types that are direct subtypes of T."""
function get_concrete_subtypes(::Type{T}) where {T}
    return [x for x in InteractiveUtils.subtypes(T) if isconcretetype(x)]
end

"""Returns an array of abstract types that are direct subtypes of T."""
function get_abstract_subtypes(::Type{T}) where {T}
    return [x for x in InteractiveUtils.subtypes(T) if isabstracttype(x)]
end

"""Returns an array of all super types of T."""
function supertypes(::Type{T}, types = []) where {T}
    super = supertype(T)
    push!(types, super)
    if super == Any
        return types
    end

    supertypes(super, types)
end

"""Converts a DataType to a Symbol, stripping off the module name(s)."""
function type_to_symbol(data_type::DataType)
    return Symbol(strip_module_name(string(data_type)))
end

"""Strips the module name off of a type."""
function strip_module_name(name::String)
    index = findfirst(".", name)
    # Account for the period being part of a parametric type.
    parametric_index = findfirst("{", name)

    if isnothing(index) ||
       (!isnothing(parametric_index) && index.start > parametric_index.start)
        basename = name
    else
        basename = name[(index.start + 1):end]
    end

    return basename
end

function strip_module_name(::Type{T}) where {T}
    return strip_module_name(string(T))
end

function strip_parametric_type(name::AbstractString)
    index = findfirst("{", name)
    if !isnothing(index)
        # Ignore the parametric type.
        name = name[1:(index.start - 1)]
    end

    return name
end

"""
Return true if all publicly exported names in mod are defined.
"""
function validate_exported_names(mod::Module)
    is_valid = true
    for name in names(mod)
        if !isdefined(mod, name)
            is_valid = false
            @error "module $mod exports $name but does not define it"
        end
    end

    return is_valid
end

"""
Recursively compares struct values. Prints all mismatched values to stdout.

# Arguments
- `x::T`: First value
- `y::T`: Second value
- `compare_uuids::Bool = false`: Compare any UUID in the object or composed objects.
"""
function compare_values(x::T, y::T; compare_uuids = false) where {T}
    match = true
    fields = fieldnames(T)
    if isempty(fields)
        match = x == y
    else
        for field_name in fields
            if T <: TimeSeriesContainer && field_name == :time_series_storage
                # This gets validated at SystemData. Don't repeat for each component.
                continue
            end
            val1 = getfield(x, field_name)
            val2 = getfield(y, field_name)
            if !isempty(fieldnames(typeof(val1)))
                if !compare_values(val1, val2, compare_uuids = compare_uuids)
                    @error "values do not match" T field_name val1 val2
                    match = false
                end
            elseif val1 isa AbstractArray
                if !compare_values(val1, val2, compare_uuids = compare_uuids)
                    @error "values do not match" T field_name val1 val2
                    match = false
                end
            else
                if val1 != val2
                    @error "values do not match" T field_name val1 val2
                    match = false
                end
            end
        end
    end

    return match
end

function compare_values(x::Vector{T}, y::Vector{T}; compare_uuids = false) where {T}
    if length(x) != length(y)
        @error "lengths do not match" T length(x) length(y)
        return false
    end

    match = true
    for i in range(1, length = length(x))
        if !compare_values(x[i], y[i], compare_uuids = compare_uuids)
            @error "values do not match" typeof(x[i]) i x[i] y[i]
            match = false
        end
    end

    return match
end

function compare_values(x::Dict, y::Dict; compare_uuids = false)
    keys_x = Set(keys(x))
    keys_y = Set(keys(y))
    if keys_x != keys_y
        @error "keys don't match" keys_x keys_y
        return false
    end

    match = true
    for key in keys_x
        if !compare_values(x[key], y[key], compare_uuids = compare_uuids)
            @error "values do not match" typeof(x[key]) key x[key] y[key]
            match = false
        end
    end

    return match
end

compare_values(::Type{T}, ::Type{T}; compare_uuids = false) where {T} = true
compare_values(::Type{T}, ::Type{U}; compare_uuids = false) where {T, U} = false

# Copied from https://discourse.julialang.org/t/encapsulating-enum-access-via-dot-syntax/11785/10
"""
Macro to wrap Enum in a module to keep the top level scope clean.

# Examples
```Julia
julia> @scoped_enum Fruit APPLE = 1 ORANGE = 2

julia> value = Fruit.APPLE
Fruit.APPLE = 1

julia> value = Fruit(1)
Fruit.APPLE = 1

julia> @scoped_enum(Fruit,
    APPLE = 1,  # comment
    ORANGE = 2,  # comment
)
```
"""
macro scoped_enum(T, args...)
    blk = esc(
        :(
            module $(Symbol("$(T)Module"))
            using JSON3
            export $T
            struct $T
                value::Int64
            end
            const NAME2VALUE =
                $(Dict(String(x.args[1]) => Int64(x.args[2]) for x in args))
            $T(str::String) = $T(NAME2VALUE[str])
            const VALUE2NAME =
                $(Dict(Int64(x.args[2]) => String(x.args[1]) for x in args))
            Base.string(e::$T) = VALUE2NAME[e.value]
            Base.getproperty(::Type{$T}, sym::Symbol) =
                haskey(NAME2VALUE, String(sym)) ? $T(String(sym)) : getfield($T, sym)
            Base.show(io::IO, e::$T) =
                print(io, string($T, ".", string(e), " = ", e.value))
            Base.propertynames(::Type{$T}) = $([x.args[1] for x in args])
            JSON3.StructType(::Type{$T}) = JSON3.StructTypes.StringType()

            Base.convert(::Type{$T}, val::Integer) = $T(val)
            Base.isless(val::$T, other::$T) = isless(val.value, other.value)
            Base.instances(::Type{$T}) = tuple($T.($(x.args[2] for x in args))...)
            end
        ),
    )
    top = Expr(:toplevel, blk)
    push!(top.args, :(using .$(Symbol("$(T)Module"))))
    return top
end

function compose_function_delegation_string(
    sender_type::String,
    sender_symbol::String,
    argid::Vector{Int},
    method::Method,
)
    s = "p" .* string.(1:(method.nargs - 1))
    s .*= "::" .* string.(fieldtype.(Ref(method.sig), 2:(method.nargs)))
    s[argid] .= "p" .* string.(argid) .* "::$sender_type"

    m = string(method.module.eval(:(parentmodule($(method.name))))) * "."
    l = "$m:(" * string(method.name) * ")(" * join(s, ", ")

    m = string(method.module) * "."
    l *= ") = $m:(" * string(method.name) * ")("
    s = "p" .* string.(1:(method.nargs - 1))

    s[argid] .= "getfield(" .* s[argid] .* ", :$sender_symbol)"
    l *= join(s, ", ") * ")"
    l = join(split(l, "#"))
    return l
end

function forward(sender::Tuple{Type, Symbol}, ::Type, method::Method)
    # Assert that function is always just one argument
    @assert method.nargs < 4 "`forward` only works for one and two argument functions"
    # Assert that function name always starts with `get_*`
    "`forward` only works for accessor methods that are defined as `get_*` or `set_*`"
    @assert startswith(string(method.name), r"set_|get_")
    sender_type = "$(parentmodule(sender[1])).$(strip_module_name(sender[1]))"
    sender_symbol = string(sender[2])
    code_array = Vector{String}()
    # Search for receiver type in method arguments
    argtype = fieldtype(method.sig, 2)
    (sender[1] == argtype) && (return code_array)
    if string(method.name)[1] == '@'
        @warn "Forwarding macros is not yet supported."
        display(method)
        println()
        return code_array
    end

    # first argument only
    push!(
        code_array,
        compose_function_delegation_string(sender_type, sender_symbol, [1], method),
    )

    tmp = split(string(method.module), ".")[1]
    code =
        "@eval " .* tmp .* " " .* code_array .* " # " .* string(method.file) .* ":" .*
        string(method.line)
    if (tmp != "Base") && (tmp != "Main")
        pushfirst!(code, "using $tmp")
    end
    code = unique(code)
    return code
end

function forward(sender::Tuple{Type, Symbol}, receiver::Type, exclusions::Vector{Symbol})
    @assert isconcretetype(sender[1])
    @assert isconcretetype(receiver)
    code = Vector{String}()
    active_methods = getfield.(InteractiveUtils.methodswith(sender[1]), :name)
    for m in InteractiveUtils.methodswith(receiver)
        m.name ∈ exclusions && continue
        m.name ∈ active_methods && continue
        if startswith(string(m.name), "get_") && m.nargs == 2
            # forwarding works for functions with 1 argument and starts with `get_`
            append!(code, forward(sender, receiver, m))
        elseif startswith(string(m.name), "set_") && m.nargs == 3
            # forwarding works for functions with 2 argument and starts with `set_`
            append!(code, forward(sender, receiver, m))
        end
    end
    return code
end

macro forward(sender, receiver, exclusions = Symbol[])
    out = quote
        list = InfrastructureSystems.forward($sender, $receiver, $exclusions)
        for line in list
            eval(Meta.parse("$line"))
        end
    end
    return esc(out)
end

"""
Return the resolution from a TimeArray.
"""
function get_resolution(ts::TimeSeries.TimeArray)
    tstamps = TimeSeries.timestamp(ts)
    timediffs = unique([tstamps[ix] - tstamps[ix - 1] for ix in 2:length(tstamps)])

    res = []

    for timediff in timediffs
        if mod(timediff, Dates.Millisecond(Dates.Day(1))) == Dates.Millisecond(0)
            push!(res, Dates.Day(timediff / Dates.Millisecond(Dates.Day(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Hour(1))) == Dates.Millisecond(0)
            push!(res, Dates.Hour(timediff / Dates.Millisecond(Dates.Hour(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Minute(1))) == Dates.Millisecond(0)
            push!(res, Dates.Minute(timediff / Dates.Millisecond(Dates.Minute(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Second(1))) == Dates.Millisecond(0)
            push!(res, Dates.Second(timediff / Dates.Millisecond(Dates.Second(1))))
        else
            throw(DataFormatError("cannot understand the resolution of the time series"))
        end
    end

    if length(res) > 1
        throw(
            DataFormatError(
                "time series has non-uniform resolution: this is currently not supported",
            ),
        )
    end

    return res[1]
end

function get_initial_timestamp(data::TimeSeries.TimeArray)
    return TimeSeries.timestamp(data)[1]
end

function get_module(module_name)
    # root_module cannot find InfrastructureSystems if it hasn't been installed by the
    # user (but has been installed as a dependency to another package).
    return module_name == "InfrastructureSystems" ? InfrastructureSystems :
           Base.root_module(Base.__toplevel__, Symbol(module_name))
end

get_type_from_strings(module_name, type) = getfield(get_module(module_name), Symbol(type))

# This function is used instead of cp given
# https://github.com/JuliaLang/julia/issues/30723
function copy_file(src::AbstractString, dst::AbstractString)
    src_path = normpath(src)
    dst_path = normpath(dst)
    if Sys.iswindows()
        run(`cmd /c copy /Y $(src_path) $(dst_path)`)
    else
        run(`cp -f $(src_path) $(dst_path)`)
    end
    return
end

function get_initial_times(
    initial_timestamp::Dates.DateTime,
    count::Int,
    interval::Dates.Period,
)
    if count == 0
        return []
    elseif interval == Dates.Second(0)
        return [initial_timestamp]
    end

    return range(initial_timestamp; length = count, step = interval)
end

function get_total_period(
    initial_timestamp::Dates.DateTime,
    count::Int,
    interval::Dates.Period,
    horizon::Int,
    resolution::Dates.Period,
)
    last_it = initial_timestamp + interval * count
    last_timestamp = last_it + resolution * (horizon - 1)
    return last_timestamp - initial_timestamp
end

function transform_array_for_hdf(data::SortedDict{Dates.DateTime, Vector{CONSTANT}})
    return hcat(values(data)...)
end

function transform_array_for_hdf(data::Vector{<:Real})
    return data
end

function transform_array_for_hdf(data::SortedDict{Dates.DateTime, Vector{POLYNOMIAL}})
    lin_cost = hcat(values(data)...)
    rows, cols = size(lin_cost)
    @assert_op length(first(lin_cost)) == 2
    t_lin_cost = Array{Float64}(undef, rows, cols, 2)
    for r in 1:rows, c in 1:cols
        tuple = lin_cost[r, c]
        for (i, value) in enumerate(tuple)
            t_lin_cost[r, c, i] = value
        end
    end
    return t_lin_cost
end

function transform_array_for_hdf(data::Vector{POLYNOMIAL})
    rows = length(data)
    @assert_op length(first(data)) == 2
    t_lin_cost = Array{Float64}(undef, rows, 1, 2)
    for r in 1:rows
        tuple = data[r]
        for (i, value) in enumerate(tuple)
            t_lin_cost[r, 1, i] = value
        end
    end
    return t_lin_cost
end

function transform_array_for_hdf(data::SortedDict{Dates.DateTime, Vector{PWL}})
    quad_cost = hcat(values(data)...)
    rows, cols = size(quad_cost)
    tuple_length = length(first(quad_cost))
    @assert_op length(first(first(quad_cost))) == 2
    t_quad_cost = Array{Float64}(undef, rows, cols, 2, tuple_length)
    for r in 1:rows, c in 1:cols
        tuple_array = quad_cost[r, c]
        for (j, tuple) in enumerate(tuple_array)
            for (i, value) in enumerate(tuple)
                t_quad_cost[r, c, i, j] = value
            end
        end
    end
    return t_quad_cost
end

function transform_array_for_hdf(data::Vector{PWL})
    rows = length(data)
    tuple_length = length(first(data))
    @assert_op length(first(first(data))) == 2
    t_quad_cost = Array{Float64}(undef, rows, 1, 2, tuple_length)
    for r in 1:rows
        tuple_array = data[r, 1]
        for (j, tuple) in enumerate(tuple_array)
            for (i, value) in enumerate(tuple)
                t_quad_cost[r, 1, i, j] = value
            end
        end
    end
    return t_quad_cost
end