text:
using DFTK
using Test
using LinearAlgebra: norm

include("testcases.jl")

function test_chi0(;symmetry=false, use_symmetry=false, temperature=0,
                   spin_polarization=:none,
                   kgrid=[3, 1, 1], fft_size=[10, 1, 10], Ecut=3)
    testcase = silicon
    n_bands  = 12
    tol      = 1e-14
    ε        = 1e-8
    testtol  = 2e-6

    collinear = spin_polarization == :collinear
    label = [
        symmetry        ? "   symm" : "no symm",
        use_symmetry    ? "   use" : "no use",
        temperature > 0 ? "temp" : "  0K",
        collinear       ? "coll" : "none",
    ]
    @testset "Computing χ0 ($(join(label, ", ")))" begin
        spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
        magnetic_moments = collinear ? [spec => 2rand(2)] : []
        model_kwargs = (temperature=temperature, symmetries=symmetry,
                        magnetic_moments=magnetic_moments,
                        spin_polarization=spin_polarization)
        basis_kwargs = (kgrid=kgrid, fft_size=fft_size, use_symmetry=use_symmetry)
        model = model_LDA(testcase.lattice, [spec => testcase.positions]; model_kwargs...)
        basis = PlaneWaveBasis(model; Ecut, basis_kwargs...)

        ρ0     = guess_density(basis, magnetic_moments)
        energies, ham0 = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)
        res = DFTK.next_density(ham0, tol=tol, eigensolver=diag_full, n_bands=n_bands)
        ρ1  = res.ρout

        # Now we make the same model, but add an artificial external potential ε * δV
        n_spin = model.n_spin_components
        δV = randn(eltype(basis), basis.fft_size..., n_spin)
        δV_sym = DFTK.symmetrize_ρ(basis, δV)
        if symmetry
            δV = δV_sym
        else
            @test δV_sym ≈ δV
        end

        εδV = ε * δV
        term_builder = basis -> DFTK.TermExternal(εδV)
        model = model_LDA(testcase.lattice, [spec => testcase.positions];
                          model_kwargs..., extra_terms=[term_builder])
        basis = PlaneWaveBasis(model; Ecut, basis_kwargs...)
        energies, ham = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)
        res = DFTK.next_density(ham, tol=tol, eigensolver=diag_full, n_bands=n_bands)
        ρ2     = res.ρout
        diff_findiff = (ρ2 - ρ1) / ε

        EVs = [eigen(Hermitian(Array(Hk))) for Hk in ham0.blocks]
        Es = [EV.values[1:n_bands] for EV in EVs]
        Vs = [EV.vectors[:, 1:n_bands] for EV in EVs]
        occ, εF = DFTK.compute_occupation(basis, Es)

        # Test apply_χ0 and compare against finite differences
        diff_applied_χ0 = apply_χ0(ham0, Vs, εF, Es, δV)
        @test norm(diff_findiff - diff_applied_χ0) < testtol

        # just to cover it here
        if temperature > 0
            N = compute_nos(εF, basis, Es)
            D = compute_dos(εF, basis, Es)
            LDOS = compute_ldos(εF, basis, Es, Vs)
        end

        if !symmetry
            # Test compute_χ0 against finite differences
            χ0 = compute_χ0(ham0)
            diff_computed_χ0 = reshape(χ0 * vec(δV), basis.fft_size..., n_spin)
            @test norm(diff_findiff - diff_computed_χ0) < testtol

            # Test that apply_χ0 is self-adjoint
            δV1 = randn(eltype(basis), basis.fft_size..., n_spin)
            δV2 = randn(eltype(basis), basis.fft_size..., n_spin)
            χ0δV1 = apply_χ0(ham0, Vs, εF, Es, δV1)
            χ0δV2 = apply_χ0(ham0, Vs, εF, Es, δV2)
            @test abs(dot(δV1, χ0δV2) - dot(δV2, χ0δV1)) < testtol
        end
    end
end

for temperature in (0, 0.03), spin_polarization in (:none, :collinear)
    for use_symmetry in (false, true), symmetry in (false, true)
        test_chi0(;symmetry=symmetry, use_symmetry=use_symmetry,
                  temperature=temperature, spin_polarization=spin_polarization)
    end
end


text:
struct CoordinateSystem{S} end

"""
    UnkownSpace
"""
const UnknownSpace = CoordinateSystem{:unknown}()

"""
    ScannerSpace

In scanner space.
"""
const ScannerSpace = CoordinateSystem{:scanner}()

"""
    AnatomicalSpace

equivalent of 'aligned' space in NIfTI standard.
"""
const AnatomicalSpace = CoordinateSystem{:anatomical}()

"""
    TailarachSpace

Tailarach space
"""
const TailarachSpace = CoordinateSystem{:tailarach}()

"""
    MNI152Space

MNI152 space
"""
const MNI152Space = CoordinateSystem{:MNI152}()

"""
    coordinate_system(x)

Return the coordinate space that `x` is in.
"""
coordinate_system(x::Any) = getter(x, "coordinatespace", CoordinateSystem, UnkownSpace)
# TODO : This probably shouldn't be mutable
function coordinate_system!(x::Any, val::CoordinateSystem)
    return setter!(x, "coordinate_system", val, CoordinateSystem)
end

ext: jl


text:
## Wing analysis case
using AeroMDAO

## Wing section setup
wing_right = HalfWing(foils     = Foil.(fill(naca4((0,0,1,2)), 3)),
                      chords    = [1.0, 0.6, 0.2],
                      twists    = [0.0, 0.0, 0.0],
                      spans     = [5.0, 0.5],
                      dihedrals = [5., 5.],
                      sweep_LEs = [5., 5.]);
wing = Wing(wing_right, wing_right)
wing_mac = mean_aerodynamic_center(wing)
S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing);
print_info(wing, "Wing")

## Assembly
ρ       = 1.225
x_w     = wing_mac[1]
ref     = [ x_w; 0.; 0. ]
V, α, β = 1.0, 5.0, 0.0
Ω       = [0.0, 0.0, 0.0]
fs      = Freestream(V, α, β, Ω)

## Evaluate case
@time nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes, Γs =
solve_case(wing, fs;
           rho_ref   = ρ,
           r_ref     = ref,
           area_ref  = S,
           span_ref  = b,
           chord_ref = c,
           span_num  = [25, 8],
           chord_num = 5,
           viscous   = true, # Only appropriate for α = β = 0, but works for other angles anyway
           x_tr      = [0.3, 0.3],
          #  spacing   = "uniform"
          );

print_coefficients(nf_coeffs, ff_coeffs, "Wing")

## Evaluate case with stability derivatives
@time nf, ff, dvs =
solve_stability_case(wing, fs;
                     rho_ref    = ρ,
                     r_ref      = ref,
                     span_num   = [25, 8],
                     chord_num  = 6,
                     name       = "My Wing",
                     viscous    = true,
                     x_tr       = [0.3, 0.3],
                     print      = false,
                    #  spacing   = ["sine", "cosine"]
                    );

#
print_coefficients(nf, ff, "Wing")
print_derivatives(dvs, "Wing")

## Plotting
using StaticArrays
using Plots
gr(size = (600, 400), dpi = 300)

## Streamlines

# Chordwise distribution
# num_points = 50
# max_z = 0.1
# y = span(wing) / 2 - 0.05
# seed = SVector.(fill(-0.1, num_points), fill(y, num_points), range(-max_z, stop = max_z, length = num_points))

# Spanwise distribution
span_points = 20
init        = chop_trailing_edge(wing, span_points)
dx, dy, dz  = 0, 0, 1e-3
seed        = [ init .+ Ref([dx, dy,  dz])  ;
                init .+ Ref([dx, dy, -dz]) ];

distance = 5
num_stream_points = 100
streams = plot_streams(fs, seed, horseshoes, Γs, distance, num_stream_points);

## Display
horseshoe_coords = plot_panels(horseshoe_panels[:])
wing_coords      = plot_wing(wing);
colpoints        = horseshoe_point.(horseshoe_panels)

# Coordinates
xs = getindex.(colpoints, 1);
ys = getindex.(colpoints, 2);
zs = getindex.(colpoints, 3);

z_limit = 5
plot(xaxis = "x", yaxis = "y", zaxis = "z",
     aspect_ratio = 1,
     camera = (30, 60),
     zlim = (-0.1, z_limit),
     size = (800, 600))
plot!.(horseshoe_coords, color = :black, label = :none)
scatter!(tupvector(colpoints)[:], marker = 1, color = :black, label = :none)
plot!.(streams, color = :green, label = :none)
plot!()

## Spanwise forces
wind_CFs = body_to_wind_axes.(CFs, fs.alpha, fs.beta)
CDis     = @. getindex(wind_CFs, 1)
CYs	     = @. getindex(wind_CFs, 2)
CLs      = @. getindex(wind_CFs, 3)

area_scale  = S ./ sum(panel_area, horseshoe_panels, dims = 1)[:]
span_CDis   = sum(CDis, dims = 1)[:] .* area_scale
span_CYs    = sum(CYs,  dims = 1)[:] .* area_scale
span_CLs    = sum(CLs,  dims = 1)[:] .* area_scale
CL_loadings = sum(Γs,   dims = 1)[:] / (0.5 * fs.V * c)

plot_CD = plot(ys[1,:], span_CDis, label = :none, ylabel = "CDi")
plot_CY = plot(ys[1,:], span_CYs, label = :none, ylabel = "CY")
plot_CL = begin
            plot(ys[1,:], span_CLs, label = :none, xlabel = "y", ylabel = "CL")
            plot!(ys[1,:], CL_loadings, label = "Normalized", xlabel = "y")
          end
plot(plot_CD, plot_CY, plot_CL, layout = (3,1))

## Lift distribution

# Exaggerated CF distribution for plot
hs_pts = bound_leg_center.(horseshoes)
hs_xs  = getindex.(hs_pts, 1)
hs_ys  = getindex.(hs_pts, 2)
hs_zs  = getindex.(hs_pts, 3)

plot(xaxis = "x", yaxis = "y", zaxis = "z",
     aspect_ratio = 1,
     camera = (60, 60),
     zlim = (-0.1, z_limit)
    )
plot!.(horseshoe_coords, color = :gray, label = :none)
# scatter!(cz_pts, zcolor = CLs[:], marker = 2, label = "CL (Exaggerated)")
quiver!(hs_xs[:], hs_ys[:], hs_zs[:], quiver=(CDis[:], CYs[:], CLs[:]) .* 100)
plot!(size = (800, 600))
plot!()

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
"""

"""
module CFTime

if VERSION >= v"0.7.0-beta.0"
    using Printf
    using Dates
    import Dates: UTInstant, Millisecond
    import Dates: year,  month,  day, hour, minute, second, millisecond
    import Dates: daysinmonth, daysinyear, yearmonthday, yearmonth
    import Dates: monthday, len, dayofyear, firstdayofyear
else
    using Missings
    import Base.Dates: UTInstant, Millisecond
    import Base.Dates: year,  month,  day, hour, minute, second, millisecond
    import Base.Dates: daysinmonth, daysinyear, yearmonthday, yearmonth
    import Base.Dates: monthday, len, dayofyear, firstdayofyear
    import Compat: occursin
end

import Base: +, -, isless, string, show, convert, reinterpret

const DEFAULT_TIME_UNITS = "days since 1900-01-01 00:00:00"

# Introduction of the Gregorian Calendar 1582-10-15
const GREGORIAN_CALENDAR = (1582,10,15)

# Time offset in days for the time origin
# if DATENUM_OFFSET = 0, then datenum_gregjulian
# corresponds to  Modified Julian Days (MJD).
# MJD is the number of days since midnight on 1858-11-17)

#const DATENUM_OFFSET = 2_400_000.5 # for Julian Days
const DATENUM_OFFSET = 0 # for Modified Julian Days

# Introduction of the Gregorian Calendar 1582-10-15
# expressed in MJD (if DATENUM_OFFSET = 0)

const DN_GREGORIAN_CALENDAR = -100840 + DATENUM_OFFSET

# DateTime(UTInstant{Millisecond}(Dates.Millisecond(0)))
# returns 0000-12-31T00:00:00
# 678576 is the output of datenum_prolepticgregorian(-1,12,31)

const DATETIME_OFFSET = Dates.Millisecond(678576 * (24*60*60*Int64(1000)))

"""
    dn = datenum_gregjulian(year,month,day,gregorian::Bool)

Days since 1858-11-17 according to the Gregorian (`gregorian` is `true`) or
Julian Calendar (`gregorian` is `false`) based on the Algorithm of
Jean Meeus [1].

The year -1, correspond to 1 BC. The year 0 does not exist in the
Gregorian or Julian Calendar.

[1] Meeus, Jean (1998) Astronomical Algorithms (2nd Edition).
Willmann-Bell,  Virginia. p. 63
"""
function datenum_gregjulian(year,month,day,gregorian::Bool)
    # turn year equal to -1 (1 BC) into year = 0
    if year < 0
        year = year+1
    end

    if gregorian
        # bring year in range of 1601 to 2000
        ncycles = (2000 - year) ÷ 400
        year = year + 400 * ncycles
        return datenum_ac(year,month,day,gregorian) - ncycles*146_097
    else
        return datenum_ac(year,month,day,gregorian)
    end

end

# Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,  Virginia. p. 63
# However, the algorithm does not work for -100:03:01 and before in
# the proleptic Gregorian Calendar

function datenum_ac(year,month,day,gregorian::Bool)

    if month <= 2
        # if the date is January or February, it is considered
        # the 13rth or 14th month of the preceeding year
        year = year - 1
        month = month + 12
    end

    B =
        if gregorian
            A = year ÷ 100
            2 - A + A ÷ 4
        else
            0
        end

    # benchmark shows that it is 40% faster replacing
    # trunc(Int64,365.25 * (year + 4716))
    # by
    # (1461 * (year + 4716)) ÷ 4
    #
    # and other floating point divisions

    # Z is the Julian Day plus 0.5
    # 1461/4 is 365.25
    # 153/5 is 30.6

    # why 153/5 (or 30.6001 ?)
    # month+1 varies between 4 (March), 5 (April), .. 14 (December),
    # 15 (January), 16 (February)

    # cm = 153 * (4:16) ÷ 5; cm[2:end]-cm[1:end-1]
    #
    # length of each month
    # --------------------
    # 31  March
    # 30  April
    # 31  May
    # 30  June
    # 31  July
    # 31  August
    # 30  September
    # 31  October
    # 30  November
    # 31  December
    # 31  January
    # 30  February (wrong, but not used, since it is the last month)

    Z = (1461 * (year + 4716)) ÷ 4 + (153 * (month+1)) ÷ 5 + day + B - 2401525
    # Modified Julan Day
    return Z + DATENUM_OFFSET
end

"""
    year, month, day = datetuple_gregjulian(Z::Integer,gregorian::Bool)

Compute year, month and day from Z which is the Modified Julian Day
for the Gregorian (true) or Julian (false) calendar.

For example:
Z = 0 for the 1858 November 17 00:00:00

Algorithm:

Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,
Virginia. p. 63
"""
function datetuple_gregjulian(Z,gregorian::Bool)
    # Z is Julian Day plus 0.5
    Z = Z + 2_400_001 - DATENUM_OFFSET

    A =
        if gregorian
            # lets magic happen
            α = trunc(Int64, (Z - 1867_216.25)/36524.25)
            #@show α,Z - 1867_216.25
            Z + 1 + α - (α ÷ 4)
        else
            Int64(Z)
        end

    # even more magic...
    B = A + 1524
    C = trunc(Int64, (B - 122.1) / 365.25)
    D = trunc(Int64, 365.25 * C)
    E = trunc(Int64, (B-D)/30.6001)

    day = B - D - trunc(Int64,30.6001 * E)
    month = (E < 14 ? E-1 : E-13)
    y = (month > 2 ? C - 4716 : C - 4715)

    # turn year 0 into year -1 (1 BC)
    if y <= 0
        y = y-1
    end
    return y,month,day
end

"""
    days,h,mi,s,ms = timetuplefrac(time::Number)

Return the number of whole days, hours (`h`), minutes (`mi`), seconds (`s`) and
millisecods (`ms`) from `time` expressed in milliseconds.
"""
function timetuplefrac(time::Number)
    # time can be negative, use fld instead of ÷
    days = fld(Int64(time), (24*60*60*1000))
    ms = Int64(time) - days * (24*60*60*1000)

    h = ms ÷ (60*60*1000)
    ms = ms - h * (60*60*1000)

    mi = ms ÷ (60*1000)
    ms = ms - mi * (60*1000)

    s = ms ÷ 1000
    ms = ms - s * 1000
    return (days,h,mi,s,ms)
end

function datenumfrac(days,h,mi,s,ms)
    ms = 60*60*1000 * h +  60*1000 * mi + 1000*s + ms
    return (24*60*60*1000) * Int64(days) + ms
end

datetuple_prolepticgregorian(Z) = datetuple_gregjulian(Z,true)
datetuple_julian(Z) = datetuple_gregjulian(Z,false)
datetuple_standard(Z) = datetuple_gregjulian(Z,Z >= DN_GREGORIAN_CALENDAR)

datenum_prolepticgregorian(y,m,d) = datenum_gregjulian(y,m,d,true)
datenum_julian(y,m,d) = datenum_gregjulian(y,m,d,false)
datenum_standard(y,m,d) = datenum_gregjulian(y,m,d,(y,m,d) >= GREGORIAN_CALENDAR)

function datenum_cal(cm, y, m, d)
    # turn year equal to -1 (1 BC) into year = 0
    if y < 0
        y = y+1
    end

    if m < 1 || m > 12
        error("invalid month $(m)")
    end

    if d < 1 || d > (cm[m+1] - cm[m])
        error("invalid day $(d) in $(@sprintf("%04d-%02d-%02d",y,m,d))")
    end

    return cm[end] * (y-1) + cm[m] + (d-1)
end

# Calendar with regular month-length

function findmonth(cm,t2)
    mo = length(cm)
    while cm[mo] > t2
        mo -= 1
    end
    return mo
end

function datetuple_cal(cm,timed_::Number)
    y = fld(Int64(timed_), cm[end])
    t2 = Int64(timed_) - cm[end]*y

    # find month
    mo = findmonth(cm,t2)
    d = t2 - cm[mo]

    # day and year start at 1 (not zero)
    d = d+1
    y = y+1

    if y <= 0
        y = y-1
    end

    return (y,mo,d)
end

for (calendar,cmm) in [
    ("allleap", (0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366)),
    ("noleap",  (0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365)),
    ("360day",    (0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360)),
]
    @eval begin
        $(Symbol(:datenum_,calendar))(y, m, d) = datenum_cal($cmm, y, m, d)
        $(Symbol(:datetuple_,calendar))(days) = datetuple_cal($cmm,days)
    end
end

abstract type AbstractCFDateTime <: Dates.TimeType
end

const RegTime = Union{Dates.Millisecond,Dates.Second,Dates.Minute,Dates.Hour,Dates.Day}

function validargs(::Type{T},arg...) where T <: AbstractCFDateTime
    return nothing
end

for (CFDateTime,calendar) in [(:DateTimeStandard,"standard"),
                              (:DateTimeJulian,"julian"),
                              (:DateTimeProlepticGregorian,"prolepticgregorian"),
                              (:DateTimeAllLeap,"allleap"),
                              (:DateTimeNoLeap,"noleap"),
                              (:DateTime360Day,"360day")]
    @eval begin
        struct $CFDateTime <: AbstractCFDateTime
            instant::UTInstant{Millisecond}
            $CFDateTime(instant::UTInstant{Millisecond}) = new(instant)
        end

        """
    $($CFDateTime)(y, [m, d, h, mi, s, ms]) -> $($CFDateTime)

Construct a `$($CFDateTime)` type by year (`y`), month (`m`, default 1),
day (`d`, default 1), hour (`h`, default 0), minute (`mi`, default 0),
second (`s`, default 0), millisecond (`ms`, default 0).
All arguments must be convertible to `Int64`.
`$($CFDateTime)` is a subtype of `AbstractCFDateTime`.

The netCDF CF calendars are defined in [the CF Standard](http://cfconventions.org/cf-conventions/cf-conventions.html#calendar).
        """
        function $CFDateTime(y::Int64, m::Int64=1, d::Int64=1,
                             h::Int64=0, mi::Int64=0, s::Int64=0, ms::Int64=0)

            days = $(Symbol(:datenum_,calendar))(y,m,d)
            totalms = datenumfrac(days,h,mi,s,ms)
            return $CFDateTime(UTInstant(Millisecond(totalms)))
        end

        # Fallback constructors
        $CFDateTime(y::Number, m=1, d=1, h=0, mi=0, s=0, ms=0) = $CFDateTime(
            Int64(y), Int64(m), Int64(d), Int64(h), Int64(mi), Int64(s),
            Int64(ms))

        if VERSION >= v"1.3-"
            function $CFDateTime(y, m, d, h, mi, s, ms, ampm)
                @assert ampm == Dates.TWENTYFOURHOUR
                return $CFDateTime(y, m, d, h, mi, s, ms)
            end
        end
        """
    $($CFDateTime)(dt::AbstractString, format::AbstractString; locale="english") -> $($CFDateTime)

Construct a $($CFDateTime) by parsing the `dt` date time string following the
pattern given in the `format` string.

!!! note
    This function is experimental and migth
    be removed in the future. It relies on some internal function of `Dates` for
    parsing the `format`.
"""
        function $CFDateTime(dt::AbstractString, format::AbstractString; locale="english")
            return parse($CFDateTime, dt, DateFormat(format, locale))
        end

        $CFDateTime(dt::AbstractString, format::DateFormat) =
            parse($CFDateTime, dt, format)

        function datetuple(dt::$CFDateTime)
            time = Dates.value(dt.instant.periods)
            days,h,mi,s,ms = timetuplefrac(time)
            y, m, d = $(Symbol(:datetuple_,calendar))(days)
            return y, m, d, h, mi, s, ms
        end

        function string(dt::$CFDateTime)
            y,mo,d,h,mi,s,ms = datetuple(dt)
            return @sprintf("%04d-%02d-%02dT%02d:%02d:%02d",y,mo,d,h,mi,s)
        end

        function show(io::IO,dt::$CFDateTime)
            write(io, string(typeof(dt)), "(",string(dt),")")
        end

        function +(dt::$CFDateTime,Δ::Dates.Year)
            y,mo,d,h,mi,s,ms = datetuple(dt)
            return $CFDateTime(y+Dates.value(Δ), mo, d, h, mi, s, ms)
        end

        function +(dt::$CFDateTime,Δ::Dates.Month)
            y,mo,d,h,mi,s,ms = datetuple(dt)
            mo = mo + Dates.value(Δ)
            mo2 = mod(mo - 1, 12) + 1
            y = y + (mo-mo2) ÷ 12
            return $CFDateTime(y, mo2, d,h, mi, s, ms)
        end

        +(dt::$CFDateTime,Δ::RegTime) = $CFDateTime(UTInstant(dt.instant.periods + Dates.Millisecond(Δ)))

        -(dt1::$CFDateTime,dt2::$CFDateTime) = dt1.instant.periods - dt2.instant.periods

        isless(dt1::$CFDateTime,dt2::$CFDateTime) = dt1.instant.periods < dt2.instant.periods

    end
end

"""
    dt2 = reinterpret(::Type{T}, dt)

Convert a variable `dt` of type `DateTime`, `DateTimeStandard`, `DateTimeJulian`,
`DateTimeProlepticGregorian`, `DateTimeAllLeap`, `DateTimeNoLeap` or
`DateTime360Day` into the date time type `T` using the same values for
year, month, day, minute, second and millisecond.
The conversion might fail if a particular date does not exist in the
target calendar.
"""
function reinterpret(::Type{T1}, dt::T2) where T1 <: Union{AbstractCFDateTime,DateTime} where T2 <: Union{AbstractCFDateTime,DateTime}
   return T1(
       Dates.year(dt),Dates.month(dt),Dates.day(dt),
       Dates.hour(dt),Dates.minute(dt),Dates.second(dt),
       Dates.millisecond(dt))
end

"""
    dt2 = convert(::Type{T}, dt)

Convert a DateTime of type `DateTimeStandard`, `DateTimeProlepticGregorian`,
`DateTimeJulian` or `DateTime` into the type `T` which can also be either
`DateTimeStandard`, `DateTimeProlepticGregorian`, `DateTimeJulian` or `DateTime`.

Converstion is done such that durations (difference of DateTime types) are
preserved. For dates on and after 1582-10-15, the year, month and days are the same for
the types `DateTimeStandard`, `DateTimeProlepticGregorian` and `DateTime`.

For dates before 1582-10-15, the year, month and days are the same for
the types `DateTimeStandard` and `DateTimeJulian`.
"""
function convert(::Type{T1}, dt::T2) where T1 <: Union{DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian} where T2 <: Union{DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian}
    return T1(dt.instant)
end

function convert(::Type{DateTime}, dt::T2) where T2 <: Union{DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian}
    DateTime(UTInstant{Millisecond}(dt.instant.periods + DATETIME_OFFSET))
end

function convert(::Type{T1}, dt::DateTime) where T1 <: Union{DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian}
    T1(UTInstant{Millisecond}(dt.instant.periods - DATETIME_OFFSET))
end

Dates.year(dt::AbstractCFDateTime) = datetuple(dt)[1]
Dates.month(dt::AbstractCFDateTime) = datetuple(dt)[2]
Dates.day(dt::AbstractCFDateTime) = datetuple(dt)[3]
Dates.hour(dt::AbstractCFDateTime)   = datetuple(dt)[4]
Dates.minute(dt::AbstractCFDateTime) = datetuple(dt)[5]
Dates.second(dt::AbstractCFDateTime) = datetuple(dt)[6]
Dates.millisecond(dt::AbstractCFDateTime) = datetuple(dt)[7]

for func in (:year, :month, :day, :hour, :minute, :second, :millisecond)
    name = string(func)
    @eval begin
        @doc """
            Dates.$($name)(dt::AbstractCFDateTime) -> Int64

        Extract the $($name)-part of a `AbstractCFDateTime` as an `Int64`.
        """ $func(dt::AbstractCFDateTime)
    end
end

-(dt::AbstractCFDateTime,Δ) = dt + (-Δ)

function parseDT(::Type{DT},str) where DT <: Union{DateTime,AbstractCFDateTime}
    str =
        if VERSION >= v"0.7.0-beta.0"
            replace(str,"T" => " ")
        else
            replace(str,"T"," ")
        end

    # remove Z time zone indicator
    # all times are assumed UTC anyway
    if endswith(str,"Z")
        str = str[1:end-1]
    end

    negativeyear = str[1] == '-'
    if negativeyear
        str = str[2:end]
    end

    y,m,d,h,mi,s,ms =
        if occursin(" ",str)
            datestr,timestr = split(str,' ')
            y,m,d = parse.(Int64,split(datestr,'-'))

            #=
            hour_min_sec = parse.(Int64,split(timestr,':'))
            s = 0
            mi = 0
            if length(hour_min_sec) == 3
                h,mi,s = hour_min_sec
            elseif length(hour_min_sec) == 2
                h,mi = hour_min_sec
            elseif length(hour_min_sec) == 1
                h, = hour_min_sec
            end

            (y,m,d,h,mi,s,Int64(0))
            =#

            time_split = split(timestr,':')

            h_str, mi_str, s_str =
                if length(time_split) == 2
                    time_split[1], time_split[2], "00"
                else
                    time_split
                end

            h = parse(Int64,h_str)
            mi = parse(Int64,mi_str)

            s,ms =
                if occursin('.',s_str)
                    # seconds contain a decimal point, e.g. 00:00:00.0
                    secfrac = parse(Float64,s_str)
                    s = floor(Int64,secfrac)
                    ms = round(Int64,1000*(secfrac - s))
                    s,ms
                else
                    (parse(Int64,s_str),Int64(0))
                end

            (y,m,d,h,mi,s,ms)
        else
            y,m,d = parse.(Int64,split(str,'-'))
            (y,m,d,Int64(0),Int64(0),Int64(0),Int64(0))
        end

    if negativeyear
        y = -y
    end

    return DT(y,m,d,h,mi,s,ms)
end

function timeunits(::Type{DT},units) where DT
    tunit_mixedcase,starttime = strip.(split(units," since "))
    tunit = lowercase(tunit_mixedcase)

    t0 = parseDT(DT,starttime)

    # make sure that plength is 64-bit on 32-bit platforms
    plength =
        if (tunit == "days") || (tunit == "day")
            24*60*60*Int64(1000)
        elseif (tunit == "hours") || (tunit == "hour")
            60*60*Int64(1000)
        elseif (tunit == "minutes") || (tunit == "minute")
            60*Int64(1000)
        elseif (tunit == "seconds") || (tunit == "second")
            Int64(1000)
        else
            error("unknown units $(tunit)")
        end

    return t0,plength
end

function timetype(calendar = "standard")
    DT =
        if (calendar == "standard") || (calendar == "gregorian")
            DateTimeStandard
        elseif calendar == "proleptic_gregorian"
            DateTimeProlepticGregorian
        elseif calendar == "julian"
            DateTimeJulian
        elseif (calendar == "noleap") || (calendar == "365_day")
            DateTimeNoLeap
        elseif (calendar == "all_leap") || (calendar == "366_day")
            DateTimeAllLeap
        elseif calendar == "360_day"
            DateTime360Day
        else
            error("Unsupported calendar: $(calendar)")
        end

    return DT
end

timetype(dt::Type{<:AbstractCFDateTime}) = dt
timetype(dt::Type{DateTime}) = dt

"""
    t0,plength = timeunits(units,calendar = "standard")

Parse time units (e.g. "days since 2000-01-01 00:00:00") and returns the start
time `t0` and the scaling factor `plength` in milliseconds.
"""
function timeunits(units, calendar = "standard")
    DT = timetype(calendar)
    return timeunits(DT,units)
end

function timedecode(::Type{DT},data::AbstractArray{Float32,N},units) where {DT,N}
    # convert to Float64
    return timedecode(DT,Float64.(data),units)
end

function timedecode(::Type{DT},data,units) where DT
    t0,plength = timeunits(DT,units)
    convert(x) = t0 + Dates.Millisecond(round(Int64,plength * x))
    return convert.(data)
end

"""
    dt = timedecode(data,units,calendar = "standard", prefer_datetime = true)

Decode the time information in data as given by the units `units` according to
the specified calendar. Valid values for `calendar` are
`"standard"`, `"gregorian"`, `"proleptic_gregorian"`, `"julian"`, `"noleap"`, `"365_day"`,
`"all_leap"`, `"366_day"` and `"360_day"`.

If `prefer_datetime` is `true` (default), dates are
converted to the `DateTime` type (for the calendars
"standard", "gregorian", "proleptic_gregorian" and "julian"). Such conversion is
not possible for the other calendars.

| Calendar            | Type (prefer_datetime=true) | Type (prefer_datetime=false) |
| ------------------- | --------------------------- | ---------------------------- |
| `standard`, `gregorian` | `DateTime`                    | `DateTimeStandard`             |
| `proleptic_gregorian` | `DateTime`                    | `DateTimeProlepticGregorian`   |
| `julian`              | `DateTime`                    | `DateTimeJulian`               |
| `noleap`, `365_day`     | `DateTimeNoLeap`              | `DateTimeNoLeap`               |
| `all_leap`, `366_day`   | `DateTimeAllLeap`             | `DateTimeAllLeap`              |
| `360_day`             | `DateTime360Day`              | `DateTime360Day`               |

## Example:

```julia
using CFTime, Dates
# standard calendar
dt = CFTime.timedecode([0,1,2,3],"days since 2000-01-01 00:00:00")
# 4-element Array{Dates.DateTime,1}:
#  2000-01-01T00:00:00
#  2000-01-02T00:00:00
#  2000-01-03T00:00:00
#  2000-01-04T00:00:00

dt = CFTime.timedecode([0,1,2,3],"days since 2000-01-01 00:00:00","360_day")
# 4-element Array{DateTime360Day,1}:
#  DateTime360Day(2000-01-01T00:00:00)
#  DateTime360Day(2000-01-02T00:00:00)
#  DateTime360Day(2000-01-03T00:00:00)
#  DateTime360Day(2000-01-04T00:00:00)
```

"""
function timedecode(data,units,calendar = "standard"; prefer_datetime = true)
    DT = timetype(calendar)
    dt = timedecode(DT,data,units)

    if prefer_datetime &&
        (DT in [DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian])

        return convert.(DateTime,dt)
    else
        return dt
    end
end

"""
    data = timeencode(dt,units,calendar = "standard")

Convert a vector or array of `DateTime` (or `DateTimeStandard`,
`DateTimeProlepticGregorian`, `DateTimeJulian`, `DateTimeNoLeap`,
`DateTimeAllLeap`, `DateTime360Day`) accoring to
the specified units (e.g. `"days since 2000-01-01 00:00:00"`) using the calendar
`calendar`.  Valid values for calendar are:
`"standard"`, `"gregorian"`, `"proleptic_gregorian"`, `"julian"`, `"noleap"`, `"365_day"`,
`"all_leap"`, `"366_day"`, `"360_day"`.
"""
function timeencode(data::AbstractArray{DT,N},units,
                    calendar = "standard") where N where DT <: Union{DateTime,AbstractCFDateTime,Union{DateTime,AbstractCFDateTime,Missing}}

    DT2 = timetype(calendar)
    t0,plength = timeunits(DT2,units)

    function encode(dt)
        if ismissing(dt)
            return missing
        end

        tmp =
            try
                convert.(DT2,dt)
            catch
                error("It is not possible to convert from $(DT) to $(DT2)")
            end

        return Dates.value(tmp - t0) / plength
    end
    return encode.(data)
end

function timeencode(data::DT,units,
                    calendar = "standard") where DT <: Union{DateTime,AbstractCFDateTime}
    return timeencode([data],units,calendar)[1]
end

# do not transform data is not a vector of DateTime
# unused, should be removed
timeencode(data,units,calendar = "standard") = data

export timeencode, timedecode, datetuple

# utility functions

"""
    monthlength = daysinmonth(::Type{DT},y,m)

Returns the number of days in a month for the year `y` and the month `m`
according to the calenar given by the type `DT`.

Example
```julia-repl
julia> daysinmonth(DateTimeAllLeap,2001,2)
29
```

"""
function daysinmonth(::Type{DT},y,m) where DT <: Union{DateTime, AbstractCFDateTime}
    t = DT(y,m,1)
    return Dates.value((t + Dates.Month(1)) - t) ÷ (24*60*60*1000)
end

"""
    monthlength = daysinmonth(t)

Returns the number of days in a month containing the date `t`

Example
```julia-repl
julia> daysinmonth(DateTimeAllLeap(2001,2,1))
29
```
"""
function daysinmonth(t::DT) where DT <: Union{DateTime, AbstractCFDateTime}
    return daysinmonth(DT,Dates.year(t),Dates.month(t))
end

"""
    yearlength = daysinyear(::Type{DT},y)

Returns the number of days in a year for the year `y`
according to the calenar given by the type `DT`.

Example
```julia-repl
julia> daysinyear(DateTimeAllLeap,2001,2)
366
```

"""
function daysinyear(::Type{DT},y) where DT <: Union{DateTime, AbstractCFDateTime}
    t = DT(y,1,1)
    return Dates.value((t + Dates.Year(1)) - t) ÷ (24*60*60*1000)
end

"""
    yearlength = daysinyear(t)

Returns the number of days in a year containing the date `t`

Example
```julia-repl
julia> daysinyear(DateTimeAllLeap(2001,2,1))
366
```
"""
function daysinyear(t::DT) where DT <: Union{DateTime, AbstractCFDateTime}
    return daysinyear(DT,Dates.year(t))
end

"""
    yearmonthday(dt::AbstractCFDateTime) -> (Int64, Int64, Int64)

Simultaneously return the year, month and day parts of `dt`.
"""
yearmonthday(dt::AbstractCFDateTime) = (Dates.year(dt),Dates.month(dt),Dates.day(dt))

"""
    yearmonth(dt::AbstractCFDateTime) -> (Int64, Int64)

Simultaneously return the year and month parts of `dt`.
"""
yearmonth(dt::AbstractCFDateTime) = (Dates.year(dt),Dates.month(dt))

"""
    monthday(dt::AbstractCFDateTime) -> (Int64, Int64)

Simultaneously return the month and day parts of `dt`.
"""
monthday(dt::AbstractCFDateTime) = (Dates.month(dt),Dates.day(dt))

"""
    firstdayofyear(dt::AbstractCFDateTime) -> Int

Return the first day of the year including the date `dt`
"""
firstdayofyear(dt::T) where T <: AbstractCFDateTime = T(Dates.year(dt),1,1,0,0,0)

"""
    dayofyear(dt::AbstractCFDateTime) -> Int

Return the day of the year for dt with January 1st being day 1.
"""
function dayofyear(dt::AbstractCFDateTime)
    t0 = firstdayofyear(dt)
    return Dates.value(dt - t0) ÷ (24*60*60*1000) + 1
end

function Dates.len(first::T, last::T, step::DT) where T <: AbstractCFDateTime where
    DT <: Union{Dates.Day,Dates.Hour,Dates.Minute,Dates.Second,Dates.Millisecond}
    return Dates.value(last-first) ÷ Dates.value(Dates.Millisecond(step))
end

function Dates.len(first::T, last::T, step) where T <: AbstractCFDateTime
    if Dates.value(step) == 0
        error("the step should not be zero")
    end
    len = 0
    next = first+step
    while next <= last
        next = next+step
        len = len+1
    end
    return len
end

export daysinmonth,
    daysinyear,
    yearmonthday,
    yearmonth,
    monthday,
    firstdayofyear,
    dayofyear,
    DateTimeStandard,
    DateTimeJulian,
    DateTimeProlepticGregorian,
    DateTimeAllLeap,
    DateTimeNoLeap,
    DateTime360Day,
    AbstractCFDateTime

function __init__()
    for CFDateTime in [DateTimeStandard,
                       DateTimeJulian,
                       DateTimeProlepticGregorian,
                       DateTimeAllLeap,
                       DateTimeNoLeap,
                       DateTime360Day]
        Dates.CONVERSION_TRANSLATIONS[CFDateTime] = Dates.CONVERSION_TRANSLATIONS[DateTime]
    end
end

end


text:
module CellularAutomata

using Graphs
using Random
using PyPlot

export CellularAutomaton
export run!
export step!
export plot

mutable struct CellularAutomaton
    initial_config::AbstractArray
    states::AbstractArray
    function CellularAutomaton(initial_config)
        new(initial_config, [initial_config])
    end
end

function run!(automaton::CellularAutomaton, rule::Int, n::Int)
    for i in 1:n
        step!(automaton, rule)
    end
    return automaton
end

function step!(automaton::CellularAutomaton, rule::Int)
    g = Graphs.grid([length(automaton.initial_config), 1], periodic=true)
    push!(automaton.states, next_iteration(last(automaton.states), g, rule))
    return automaton
end

function next_iteration(iter, g, rule)
    return [next_state(get_bit_rep(iter, i, g), rule) for i in 1:length(iter)]
end

function next_state(bit_rep, rule)
    bit_rule = last(bitstring(rule), 8)
    rule_index = 8 - parse(Int, bit_rep, base=2)
    next_state = bit_rule[rule_index]
    return parse(Int, next_state)
end

function get_bit_rep(iter, idx, g)
    n = copy(neighbors(g, idx))
    push!(n, idx)
    n_with_mod = [(i, i % 10) for i in n]
    sort!(n_with_mod, by=x->x[2])
    n = [i[1] for i in n_with_mod]
    cell_bits = [iter[i] for i in n]
    bit_rep = reduce(*, string.(Int.(cell_bits)))
    return bit_rep
end

function plot(automaton::CellularAutomaton)
    imshow(automaton.states, cmap="gray")
end

end  # end module


text:
using RadeonProRender, GeometryTypes, GLAbstraction, GLVisualize, Colors, ModernGL, FileIO, Reactive
using GLWindow
const FR = RadeonProRender
w = glscreen()
# Create OpenCL context using a single GPU, which is the default
context, glframebuffer = interactive_context(w)

# Create a scene
scene = FR.Scene(context)
N = 256
const xrange = linspace(-5.0f0, 5.0f0, N)
const yrange = linspace(-5.0f0, 5.0f0, N)

z = Float32[sin(1.3 * x) * cos(0.9 * y) + cos(0.8 * x) * sin(1.9 * y) + cos(y * 0.2 * x)
            for x in xrange, y in yrange]
mini = minimum(z)
maxi = maximum(z)
const cmap = map(x -> RGBA{U8}(x, 1.0), colormap("Blues"))
to_color(val, mini, maxi) = cmap[floor(Int, (((val - mini) / (maxi - mini)) * (length(cmap) - 1))) + 1]
image = map(z) do height
    return to_color(height, mini, maxi)
end
plane = FR.Shape(context, SimpleRectangle(-5, -5, 10, 10), (N, N))
push!(scene, plane)
const t = Signal(2.0f0 * pi)
matsys = FR.MaterialSystem(context, 0)
reflective = FR.MaterialNode(matsys, FR.MATERIAL_NODE_MICROFACET)
tex = FR.MaterialNode(matsys, FR.MATERIAL_NODE_IMAGE_TEXTURE)
displace = FR.Image(context, z)
color = FR.Image(context, image)
set!(tex, "data", color)
set!(reflective, "color", tex)
set!(plane, reflective)
set!(plane, displace)
setdisplacementscale!(plane, 1.0)
setsubdivisions!(plane, 1)

preserve(map(t) do i
             z = Float32[sin(1.3 * x * i) * cos(0.9 * y * i) +
                         cos(0.8 * x * i) * sin(1.9 * y * i) +
                         cos(y * 0.2 * x * i) for x in xrange, y in yrange]
             mini = minimum(z)
             maxi = maximum(z)
             image = map(z) do height
                 return to_color(height, mini, maxi)
             end
             displace = FR.Image(context, z)
             color = FR.Image(context, image)
             set!(tex, "data", color)
             return set!(plane, displace)
         end)

camera = FR.Camera(context)
lookat!(camera, Vec3f(0, 9.5, 11), Vec3f(0), Vec3f(0, 0, 1))
FR.CameraSetFocusDistance(camera.x, 7.0f0)
FR.CameraSetFStop(camera.x, 1.8f0)

set!(scene, camera)
set!(context, scene)
ibl = FR.EnvironmentLight(context)
imgpath = joinpath("C:\\", "Program Files", "KeyShot6", "bin", "Materials 2k.hdr")
set!(ibl, context, imgpath)
set!(scene, ibl)
pl = FR.PointLight(context);
setradiantpower!(pl, fill(10.0f0^3, 3)...)
transform!(pl, translationmatrix(Vec3f(0, 0, 7)))

push!(scene, pl)

set_standard_tonemapping!(context)

frame = 1
for i in (0.5f0:0.01f0:(pi * 2.0f0))
    push!(t, i)
    yield()
    isopen(w) || break
    clear!(glframebuffer)
    for i in 1:10
        glBindTexture(GL_TEXTURE_2D, 0)
        @time render(context)
        isopen(w) || break
        GLWindow.render_frame(w)
    end
    isopen(w) || break
    try
        screenshot(w; path="test$frame.png")
    catch e
        println(e)
    end
    frame += 1
end


text:
using WaterLily
using LinearAlgebra: norm2
include("TwoD_plots.jl")

function TwoD_julia_video(;p=6,Re=250,stop=60.)
    # Set simulation size & physical parameters
    n,m = 2^p,2^p
    U,R,r = [0.,-1.], m/16, m/16/0.75
    ν=norm2(U)*R/Re
    @show R,ν

    # Immerse three well placed circles (change for other shapes)
    z = [R*exp(im*θ) for θ ∈ range(0,2π,length=33)]
    centers = [n/2+im*3*m/4+r*exp(im*ϕ) for ϕ ∈ [π/2,π/2+2π/3,π/2-2π/3]]
    colors = [:forestgreen,:brown3,:mediumorchid3]
    body = AutoBody() do x,t  # signed distance function
            minimum(centers) do center
                norm2(complex(x...) - center) - R
            end
    end

    # Initialize simulation
    sim = Simulation((n+2,m+2),U,R;ν,body)

    # Solve flow and make nice gif
    t = range(0.,stop;step=0.3)
    gr(show = false, size=(700,600))
    @time @gif for tᵢ in t
        println("tU/L=",round(tᵢ,digits=4))
        sim_step!(sim,tᵢ)
        @inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
        flood(sim.flow.σ,shift=(-0.5,-0.5),clims=(-5,5),
            cfill=:Blues,legend=false,border=:none)
        for (center,color) ∈ zip(centers,colors)
            addbody(real(z.+center),imag(z.+center),c=color)
        end
    end
    return sim
end


text:
@inline function MajorMinorVersion(version::VersionNumber)
    major = version.major::Base.VInt
    minor = version.minor::Base.VInt
    result = MajorMinorVersion(major, minor)
    return result
end

@inline MajorMinorVersion(::Nothing) = nothing

@inline function VersionNumber(x::MajorMinorVersion)
    major = x.major::Base.VInt
    minor = x.minor::Base.VInt
    result = VersionNumber(x, y)::VersionNumber
    return result
end

@inline function Base.isless(a::MajorMinorVersion, b::MajorMinorVersion)
    a_versionnumber = VersionNumber(a)::VersionNumber
    b_versionnumber = VersionNumber(b)::VersionNumber
    result = Base.isless(a_versionnumber, b_versionnumber)
    return result
end


text:
# This file is part of RLinearAlgebra.jl

# using LinearAlgebra

"""
    LSLogOracle <: LinSysSolverLog

A mutable structure that stores information about a randomized linear solver's behavior.
    The log assumes that the true solution of the system is known and will be supplied.
    The goal of this log is usually for research, development, or testing.

# Fields
- `solution::AbstractVector`, a solution to the given linear system.
- `collection_rate::Int64`, the frequency with which to record information to append to the
    remaining fields, starting with the initialization (i.e., iteration 0).
- `error_hist::Vector{Float64}`, retains a vector of numbers corresponding to the error
    between the iterates of the solver and the `solution`
- `error_norm::Function`, a function that accepts a single vector argument and returns a
    scalar. Used to compute the error size.
- `resid_hist::Vector{Float64}`, retains a vector of numbers corresponding to the residual
    (uses the whole system to compute the residual)
- `resid_norm::Function`, a function that accepts a single vector argument and returns a
    scalar. Used to compute the residual size.
- `iterations::Int64`, the number of iterations of the solver.
- `converged::Bool`, a flag to indicate whether the system has converged by some measure

# Constructors
- Calling `LSLogOracle(x_star::Vector{Float64})` sets `solution = x_star`,
    `collection_rate = 1`, `error_hist = Float64[]`, `error_norm = norm` (Euclidean norm),
    `resid_hist = Float64[]`, `resid_norm = norm` (Euclidean norm), `iterations=-1`,
    and `converged = false`.
- Calling `LSLogOracle(x_star::Vector{Float64}, cr::Int64)` sets the structure with the
    same parameters as for `LSLogOracle(x_star)` except `collection_rate = cr`.
"""
mutable struct LSLogOracle <: LinSysSolverLog
    solution::AbstractVector
    collection_rate::Int64
    error_hist::Vector{Float64}
    error_norm::Function
    resid_hist::Vector{Float64}
    resid_norm::Function
    iterations::Int64
    converged::Bool
end
LSLogOracle(x_star::Vector{Float64}) = LSLogOracle(x_star, 1, Float64[], norm, Float64[],
    norm, -1, false)
LSLogOracle(x_star::Vector{Float64}, cr::Int64) = LSLogOracle(x_star, cr, Float64[], norm,
    Float64[],norm, -1, false)

# Common interface for update
function log_update!(
    log::LSLogOracle,
    sampler::LinSysSampler,
    x::AbstractVector,
    samp::Tuple,
    iter::Int64,
    A::AbstractArray,
    b::AbstractVector
)
    # Update iteration counter
    log.iterations = iter

    # Check whether to record information
    if mod(iter, 1:log.collection_rate) == log.collection_rate
        # Push absolute error between iterate and solution
        push!(log.error_hist, log.error_norm(x - log.solution))

        # Push residual norm of iterate
        push!(log.resid_hist, log.resid_norm(A * x - b))
    end

    return nothing
end


text:
ts = TS(df_timetype_index)

function test_types(obj::TS)
    @test typeof(obj.coredata) == DataFrame
end

# getindex(ts, i::Int)
ind = 1
test_types(ts[ind])
@test ts[ind].coredata == df_timetype_index[[ind], :]

# getindex(ts, r::UnitRange)
ind = 4:10
test_types(ts[ind])
@test ts[ind].coredata == df_timetype_index[ind, :]

# getindex(ts, a::AbstractVector{Int64})
ind = collect(4:10)
test_types(ts[ind])
@test ts[ind].coredata == df_timetype_index[ind, :]

# getindex(ts, d::Date)
ind = Date(2007, 2, 1)
test_types(ts[ind])
@test TSx.index(ts[ind]) == [ind]

# getindex(ts, y::Year)
ind = Year(2007)
test_types(ts[ind])
@test unique(Dates.year.(TSx.index(ts[ind]))) == [2007]

# getindex(ts, y::Year, m::Month)
y = Year(2007)
m = Month(3)
test_types(ts[y, m])
@test unique(Dates.yearmonth.(TSx.index(ts[y, m]))) == [(2007, 3)]

# getindex(ts, i::String)
ind = "2007-10-01"
d = Date(2007, 10, 1)
test_types(ts[ind])
@test TSx.index(ts[ind]) == [d]

# getindex(ts, i::Int, j::Int)
i = 1; j = 1
t = ts[i, j]
test_types(t)
@test t.coredata[!, :Index] == [df_timetype_index[i, 1]]
@test t.coredata[!, :data] == [df_timetype_index[i, j+1]]

# getindex(ts, i::UnitRange, j::Int)
i = 1:10; j = 1
t = ts[i, j]
test_types(t)
@test DataFrames.nrow(t.coredata) == length(i)
@test DataFrames.ncol(t.coredata) == length(j)+1 # Index + Ncols
@test t.coredata[!, :Index] == df_timetype_index[i, :Index]
@test t.coredata[!, :data] == df_timetype_index[i, :data]

# getindex(ts::TS, i::Int, j::UnitRange)
i = 2; j = 1:1
t = ts[i, j]
test_types(t)
@test DataFrames.nrow(t.coredata) == length(i)
@test DataFrames.ncol(t.coredata) == length(j)+1
@test t.coredata[!, :Index] == [df_timetype_index[i, :Index]]
@test t.coredata[!, :data] == df_timetype_index[[i], j]

# getindex(ts::TS, i::Int, j::Symbol)
i = 1; j = :data
t = ts[i, j]
test_types(t)
@test t.coredata[!, j] == df_timetype_index[i, j]

# getindex(ts::TS, i::UnitRange, j::Symbol)
i = 1:10; j = :data
t = ts[i, j]
test_types(t)
@test t.coredata[!, j] == df_timetype_index[i, j]

# getindex(ts::TS, i::Int, j::String)
i = 1; j = "data"
t = ts[i, j]
test_types(t)
@test t.coredata[!, j] == df_timetype_index[i, j]

# getindex(ts::TS, i::UnitRange, j::String)
i = 1:10; j = "data"
t = ts[i, j]
test_types(t)
@test t.coredata[!, j] == df_timetype_index[i, j]

# getindex(ts::TS, r::StepRange{T, V}) where {T<:TimeType, V<:Period}
ind = Date(2007, 1, 1):Day(1):Date(2007, 2, 1)
t = ts[ind]
test_types(ts[ind])
@test t.coredata[!, :Index] == collect(ind)


text:
include("./lib/probs.jl")
include("./lib/utils.jl")
include("./lib/pccf.jl")

module Sim2
# ;; SIMULATION TO DEMONSTRATE IMPROVEMENT DEPENDING ON THE THRESHOLD
# ;; PLOT: CALC FRACTION OF REMOVED FAs and PREDICTED TPs
using Pccf
using Probs

# 1
function generateErrors(ratePerMu, mu, delayConfirm, maxDelay)
	# Generate locations of errors within _cinfirmation_ delay interval uniformly on average rate-per-my errors per mu points
	# maxDelay ~ length of the signal
	avgNumOfErrors = div(maxDelay, mu) * ratePerMu
	filter(x -> x < delayConfirm, rand(1:maxDelay, Int(avgNumOfErrors)))
end

# 2
function fracOfRemovedFa(h, probs, delayConfirm, allErrors)
    # Given PCCF and threshold values function counts fraction of
    # the candidates which fall in the regions where probs are
    # lower than threshold and therefore we predict that they are
    # FA. If probs are higher than a threshold then we expect that
    # this is a change.
    idxsBelowH = Set(Probs.belowH(h, probs, delayConfirm))
    errors     = Set(filter(x -> x <= delayConfirm, allErrors))
    inter      = intersect(errors, idxsBelowH)
    n          = length(errors)
    if 1 == 0
            println("...")
            println("Errors within a confirmation delay:", errors)
            println("Pccf vals below h:", idxsBelowH)
            println("Intersection:", inter)
    end
    n > 0 ? length(inter) / n : 0
end

# 3
function fracOfPredictedTp(h, probs, delayConfirm, allChanges)
    idxsAboveH = Set(Probs.aboveH(h, probs, delayConfirm))
    chps       = Set(filter(x -> x <= delayConfirm, allChanges))
    inter      = intersect(chps, idxsAboveH)
    n          = length(chps)
    n > 0 ? length(inter)/n : 0
end

# 4
function genAndCountFa(h, probs, errRate, mu, delayConfirm, maxDelay)
    # Single simulation: given probs and threshold generate errors
    # locations (given erro rate per mu) count fraction of those who
    # are below threshold
    errorLocations=filter(x->x<=delayConfirm, generateErrors(errRate, mu, delayConfirm, maxDelay))
    length(errorLocations)==0 ? 0 : fracOfRemovedFa(h, probs, delayConfirm, errorLocations)
end

# 5
function genAndCountTp(h, probs, mu, sigma, n, delayConfirm, width)
    # Generate locations of true changes and count TP which are
    # coorectly predicted by the prediction intervals (probs > h)
    chpsLocations = Probs.generateChanges(n, mu, sigma, width)
    length(chpsLocations)==0 ? 0 : fracOfPredictedTp(h,probs, delayConfirm, chpsLocations)
end

# 6
function improvementVsDelay(mu, sigma, width, h, chpsnum, errRate, probs,
                            delayConfirm, maxDelay, numSim)
    # N simulations for a given delay; function returns averaged fractions
    removedFractionFA = Float64[]
    predictedTP       = Float64[]
    for i in 1:numSim
        frfa = genAndCountFa(h, probs, errRate, mu, delayConfirm, maxDelay)
        frtp = genAndCountTp(h, probs, mu, sigma, chpsnum, delayConfirm, width)
        push!(removedFractionFA, frfa)
        push!(predictedTP, frtp)
    end
    return mean(removedFractionFA), mean(predictedTP)
end

# 7
function improvementVsDelayMulti(mu, sigma, width, chpsnum, errRate, delayFrom,
                                 delayTo, numSim, outfile, mainpath)
    # Pairs [delay vs improvement in fa vs improvement in tp ]
    threshold = errRate/mu
    probs = Pccf.pccf(width*chpsnum, mu, sigma) # ?
    delVsImpr = zeros(delayTo-delayFrom+1, 3)
    for d in delayFrom:delayTo
        nextAvgImpr = improvementVsDelay(mu, sigma, width, threshold, chpsnum, errRate, probs, d, delayTo, numSim)
        delVsImpr[d-delayFrom+1, :] = [d, nextAvgImpr[1], nextAvgImpr[2]]
    end
    fullpath = mainpath * outfile
    writedlm_header(fullpath, delVsImpr, "Delay, CorrectFA, CorrectTP")
    println("... Results are written into: ", fullpath)
end

function run()
    s = Dict("mu"=>10.0, "sigma"=>0.7, "width"=>20, "chpsnum"=>50,
             "errrate"=>0.5, "numSim"=>500, "minDealay"=>5, "maxDelay"=>300)
    mpath="D:/MyTemp/gitlocal/Pccf/r-code/out/"

    improvementVsDelayMulti(s["mu"], s["sigma"], s["width"], s["chpsnum"],
                            s["errrate"], s["minDealay"], s["maxDelay"],
                            s["numSim"], "sim2_a.txt", mpath)

    improvementVsDelayMulti(s["mu"], s["sigma"], s["width"], s["chpsnum"],
                            2           , s["minDealay"], s["maxDelay"],
                            s["numSim"], "sim2_b.txt", mpath)
end

end # module

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
# Julia script to compute robsutness coefficients -
include("SolveBalances.jl")
using PyPlot
using PyCall
@pyimport numpy as np

# main: calculates the ratio of the AUC in the perturbed state, to the base state
function main(readout_index,perturbed_index_array,delta_array) # what we want to look at, how we want to perturb, difference

  # Load the ensemble files from disk -
  pc_array_full = readdlm("./data/pc_array_O1_O2.dat")
  rank_array = readdlm("./data/rank_array_O1_O2.dat")

  # Select the desired rank -
  idx_rank = find(rank_array .<= 5.0)

  # initialized the robsutness array -
  alpha_array = Float64[]

  # Run the simulation for this rank -
  number_of_samples = 1
  for (index,rank_index_value) in enumerate(idx_rank)

    if (mod(index,1)==0)

      # Grab the parameter set from the cache -
      parameter_array = pc_array_full[:,rank_index_value]

      # Calculate the AUC for the base state -
      AUC_BASE = calculate_auc_base(parameter_array,readout_index)

      # Calculate the AUC for the perturbed case -
      AUC_PERTURBED = calculate_auc_perturbed(parameter_array,readout_index,perturbed_index_array,delta_array)

      @show AUC_BASE,AUC_PERTURBED

      alpha_value = log10(abs(AUC_PERTURBED/AUC_BASE))
      push!(alpha_array,alpha_value)

    end
  end

  return alpha_array
end

function calculate_auc_base(parameter_array,readout_index)

  tStart = 0.0;
  tStop = 25.0;
  tStep = 0.1;

  # Load base experiment
  # will need to rerun to steady state and load a blank/default data dictionary
  # need to initialize each model using a parameter array from POETs

  data_dictionary = Dict()
  initial_condition_array = zeros(19)
  initial_condition_array[1] = 1.0      # 1 zymosan
  initial_condition_array[2] = 1.9      # 2 C4
  initial_condition_array[3] = 0.322    # 3 C2
  initial_condition_array[8] = 7.57     # 4 C3
  initial_condition_array[14] = 0.195   # 5 C5
  initial_condition_array[17] = 2.23    # 6 Factor H
  initial_condition_array[18] = 0.417   # 7 C4BP
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  data_dictionary["PARAMETER_ARRAY"] = parameter_array

  # Run the model -
  (t,x) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # Calculate the AUC for readout_index -
  return trapz(t,x[:,readout_index])
end

function calculate_auc_perturbed(parameter_array,readout_index,perturbed_index_array,delta_array)
  # perturbed_index_array -- here we are selecting a species initial condition to perturb but I will
  # need to send the simulation files
  # Need to do calculate the AOC from the simulations twice
  # Because we need both ecadherin robustness and vimentin robustness
  # Sample 100 sets (mod 100)
  # at end just plot a satter (log)
  # keep in mind: definition of robustness - can cut off the smalls or run it with a 1

  tStart = 0.0;
  tStop = 25.0;
  tStep = 0.1;

  data_dictionary = Dict()
  initial_condition_array = zeros(19)
  initial_condition_array[1] = 1.0    # 1 zymosan
  initial_condition_array[2] = 1.9      # 2 C4
  initial_condition_array[3] = 0.322    # 3 C2
  initial_condition_array[8] = 7.57     # 4 C3
  initial_condition_array[14] = 0.195   # 5 C5
  initial_condition_array[17] = 2.23    # 6 Factor H
  initial_condition_array[18] = 0.417   # 7 C4BP

  # Issue the perturbation -
  for (index,perturbed_index) in enumerate(perturbed_index_array)
    delta = delta_array[index]
    initial_condition_array[perturbed_index] = delta*initial_condition_array[perturbed_index]
  end

  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  data_dictionary["PARAMETER_ARRAY"] = parameter_array

  # Run the model -
  (t,x) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # Calculate the AUC for readout_index -
  return trapz(t,x[:,readout_index])
end

function trapz{Tx<:Number, Ty<:Number}(x::Vector{Tx}, y::Vector{Ty})
    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    #= correction -h^2/12 * (f'(b) - f'(a))
    ha = x[2] - x[1]
    he = x[end] - x[end-1]
    ra = (y[2] - y[1]) / ha
    re = (y[end] - y[end-1]) / he
    r/2 - ha*he/12 * (re - ra)
    =#
    return r/2
end


text:
testfile_specs = Dict(EvenTimeSeries =>
                      [Dict("fname" => "./test-files/delta-etime.sac",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1.0; x; end,
                            "ascii" => false),
                       Dict("fname" => "./test-files/be-delta-etime.sac",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1.0; x; end,
                            "ascii" => false),
                       Dict("fname" => "./test-files/delta-etime.txt",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1.0; x; end,
                            "ascii" => true)],
                      UnevenTimeSeries =>
                      [Dict("fname" => "./test-files/delta-utime.sac",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1; x; end,
                            "data2" => log(Float32[1:100;]),
                            "ascii" => false),
                       Dict("fname" => "./test-files/delta-utime.txt",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1; x; end,
                            "data2" => log(Float32[1:100;]),
                            "ascii" => true)],
                      AmplitudeSpectrum =>
                      [Dict("fname" => "./test-files/delta-amph.sac",
                            "data1" => ones(Float32, 128),
                            "data2" => zeros(Float32, 128),
                            "ascii" => false),
                       Dict("fname" => "./test-files/delta-amph.txt",
                            "data1" => ones(Float32, 128),
                            "data2" => zeros(Float32, 128),
                            "ascii" => true)],
                      ComplexSpectrum =>
                      [Dict("fname" => "./test-files/delta-rlim.sac",
                            "data1" => fill(1+0im, 128),
                            "ascii" => false),
                       Dict("fname" => "./test-files/delta-rlim.txt",
                            "data1" => fill(1+0im, 128),
                            "ascii" => true)],
                      GeneralXY =>
                      [Dict("fname" => "./test-files/delta-xy.sac",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1.0; x; end,
                            "data2" => log(Float32[1:100;]),
                            "ascii" => false),
                       Dict("fname" => "./test-files/delta-xy.txt",
                            "data1" => begin x = zeros(Float32, 100); x[50] = 1.0; x; end,
                            "data2" => log(Float32[1:100;]),
                            "ascii" => true)])


text:
abstract type AtomicStructure end

struct BCC <: AtomicStructure end
struct FCC <: AtomicStructure end

# TODO: introduce type Slipsystem instead of using NTuple

#slip planes and slip directions for BCC
function _get_slipsystems(::FCC)
    slip_planes = [ [ 1., 1., 1.],
                    [ 1., 1.,-1.],
                    [ 1.,-1.,-1.],
                    [ 1.,-1., 1.]]
    slip_directions = [ ([ 1., 0.,-1.], [ 0., 1.,-1.], [ 1.,-1., 0.]);
                        ([ 0., 1., 1.], [ 1., 0., 1.], [ 1.,-1., 0.]);
                        ([ 1., 0., 1.], [ 0., 1.,-1.], [ 1., 1., 0.]);
                        ([ 0., 1., 1.], [ 1., 0.,-1.], [ 1., 1., 0.]) ]
    return slip_planes, slip_directions
end

function _get_slipsystems(::BCC)
    slip_planes = [ [ 1., 1., 0.],
                    [ 1.,-1., 0.],
                    [ 1., 0., 1.],
                    [-1., 0.,-1.],
                    [ 0., 1., 1.],
                    [ 0., 1.,-1.]]
    slip_directions = [ ([ 1.,-1., 1.], [ 1.,-1.,-1.]);
                        ([ 1., 1., 1.], [ 1., 1.,-1.]);
                        ([-1., 1., 1.], [-1.,-1., 1.]);
                        ([ 1., 1., 1.], [ 1.,-1., 1.]);
                        ([ 1., 1.,-1.], [-1., 1.,-1.]);
                        ([-1., 1., 1.], [ 1., 1., 1.]) ]
    return slip_planes, slip_directions
end

function slipsystems(as::AtomicStructure, R::Rotations.Rotation)
    # generate rotated slip systems
    slip_systems = NTuple{2, Tensor{1,3}}[]
    slip_planes, slip_directions = _get_slipsystems(as)
    for (i, sp) in enumerate(slip_planes)
        for sd in slip_directions[i]
            push!(slip_systems, (Tensor{1,3}(R*sp), Tensor{1,3}(R*sd)))
        end
    end

    normalized_slip_systems = Vector{NTuple{2, Tensor{1,3,Float64,3}}}(undef, length(slip_systems))
    for i=1:length(slip_systems)
        nss = (slip_systems[i][1] / norm(slip_systems[i][1]), slip_systems[i][2] / norm(slip_systems[i][2]))
        normalized_slip_systems[i] = nss
    end
    return normalized_slip_systems
end


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("lattice", "singer"), x=:VoicePart, y=:Height, color=:VoicePart, Geom.violin)


text:
using LinearAlgebra, ForwardDiff, IterativeSolvers

function Pκ(z)
    z0 = z[1]
    z̄ = z[2:end]
    n = length(z)

    η1 = z0 - norm(z̄)
    η2 = z0 + norm(z̄)

    w̄ = ones(n - 1) ./ norm(ones(n - 1))

    if all(z̄ .!= 0.0)
        u1 = 0.5 * [1.0; -z̄ ./ norm(z̄)]

        u2 = 0.5 * [1.0; z̄ ./ norm(z̄)]
    else
        u1 = 0.5 * [1.0; -w̄]
        u2 = 0.5 * [1.0; w̄]
    end

    return max(0.0, η1) * u1 + max(0.0, η2) * u2
end

function ∇Pκ(z)
    z0 = z[1]
    z̄ = z[2:end]
    n = length(z)

    if z0 != norm(z̄) && z0 != norm(z̄)
        if z0 < -norm(z̄)
            @warn "projection weirdness"
            return zeros(n, n)
        elseif z0 > norm(z̄)
            return Diagonal(ones(n))
        elseif -norm(z̄) < z0 && z0 < norm(z̄)
            w̄ = z̄ ./ norm(z̄)
            H = ((1.0 + z0 / norm(z̄)) * Diagonal(ones(n-1))
                - z0 / norm(z̄) * w̄ * w̄')
            return 0.5 * [1.0 w̄'; w̄ H]
        end

    else
        if any(z̄ .!= 0.0) && z0 == norm(z̄)
            V1 = Diagonal(ones(n))

            w̄ = z̄ ./ norm(z̄)
            H = (2.0 * Diagonal(ones(n-1))
                - w̄ * w̄')
            V2 = 0.5 * [1.0 w̄'; w̄ H]
            return V2
        elseif any(z̄ .!= 0.0) && z0 == -norm(z̄)
            V1 = zeros(n, n)

            w̄ = z̄ ./ norm(z̄)
            H = - w̄ * w̄'
            V2 = 0.5 * [1.0 w̄'; w̄ H]
            return V2
        elseif all(z̄ .== 0.0) && z0 .== 0.0
            V1 = zeros(n, n)
            V2 = Diagonal(ones(n))

            w̄ = ones(n - 1) ./ norm(ones(n - 1))
            ρ = 0.5
            H = ((1.0 + ρ) * Diagonal(ones(n-1))
                - ρ * w̄ * w̄')
            V3 = 0.5 * [1.0 w̄'; w̄ H]

            return V3
        else
            @error "projection jacobian weirdness"
        end
    end
end

z = zeros(3)

Pκ(z)
norm(∇Pκ(z) - ForwardDiff.jacobian(Pκ, z))

# example 3.6
v = ones(2)
y = 1.0
c = [0.0; v]
A = [1.0 0.0 0.0]
b = [y]
f(x) = c' * x

m = 1
# A = zeros(m, n)
# b = zeros(m)

function M(z)
    x = z[1:n]
    λ = z[n .+ (1:n)]
    μ = z[2n .+ (1:m)]

    return [ForwardDiff.gradient(f, x) - A' * μ - λ;
            A * x - b;
            x - Pκ(x - λ)]
end

function ∇M(z)
    x = z[1:n]
    λ = z[n .+ (1:n)]
    μ = z[2n .+ (1:m)]

    return [ForwardDiff.hessian(f, x) -A' -Diagonal(ones(n));
            A zeros(m, m) zeros(m, n)
            (Diagonal(ones(n)) - ∇Pκ(x - λ)) zeros(n, m) ∇Pκ(x - λ)]
end

z = rand(2n + m)
M(z)
rank(∇M(z))

function solve()
    x = zeros(n)
    λ = zeros(n)
    μ = zeros(m)
    z = [x; λ; μ]

    ρ = 1.0e-5

    for i = 1:10
        r = M(z)
        println("norm: $(norm(r))")

        if norm(r) < 1.0e-8
            return z
        end

        ∇r = ∇M(z)
        Δ = -1.0 * (∇r' * ∇r + ρ * I) \ (∇r' * r)
        # Δ = -1.0 * ∇r \ r
        # Δ = gmres(∇r, -1.0 * r)

        α = 1.0
        iter = 0
        while norm(M(z + α * Δ)) > (1.0 - 0.001 * α) * norm(r)
            α = 0.5 * α
            iter += 1

            if iter > 25
                @error "line search failed"
            end
        end

        z .+= α * Δ

    end

    return z
end

z_sol = solve()


text:
# Least-squares RTM of the BP synthetic 2007 data set
# Author: Philipp Witte
# Date: March 2018
#

# TO DO
# Set up path where data will be saved
data_path = "/path/to/data/"

using JUDI, SegyIO, JLD, PyPlot, JOLI, Random, LinearAlgebra

# Load velocity model(replace with correct paths)
if !isfile("bp_synthetic_2004_migration_velocity.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_migration_velocity.jld`)
end
vp = load(join([pwd(), "/bp_synthetic_2004_migration_velocity.jld"]))["vp"] / 1f3

if !isfile("bp_synthetic_2004_water_bottom.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_water_bottom.jld`)
end
water_bottom = load(join([pwd(), "/bp_synthetic_2004_water_bottom.jld"]))["wb"]

# Set up model structure
d = (6.25, 6.25)
o = (0., 0.)
m0 = (1f0 ./ vp).^2
n = size(m0)
model0 = Model(n, d, o, m0)

# Scan directory for segy files and create out-of-core data container
container = segy_scan(data_path, "bp_observed_data", ["GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container; segy_depth_key = "SourceDepth")

# Set up source
src_geometry = Geometry(container; key = "source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.020)  # 27 Hz peak frequency
q = judiVector(src_geometry, wavelet)

# Info structure for linear operators
ntComp = get_computational_nt(src_geometry, d_obs.geometry, model0)    # no. of computational time steps
info = Info(prod(model0.n), d_obs.nsrc, ntComp)

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

# Set options
opt = Options(limit_m = true,
              buffer_size = 3000f0,
              isic = true,
			  dft_subsampling_factor=8)

# Setup operators
F = judiModeling(info, model0, q.geometry, d_obs.geometry; options=opt)
J = judiJacobian(F, q)

# Right-hand preconditioners
D = judiDepthScaling(model0)
T = judiTopmute(model0.n, (1 .- water_bottom), [])
Mr = D*T

# Linearized Bregman parameters
x = zeros(Float32, info.n)
z = zeros(Float32, info.n)
batchsize = 200
niter = 20
nfreq = 20
fval = zeros(Float32, niter)
q_dist = generate_distribution(q)
J.options.frequencies = Array{Any}(undef, d_obs.nsrc)

# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64,lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32,lambda), 0f0)
C = joCurvelet2D(model0.n[1], model0.n[2]; zero_finest=true, DDT=Float32, RDT=Float64)
lambda = []
t = []  # t = 1f-4

# Main loop
for j=1:niter
    println("Iteration: ", j)

   # Set randomized frequencies
    for k=1:d_obs.nsrc
        J.options.frequencies[k] = select_frequencies(q_dist; fmin=0.003, fmax=0.04, nf=nfreq)
    end

    # Select batch and set up left-hand preconditioner
    i = randperm(d_obs.nsrc)[1:batchsize]
    d_sub = get_data(d_obs[i])
    Ml = judiMarineTopmute2D(35, d_sub.geometry)

    # Compute residual and estimate source
    if j > 1
        d_pred = J[i]*Mr*x
        r = Ml*d_pred - Ml*d_sub
    else
        r = Ml*d_sub*(-1f0)    # skip forward modeling in first iteration
    end

    # Residual and gradient
    g = adjoint(Mr)*adjoint(J[i])*adjoint(Ml)*r

    # Step size and update variable
    fval[j] = .5*norm(r)^2
    t = norm(r)^2/norm(g)^2 # divide by 10
    println("    Stepsize: ", t)

    j==1 && (global lambda = 0.03*norm(C*t*g, Inf))   # estimate thresholding parameter in 1st iteration

    # Update variables
    global z -= t*g
    global x = adjoint(C)*soft_thresholding(C*z, lambda)

    # Save snapshot
    save(join([path, "/results/splsrtm_freq_iteration_", string(j), ".jld"]), "x", reshape(x, model0.n), "z", reshape(z, model0.n), "t", t, "lambda", lambda, "fval", fval[j])
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using ModelingToolkit, LightGraphs

# use a ReactionSystem to generate systems for testing
@parameters k1 k2 t
@variables S(t) I(t) R(t)

rxs = [Reaction(k1, nothing, [S]),
       Reaction(k1, [S], nothing),
       Reaction(k2, [S,I], [I], [1,1], [2]),
       Reaction(k2, [S,R], [S], [2,1], [2]),
       Reaction(k1*I, nothing, [R]),
       Reaction(k1*k2/(1+t), [S], [R])]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])

#################################
#  testing for Jumps / all dgs
#################################
js = convert(JumpSystem, rs)
S  = convert(Variable,S); I = convert(Variable,I); R = convert(Variable,R)
k1 = convert(Variable,k1); k2 = convert(Variable,k2)

# eq to vars they depend on
eq_sdeps  = [Variable[], [S], [S,I], [S,R], [I], [S]]
eq_sidepsf = [Int[], [1], [1,2], [1,3], [2], [1]]
eq_sidepsb = [[2,3,4,6], [3,5],[4]]
deps = equation_dependencies(js)
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(rxs))
depsbg = asgraph(js)
@test depsbg.fadjlist == eq_sidepsf
@test depsbg.badjlist == eq_sidepsb

# eq to params they depend on
eq_pdeps    = [[k1],[k1],[k2],[k2],[k1],[k1,k2]]
eq_pidepsf = [[1],[1],[2],[2],[1],[1,2]]
eq_pidepsb = [[1,2,5,6],[3,4,6]]
deps = equation_dependencies(js, variables=parameters(js))
@test all(i -> isequal(Set(eq_pdeps[i]),Set(deps[i])), 1:length(rxs))
depsbg2 = asgraph(js, variables=parameters(js))
@test depsbg2.fadjlist == eq_pidepsf
@test depsbg2.badjlist == eq_pidepsb

# var to eqs that modify them
s_eqdepsf = [[1,2,3,6],[3],[4,5,6]]
s_eqdepsb = [[1],[1],[1,2],[3],[3],[1,3]]
ne        = 8
bg        = BipartiteGraph(ne, s_eqdepsf, s_eqdepsb)
deps2     = variable_dependencies(js)
@test isequal(bg,deps2)

# eq to eqs that depend on them
eq_eqdeps = [[2,3,4,6],[2,3,4,6],[2,3,4,5,6],[4],[4],[2,3,4,6]]
dg = SimpleDiGraph(6)
for (eqidx,eqdeps) in enumerate(eq_eqdeps)
    for eqdepidx in eqdeps
       add_edge!(dg, eqidx, eqdepidx)
    end
end
dg3 = eqeq_dependencies(depsbg,deps2)
@test dg == dg3

# var to vars that depend on them
var_vardeps = [[1,2,3],[1,2,3],[3]]
ne = 7
dg = SimpleDiGraph(3)
for (vidx,vdeps) in enumerate(var_vardeps)
    for vdepidx in vdeps
       add_edge!(dg, vidx, vdepidx)
    end
end
dg4 = varvar_dependencies(depsbg,deps2)
@test dg == dg4

#####################################
#       testing for ODE/SDEs
#####################################
os = convert(ODESystem, rs)
deps = equation_dependencies(os)
eq_sdeps  = [[S,I], [S,I], [S,I,R]]
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(deps))

sdes = convert(SDESystem, rs)
deps = equation_dependencies(sdes)
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(deps))

deps = variable_dependencies(os)
s_eqdeps = [[1],[2],[3]]
@test deps.fadjlist == s_eqdeps

#####################################
#       testing for nonlin sys
#####################################
@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ ρ-y,
       0 ~ y - β*z]
ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
deps = equation_dependencies(ns)
eq_sdeps = [[x,y],[y],[y,z]]
@test all(i -> isequal(Set(deps[i]),Set(convert.(Variable,eq_sdeps[i]))), 1:length(deps))


text:
@testset "medianvariance.jl" begin
    @testset "Constructors" begin
        @test_throws ErrorException MedianVarianceBinning(-1)
        @test_throws ErrorException MedianVarianceBinning(0)
        @test_throws ErrorException MedianVarianceBinning(10, -1)
        @test_throws ErrorException MedianVarianceBinning(10, 0)
    end

    @testset "Basic tests ($nclasses classes)" for nclasses in (2, 10, 100)
        nsamples = 1_000
        dist = Dirichlet(nclasses, 1.0)
        predictions = [rand(dist) for _ in 1:nsamples]
        targets = rand(1:nclasses, nsamples)

        # set minimum number of samples
        for minsize in (1, 10, 100, 500, 1_000)
            bins = @inferred(
                CalibrationErrors.perform(
                    MedianVarianceBinning(minsize), predictions, targets
                )
            )

            @test all(bin -> bin.nsamples ≥ minsize, bins)
            @test all(bin -> sum(bin.mean_predictions) ≈ 1, bins)
            @test all(bin -> sum(bin.proportions_targets) ≈ 1, bins)

            @test sum(bin -> bin.nsamples, bins) == nsamples
            @test sum(bin -> bin.nsamples .* bin.mean_predictions, bins) ≈ sum(predictions)
            @test sum(bin -> bin.nsamples .* bin.proportions_targets, bins) ≈
                  counts(targets, nclasses)
        end

        # set maximum number of bins
        for maxbins in (1, 10, 100, 500, 1_000)
            bins = @inferred(
                CalibrationErrors.perform(
                    MedianVarianceBinning(1, maxbins), predictions, targets
                )
            )

            @test length(bins) ≤ maxbins

            @test all(bin -> bin.nsamples ≥ 1, bins)
            @test all(bin -> sum(bin.mean_predictions) ≈ 1, bins)
            @test all(bin -> sum(bin.proportions_targets) ≈ 1, bins)

            @test sum(bin -> bin.nsamples, bins) == nsamples
            @test sum(bin -> bin.nsamples .* bin.mean_predictions, bins) ≈ sum(predictions)
            @test sum(bin -> bin.nsamples .* bin.proportions_targets, bins) ≈
                  counts(targets, nclasses)
        end
    end

    @testset "Simple example" begin
        predictions = [[0.4, 0.1, 0.5], [0.5, 0.3, 0.2], [0.3, 0.7, 0.0]]
        targets = [1, 2, 3]
        # maximum possible steps in the order they might occur:
        # first step: [1], [2, 3] -> create bin with [1]
        # second step: [2], [3] -> create bins with [2] and [3]

        bins = CalibrationErrors.perform(MedianVarianceBinning(1), predictions, targets)
        @test length(bins) == 3
        @test all(bin -> bin.nsamples == 1, bins)
        for (i, idx) in enumerate((1, 2, 3))
            @test bins[i].mean_predictions == predictions[idx]
            @test bins[i].proportions_targets == Matrix{Float64}(I, 3, 3)[:, targets[idx]]
        end

        bins = CalibrationErrors.perform(MedianVarianceBinning(2), predictions, targets)
        @test length(bins) == 1
        @test bins[1].nsamples == 3
        @test bins[1].mean_predictions == mean(predictions)
        @test bins[1].proportions_targets == [1 / 3, 1 / 3, 1 / 3]

        predictions = [
            [0.4, 0.1, 0.5],
            [0.5, 0.3, 0.2],
            [0.3, 0.7, 0.0],
            [0.1, 0.0, 0.9],
            [0.8, 0.1, 0.1],
        ]
        targets = [1, 2, 3, 1, 2]
        # maximum possible steps in the order they might occur:
        # first step: [3, 5], [1, 2, 4]
        # second step: [5], [3], [1, 2, 4] -> create bins with [5] and [3]
        # third step:  [2], [1, 4] -> create bin with [2]
        # fourth step: [1], [4] -> create bins with [1] and [4]

        bins = CalibrationErrors.perform(MedianVarianceBinning(1), predictions, targets)
        @test length(bins) == 5
        @test all(bin -> bin.nsamples == 1, bins)
        for (i, idx) in enumerate((5, 3, 2, 1, 4))
            @test bins[i].mean_predictions == predictions[idx]
            @test bins[i].proportions_targets == Matrix{Float64}(I, 3, 3)[:, targets[idx]]
        end

        bins = CalibrationErrors.perform(MedianVarianceBinning(2), predictions, targets)
        @test length(bins) == 2
        @test all(bin -> sum(bin.mean_predictions) ≈ 1, bins)
        @test all(bin -> sum(bin.proportions_targets) ≈ 1, bins)
        for (i, idxs) in enumerate(([3, 5], [1, 2, 4]))
            @test bins[i].nsamples == length(idxs)
            @test bins[i].mean_predictions ≈ mean(predictions[idxs])
            @test bins[i].proportions_targets ==
                  vec(mean(Matrix{Float64}(I, 3, 3)[:, targets[idxs]]; dims=2))
        end

        bins = CalibrationErrors.perform(MedianVarianceBinning(3), predictions, targets)
        @test length(bins) == 1
        @test bins[1].nsamples == 5
        @test bins[1].mean_predictions == mean(predictions)
        @test bins[1].proportions_targets == [0.4, 0.4, 0.2]
    end
end


text:
using Fed
using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Momentum

struct Client
    model::Chain
    dataset::DataLoader
    optimizer::Momentum
    epochs::Int

    Client(dataset) = new(
        newmodel(),
        dataset,
        Momentum(0.01, 0.5),
        5 # epochs
    )
end

"""
    start_client(dataset, port)

Starts one client with the given dataset on the given localhost port.
"""
function start_client(dataset, port::Int, store_port::Int)
    client = Client(dataset)
    config = newconfig(store_port, true)

    host = "127.0.0.1"
    train_hook = (weights::Vector{Float32}) -> fit(client, weights)

    node = Fed.Client.Node(host, port, train_hook, config)

    @info "Client started on [http://$host:$port]"
    Fed.Client.start(node)
end

"""
    start_client()

Starts 100 clients asynchronously adn return the tasks running them.
"""
function start_clients(num_clients::Int)
    num_sample = 60000
    step = Int(floor(num_sample / num_clients))

    tasks = [
        @async start_client(get_traindata((i * step + 1):((i+1) * step), isflat=true), 8081+i, 9091+i)
        for i in 0:num_clients-1
    ]

    return tasks
end

"""
    stop_clients(tasks)

Stops the given `tasks` which are running the clients.
"""
function stop_clients(tasks)
    for task in tasks
        @async Base.throwto(task, InterruptException())
    end
end


text:
using BenchmarkTools
using Statistics
using StatsBase: sample, quantile
using Distributions
using Random
using Plots
using Revise
using EvoTrees
using EvoTrees: sigmoid, logit
# using ProfileView

# prepare a dataset
Random.seed!(123)
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(T=Float64,
    loss=:linear, metric=:mse,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25);
# 67.159 ms (77252 allocations: 28.06 MiB)

# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 25, metric=:mae)
@time pred_train_linear = predict(model, X_train)
mean(abs.(pred_train_linear .- Y_train))
sqrt(mean((pred_train_linear .- Y_train) .^ 2))

# logistic / cross-entropy
params1 = EvoTreeRegressor(T=Float64,
    loss=:logistic, metric = :logloss,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# 218.040 ms (123372 allocations: 34.71 MiB)
# @btime model = fit_evotree($params1, $X_train, $Y_train, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train_logistic = predict(model, Float32.(X_train))
sqrt(mean((pred_train_logistic .- Y_train) .^ 2))

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, msize = 1, mcolor = "gray", mswidth=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
plot!(X_train[:,1][x_perm], pred_train_logistic[x_perm], color = "darkred", linewidth = 1.5, label = "Logistic")
# plot!(X_train[:,1][x_perm], pred_train_poisson[x_perm], color = "green", linewidth = 1.5, label = "Poisson")
# plot!(X_train[:,1][x_perm], pred_train_L1[x_perm], color = "pink", linewidth = 1.5, label = "L1")
savefig("figures/regression_sinus_gpu.png")

###############################
## gaussian
###############################
params1 = EvoTreeGaussian(T=Float64,
    loss=:gaussian, metric=:gaussian,
    nrounds=200, nbins=64,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0, rng=123,
    device = "gpu")

@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n = 10);
# @time model = fit_evotree(params1, X_train, Y_train, print_every_n = 10);
@time pred_train_gaussian = EvoTrees.predict(model, X_train)

pred_gauss = [Distributions.Normal(pred_train_gaussian[i,1], pred_train_gaussian[i,2]) for i in 1:size(pred_train_gaussian,1)]
pred_q80 = quantile.(pred_gauss, 0.8)
pred_q20 = quantile.(pred_gauss, 0.2)

mean(Y_train .< pred_q80)
mean(Y_train .< pred_q20)

x_perm = sortperm(X_train[:,1])
plot(X_train[:, 1], Y_train, ms = 1, mcolor = "gray", mswidth=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_gaussian[x_perm, 1], color = "navy", linewidth = 1.5, label = "mu")
plot!(X_train[:,1][x_perm], pred_train_gaussian[x_perm, 2], color = "darkred", linewidth = 1.5, label = "sigma")
plot!(X_train[:,1][x_perm], pred_q20[x_perm, 1], color = "green", linewidth = 1.5, label = "q20")
plot!(X_train[:,1][x_perm], pred_q80[x_perm, 1], color = "green", linewidth = 1.5, label = "q80")
savefig("figures/gaussian_sinus_gpu.png")


text:
module Prim

using ..Omega
using ..Omega: Ω, RandVar, URandVar, MaybeRV, ID, lift, uid, elemtype, isconstant
import ..Omega: params, name, ppapl, apl, reify
import Statistics: mean, var, quantile
import ..Causal: ReplaceRandVar
using ..Util
using Spec
import Distributions
const Djl = Distributions
import Base: minimum, maximum
import Random
using DocStringExtensions: SIGNATURES

export  bernoulli,
        betarv,
        β,
        categorical,
        # dirichlet,
        exponential,
        gammarv,
        Γ,
        invgamma,
        kumaraswamy,
        logistic,
        # mvnormal,
        normal,
        poisson,
        rademacher,
        uniform,

        mean

"Primitive random variable of known distribution"
abstract type PrimRandVar <: RandVar end

"Name of a distribution"
function name end

name(t::T) where {T <: PrimRandVar} = T.name.name

"Parameters of `rv`"
@generated function params(rv::PrimRandVar)
  fields = [Expr(:., :rv, QuoteNode(f)) for f in fieldnames(rv) if f !== :id]
  Expr(:tuple, fields...)
end

ppapl(rv::PrimRandVar, ωπ) = rvtransform(rv)(ωπ, reify(ωπ, params(rv))...)

@generated function anysize(args::Union{<:AbstractArray, Real}...)
  isarr = (arg -> arg <: AbstractArray).([args...])
  firstarr = findfirst(isarr)
  if isempty(isarr)
    :(())
  else
    :(size(args[$firstarr]))
  end
end
@spec same(size.(filter(a -> a isa AbstractArray, args)))

include("univariate.jl")      # Univariate Distributions
include("multivariate.jl")    # Multivariate Distributions
include("statistics.jl")      # Distributional properties: mean, variance, etc
include("randx.jl")      # Distributional properties: mean, variance, etc
export  succprob,
        failprob,
        maximum,
        minimum,
        islowerbounded,
        isupperbounded,
        isbounded,
        std,
        median,
        mode,
        modes,

        skewness,
        kurtosis,
        isplatykurtic,
        ismesokurtic,

        isleptokurtic,
        entropy,
        mean,
        samplemean,
        samplemeanᵣ,
        sampleprob,
        sampleprobᵣ,
        prob,
        lprob

# Lifted distributional functions
export  lsuccprob,
        lfailprob,
        lmaximum,
        lminimum,
        lislowerbounded,
        lisupperbounded,
        lisbounded,
        lstd,
        lmedian,
        lmode,
        lmodes,

        lskewness,
        lkurtosis,
        lisplatykurtic,
        lismesokurtic,

        lisleptokurtic,
        lentropy,
        lmean
include("djl.jl")             # Distributions.jl interop

end


text:
"""
Module to access nflfastR-data files

See https://github.com/nflverse/nflfastR-data for additional data documentation.
"""
module NFLFastR

using  CSV: CSV
using  DataFrames: DataFrame
using  Downloads: Downloads
using  Logging
using  Scratch: @get_scratch!

export getdata

"""
seasons with nflfastR data
"""
const SEASONS = tuple(1999:2021...)

"""
The download root of the nflfastR-data GitHub repo.
"""
const REPOROOT = "https://github.com/nflverse/nflfastR-data/raw/master/"

function getfilepath(season::Integer)::String
    season in SEASONS || error("Invalid season $season")
    return "play_by_play_$(season).csv.gz"
end
function getfilepath(symbol::Symbol)::String
    symbol === :player_stats || symbol === :player_stats_kicking || error("invalid symbol: $symbol")
    return "$(symbol).csv.gz"
end

function load_data_from_disk(filepath; redownload=false, reporoot=REPOROOT)
    dirpath = @get_scratch!("nflfastR")
    localpath = joinpath(dirpath, filepath)
    if redownload || !isfile(localpath)
        remotepath = joinpath(reporoot, "data", filepath)
        @info "downloading..." remotepath localpath
        Downloads.download(remotepath, localpath)
    end
    return DataFrame(CSV.File(localpath, missingstring="NA"))
end

"""
    getdata(season::Integer)
    getdata(label::Symbol)

Create a dataframe of play-by-play data for `season` or `label` data for all seasons`.

# Examples
```jldoctest
julia> df = NFLFastR.getdata(2010);

julia> show(first(df[:, [:home_team, :away_team, :yardline_100, :half_seconds_remaining]], 5), eltypes=false)
5×4 DataFrame
 Row │ home_team  away_team  yardline_100  half_seconds_remaining
─────┼────────────────────────────────────────────────────────────
   1 │ LA         ARI             missing                    1800
   2 │ LA         ARI                  30                    1800
   3 │ LA         ARI                  78                    1795
   4 │ LA         ARI                  78                    1764
   5 │ LA         ARI                  73                    1723
```
"""
getdata(label; redownload=false, reporoot=REPOROOT) = load_data_from_disk(getfilepath(label); redownload, reporoot)

end  # module NFLFastR


text:
module REST

using HTTP
using MbedTLS
using Dates
using Base64

const API_VER = "2016-05-31"

const DEFAULT_KWARGS = Dict(:retries=>0, :readtimeout=>300, :status_exception=>false)

mutable struct StandardHeaders
    content_encoding::Union{String,Nothing}
    content_language::Union{String,Nothing}
    content_length::Union{String,Nothing}
    content_md5::Union{String,Nothing}
    content_type::Union{String,Nothing}
    date::Union{String,Nothing}
    if_modified_since::Union{String,Nothing}
    if_match::Union{String,Nothing}
    if_none_match::Union{String,Nothing}
    if_unmodified_since::Union{String,Nothing}
    range::Union{String,Nothing}
end

function StandardHeaders(; kwargs...)
    obj = StandardHeaders(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
    for (n,v) in kwargs
        setfield!(obj, n, v)
    end
    obj
end

struct ServiceRequest
    account::String
    verb::String
    resource::String
    standard_headers::StandardHeaders
    headers::Dict{String,String}

    function ServiceRequest(account::String, verb::String, resource::String; headers...)
        hdrdate = Dates.format(now(Dates.UTC), Dates.RFC1123Format)
        endswith(hdrdate, "GMT") || (hdrdate *= " GMT")

        reqhdrs = Dict{String,String}("x-ms-date"=>hdrdate, "x-ms-version"=>API_VER)
        stdhdrparams = Dict{Symbol,String}()
        for (n,v) in headers
            if (n in fieldnames(StandardHeaders))
                stdhdrparams[n] = v
                reqhdrs[to_http_header_name(n)] = v
            else
                reqhdrs[string(n)] = v
            end
        end
        (verb == "PUT") && (reqhdrs["Content-Length"] = "0")
        stdhdr = StandardHeaders(; stdhdrparams...)

        new(account, verb, resource, stdhdr, reqhdrs)
    end
end

to_http_header_name(n::Symbol) = join(map(ucfirst, split(string(n), '_')), '-')

sign_strip(str::String) = replace(strip(str), "\n"=>" ")
sign_hdr(nv::Pair{String,String}) = lowercase(sign_strip(nv[1])) * ":" * sign_strip(nv[2])
sign_hdr(nv::Pair{T1,T2}) where {T1<:AbstractString, T2<:AbstractString} = sign_hdr(String(nv[1])=>String(nv[2]))

function sign_sharedkey(req::ServiceRequest, key::String)
    iob = IOBuffer()

    println(iob, req.verb)

    for n in fieldnames(StandardHeaders)
        v = getfield(req.standard_headers, n)
        # special case rule for Content-Length header in api version 2015-02-21 and later
        (n == :content_length) && (v !== nothing) && !isempty(v) && (parse(Int, v) == 0) && (v = nothing)
        println(iob, v === nothing ? "" : v)
    end

    x_ms_headers = sort(map(sign_hdr, collect(filter((p)->startswith(lowercase(strip(p[1])), "x-ms-"), req.headers))))
    for v in x_ms_headers
        println(iob, v)
    end

    print(iob, canonicalize_resource(req.account, req.resource))

    signingstr = String(take!(iob))
    hmacsign = base64encode(digest(MD_SHA256, signingstr, base64decode(key)))
    req.headers["Authorization"] = "SharedKey $(req.account):$(hmacsign)"
    nothing
end

function canonicalize_resource(resource_account::String, uri::String)
    url = HTTP.URIs.URI(uri)
    container = url.path
    query = url.query

    iob = IOBuffer()
    print(iob, "/", resource_account, container)
    if !isempty(query)
        qparts = sort(map(sign_hdr, collect(HTTP.URIs.queryparams(url))))
        for q in qparts
            print(iob, "\n", q)
        end
    end
    String(take!(iob))
end

function execute(req::ServiceRequest, key::String; retry_count::Int=0, retry_interval::Int=60)
    success = false
    resp = HTTP.Response()
    count = 0
    sign_sharedkey(req, key)
    while !success && count <= retry_count
        count += 1
        (count == 1) || sleep(retry_interval)
        resp = HTTP.request(uppercase(req.verb), HTTP.URIs.URI(req.resource), req.headers; DEFAULT_KWARGS...)
        success = (200 <= resp.status <= 206)
        success && (return resp)
        @warn("Storage service request failed. ", resp)
    end
    resp
end

issuccess(resp::HTTP.Response) = (200 <= resp.status <= 206)

end # module REST


text:
@testset "Node pop report" begin

fName = "report/nodepop"

excelPopReport( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career", filename=fName )
excelPopEvolutionReport( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career", filename=fName,
    overwrite=false )

end  # @testset "Node pop report"


text:
using Compat; import Compat.String

_winston_config = Inifile()

begin
    local fn
    for dir in [".";Pkg.dir();LOAD_PATH]
        fn = joinpath(dir, "Winston.ini")
        if isfile(fn) break end
        fn = joinpath(dir, "Winston", "src", "Winston.ini")
        if isfile(fn) break end
    end
    read(_winston_config, fn)
end

if VERSION < v"0.4-"
    split_keep_false(a, b) = split(a, b, false)
else
    split_keep_false(a, b) = split(a, b, keep=false)
end

function _atox(s::AbstractString)
    x = strip(s)
    if x == "nothing"
        return nothing
    elseif x == "true"
        return true
    elseif x == "false"
        return false
    elseif length(x) > 2 && lowercase(x[1:2]) == "0x"
        try
            h = parse(Int, x[3:end], 16)
            return h
        end
    elseif x[1] == '{' && x[end] == '}'
        style = Dict{Symbol,Any}()
        pairs = map(strip, split_keep_false(x[2:end-1], ','))
        for pair in pairs
            kv = split_keep_false(pair, ':')
            style[ Symbol(strip(kv[1])) ] = _atox(strip(kv[2]))
        end
        return style
    elseif x[1] == '"' && x[end] == '"'
        return x[2:end-1]
    end
    if ismatch(r"^[+-]?\d+$",x)
        return parse(Int,x)
    end
    r  = tryparse(Float64, x)
    isnull(r) ? x : get(r)
end

function config_value(section, option)
    strval = get(_winston_config, section, option, "nothing")
    _atox(strval)
end

function config_options(sec::AbstractString)
    opts = Dict{Symbol,Any}()
    if sec == "defaults"
        for (k,v) in _winston_config.defaults
            opts[Symbol(k)] = _atox(v)
        end
    elseif has_section(_winston_config, sec)
        for (k,v) in section(_winston_config, sec)
            opts[Symbol(k)] = _atox(v)
        end
    end
    opts
end