text:
"""
    yyyymmdd2yyyy(d)

Takes an Int64-formatted date and returns the year in Int64 format.
```@example
dt = 20100811
@show yyyymmdd2yyyy(dt)
```
"""
yyyymmdd2yyyy(d) = div(d,10000)

"""
    yyyymmdd2yyyymm(d)

Takes an Int64-formatted date and returns the year and month in Int64 format.
```@example
dt = 20110909
@show yyyymmdd2yyyymm(dt)
```
"""
yyyymmdd2yyyymm(d) = div(d,100)

"""
    yyyymmdd2date(yyyymmdd::Number)

Takes in an Int-64 formatted date and returns a Date object in y-m-d format.
```@example
dt = 19980430
@show yyyymmdd2date(dt)
```
"""
function yyyymmdd2date(yyyymmdd::Number)
    yyyy=div(yyyymmdd,10000)
    mm=div(mod(yyyymmdd,10000),100)
    dd=mod(yyyymmdd,100)
    Date(yyyy,mm,dd)
end

"""
    yyyymm2date(yyyymm::Number)

Takes in a Int64-formatted year and date and returns a Date object,
with the last day of the month attached.

# Example
```@example
dt = 199804
@show yyyymm2date(dt)
```
"""
function yyyymm2date(yyyymm::Number)
    yyyy=div(yyyymm,100)
    mm=mod(yyyymm,100)
    Dates.lastdayofmonth(Date(yyyy,mm,1))
end

const dduuuyyyydateformat = Dates.DateFormat("dduuuyyyy")
"""
    dduuuyyyy2date(dduuuyyyy)

Takes in a string of the form "31Jan1997" and returns a Date
object in the format y-m-d.
```@example
d = "31Jan1997"
@show dduuuyyyy2date(d)
```
"""
function dduuuyyyy2date(dduuuyyyy)
    Date(dduuuyyyy,dduuuyyyydateformat)
end

const uuuddyyyydateformat = Dates.DateFormat("uuu dd, yyyy")
"""
    uuuddyyyy2date(uuuddyyyy)

Takes in a string of the form "Month Day, Year" and returns a Date in the
format y-m-d
```@example
u = "Jan 31, 1891"
@show uuuddyyyy2date(u)
```
"""
function uuuddyyyy2date(uuuddyyyy)
    Date(uuuddyyyy,uuuddyyyydateformat)
end

"""
    datetimestring2datetime(str,datetimeformat)

Takes in a string of the form "" and converts it to the user-specified
DateTime format.
```@example
using Dates #hide
str = "1989-07-07"
dt_fmt = DateFormat("y-m-d")
@show datetimestring2datetime(str, dt_fmt)
```
"""
function datetimestring2datetime(str,datetimeformat)
  try
    dt=DateTime(str,datetimeformat)
    y=Dates.year(dt)
    if 70<=y<=99
      y+=1900
    elseif 0<=y<70
      y+=2000
    end
    dt=Dates.DateTime(y,Dates.month(dt),Dates.day(dt),Dates.hour(dt),Dates.minute(dt),Dates.second(dt))
  catch
    missing
  end
end

"""
    date2yyyymmdd(date)

Takes in a Date object and converts it back to Int64.
```@example
using Dates #hide
dt = Date(2014,1,29)
@show date2yyyymmdd(dt)
```
"""
function date2yyyymmdd(date)
    try
        Dates.year(date)*10000+Dates.month(date)*100+Dates.day(date)
    catch
        missing
    end
end

"""
    date2yyyymm(date)

Takes in a Date object and transforms it back to Int64.
```@example
using Dates #hide
dt = Date(2010,8,11)
@show date2yyyymm(dt)
```
"""
function date2yyyymm(date)
    try
        Dates.year(date)*100+Dates.month(date)
    catch
        missing
    end
end

"""
    date2yyyy(date)

Takes in a Date object and returns the year as an Int64 type.
```@example
using Dates #hide
dt = Date(2011,7,24)
@show date2yyyy(dt)
```
"""
function date2yyyy(date)
    try
        Dates.year(date)
    catch
        missing
    end
end

"""
    datetimestring2yyyymmdd(str,datetimeformat)

Takes in a string and a DateTime format, converts the date back into Int64 type.
```@example
using Dates #hide
str = "1999-07-24"
dt_fmt = DateFormat("y-m-d")
@show datetimestring2yyyymmdd(str,dt_fmt)
```
"""
function datetimestring2yyyymmdd(str,datetimeformat)
    d = datetimestring2datetime(str,datetimeformat)
    if ismissing(d)
        missing
    else
        date2yyyymmdd(d)
    end
end

"""
    date2yyyyqq(d)

Takes in a Date object, formatted as y-m-d, and returns the year and quarter in Int64 type.
```@example
using Dates #hide
dt = Date(2009,6,22)
@show date2yyyyqq(dt)
```
"""
function date2yyyyqq(d)
  try
    date2yyyy(d)*100 + Dates.quarterofyear(d)
  catch
    missing
  end
end

"""
    date2yyyyww(d)

Takes in a Date object, formatted as y-m-d, and returns the year and
week in Int64 format.
```@example
using Dates #hide
dt = Date(2011,9,8)
@show dt
@show date2yyyyww(dt)
```
"""
function date2yyyyww(d)
  try
    date2yyyy(d)*100 + Dates.week(d)
  catch
    missing
  end
end

"""
    yearqtr2date(s, dlm="-Q")

convert something like 2018-Q3 into end of quarter date
`dlm` defines the string delimiter characters passed to split()
"""
function yearqtr2date(s, dlm="-Q")
  sp = split(s, dlm)
  y = parse(Int,sp[1])
  q = parse(Int,sp[2])
  Dates.lastdayofquarter(Date(y) + Dates.Month(3*(q-1)))
end

"""
    yearmonth2date(s, dlm="-")

convert something like 2018-11 into year month
`dlm` defines the string delimiter characters passed to split()
"""
function yearmonth2date(s, dlm="-")
  sp = split(s, dlm)
  y = parse(Int,sp[1])
  m = parse(Int,sp[2])
  Dates.lastdayofmonth(Date(y,m,1))
end


text:
# Script for processing BATSRUS explicit shocktube test weak scaling runs.
#
# The runlogs are stored in Google Drive, ScalingTests.
#
# Hongyang Zhou, hyzhou@umich.edu 05/09/2019, last modified 09/03/2020
include("Utility.jl")

using PyPlot
#using Plots
#pyplot()

Machine = "Frontera"

if Machine == "Bluewaters"

   # New results after modifying the message passing
   res = postprocess("results/gnu_BATL_update")
   #res = postprocess("results/gnu")

   result_0 = sort( res[res.Thread.==0,:],(:Proc))
   result_1 = sort( res[res.Thread.==1,:],(:Proc))
   result_2 = sort( res[res.Thread.==2,:],(:Proc))
   result_4 = sort( res[res.Thread.==4,:],(:Proc))
   result_8 = sort( res[res.Thread.==8,:],(:Proc))
   result_16 = sort(res[res.Thread.==16,:],(:Proc))
   result_32 = sort(res[res.Thread.==32,:],(:Proc))

   time_0 = copy(result_0.Time)
   time_1 = copy(result_1.Time)
   time_2 = copy(result_2.Time)
   time_4 = copy(result_4.Time)
   time_8 = copy(result_8.Time)
   time_16 = copy(result_16.Time)
   time_32 = copy(result_32.Time)

   append!(time_0,fill(NaN,size(time_32)[1]-size(time_0)[1]))
   append!(time_1,fill(NaN,size(time_32)[1]-size(time_1)[1]))
   append!(time_2,fill(NaN,size(time_32)[1]-size(time_2)[1]))
   append!(time_4,fill(NaN,size(time_32)[1]-size(time_4)[1]))
   append!(time_8,fill(NaN,size(time_32)[1]-size(time_8)[1]))
   append!(time_16,fill(NaN,size(time_32)[1]-size(time_16)[1]))

   # Calculate the number of cell updates per second for 2nd order scheme
   cell_size = sort(8^3 .* 256 .* result_32.Thread .* result_32.Proc .* 2)

   speed_0 = cell_size ./ time_0
   speed_1 = cell_size ./ time_1
   speed_2 = cell_size ./ time_2
   speed_4 = cell_size ./ time_4
   speed_8 = cell_size ./ time_8
   speed_16 = cell_size ./ time_16
   speed_32 = cell_size ./ time_32

   proc = result_32.Proc .* result_32.Thread

   speed = hcat(speed_0, speed_1, speed_2, speed_4, speed_8, speed_16, speed_32)
   speed_ref = speed_0[1] ./ proc[1] .* proc

   fig = plt.figure(figsize=(18,5))
   # For 1,2,4 threads
   ax1 = plt.subplot(121)

   plt.plot(proc[1:7], speed[1:7,:]/1e9, "-o", lw=1, markersize=3, alpha=0.8)
   plt.plot(proc[1:7], speed_ref[1:7]/1e9, ":")
   plt.legend(loc="upper left",
         ("pure MPI", "1 thread", "2 threads", "4 threads", "8 threads",
         "16 threads", "32 threads", "linear ref."))
   plt.grid(true)
   ax1.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Explicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)

   # For 1,2,4,8,16,32 threads
   ax2 = plt.subplot(122)
   plt.plot(proc, speed/1e9, "-o", markersize=3, alpha=0.8)
   plt.plot(proc, speed_ref/1e9, ":")
   plt.legend(loc="upper left",
         ("pure MPI", "1 thread", "2 threads", "4 threads","8 threads",
         "16 threads", "32 threads", "linear ref."))
   plt.grid(true)
   ax2.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Explicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)
   #plt.text(0.9,-0.08, "number of cores", transform=ax2.transAxes)

   #=
   mylabel = ["1 thread" "16 threads" "32 threads" "" "" ""]

   plot(proc[1:end],speed[1:end,:], seriestype=[:line, :scatter],
      markersize=4, lw=2, markercolor=[:red :green :purple],
      markerstrokewidth = 0.1,
      linecolor=[:red :green :purple],
      xscale=:identity, yscale=:identity, legend=:best,
      xlim=[10^2, 1.4e5], ylim=[10^6,1.8e9],
      label=mylabel,
      xlabel="number of cores",
      ylabel="cell updates per second",
      title="BATSRUS Hybrid Weak Scaling",
      xtickfontsize=16, ytickfontsize=16,
      legendfontsize=16, axis=:equal,
      size=(900,600))

   plot!(proc[1:end],speed_ref[1:end], linestyle=:dashdot, linecolor=:black,
         label="linear ref.")
   =#

   #=
   speed = hcat(speed_1,speed_2,speed_4,speed_8)
   mylabel = ["1 thread" "2 threads" "4 threads" "8 threads" "" "" "" ""]

   plot(proc[1:6],speed[1:6,:], seriestype=[:line, :scatter],
      markersize=4, lw=2, markercolor=[1 2 3 4],
      markerstrokewidth = 0.1,
      linecolor=[1 2 3 4],
      xscale=:identity, yscale=:identity, legend=:best,
      label=mylabel,
      xlabel="number of cores",
      ylabel="cell updates per second",
      xtickfontsize=16, ytickfontsize=16,
      legendfontsize=16, axis=:equal,
      reuse=false)
   =#

   #using StatsPlots

   #=
   @df results plot(:Thread, :Time, xticks=[1,2,4,8,16,32],
      seriestype=:scatter,
      title="Timing")

   @df results plot!(:Thread, :Time, xticks=[1,2,4,8,16,32],
      seriestype=:line,
      title="Timing")
   =#

   #plotly()

   #=
   @df results plot(:Thread, :Time, xticks=[1,2,4,8,16,32], group=(:Compiler),
      seriestype=[:line,:scatter],
      markersize=0.5 .* log2.(:Proc),
      xlabel="thread",
      ylabel="execution time [s]",
      title="Timing")
   =#

elseif Machine == "Stampede"

   # New results after modifying the message passing
   res = postprocess("results/intel_stampede")

   result_1 = sort( res[res[:Thread].==1,:],(:Proc))
   result_2 = sort( res[res[:Thread].==2,:],(:Proc))
   result_4 = sort( res[res[:Thread].==4,:],(:Proc))
   result_8 = sort( res[res[:Thread].==8,:],(:Proc))
   result_12 = sort( res[res[:Thread].==12,:],(:Proc))
   result_16 = sort(res[res[:Thread].==16,:],(:Proc))

   time_1 = result_1[:Time]
   time_2 = result_2[:Time]
   time_4 = result_4[:Time]
   time_8 = result_8[:Time]
   time_12 = result_12[:Time]
   time_16 = result_16[:Time]

   append!(time_1,fill(NaN,size(time_16)[1]-size(time_1)[1]))
   append!(time_2,fill(NaN,size(time_16)[1]-size(time_2)[1]))
   append!(time_4,fill(NaN,size(time_16)[1]-size(time_4)[1]))
   append!(time_8,fill(NaN,size(time_16)[1]-size(time_8)[1]))
   append!(time_12,fill(NaN,size(time_16)[1]-size(time_12)[1]))
   append!(time_16,fill(NaN,size(time_16)[1]-size(time_16)[1]))

   # Calculate the number of cell updates per second for 2nd order scheme
   cell_size = sort(8^3 .* 256 .* result_16[:Thread] .* result_16[:Proc] .* 2)

   speed_1 = cell_size ./ time_1
   speed_2 = cell_size ./ time_2
   speed_4 = cell_size ./ time_4
   speed_8 = cell_size ./ time_8
   speed_12 = cell_size ./ time_12
   speed_16 = cell_size ./ time_16

   proc = result_16[:Proc] .* result_16[:Thread]
   #=
   speed = hcat(speed_1,speed_16,speed_32)

   mylabel = ["1 thread" "16 threads" "32 threads" "" "" ""]

   plot(proc,speed, seriestype=[:line, :scatter],
      markersize=4, lw=2, markercolor=[:red :green :purple],
      markerstrokewidth = 0.1,
      linecolor=[:red :green :purple],
      xscale=:identity, yscale=:identity, legend=:best,
      xlim=[10^2, 5.5e5], ylim=[10^6,4.5e9],
      label=mylabel,
      xlabel="number of cores",
      ylabel="cell updates per second",
      title="BATSRUS Hybrid Weak Scaling",
      xtickfontsize=16, ytickfontsize=16,
      legendfontsize=16, axis=:equal)

   speed_ref = speed_1[1] ./ proc[1] .* proc
   plot!(proc,speed_ref, linestyle=:dashdot, linecolor=:black,
         label="linear ref.")
   =#
   speed = hcat(speed_1,speed_2,speed_4,speed_8,speed_12,speed_16)
   speed_ref = speed_1[1] ./ proc[1] .* proc

   plt.figure(figsize=(5,5))
   plt.plot(proc,speed/1e9,"-o",markersize=3)
   plt.plot(proc,speed_ref/1e9,"--")
   plt.legend(loc="upper left",
         ("1 thread", "2 threads", "4 threads", "8 threads", "12 threads",
         "16 threads","linear ref."))
   plt.grid(true)
   plt.title("BATSRUS Weak Scaling, Explicit Scheme",y=1.03)
   plt.xlabel("number of cores")
   plt.ylabel(L"giga cell ($10^9$) updates per second")

   #=
   mylabel = ["1 thread", "2 threads", "4 threads", "8 threads", "16 threads",
         "", "", "", "", ""]
   mylabel = reshape(mylabel,1,:)

   plot(proc,speed[:,:], seriestype=[:line, :scatter],
      markersize=4, lw=2, markercolor=[1 2 3 4 5],
      markerstrokewidth = 0.1,
      linecolor=[1 2 3 4 5],
      xscale=:identity, yscale=:identity, legend=:best,
      label=mylabel,
      xlabel="number of cores",
      ylabel="cell updates per second",
      xtickfontsize=16, ytickfontsize=16,
      legendfontsize=16, axis=:equal,
      reuse=false)

   plot!(proc,speed_ref, linestyle=:dashdot, linecolor=:black,
         label="linear ref.")
   =#
elseif Machine == "Frontera"

   res = postprocess("RESULT_gnu_Expl")

   result_1 = sort( res[res.Thread.==1,:],(:Proc))
   result_2 = sort( res[res.Thread.==2,:],(:Proc))
   result_4 = sort( res[res.Thread.==4,:],(:Proc))
   result_7 = sort( res[res.Thread.==7,:],(:Proc))
   result_14 = sort( res[res.Thread.==14,:],(:Proc))
   result_28 = sort( res[res.Thread.==28,:],(:Proc))
   result_56 = sort(res[res.Thread.==56,:],(:Proc))

   time_1 = copy(result_1.Time)
   time_2 = copy(result_2.Time)
   time_4 = copy(result_4.Time)
   time_7 = copy(result_7.Time)
   time_14 = copy(result_14.Time)
   time_28 = copy(result_28.Time)
   time_56 = copy(result_56.Time)

   append!(time_1,fill(NaN,size(time_56)[1]-size(time_1)[1]))
   append!(time_2,fill(NaN,size(time_56)[1]-size(time_2)[1]))
   append!(time_4,fill(NaN,size(time_56)[1]-size(time_4)[1]))
   append!(time_7,fill(NaN,size(time_56)[1]-size(time_7)[1]))
   append!(time_14,fill(NaN,size(time_56)[1]-size(time_14)[1]))
   append!(time_28,fill(NaN,size(time_56)[1]-size(time_28)[1]))
   append!(time_56,fill(NaN,size(time_56)[1]-size(time_56)[1]))

   # Calculate the number of cell updates per second for 2nd order scheme
   cell_size = sort(8^3 .* 224 .* result_56.Thread .* result_56.Proc .* 2)

   speed_1 = cell_size ./ time_1
   speed_2 = cell_size ./ time_2
   speed_4 = cell_size ./ time_4
   speed_7 = cell_size ./ time_7
   speed_14 = cell_size ./ time_14
   speed_28 = cell_size ./ time_28
   speed_56 = cell_size ./ time_56

   proc = result_56.Proc .* result_56.Thread

   speed = hcat(speed_1, speed_2, speed_4, speed_7, speed_14, speed_28, speed_56)
   speed_ref = speed_1[1] ./ proc[1] .* proc

   fig = plt.figure(figsize=(6,5))

   ax1 = plt.subplot(111)

   plt.plot(proc, speed/1e9, "-o", lw=1, markersize=3, alpha=0.8)
   plt.plot(proc, speed_ref/1e9, ":")
   plt.legend(loc="upper left",
         ("1 thread", "2 thread", "4 thread", "7 threads", "14 threads",
         "28 threads", "56 threads", "linear ref."))
   plt.grid(true)
   ax1.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Explicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)

end

ext: jl


text:
using Test, Dates, CSV, MonthlyDates

replstr(x, kv::Pair...) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80), kv...), MIME("text/plain"), x), x)

##############################################################################
##
## Quarter
##
##############################################################################

@test Quarter(Date(1990, 1)) == Quarter(1)
@test quarterofyear(Date(1990, 1)) == 1
@test Month(Quarter(3)) == Month(9)
@test_throws InexactError Quarter(Month(4))

##############################################################################
##
## MonthlyDate
##
##############################################################################
# constructur
@test_throws ArgumentError MonthlyDate(1, 0)
@test_throws ArgumentError MonthlyDate(1, 13)

# since Dates.value(Date(1, 1, 1)) == 1
@test Dates.value(MonthlyDate(1, 1)) == 1
@test Dates.value(MonthlyDate(1.0, 1.0)) == 1

@test MonthlyDate(1990, 1) >= MonthlyDate(1989, 1)
@test Date(MonthlyDate(Date(1990, 1))) == Date(1990, 1)
@test Date(MonthlyDate(1990, 1)) == Date(1990, 1)
# similar definition for convert(Day, Date(1990, 1, 1))
@test convert(Month, MonthlyDate(1990, 1)) == Month(23869)
@test convert(MonthlyDate, Month(23869)) == MonthlyDate(1990, 1)

# promotion rules
@test MonthlyDate(1990, 1) == Date(1990, 1, 1)
@test MonthlyDate(1990, 1) != Date(1990, 1, 2)

# accessor
@test year(MonthlyDate(1990, 12)) == 1990
@test quarterofyear(MonthlyDate(1990, 1)) == 1
@test month(MonthlyDate(1990, 12)) == 12

@test Year(MonthlyDate(1990, 1)) == Year(1990)
@test Quarter(MonthlyDate(1990, 1)) == Quarter(1)
@test Month(MonthlyDate(1990, 1)) == Month(1)

# adjusters
@test firstdayofmonth(MonthlyDate(1990, 3)) == Date(1990, 3, 1)
@test lastdayofmonth(MonthlyDate(2000, 2)) == Date(2000, 2, 29)
@test firstdayofquarter(MonthlyDate(2000, 2)) == Date(2000, 1, 1)
@test lastdayofquarter(MonthlyDate(2000, 2)) == Date(2000, 3, 31)

# arithmetic
@test MonthlyDate(1990, 1) + Year(3) == MonthlyDate(1993, 1)
@test MonthlyDate(1990, 1) - Year(3) == MonthlyDate(1987, 1)
@test MonthlyDate(1990, 1) - Month(3) == MonthlyDate(1989, 10)
@test MonthlyDate(1990, 1) + Month(12) == MonthlyDate(1991, 1)
@test MonthlyDate(1990, 1) - Month(12) == MonthlyDate(1989, 1)
@test MonthlyDate(1990, 1) + Month(13) == MonthlyDate(1991, 2)
@test MonthlyDate(1990, 1) - Month(13) == MonthlyDate(1988, 12)
@test MonthlyDate(1990, 1) + Quarter(1) == MonthlyDate(1990, 4)
@test MonthlyDate(1990, 1) - Quarter(1) == MonthlyDate(1989, 10)
@test MonthlyDate(1990, 1) - MonthlyDate(1989, 1) == Month(12)
@test MonthlyDate(1990, 1) + Quarter(2) + Month(2) == MonthlyDate(1990, 9)

# ranges
@test length(range(MonthlyDate(1990, 1), MonthlyDate(1991, 1), step = Month(1))) == 13
@test length(range(MonthlyDate(1990, 1), MonthlyDate(1991, 1), step = Month(3))) == 5
@test length(range(MonthlyDate(1990, 1), MonthlyDate(1991, 1), step = Quarter(1))) == 5
@test length(range(MonthlyDate(1990, 1), MonthlyDate(1992, 1), step = Year(1))) == 3

# print
@test replstr([MonthlyDate(1990, 1)]) == "1-element $(string(Array{MonthlyDate,1})):\n 1990-01"
@test Dates.format(MonthlyDate(1990, 3), dateformat"Y-mm") == "1990-03"

# parse
@test MonthlyDate("1990-01") == MonthlyDate(1990, 1)
@test MonthlyDate("1990/01", "y/m") == MonthlyDate(1990, 1)
@test MonthlyDate("1990m01", dateformat"y\mm") == MonthlyDate(1990, 1)

# csv
io = IOBuffer()
df = (x = [MonthlyDate(1990, 1), MonthlyDate(1990, 2)], )
CSV.write(io, df)
@test String(take!(io)) == "x\n1990-01\n1990-02\n"
##############################################################################
##
## QuarterlyDate
##
##############################################################################
# constructur
@test_throws ArgumentError QuarterlyDate(1, 0)
@test_throws ArgumentError QuarterlyDate(1, 5)

@test Dates.value(QuarterlyDate(1, 1)) == 1
@test QuarterlyDate(1990, 1) - QuarterlyDate(1989, 1) == Quarter(4)
@test QuarterlyDate(1990, 1) >= QuarterlyDate(1989, 1)
@test Date(QuarterlyDate(1990, 1)) == Date(1990, 1)
@test QuarterlyDate(Date(1990, 1)) == QuarterlyDate(1990, 1)
@test QuarterlyDate(Date(1990, 2)) == QuarterlyDate(1990, 1)
@test DateTime(QuarterlyDate(1990, 1)) == DateTime(1990, 1)
@test QuarterlyDate(DateTime(1990, 1)) == QuarterlyDate(1990, 1)
@test QuarterlyDate(DateTime(1990, 2)) == QuarterlyDate(1990, 1)
@test MonthlyDate(QuarterlyDate(1990, 1)) == MonthlyDate(1990, 1)
@test QuarterlyDate(MonthlyDate(1990, 1)) == QuarterlyDate(1990, 1)
# similar definition for convert(Day, Date(1990, 1, 1))
@test convert(Quarter, QuarterlyDate(1990, 1)) == Quarter(7957)
@test convert(QuarterlyDate, Quarter(7957)) == QuarterlyDate(1990, 1)

# promotion rules
@test QuarterlyDate(1990, 1) != Date(1990, 1, 2)
@test QuarterlyDate(1990, 1) == MonthlyDate(1990, 1)
@test QuarterlyDate(1990, 1) != MonthlyDate(1990, 2)

# accessor
@test year(QuarterlyDate(1990, 4)) == 1990
@test quarterofyear(QuarterlyDate(1990, 4)) == 4

@test Year(QuarterlyDate(1990, 1)) == Year(1990)
@test Quarter(QuarterlyDate(1990, 1)) == Quarter(1)

# adjusters
@test firstdayofquarter(QuarterlyDate(1990, 3)) == Date(1990, 7, 1)
@test lastdayofquarter(QuarterlyDate(1990, 3)) == Date(1990, 9, 30)

# arithmetic
@test QuarterlyDate(1990, 1) + Year(1) == QuarterlyDate(1991, 1)
@test QuarterlyDate(1990, 1) - Year(1) == QuarterlyDate(1989, 1)
@test QuarterlyDate(1990, 1) + Quarter(1) == QuarterlyDate(1990, 2)
@test QuarterlyDate(1990, 1) - Quarter(1) == QuarterlyDate(1989, 4)
@test QuarterlyDate(1990, 1) + Year(1) + Quarter(3) == QuarterlyDate(1991, 4)

# ranges
@test length(range(QuarterlyDate(1990, 1), QuarterlyDate(1991, 1), step = Quarter(1))) == 5
@test length(range(QuarterlyDate(1990, 1), QuarterlyDate(1991, 1), step = Quarter(2))) == 3
@test length(range(QuarterlyDate(1990, 1), QuarterlyDate(1992, 1), step = Year(1))) == 3

# io
@test replstr([QuarterlyDate(1990, 1)]) == "1-element $(string(Array{QuarterlyDate,1})):\n 1990-Q1"

# parse
@test QuarterlyDate("1990-07") == QuarterlyDate(1990, 3)
@test QuarterlyDate("1990-Q2") == QuarterlyDate(1990, 2)
@test QuarterlyDate("1990/01", "y/m") == QuarterlyDate(1990, 1)
@test QuarterlyDate("1990m01", dateformat"y\mm") == QuarterlyDate(1990, 1)
@test QuarterlyDate("1990m07", dateformat"y\mm") == QuarterlyDate(1990, 3)
@test QuarterlyDate("1990-03", dateformat"y-q") == QuarterlyDate(1990, 3)

# formater
@test Dates.format(QuarterlyDate(1990, 3), dateformat"Y-\Qq") == "1990-Q3"

# csv
io = IOBuffer()
df = (x = [QuarterlyDate(1990, 1), QuarterlyDate(1990, 2)], )
CSV.write(io, df)
@test String(take!(io)) == "x\n1990-Q1\n1990-Q2\n"


text:
using SparseArrays
using Mitosis
using LinearAlgebra, Statistics, Random

"""
Graph Laplacian of a `m×n` lattice.
"""
function laplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1
                    S[linear[i2, j2], linear[i, j]] -= 1

                    S[linear[i, j], linear[i, j]] += 1
                    S[linear[i2, j2], linear[i2, j2]] += 1
                end
            end
        end
    end
    S
end
# indices
m, n = 20, 20
C = CartesianIndices((m, n))
L = LinearIndices((m, n))

# tuning parameters
κ = 4 # κ-times integrated white noise, so κ is a smoothing parameter
ϵ = 0.1 # ϵ regularises
σ = 0.02 # σ the amount of white noise in the image, in absence of pixel noise regularises
λ = .2 # λ is tradeoff between exploration and exploitation in acquisition (this is just the upper confidence level acquisition function)

# data
W = [exp(-norm(([i,j]-[m÷3,n÷3]))^2/25.0) + σ*randn() for i in 1:m, j in 1:n]

# prior (power of the graph laplacian)
Λ = laplacian(Float32, m, n)
u = Gaussian{(:F,:Γ)}(zeros(m*n), (Λ + ϵ*I)^κ)

# Bayesian optimization
i = CartesianIndex(m÷2, n÷2)
is = [Tuple(i)]
for k in 1:11
    global i, u
    # observation scheme
    H = sparse(zeros(1, m*n))
    H[L[i]] = 1
    k = kernel(Gaussian; μ=AffineMap(H, [0.]), Σ=ConstantMap(Matrix(σ^2*I(1))))
    # update
    _, p = Mitosis.backward(BF(), k, [W[i]])
    # fusion
    _, u = fuse(u, p)
    # aquisition
    if m + n > 50
        _, i_ = findmax(mean(u) + 2λ./sqrt.(diag(u.Γ))) # proxy for inverse
    else
        _, i_ = findmax(mean(u) + λ*sqrt.(diag(inv(Matrix(u.Γ)))))
    end
    i = C[i_]
    push!(is, Tuple(i))
    println(i)
end

imax = C[findmax(mean(u))[2]] # estimate of the max
i = findmax(W)[2] # truth

# plot
using Makie
img(x) = heatmap(reshape(x, m, n))
p1 = img(mean(u))
scatter!(p1, first.(is), last.(is))
scatter!(p1, [i[1]], [i[2]], color=:green)
scatter!(p1, [imax[1]], [imax[2]], color=:red) # red estimate

p2 = img(W)
scatter!(p2, [i[1]], [i[2]], color=:green)
scatter!(p2, [imax[1]], [imax[2]], color=:red) # red estimate
hbox(p1, p2)


text:
using DataFrames, CSV
# using Pkg
# Pkg.activate(ENV["JULIA_DQMC"])
push!(LOAD_PATH, ENV["JULIA_DQMC"])

include(joinpath(ENV["JULIA_DQMC"], "src/dqmc_framework.jl"))

SAFE_MULT = 10

function estimate_memory_usage_tdgfs()
    df = DataFrame(L=Int[], B=Int[], mem=Float64[])

    for L in [4,8,10,12,14]
        for beta in [5,10,20,40]
            mem = estimate_memory_usage_tdgfs(L, beta; safe_mult=SAFE_MULT)
            push!(df,[L, beta, mem])
        end
    end

    return df
end

df = estimate_memory_usage_tdgfs()

CSV.write("tdgfs_memory_safe_mult_$(SAFE_MULT).csv", df)

using PyPlot, Polynomials
plt.style.use("publication_tex")
# plt.style.use("default")

function plot_memusage(df)
    g = groupby(df, :L);

    for (i, gl) in enumerate(g)
        x = gl[:B]
        y = gl[:mem] ./ 1024 # MB -> GB
        p = polyfit(x, y, 1)
        m = round(p.a[2], digits=3)
        b = round(p.a[1], digits=3)
        plot(x, y, "o-", color="C$i", label="L = $(gl[:L][1]) ($(m)x + $b)")
    end
    handles, labels = gca().get_legend_handles_labels()
    legend(reverse(handles), reverse(labels), frameon=false, prop=Dict("size" => 15))
    xlim([3,42])
    # ylim([-0.7, 11])
    xlabel(L"inverse temperature $\beta$")
    ylabel("tdgfs memory usage in GB")
    tight_layout()
    savefig("tdgfs_memory_safe_mult_$(SAFE_MULT).pdf")
    nothing
end

plot_memusage(df)

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using FProfile
using Base.Test
using NBInclude

nbinclude("../Manual.ipynb")


text:
export DecayCutPruningAlgo, DecayCutPruner

"""
    DecayCutPruningAlgo <: AbstractCutPruningAlgo

Removes the cuts with lower trust where the trust is initially
`newcuttrust + bonus` and is updated using `trust -> λ * trust + used`
after each optimization done with it.
The value `used` is 1 if the cut was used and 0 otherwise.
It has a bonus equal to `mycutbonus` if the cut was generated using a trial
given by the problem using this cut.
We say that the cut was used if its dual value is nonzero.
"""
struct DecayCutPruningAlgo <: AbstractCutPruningAlgo
    maxncuts::Int
    λ::Float64
    newcuttrust::Float64
    mycutbonus::Float64
    function DecayCutPruningAlgo(maxncuts::Int, λ=0.9, newcuttrust=0.8, mycutbonus=1)#newcuttrust=(1/(1/0.9-1))/2, mycutbonus=(1/(1/0.9-1))/2)
        new(maxncuts, λ, newcuttrust, mycutbonus)
    end
end

mutable struct DecayCutPruner{N, T} <: AbstractCutPruner{N, T}
    # used to generate cuts
    isfun::Bool
    islb::Bool
    lazy_minus::Bool
    A::AbstractMatrix{T}
    b::AbstractVector{T}

    maxncuts::Int

    trust::Vector{Float64}
    ids::Vector{Int}
    id::Int

    λ::Float64
    newcuttrust::Float64
    mycutbonus::Float64

    # check redundancy with old cuts when adding new cuts
    excheck::Bool

    # tolerance to check redundancy between two cuts
    TOL_EPS::Float64

    function (::Type{DecayCutPruner{N, T}})(sense::Symbol, maxncuts::Int, λ=0.9,
            newcuttrust=0.8, mycutbonus=1, lazy_minus::Bool=false, tol=1e-6, excheck::Bool=false) where {N, T} #newcuttrust=(1/(1/0.9-1))/2, mycutbonus=(1/(1/0.9-1))/2)
        isfun, islb = gettype(sense)
        new{N, T}(isfun, islb, lazy_minus, spzeros(T, 0, N), T[], maxncuts, Float64[], Int[], 0, λ, newcuttrust, mycutbonus, excheck, tol)
    end
end

(::Type{CutPruner{N, T}})(algo::DecayCutPruningAlgo, sense::Symbol, lazy_minus::Bool=false) where {N, T} = DecayCutPruner{N, T}(sense, algo.maxncuts, algo.λ, algo.newcuttrust, algo.mycutbonus, lazy_minus)

# COMPARISON

function addusage!(man::DecayCutPruner, σρ)
    if ncuts(man) > 0
        man.trust .*= man.λ
        man.trust[σρ .> 1e-6] .+= 1
    end
end

function initialtrust(man::DecayCutPruner, mycut)
    if mycut
        man.newcuttrust + man.mycutbonus
    else
        man.newcuttrust
    end
end

function isbetter(man::DecayCutPruner, i::Int, mycut::Bool)
    if mycut
        # If the cut has been generated, that means it is useful
        false
    else
        # The new cut has initial trust initialtrust(man, false)
        # but it is a bit disadvantaged since it is new so
        # as we advantage the new cut if mycut == true,
        # we advantage this cut by taking initialtrust(man, true)
        # with true instead of false
        man.trust[i] > initialtrust(man, mycut)
    end
end


text:
using DocStringExtensions
"""
Select the largest peak among a set of peak candidates
(which must be passed as an argument to the function along with their time of occurrence),
that are within a minimum and maximum latency window. If no suitable
peak is found the function returns `missing`.

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the waveform for which the peak is sought.
* `peakPnts::AbstractVector{Real}`: vector with candidate peak points. Peak candidates can be found with the `findPeaks` function.
* `peakTimes::AbstractVector{Real}`: vector with candidate peak times. Peak candidates can be found with the `findPeaks` function.
* `minLat::Real`: minimum latency of the peak.
* `maxLat::Real`: maximum latency of the  peak.

##### Returns

* `peakPoint::Real`: index of the point at which the highest peak is found. If no peak is found `missing` is returned.
* `peakTime`: time at which the highest peak is found, in seconds. If no peak is found `missing` is returned.

##### Examples

```julia
    using BrainWave
    sampRate = 256; nTaps=64; minLat=0.5; maxLat=2.5
    rec, evtTab = simulateRecording(nChans=1, dur=4, sampRate=sampRate)
    filterContinuous!(rec, sampRate, "lowpass", nTaps, [2], channels=[1], transitionWidth=0.2)
    pkPnts, pkTimes = findPeaks(rec[1,:], sampRate)
    peakPnt, peakTime = selectLargestPeakInWindow(rec[1,:], pkPnts, pkTimes, minLat, maxLat)
    ## not run
    ## using PlotlyJS
    ## tArr = collect(1:length(rec[1,:]))/sampRate
    ## s1 = scatter(;x=tArr, y=rec[1,:], name="Waveform")
    ## s2 = scatter(;x=tArr[pkPnts], y=rec[1, pkPnts], mode="markers", name="Peak Candidates")
    ## s3 = scatter(;x=[tArr[peakPnt]], y=[rec[1, peakPnt]], mode="markers", marker_size=10, name="Largest Peak")
    ## shapes = rect([minLat], [maxLat], [0], [1], fillcolor="gray", opacity=0.2, line_width=0, xref="x", yref="paper")
    ## plot([s1, s2, s3], Layout(shapes=shapes, xaxis_title="Time (s)", yaxis_title="Amplitude (a.u.)"))

```

"""
function selectLargestPeakInWindow(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, peakPnts::AbstractVector{P}, peakTimes::AbstractVector{S}, minLat::Real, maxLat::Real) where {T<:Real, P<:Real, S<:Real}

    peakCandidates = findall((peakTimes .<= maxLat ) .& (peakTimes .>= minLat))
    nCandidates = length(peakCandidates)
    if nCandidates == 0
        peakPnt = missing
        peakTime = missing
    elseif nCandidates == 1
        peakPnt = peakPnts[peakCandidates[1]]
        peakTime = peakTimes[peakCandidates[1]]
    else
        idx = findall(sig[peakPnts[peakCandidates]] .== maximum(sig[peakPnts[peakCandidates]]))[1]
        peakPnt = peakPnts[peakCandidates[idx]]
        peakTime = peakTimes[peakCandidates[idx]]
    end

    return peakPnt, peakTime

end

"""
Select the largest trough among a set of trough candidates
(which must be passed as an argument to the function along with their time of occurrence),
that are within a minimum and maximum latency window. If no suitable
trough is found the function returns `missing`.

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the waveform for which the trough is sought.
* `troughPnts::AbstractVector{Real}`: vector with candidate trough points. Trough candidates can be found with the `findTroughs` function.
* `troughTimes::AbstractVector{Real}`: vector with candidate trough times. Trough candidates can be found with the `findTroughs` function.
* `minLat::Real`: minimum latency of the trough.
* `maxLat::Real`: maximum latency of the  trough.

##### Returns

* `troughPoint::Real`: index of the point at which the highest trough is found. If no trough is found `missing` is returned.
* `troughTime`: time at which the highest trough is found, in seconds. If no trough is found `missing` is returned.

##### Examples

```julia
    sampRate = 256; nTaps=64; minLat=0.5; maxLat=2.5
    rec, evtTab = simulateRecording(nChans=1, dur=4, sampRate=sampRate)
    filterContinuous!(rec, sampRate, "lowpass", nTaps, [2], channels=[1], transitionWidth=0.2)
    trPnts, trTimes = findTroughs(rec[1,:], sampRate)
    troughPnt, troughTime = selectLargestTroughInWindow(rec[1,:], trPnts, trTimes, minLat, maxLat)
    ## not run
    ## using PlotlyJS
    ## tArr = collect(1:length(rec[1,:]))/sampRate
    ## s1 = scatter(;x=tArr, y=rec[1,:], name="Waveform")
    ## s2 = scatter(;x=tArr[pkPnts], y=rec[1, pkPnts], mode="markers", name="Trough Candidates")
    ## s3 = scatter(;x=[tArr[troughPnt]], y=[rec[1, troughPnt]], mode="markers", marker_size=10, name="Largest Trough")
    ## shapes = rect([minLat], [maxLat], [0], [1], fillcolor="gray", opacity=0.2, line_width=0, xref="x", yref="paper")
    ## plot([s1, s2, s3], Layout(shapes=shapes, xaxis_title="Time (s)", yaxis_title="Amplitude (a.u.)"))

```

"""
function selectLargestTroughInWindow(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, troughPnts::AbstractVector{P}, troughTimes::AbstractVector{S}, minLat::Real, maxLat::Real, maxAmp::Real=Inf) where {T<:Real, P<:Real, S<:Real}

    troughCandidates = findall((troughTimes .<= maxLat ) .& (troughTimes .>= minLat))
    #make sure trough amplitude is smaller than peak amplitude
    troughCandidates = troughCandidates[findall(sig[troughPnts[troughCandidates]] .< maxAmp)]
    nCandidates = length(troughCandidates)
    if nCandidates == 0
        troughPnt = missing
        troughTime = missing
    elseif nCandidates == 1
        troughPnt = troughPnts[troughCandidates[1]]
        troughTime = troughTimes[troughCandidates[1]]
    else
        idx = findall(sig[troughPnts[troughCandidates]] .== minimum(sig[troughPnts[troughCandidates]]))[1]
        troughPnt = troughPnts[troughCandidates[idx]]
        troughTime = troughTimes[troughCandidates[idx]]
    end

    return troughPnt, troughTime

end

"""
Select the strongest inflection point among a set of inflection point candidates
(which must be passed as an argument to the function along with their time of occurrence),
that are within a minimum and maximum latency window. If no suitable
inflection point is found the function returns `missing`.

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the waveform for which the inflection point is sought.
* `dy::AbstractVector{Real}`: first derivative of `sig`.
* `inflectionPnts::AbstractVector{Real}`: vector with candidate inflection points. Inflection point candidates can be found with the `findInflections` function.
* `inflectionTimes::AbstractVector{Real}`: vector with candidate inflection times. Inflection point candidates can be found with the `findInflections` function.
* `minLat::Real`: minimum latency of the inflection.
* `maxLat::Real`: maximum latency of the  inflection.

##### Returns

* `inflectionPoint::Real`: index of the point at which the highest inflection point is found. If no inflection point is found `missing` is returned.
* `inflectionTime`: time at which the highest inflection point is found, in seconds. If no inflection point is found `missing` is returned.

##### Examples

```julia
    sampRate = 256; nTaps=64; minLat=1; maxLat=2
    rec, evtTab = simulateRecording(nChans=1, dur=4, sampRate=sampRate)
    filterContinuous!(rec, sampRate, "lowpass", nTaps, [2], channels=[1], transitionWidth=0.2)
    inflPnts, inflTimes = findInflections(rec[1,:], sampRate)
    dy = zeros(size(rec[1,:])); dy[2:end] = diff(rec[1,:])
    inflectionPnt, inflectionTime = selectStrongestInflectionInWindow(rec[1,:], dy, inflPnts, inflTimes, minLat, maxLat)
    ## not run
    ## using PlotlyJS
    ## tArr = collect(1:length(rec[1,:]))/sampRate
    ## s1 = scatter(;x=tArr, y=rec[1,:], name="Waveform")
    ## s2 = scatter(;x=tArr[pkPnts], y=rec[1, pkPnts], mode="markers", name="Inflection Candidates")
    ## s3 = scatter(;x=[tArr[inflectionPnt]], y=[rec[1, inflectionPnt]], mode="markers", marker_size=10, name="Strongest Inflection")
    ## shapes = rect([minLat], [maxLat], [0], [1], fillcolor="gray", opacity=0.2, line_width=0, xref="x", yref="paper")
    ## plot([s1, s2, s3], Layout(shapes=shapes, xaxis_title="Time (s)", yaxis_title="Amplitude (a.u.)"))

```

"""
function selectStrongestInflectionInWindow(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, dy::Union{AbstractMatrix{Q}, AbstractVector{Q}}, inflPnts::AbstractVector{P}, inflTimes::AbstractVector{S}, minLat::Real, maxLat::Real; maxAmp::Real=Inf) where {T<:Real, Q<:Real, P<:Real, S<:Real}

    inflectionCandidates = findall((inflTimes .<= maxLat ) .& (inflTimes .>= minLat))
    #make sure inflection amplitude is smaller than peak amplitude
    inflectionCandidates = inflectionCandidates[findall(sig[inflPnts[inflectionCandidates]] .< maxAmp)]
    nCandidates = length(inflectionCandidates)
    if nCandidates == 0
        inflectionPnt = missing
        inflectionTime = missing
    elseif nCandidates == 1
        inflectionPnt = inflPnts[inflectionCandidates[1]]
        inflectionTime = inflTimes[inflectionCandidates[1]]
    else
        #strategy 1: find the minimum of signal (useful as surrogate of finding troughs)
        #idx = find(sig[inflPnts[inflectionCandidates]] .== minimum(sig[inflPnts[inflectionCandidates]]))[1]
        # strategy 2: find point at which first derivative is closer to zero (closer to a trough)
        idx = findall(abs.(dy[inflPnts[inflectionCandidates]]) .== minimum(abs.(dy[inflPnts[inflectionCandidates]])))[1]
        inflectionPnt = inflPnts[inflectionCandidates[idx]]
        inflectionTime = inflTimes[inflectionCandidates[idx]]
    end

    return inflectionPnt, inflectionTime

end

"""
Attempt to find peaks and troughs of the ABR response for a click of a given level. The algorithm is
largely based on Bradley and Wilson (2004).

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the ABR waveform for which the peaks and troughs are sought.
* `stimLevel::Real`: the level of the click used to evoke the ABR response.
* `sampRate::Real`: the sampling rate of the ABR recording.
* `epochStart::Real`: the time, in seconds, at which the epoch starts.
* `waveVLatencyOffset::Real`: additional wave V latency delay (e.g. if stimulus is different than a standard click you may want to specify an additional delay on top of the delay computed by the formula)
* `dBRef::String`: whether the stimulus level is specified in dB SPL `SPL`, or in dB normal hearing level `nHL`.

##### Returns

A dataframe with the following columns:

* `wave::String`: wave label.
* `peakPoint::Real`: index of the point at which the peak is detected in the waveform.
* `troughPoint::Real`: index of the point at which the peak is detected in the waveform.
* `peakLatency::Real`: latency of the peak, in seconds.
* `troughLatency::Real`: latency of the trough, in seconds.
* `peakAmp::Real`: amplitude of the peak, in microvolts.
* `troughAmp::Real`: amplitude of the trough, in microvolts.
* `peakTroughAmp::Real`: peak-to-trough amplitude, in microvolts.
* `minPeakLat::Real`: minimum peak latency, in seconds, used by the algorithm to find the peak.
* `maxPeakLat::Real`: maximum peak latency, in seconds, used by the algorithm to find the peak.
* `minTroughLat::Real`: minimum trough latency, in seconds, used by the algorithm to find the trough.
* `maxTroughLat::Real`: maximum trough latency, in seconds, used by the algorithm to find the trough.

##### References

* Bradley, A. P., & Wilson, W. J. (2004). Automated Analysis of the Auditory Brainstem Response. Proceedings of the 2004 Intelligent Sensors, Sensor Networks and Information Processing Conference, 2004., 541–546. http://doi.org/10.1109/ISSNIP.2004.1417519

* Prosser, S., & Arslan, E. (1987). Prediction of auditory brainstem wave V latency as a diagnostic tool of sensorineural hearing loss. Audiology : Official Organ of the International Society of Audiology, 26(3), 179–87.

##### Examples

```julia

```

"""
function findABRWaves(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, stimLevel::Real, sampRate::Real; epochStart::Real=0, waveVLatencyOffset::Real=0, dBRef::String="SPL") where {T<:Real}
    sig = vec(sig)
    if dBRef == "SPL"
        stimLevel = stimLevel - 31
    elseif dBRef == "nHL"
        stimLevel = stimLevel
    else
        error("`stimLevel` needs to be specified in `SPL` or `nHL`")
    end

    #equation derived from Prosser, S., & Arslan, E. (1987). Prediction of auditory brainstem wave V latency as a diagnostic tool of sensorineural hearing loss. Audiology : Official Organ of the International Society of Audiology, 26(3), 179–87. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/3662941
    avPeakVLat = (-5.859e-08*stimLevel^4 + 1.274e-05*stimLevel^3 - 0.0005424*stimLevel^2 - 0.05297*stimLevel + 9.3)/1000
    avPeakVLat = avPeakVLat + waveVLatencyOffset

    minPeakVLat = avPeakVLat-(0.00024*3)
    maxPeakVLat = avPeakVLat+(0.00024*3)
    #todo
    #shorten peak II trough II latency

    ##data from Picton book
    avPeakVpeakI_IPL = 0.00395
    avPeakVpeakIII_IPL = 0.00207

    minPeakVpeakI_IPL = avPeakVpeakI_IPL - (0.00023*2)
    maxPeakVpeakI_IPL = avPeakVpeakI_IPL + (0.00023*2)
    minPeakVpeakIII_IPL = avPeakVpeakIII_IPL - (0.00019*2)
    maxPeakVpeakIII_IPL = avPeakVpeakIII_IPL + (0.00019*2)

    ## guesstimates
    minPeakITroughILat = 0.25/1000
    maxPeakITroughILat = minPeakITroughILat+0.75/1000

    minPeakIIITroughIIILat = 0.25/1000
    maxPeakIIITroughIIILat = minPeakIIITroughIIILat+0.75/1000

    minPeakVTroughVLat = 0.25/1000
    maxPeakVTroughVLat = 1.75/1000

    minPeakIITroughIILat = 0.125/1000
    minPeakIVTroughIVLat = 0.125/1000

    # find all peaks, troughs and inflection points in the waveform
    peakPnts, peakTimes = findPeaks(sig, sampRate, epochStart=epochStart)
    troughPnts, troughTimes = findTroughs(sig, sampRate, epochStart=epochStart)
    inflPnts, inflTimes = findInflections(sig, sampRate, epochStart=epochStart)
    #dy = zeros(size(sig)); dy[:,2:size(sig)[2]] = diff(sig, 2)
    dy = zeros(size(sig)); dy[2:end] = diff(sig)

    peakVPnt, peakVTime = selectLargestPeakInWindow(sig, peakPnts, peakTimes, minPeakVLat, maxPeakVLat)
    if ismissing(peakVPnt) == true
        peakIPnt = missing
        peakITime = missing
        peakIIPnt = missing
        peakIITime = missing
        peakIIIPnt = missing
        peakIIITime = missing
        peakIVPnt = missing
        peakIVTime = missing

        troughIPnt = missing
        troughITime = missing
        troughIIPnt = missing
        troughIITime = missing
        troughIIIPnt = missing
        troughIIITime = missing
        troughIVPnt = missing
        troughIVTime = missing
        troughVPnt = missing
        troughVTime = missing
        minPeakILat = missing
        maxPeakILat = missing
        minPeakIILat = missing
        maxPeakIILat = missing
        minPeakIIILat = missing
        maxPeakIIILat = missing
        minPeakIVLat = missing
        maxPeakIVLat = missing

        minTroughILat = missing
        maxTroughILat = missing
        minTroughIILat = missing
        maxTroughIILat = missing
        minTroughIIILat = missing
        maxTroughIIILat = missing
        minTroughIVLat = missing
        maxTroughIVLat = missing
        minTroughVLat = missing
        maxTroughVLat = missing
    else
        minPeakILat = peakVTime - maxPeakVpeakI_IPL
        maxPeakILat = peakVTime - minPeakVpeakI_IPL
        minPeakIIILat = peakVTime - maxPeakVpeakIII_IPL
        maxPeakIIILat = peakVTime - minPeakVpeakIII_IPL

        minTroughVLat = peakVTime + minPeakVTroughVLat
        maxTroughVLat = peakVTime + maxPeakVTroughVLat

        #Peak I
        peakIPnt, peakITime = selectLargestPeakInWindow(sig, peakPnts, peakTimes, minPeakILat, maxPeakILat)
        #Peak III
        peakIIIPnt, peakIIITime = selectLargestPeakInWindow(sig, peakPnts, peakTimes, minPeakIIILat, maxPeakIIILat)

        #Trough V
        troughVPnt, troughVTime = selectLargestTroughInWindow(sig, troughPnts, troughTimes, minTroughVLat, maxTroughVLat, sig[peakVPnt])
        if ismissing(troughVPnt)
            troughVPnt, troughVTime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minTroughVLat, maxTroughVLat, maxAmp=sig[peakVPnt])
        end

        #Trough I
        if ismissing(peakIPnt) == true
            troughIPnt = missing
            troughITime = missing
            minTroughILat = missing
            maxTroughILat = missing
        else
            minTroughILat = peakITime + minPeakITroughILat
            maxTroughILat = peakITime + maxPeakITroughILat
            troughIPnt, troughITime = selectLargestTroughInWindow(sig, troughPnts, troughTimes, minTroughILat, maxTroughILat, sig[peakIPnt])
            if ismissing(troughIPnt)
                troughIPnt, troughITime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minTroughILat, maxTroughILat, maxAmp=sig[peakIPnt])
            end
        end

        #Trough III
        if ismissing(peakIIIPnt) == true
            troughIIIPnt = missing
            troughIIITime = missing
            minTroughIIILat = missing
            maxTroughIIILat = missing
        else
            minTroughIIILat = peakIIITime + minPeakIIITroughIIILat
            maxTroughIIILat = peakIIITime + maxPeakIIITroughIIILat
            troughIIIPnt, troughIIITime = selectLargestTroughInWindow(sig, troughPnts, troughTimes, minTroughIIILat, maxTroughIIILat, sig[peakIIIPnt])
            if ismissing(troughIIIPnt)
                troughIIIPnt, troughIIITime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minTroughIIILat, maxTroughIIILat, maxAmp=sig[peakIIIPnt])
            end
        end

        #Peak II
        minPeakIILat = missing
        maxPeakIILat = missing
        if ismissing(peakIIITime) == false
            maxPeakIILat = peakIIITime - 0.25/1000
        end
        if ismissing(troughITime) == false
            minPeakIILat = troughITime + 0.25/1000
        end
        if ismissing(maxPeakIILat) == true
            if ismissing(peakVTime) == false
                if ismissing(minPeakIILat) == false
                    maxPeakIILat = (peakVTime + minPeakIILat)/2.5
                end
            end
        end

        peakIIPnt=missing; peakIITime=missing
        if ((ismissing(minPeakIILat) == false) .& (ismissing(maxPeakIILat) == false))
            peakIIPnt, peakIITime = selectLargestPeakInWindow(sig, peakPnts, peakTimes, minPeakIILat, maxPeakIILat)
            if ismissing(peakIIPnt)
                peakIIPnt, peakIITime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minPeakIILat, maxPeakIILat)
            end
        end

        #Trough II
        minTroughIILat = missing
        maxTroughIILat = missing
        if ismissing(peakIIITime) == false
            maxTroughIILat = peakIIITime - 0.25/1000
        end
        if ismissing(peakIITime) == false
            minTroughIILat = peakIITime + minPeakIITroughIILat
        end

        troughIIPnt=missing; troughIITime=missing
        if ((ismissing(minTroughIILat) == false) .& (ismissing(maxTroughIILat) == false))
            troughIIPnt, troughIITime = selectLargestTroughInWindow(sig, troughPnts, troughTimes, minTroughIILat, maxTroughIILat, sig[peakIIPnt])
            if ismissing(troughIIPnt)
                troughIIPnt, troughIITime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minTroughIILat, maxTroughIILat, maxAmp=sig[peakIIPnt])
            end
        end

        #peak IV
        minPeakIVLat = missing
        maxPeakIVLat = missing
        if ismissing(peakVTime) == false
            maxPeakIVLat = peakVTime - 0.25/1000
        end
        if ismissing(troughIIITime) == false
            minPeakIVLat = troughIIITime + 0.25/1000
        end

        peakIVPnt=missing; peakIVTime=missing
        if ((ismissing(minPeakIVLat) == false) .& (ismissing(maxPeakIVLat) == false))
            peakIVPnt, peakIVTime = selectLargestPeakInWindow(sig, peakPnts, peakTimes, minPeakIVLat, maxPeakIVLat)
            if ismissing(peakIVPnt)
                peakIVPnt, peakIVTime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minPeakIVLat, maxPeakIVLat)
            end
        end

        #trough IV
        minTroughIVLat = missing
        maxTroughIVLat = missing
        if ismissing(peakVTime) == false
            maxTroughIVLat = peakVTime - 0.25/1000
        end
        if ismissing(peakIVTime) == false
            minTroughIVLat = peakIVTime + minPeakIVTroughIVLat
        end

        troughIVPnt=missing; troughIVTime=missing
        if ((ismissing(minTroughIVLat) == false) .& (ismissing(maxTroughIVLat) == false))
            troughIVPnt, troughIVTime = selectLargestTroughInWindow(sig, troughPnts, troughTimes, minTroughIVLat, maxTroughIVLat, sig[peakIVPnt])
            if ismissing(troughIVPnt)
                troughIVPnt, troughIVTime = selectStrongestInflectionInWindow(sig, dy, inflPnts, inflTimes, minTroughIVLat, maxTroughIVLat, maxAmp=sig[peakIVPnt])
            end
        end
    end

    peakPoints = [peakIPnt; peakIIPnt; peakIIIPnt; peakIVPnt; peakVPnt]
    peakLatencies = [peakITime; peakIITime; peakIIITime; peakIVTime; peakVTime]
    troughPoints = [troughIPnt; troughIIPnt; troughIIIPnt; troughIVPnt; troughVPnt]
    troughLatencies = [troughITime; troughIITime; troughIIITime; troughIVTime; troughVTime]

    peakAmps = missings(Float64, length(peakPoints))
    troughAmps = missings(Float64, length(peakPoints))
    #peakAmpsEstMax = zeros(length(peakPoints))
    #troughAmpsEstMax = zeros(length(peakPoints))
    for i=1:length(peakAmps)
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            #peakAmpsEstMax[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            #peakAmpsEstMax[i] = maximum(sig[,)
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
        end
    end

    peakTroughAmps = peakAmps-troughAmps

    waveLabels = ["I", "II", "III", "IV", "V"]

    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps,
                   minPeakLat=[minPeakILat, minPeakIILat, minPeakIIILat, minPeakIVLat, minPeakVLat],
                   maxPeakLat=[maxPeakILat, maxPeakIILat, maxPeakIIILat, maxPeakIVLat, maxPeakVLat],
                   minTroughLat=[minTroughILat, minTroughIILat, minTroughIIILat, minTroughIVLat, minTroughVLat],
                   maxTroughLat=[maxTroughILat, maxTroughIILat, maxTroughIIILat, maxTroughIVLat, maxTroughVLat])

    return df
end


text:
push!(LOAD_PATH, "./")
using Weno
using Burghers
using Convection
using RK3TVD
function Init!(X,L)
    size=length(X)
    h=L/size
    for i=0:size-1
        if i>floor(size//8) && i<floor(size//2)+floor(size//8)
            X[i+1]=1.-2*(i-floor(size/8))*h/L;
        else
            X[i+1]=0.0
        end
    end
end

@inline NumfluxGodunov(Fminf,Fmaxf,a::Float64,b::Float64)=
                     ifelse(a<=b,Fminf(a,b),Fmaxf(b,a))
@inline NumfluxLaxFriedrichs(F,a::Float64,b::Float64,alpha::Float64)=
                      0.5*(F(a)+F(b) - alpha*(b-a))

const size=1000
const L=1.
const T=1.0
const dt=0.8/size
println("size= ",size," dt= ",dt, " nteps= ",floor(T/dt))

In=Array{Float64}(size)
Init!(In,L)
Out=Array{Float64}(size)

W=WenoData(size)
R=RK3TVDData(size)

# Choose equation:
#const EQ=Convection
const EQ=Burghers

# Choose numerical flux:
const NF=NumfluxGodunov
#const NF=NumfluxLaxFriedrichs

if NF==NumfluxGodunov
    @inline NumFlux(X::Float64,Y::Float64)=NumfluxGodunov(EQ.minf,EQ.maxf,X,Y)
else
    @inline NumFlux(X::Float64,Y::Float64)=NumfluxLaxFriedrichs(EQ.flux,X,Y,1.0)
end

S(X::Array{Float64,1},Y::Array{Float64,1})=weno!(W,NumFlux,L,X,Y)

# be sure to run once before actually running the benchmark!
Init!(In,L)
Rk3tvd!(R,S,dt,In,Out)
Out,In=In,Out
#
Init!(In,L)

f=open("gp0","w")
writedlm(f, In)
close(f)

println("start computation")
t=0.
t1 = time_ns()

#@profile
while t<T
    Rk3tvd!(R,S,dt,In,Out)
    Out,In=In,Out
    t+=dt
end
#Profile.print(format=:flat)
tcomp = (time_ns() - t1)*1.e-9

println("computing time: ",tcomp," seconds.")

f=open("gp","w")
writedlm(f, In)
close(f)
println("A file 'gp' with the final solution was created.")
open("RunningOn"*gethostname(),"w") do f
    write(f,string(tcomp),"\n")
end

max_stars_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_issues_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_forks_repo_name: wolfv/BenchmarksPythonJuliaAndCo


text:
@testset "POMCPOW" begin
    pomdp = SingleCovidPOMDP(actions = CovidPOMDPs.CovidActionSpace(7;zero=true))
    b0 = initialstate(pomdp)

    sol = POMCPOWSolver(
        criterion = MaxUCB(2.0),
        max_depth = 30,
        tree_queries = 1_000,
        check_repeat_act = true,
        check_repeat_obs = false,
        enable_action_pw = false,
        k_observation = 3,
        alpha_observation = 0.05,
        estimate_value = CovidPOMDPs.ProportionalControlSolver(10.0)
    )

    planner = solve(sol, pomdp)

    a = action(planner, b0)
    @test a isa CovidAction
    @test 0. ≤ a.testing_prop ≤ 1.

    hist = simulate(
        pomdp,
        b0,
        planner;
        T = 50,
        upd = CovidPOMDPs.NoisyCovidFilter,
        n_p = 1_000,
        progress = false
    )

    validate_history(hist)
end


text:
using ReverseDiff: TrackedArray, @grad
const data = ReverseDiff.value
const revtrack = ReverseDiff.track

nobacksies(f, x) = revtrack(nobacksies, f, x)
nobacksies(f, xs::Tuple) = map(x -> nobacksies(f, x), xs)
@grad nobacksies(f::Symbol, x) = data(x), Δ -> error("Nested AD not defined for $f")
@grad nobacksies(f::String, x) = data(x), Δ -> error(f)

softmax(xs::TrackedArray) = revtrack(softmax, xs)

@grad softmax(xs) = softmax(data(xs)), Δ -> (nobacksies(:softmax, ∇softmax(data(Δ), data(xs))),)

logsoftmax(xs::TrackedArray) = revtrack(logsoftmax, xs)

@grad logsoftmax(xs) = logsoftmax(data(xs)), Δ -> (nobacksies(:logsoftmax, ∇logsoftmax(data(Δ), data(xs))),)

depthwiseconv(x::TrackedArray, w::TrackedArray, cdims::DepthwiseConvDims; kw...) = revtrack(depthwiseconv, x, w, cdims; kw...)
depthwiseconv(x::AbstractArray, w::TrackedArray, cdims::DepthwiseConvDims; kw...) = revtrack(depthwiseconv, x, w, cdims; kw...)
depthwiseconv(x::TrackedArray, w::AbstractArray, cdims::DepthwiseConvDims; kw...) = revtrack(depthwiseconv, x, w, cdims; kw...)

@grad depthwiseconv(x, w, cdims::DepthwiseConvDims; kw...) =
  depthwiseconv(data(x), data(w), cdims; kw...),
    Δ -> nobacksies(:depthwiseconv,
      (NNlib.∇depthwiseconv_data(data.((Δ, w))..., cdims; kw...),
       NNlib.∇depthwiseconv_filter(data.((x, Δ))..., cdims; kw...),
       nothing))

conv(x::TrackedArray,  w::TrackedArray, cdims::DenseConvDims;  kw...) = revtrack(conv, x, w, cdims; kw...)
conv(x::AbstractArray, w::TrackedArray, cdims::DenseConvDims;  kw...) = revtrack(conv, x, w, cdims; kw...)
conv(x::TrackedArray,  w::AbstractArray, cdims::DenseConvDims; kw...) = revtrack(conv, x, w, cdims; kw...)

@grad conv(x, w, cdims::DenseConvDims; kw...) =
  conv(data(x), data(w), cdims; kw...),
    Δ -> nobacksies(:conv,
      (NNlib.∇conv_data(data.((Δ, w))..., cdims; kw...),
       NNlib.∇conv_filter(data.((x, Δ))..., cdims; kw...),
       nothing))

∇conv_data(x::TrackedArray,  w::TrackedArray, cdims::DenseConvDims;  kw...) = revtrack(∇conv_data, x, w, cdims; kw...)
∇conv_data(x::AbstractArray, w::TrackedArray, cdims::DenseConvDims;  kw...) = revtrack(∇conv_data, x, w, cdims; kw...)
∇conv_data(x::TrackedArray,  w::AbstractArray, cdims::DenseConvDims; kw...) = revtrack(∇conv_data, x, w, cdims; kw...)

@grad function ∇conv_data(y, w, cdims::DenseConvDims; kw...)
  return (
    ∇conv_data(data(y), data(w), cdims; kw...),
    Δ -> begin
      return nobacksies(:conv,
        (NNlib.conv(data.((Δ, w))..., cdims; kw...),
         NNlib.∇conv_filter(data.((Δ, y))..., cdims; kw...),
         nothing)
      )
    end
  )
end

maxpool(x::TrackedArray, pdims::PoolDims; kw...) = revtrack(maxpool, x, pdims; kw...)

@grad function maxpool(x, pdims::PoolDims; kw...)
  y = maxpool(data(x), pdims; kw...)
  y, Δ -> (nobacksies(:maxpool, NNlib.∇maxpool(data.((Δ, y, x))..., pdims; kw...)), nothing)
end

meanpool(x::TrackedArray, pdims::PoolDims; kw...) = revtrack(meanpool, x, pdims; kw...)

@grad function meanpool(x, pdims::PoolDims; kw...)
  y = meanpool(data(x), pdims; kw...)
  y, Δ -> (nobacksies(:meanpool, NNlib.∇meanpool(data.((Δ, y, x))..., pdims; kw...)), nothing)
end


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Compressible Euler Multicomponent" begin
  # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the
  #       simulations are set up with dimensional states. For example, the reference pressure in SI
  #       units is 101325 Pa, i.e., pressure has values of O(10^5)

  @trixi_testset "elixir_eulermulti_shock_bubble.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"),
      l2   = [7.38219722e+01, 8.80169627e-01, 5.79496021e+04, 1.83101743e-01, 1.11794879e-02],
      linf = [1.96834444e+02, 7.19013908e+00, 1.58095793e+05, 8.19793567e-01, 8.00831673e-02],
      tspan = (0.0, 0.001))
  end

  @trixi_testset "elixir_eulermulti_shock_bubble_3components.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble_3components.jl"),
      l2   = [6.69469420e+01, 1.24137635e+00, 5.19736818e+04, 1.67537081e-01, 6.94144768e-03, 7.33198795e-02],
      linf = [2.14068118e+02, 1.82870581e+01, 1.93352191e+05, 1.02911816e+00, 7.00111358e-02, 9.34275213e-01],
      tspan = (0.0, 0.001))
  end

  @trixi_testset "elixir_eulermulti_shock_bubble_7components.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble_7components.jl"),
    l2   = [6.70895845e+01, 1.50857110e+00, 5.35958153e+04, 1.88695032e-01, 5.30253136e-03, 3.89776034e-02, 7.99227155e-02, 3.50734339e-03, 3.06391963e-02, 1.17354254e-01],
    linf = [2.27410552e+02, 1.95446320e+01, 2.96650723e+05, 1.09905431e+00, 4.40228926e-02, 5.08315400e-01, 1.06801398e+00, 1.76102931e-02, 4.01158590e-01, 1.47984635e+00],
    tspan = (0.0, 0.001))
  end

  @trixi_testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [5.01948074e-02, 5.02023248e-02, 2.25886833e-01, 6.17286464e-02],
      linf = [3.06937711e-01, 3.06807092e-01, 1.06295287e+00, 2.98135725e-01])
  end

  @trixi_testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [4.96608585e-02, 4.96617550e-02, 2.24265059e-01, 4.08685600e-03, 8.17371200e-03, 1.63474240e-02, 3.26948480e-02],
      linf = [2.48820395e-01, 2.48311640e-01, 9.35790843e-01, 1.72574990e-02, 3.45149980e-02, 6.90299959e-02, 1.38059992e-01])
  end

  @trixi_testset "elixir_eulermulti_eoc_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_eoc_ec.jl"),
      l2   = [1.15269058e-04, 1.15269058e-04, 1.66578747e-04, 4.61097337e-05, 9.22194673e-05],
      linf = [5.05618608e-04, 5.05618608e-04, 7.87591809e-04, 2.06476960e-04, 4.12953919e-04])
  end

  @trixi_testset "elixir_eulermulti_eoc_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_eoc_es.jl"),
      l2   = [1.74779183e-06, 1.74779183e-06, 5.62550396e-06, 6.15682712e-07, 1.23136542e-06],
      linf = [1.51966473e-05, 1.51966473e-05, 5.51065526e-05, 5.58096521e-06, 1.11619304e-05])
  end
end

end # module


text:
"""
    const FASTAI_DATA_REGISTRY

The default `DataRegistry` containing every dataset in
the fastai dataset collection.
"""
const FASTAI_DATA_REGISTRY = DatasetRegistry(
    Dict(d => () -> datasetpath(d) for d in DATASETS),
    Dict{String,Vector{DatasetRecipe}}(),
)


text:
## target machine

export TargetMachine, dispose,
       target, triple, cpu, features, asm_verbosity!,
       emit, add_transform_info!, add_library_info!
export JITTargetMachine

@checked struct TargetMachine
    ref::API.LLVMTargetMachineRef
end

Base.unsafe_convert(::Type{API.LLVMTargetMachineRef}, tm::TargetMachine) = tm.ref

TargetMachine(t::Target, triple::String, cpu::String="", features::String="";
              optlevel::API.LLVMCodeGenOptLevel=API.LLVMCodeGenLevelDefault,
              reloc::API.LLVMRelocMode=API.LLVMRelocDefault,
              code::API.LLVMCodeModel=API.LLVMCodeModelDefault) =
    TargetMachine(API.LLVMCreateTargetMachine(t, triple, cpu, features, optlevel,
                                              reloc, code))

dispose(tm::TargetMachine) = API.LLVMDisposeTargetMachine(tm)

function TargetMachine(f::Core.Function, args...; kwargs...)
    tm = TargetMachine(args...; kwargs...)
    try
        f(tm)
    finally
        dispose(tm)
    end
end

target(tm::TargetMachine) = Target(API.LLVMGetTargetMachineTarget(tm))
triple(tm::TargetMachine) = unsafe_message(API.LLVMGetTargetMachineTriple(tm))
triple() = unsafe_message(API.LLVMGetDefaultTargetTriple())
normalize(triple) = unsafe_message(API.LLVMNormalizeTargetTriple(triple))
cpu(tm::TargetMachine) = unsafe_message(API.LLVMGetTargetMachineCPU(tm))
features(tm::TargetMachine) = unsafe_message(API.LLVMGetTargetMachineFeatureString(tm))

asm_verbosity!(tm::TargetMachine, verbose::Core.Bool) =
    API.LLVMSetTargetMachineAsmVerbosity(tm, convert(Bool, verbose))

function emit(tm::TargetMachine, mod::Module, filetype::API.LLVMCodeGenFileType)
    out_error = Ref{Cstring}()
    out_membuf = Ref{API.LLVMMemoryBufferRef}()
    status = convert(Core.Bool,
        API.LLVMTargetMachineEmitToMemoryBuffer(tm, mod, filetype,
                                                out_error, out_membuf))

    if status
        error = unsafe_message(out_error[])
        throw(LLVMException(error))
    end

    return convert(Vector{UInt8}, MemoryBuffer(out_membuf[]))
end

function emit(tm::TargetMachine, mod::Module, filetype::API.LLVMCodeGenFileType, path::String)
    out_error = Ref{Cstring}()
    status = convert(Core.Bool,
        API.LLVMTargetMachineEmitToFile(tm, mod, path, filetype, out_error))

    if status
        error = unsafe_message(out_error[])
        throw(LLVMException(error))
    end

    return nothing
end

function add_transform_info!(pm::PassManager, tm::Union{Nothing,TargetMachine}=nothing)
    if tm !== nothing
        API.LLVMAddAnalysisPasses(tm, pm)
    else
        API.LLVMExtraAddGenericAnalysisPasses(pm)
    end
end
add_library_info!(pm::PassManager, triple::String) =
    API.LLVMAddTargetLibraryInfoByTriple(triple, pm)

function JITTargetMachine(triple = LLVM.triple(),
                          cpu = "", features = "";
                          optlevel = LLVM.API.LLVMCodeGenLevelDefault)

    # Force ELF on windows,
    # Note: Without this call to normalize Orc get's confused
    #       and chooses the x86_64 SysV ABI on Win x64
    triple = LLVM.normalize(triple)
    if Sys.iswindows()
        triple *= "-elf"
    end
    target = LLVM.Target(triple=triple)
    @debug "Configuring OrcJIT with" triple cpu features optlevel

    tm = TargetMachine(target, triple, cpu, features;
                       optlevel,
                       reloc = LLVM.API.LLVMRelocStatic, # Generate simpler code for JIT
                       code = LLVM.API.LLVMCodeModelJITDefault, # Required to init TM as JIT
                       )
    return tm
end


text:
function build_Ac(
  num_nodes,  # number of nodes.
  branches,  # array of branches.
  num_branches  # number of branches (sum of phases).
  )  #=

builds the branch-node connectivity matrix Ac from an array of Branch and returns it.

Author: Taku Noda
Started on Nov. 5, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

Ac = spzeros( num_branches, num_nodes )

irow = 1
for branch in branches
  for iph = 1:( branch.nph )
    jcol = branch.lhs[ iph ]
    if jcol > 0
      Ac[ irow, jcol ] = 1
    end
    jcol = branch.rhs[ iph ]
    if jcol > 0
      Ac[ irow, jcol ] = -1
    end
    irow = irow + 1
  end
end

return Ac

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
using Test
using Random
#using LinearAlgebra
#Random.seed!(9874984737484)

for tests in [
			"runtestString.jl",		#
            #"runtestAdd.jl", 			#
			#"runtestLogicalOr.jl",		#
			#"runtestLogicalXor.jl", 	#
			#"runtestReverseString.jl",	#
			#"runtestLengthString.jl",	#
			#"runtestTimesTwo.jl",		#
			#"runtestTimesThree.jl"
			#"runtestLogicalAnd.jl",	#
			#"runtestRepeat.jl",		#
			#"runtestWarningCountdown.jl" # BROKEN ! Doesn't work, NENI
			#"runtestCountdown.jl"
			#"runtestExtractInQuote.jl",
			#"runtestExtractInQuoteInside.jl" # Never performed
			#"runtestFibonacci.jl"
			#"runtestTrimLefToFQuote.jl"
			#"runtestSubtract.jl"
			]
    include(tests)
end


text:
using ModernGL, GLAbstraction, GLFW, GeometryTypes

const GLA = GLAbstraction
const window = GLFW.Window(name="Example")
GLFW.MakeContextCurrent(window)
GLA.set_context!(window)

const vsh = GLA.vert"""
#version 150
in vec2 position;

void main(){
    gl_Position = vec4(position, 0, 1.0);
}
"""

const fsh = GLA.frag"""
#version 150
out vec4 outColor;

void main() {
    outColor = vec4(1.0, 1.0, 1.0, 1.0);
}
"""
prog = GLA.Program(vsh, fsh)
const triangle = GLA.VertexArray(GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position=Point2f0[(0.0, 0.5), (0.5, -0.5), (-0.5,-0.5)]))

glClearColor(0, 0, 0, 1)
GLA.bind(prog)
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    GLA.bind(triangle)
    GLA.draw(triangle)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end
GLFW.DestroyWindow(window)


text:
@testset "Testing the HMatrix builders" begin
    # builder 'single value'
    H = HMatrix(10,20,1e-3,1.0)
    @test norm(data(H) - ones(10,20)) ≤ 1e-12
end

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
"""
	corrp(x,y)
Compute Pearson's correlation coefficient (`r`) and corresponding `p` value with
null hypothesis of no correlation (against alternative of a non-zero correlation).
`p` values (two-tailed) below 0.05 typically indicate that `r` is statistically
significant, see e.g., https://www.medcalc.org/manual/correlation.php

**Example**
```
x,y = collect(1.:1:10),collect(0.:10:90);
y[end] = 80;
r,p = corrp(x,y);
```
"""
function corrp(x::Vector{Float64},y::Vector{Float64})
	xi,yi = HypoTest.prepdata(x,y);
	n = length(xi);
	r = Statistics.cor(xi,yi);
	t = HypoTest.tpearval(r,n);
	return r, HypoTest.pvalcorr(t,n)
end

# t-value for Pearson's correlation coefficient: https://www.danielsoper.com/statcalc/formulas.aspx?id=44
tpearval(r,n) = r/sqrt((1-r^2)/(n-2));

# two-tailed probability for correlation (=t-distribution)
pvalcorr(t,n) = 2*(1-Distributions.cdf(Distributions.TDist(n-2),t));

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
module CoffeeMachine

export Machine, Cup, makecoffee!, fillwith!, readdisplay

mutable struct Machine
    coffee::Float64
    milk::Float64
    message::String

    Machine(; coffee=0.0, milk=0.0) = new(coffee, milk, "")
end

struct Cup
    coffee::Float64
    milk::Float64
end

function fillwith!(m::Machine; coffee=0.0, milk=0.0)
    m.coffee += coffee
    m.milk += milk
end

function makecoffee!(m::Machine; withmilk=false) :: Union{Cup, Nothing}
    if m.coffee <= 0.0
        display!(m, "Out of coffee")
        return nothing
    end

    if withmilk && m.milk <= 0.0
        display!(m, "Out of milk")
        return nothing
    end

    milkincup = if withmilk
        1.0
    else
        0.0
    end

    m.coffee -= 1.0
    m.milk -= milkincup

    display!(m, "Enjoy")

    Cup(1.0, milkincup)
end

readdisplay(m::Machine) = m.message
display!(m::Machine, msg::String) = m.message = msg

end # module


text:
immutable NoncentralChisq <: ContinuousUnivariateDistribution
    df::Float64
    ncp::Float64
    function NoncentralChisq(d::Real, nc::Real)
    	d >= zero(d) && nc >= zero(nc) || error("df and ncp must be non-negative")
    	new(float64(d), float64(nc))
    end
end

mean(d::NoncentralChisq) = d.df + d.ncp
var(d::NoncentralChisq) = 2.0*(d.df + 2.0*d.ncp)
skewness(d::NoncentralChisq) = 2.0*√2*(d.df + 3.0*d.ncp)/sqrt(d.df + 2.0*d.ncp)^3
kurtosis(d::NoncentralChisq) = 12.0*(d.df + 4.0*d.ncp)/(d.df + 2.0*d.ncp)^2

@_jl_dist_2p NoncentralChisq nchisq

insupport(::NoncentralChisq, x::Real) = zero(x) < x < Inf
insupport(::Type{NoncentralChisq}, x::Real) = zero(x) < x < Inf


text:
# Example of vOptSpecific, for the 2LAP
# July 2017

n  =  3

C1 = [  3    9   7   ;
       16   10   6   ;
        2    7  11     ]

C2 = [ 16   15   6   ;
        5    7  13   ;
        1    2  13     ]

using vOptSpecific
m = set2LAP( n , C1 , C2 )
z1, z2, S = vSolve( m )

max_stars_repo_name: faithyyyy/vOptSolver

max_issues_repo_name: faithyyyy/vOptSolver

max_forks_repo_name: faithyyyy/vOptSolver


text:
using Functors, Test

struct Foo
  x
  y
end
@functor Foo

struct Bar
  x
end
@functor Bar

struct Baz
  x
  y
  z
end
@functor Baz (y,)

struct NoChildren
  x
  y
end

@testset "Nested" begin
  model = Bar(Foo(1, [1, 2, 3]))

  model′ = fmap(float, model)

  @test model.x.y == model′.x.y
  @test model′.x.y isa Vector{Float64}
end

@testset "Exclude" begin
  f(x::AbstractArray) = x
  f(x::Char) = 'z'

  x = ['a', 'b', 'c']
  @test fmap(f, x)  == ['z', 'z', 'z']
  @test fmap(f, x; exclude = x -> x isa AbstractArray) == x

  x = (['a', 'b', 'c'], ['d', 'e', 'f'])
  @test fmap(f, x)  == (['z', 'z', 'z'], ['z', 'z', 'z'])
  @test fmap(f, x; exclude = x -> x isa AbstractArray) == x
end

@testset "Walk" begin
  model = Foo((0, Bar([1, 2, 3])), [4, 5])

  model′ = fmapstructure(identity, model)
  @test model′ == (; x=(0, (; x=[1, 2, 3])), y=[4, 5])
end

@testset "Property list" begin
  model = Baz(1, 2, 3)
  model′ = fmap(x -> 2x, model)

  @test (model′.x, model′.y, model′.z) == (1, 4, 3)
end

@testset "fcollect" begin
  m1 = [1, 2, 3]
  m2 = 1
  m3 = Foo(m1, m2)
  m4 = Bar(m3)
  @test all(fcollect(m4) .=== [m4, m3, m1, m2])
  @test all(fcollect(m4, exclude = x -> x isa Array) .=== [m4, m3, m2])
  @test all(fcollect(m4, exclude = x -> x isa Foo) .=== [m4])

  m1 = [1, 2, 3]
  m2 = Bar(m1)
  m0 = NoChildren(:a, :b)
  m3 = Foo(m2, m0)
  m4 = Bar(m3)
  println(fcollect(m4))
  @test all(fcollect(m4) .=== [m4, m3, m2, m1, m0])
end


text:
using Base.Test
using Vinyl: @overdub, @hook

buf = IOBuffer()

struct TraceCtx end
@hook    TraceCtx   (f::Any)(xs...) = println(buf, "Called $(:($f($(xs...))))")
@overdub TraceCtx() 1+2

@test String(take!(buf)) == """
Called (+)(1, 2)
Called (add_int)(1, 2)
"""


text:
"""
test_station(model_path, stn_df, ycol, stn_code, stn_name,
    sample_size, output_size,
    output_dir, output_prefix;
    test_fdate=ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul"),
    test_tdate==ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul"),
    eltype=Float32,
    is_zscored=false)

Dirty codes for post processing
"""
function test_station(model_path::String, stn_df::DataFrame,
    ycol::Symbol, stn_code::Integer, stn_name::String,
    sample_size::Integer, output_size::Integer,
    output_dir::String, output_prefix::String,
    test_fdate::ZonedDateTime=ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul"),
    test_tdate::ZonedDateTime=ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul"),
    eltype::DataType=Float32, is_normalize::Bool=false)

    filepath = model_path * string(ycol) * ".bson"

    BSON.@load filepath cpu_model weights μ σ
    Flux.loadparams!(cpu_model, weights)
    model = cpu_model |> gpu

    # norm_features and features
    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    norm_prefix = "norm_"
    _norm_feas = [Symbol(eval(norm_prefix * String(f))) for f in features]

    if (is_normalize == false)
        min_max_scaling!(stn_df, features, _norm_feas)
    end

    # type conversion
    for fea in features
        stn_df[!, fea] = eltype.(stn_df[!, fea])
    end

    for nfea in _norm_feas
        stn_df[!, nfea] = eltype.(stn_df[!, nfea])
    end

    # remove target (ycol)
    norm_ycol = Symbol(norm_prefix, ycol)
    norm_feas = copy(_norm_feas)
    # remove ycol itself
    # deleteat!(norm_feas, findall(x -> x == norm_ycol, norm_feas))
    test_wd = window_df(stn_df, sample_size, output_size, test_fdate, test_tdate)
    test_wd_size = length(test_wd)

    μσs = mean_and_std_cols(stn_df, features)
    stn_μ = μσs[String(ycol), "μ"].value
    stn_σ = μσs[String(ycol), "σ"].value
    total_min = float(μσs[String(ycol), "minimum"].value)
    total_max = float(μσs[String(ycol), "maximum"].value)

    μσ = ndsparse((
        dataset = ["total", "total"],
        type = ["μ", "σ"]),
        (value = [μ, σ],))

    @info "    Construct Test Set..."
    p = Progress(length(test_wd), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    test_set = [(ProgressMeter.next!(p);
        make_pair_DNN(df, norm_ycol, norm_feas, sample_size, output_size, eltype)) for df in test_wd]

    _RMSE = RMSE(test_set, model, μσ)
    _MAE = MAE(test_set, model, μσ)
    _MSPE = MSPE(test_set, model, μσ)
    _MAPE = MAPE(test_set, model, μσ)

    #test_dates = collect(test_fdate:Hour(1):test_tdate)
    test_dates = collect(test_fdate + Hour(sample_size - 1):Hour(1):test_tdate - Hour(output_size))

    #dnn_table = predict_model(test_set, model, ycol, μ, σ, output_size, output_dir)
    dnn_table = predict_model_minmax(test_set, model, ycol, total_min, total_max,
        output_size, "/mnt/")
    dfs_out = export_CSV(DateTime.(test_dates), dnn_table, ycol, output_size, output_dir, output_prefix)
    df_corr = compute_corr(dnn_table, output_size, output_dir, output_prefix)

    plot_datefmt = @dateformat_str "yyyymmddHH"
    test_fdate_str = Dates.format(test_fdate, plot_datefmt)
    test_tdate_str = Dates.format(test_tdate, plot_datefmt)
    plot_DNN_lineplot(DateTime.(test_dates), dnn_table,
        DateTime(test_fdate), DateTime(test_tdate), ycol, output_size, output_dir,
        output_prefix * "_$(test_fdate_str)_$(test_tdate_str)")

    @info " $(string(ycol)) RMSE        : ", _RMSE
    @info " $(string(ycol)) MAE         : ", _MAE
    @info " $(string(ycol)) MSPE        : ", _MSPE
    @info " $(string(ycol)) MAPE        : ", _MAPE

    stn_lat = stn_df[1, :lat]
    stn_lon = stn_df[1, :lon]
    return [stn_code, stn_name, stn_lat, stn_lon, string(ycol), _RMSE, _MAE, _MSPE, _MAPE]
end

function test_station(model_path::String, stn_df_in::DataFrame, stn_df_out::DataFrame,
    ycol::Symbol, stn_code::Integer, stn_name::String,
    sample_size::Integer, output_size::Integer,
    output_dir::String, output_prefix::String,
    test_fdate::ZonedDateTime=ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul"),
    test_tdate::ZonedDateTime=ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul"),
    eltype::DataType=Float32, is_normalize::Bool=false)

    filepath = model_path * string(ycol) * ".bson"

    BSON.@load filepath cpu_model weights μ σ
    Flux.loadparams!(cpu_model, weights)
    model = cpu_model |> gpu

    # create dataset
    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    norm_prefix = "norm_"
    _norm_feas = [Symbol(eval(norm_prefix * String(f))) for f in features]

    if (is_normalize == false)
        #zscore!(stn_df_in, features, _norm_feas)
        min_max_scaling!(stn_df_in, features, _norm_feas)
    end

    # type conversion
    for fea in features
        stn_df_in[!, fea] = eltype.(stn_df_in[!, fea])
    end

    for nfea in _norm_feas
        stn_df_in[!, nfea] = eltype.(stn_df_in[!, nfea])
    end

    # remove target (ycol)
    norm_ycol = Symbol(norm_prefix, ycol)
    norm_feas = copy(_norm_feas)
    # remove ycol itself
    # deleteat!(norm_feas, findall(x -> x == norm_ycol, norm_feas))
    test_wd = window_df(stn_df_in, sample_size, output_size, test_fdate, test_tdate)
    test_wd_size = length(test_wd)

    μσs = mean_and_std_cols(stn_df_out, features)
    stn_μ = μσs[String(ycol), "μ"].value
    stn_σ = μσs[String(ycol), "σ"].value
    total_min = float(μσs[String(ycol), "minimum"].value)
    total_max = float(μσs[String(ycol), "maximum"].value)

    μσ = ndsparse((
        dataset = ["total", "total"],
        type = ["μ", "σ"]),
        (value = [μ, σ],))

    p = Progress(length(test_wd), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    test_set = [(ProgressMeter.next!(p);
        @show length(df), output_size;
        make_pair_DNN(df, norm_ycol, norm_feas, sample_size, output_size, eltype)) for df in test_wd]

    _RMSE = RMSE(test_set, model, μσ)
    _MAE = MAE(test_set, model, μσ)
    _MSPE = MSPE(test_set, model, μσ)
    _MAPE = MAPE(test_set, model, μσ)

    #test_dates = collect(test_fdate:Hour(1):test_tdate)
    test_dates = collect(test_fdate + Hour(sample_size - 1):Hour(1):test_tdate - Hour(output_size))

    #dnn_table = predict_model(test_set, model, ycol, μ, σ, output_size, output_dir)
    dnn_table = predict_model_minmax(test_set, model, ycol, total_min, total_max,
        output_size, "/mnt/")
    dfs_out = export_CSV(DateTime.(test_dates), dnn_table, ycol, output_size, output_dir, output_prefix)
    df_corr = compute_corr(dnn_table, output_size, output_dir, output_prefix)

    plot_datefmt = @dateformat_str "yyyymmddHH"
    test_fdate_str = Dates.format(test_fdate, plot_datefmt)
    test_tdate_str = Dates.format(test_tdate, plot_datefmt)
    plot_DNN_lineplot(DateTime.(test_dates), dnn_table,
        DateTime(test_fdate), DateTime(test_tdate), ycol, output_size, output_dir,
        output_prefix * "_$(test_fdate_str)_$(test_tdate_str)")

    @info " $(string(ycol)) RMSE        : ", _RMSE
    @info " $(string(ycol)) MAE         : ", _MAE
    @info " $(string(ycol)) MSPE        : ", _MSPE
    @info " $(string(ycol)) MAPE        : ", _MAPE

    stn_lat = stn_df_in[1, :lat]
    stn_lon = stn_df_in[1, :lon]
    return [stn_code, stn_name, stn_lat, stn_lon, string(ycol), _RMSE, _MAE, _MSPE, _MAPE]
end

function test_classification(model_path::String, stn_df::DataFrame, ycol::Symbol, stn_code::Integer, stn_name::String,
    sample_size::Integer, output_size::Integer,
    output_dir::String, output_prefix::String,
    test_fdate::ZonedDateTime=ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul"),
    test_tdate::ZonedDateTime=ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul"),
    eltype::DataType=Float32, is_zscored::Bool=false)

    filepath = model_path * string(ycol) * ".bson"

    BSON.@load filepath cpu_model weights μ σ
    Flux.loadparams!(cpu_model, weights)
    model = cpu_model |> gpu

    # create dataset
    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    norm_prefix = "norm_"
    _norm_feas = [Symbol(eval(norm_prefix * String(f))) for f in features]

    # type conversion
    for fea in features
        stn_df[!, fea] = eltype.(stn_df[!, fea])
    end

    for nfea in _norm_feas
        stn_df[!, nfea] = eltype.(stn_df[!, nfea])
    end

    if (is_zscored == false)
        zscore!(stn_df, features, _norm_feas)
    end
    # remove target (ycol)
    norm_ycol = Symbol(norm_prefix, ycol)
    norm_feas = copy(_norm_feas)
    # remove ycol itself
    # deleteat!(norm_feas, findall(x -> x == norm_ycol, norm_feas))
    test_wd_idxs = window_df(stn_df, sample_size, output_size, test_fdate, test_tdate)
    sidx = findrow(stn_df, :date, test_fdate)
    fidx = findrow(stn_df, :date, test_tdate)
    if (sidx == 0 || fidx == 0)
        @error "Error; cannot found date: ", sidx, fidx
        return
    end
    test_wd_size = length(test_wd_idxs)
    test_size = length(test_wd_idxs)
    test_idxs = create_idxs(test_wd_idxs, test_size)

    μσs = mean_and_std_cols(stn_df, features)
    stn_μ = μσs[String(ycol), "μ"].value
    stn_σ = μσs[String(ycol), "σ"].value

    μσ = ndsparse((
        dataset = ["total", "total"],
        type = ["μ", "σ"]),
        (value = [μ, σ],))

    @info "    Construct Test Set..."
    test_set = [make_pair_DNN(stn_df, norm_ycol, idx, norm_feas, sample_size, output_size) for idx in test_wd_idxs]

    forecast_all, forecast_high = classification(test_set, model, ycol)

    stn_lat = stn_df[1, :lat]
    stn_lon = stn_df[1, :lon]
    return [stn_code, stn_name, stn_lat, stn_lon, string(ycol), forecast_all, forecast_high]
end


text:
using Test
import ArchGDAL;
const AG = ArchGDAL;

@testset "test_utils.jl" begin
    @testset "metadataitem" begin
        driver = AG.getdriver("DERIVED")
        @test AG.metadataitem(driver, "DMD_EXTENSIONS") == ""
        driver = AG.getdriver("GTiff")
        @test AG.metadataitem(driver, "DMD_EXTENSIONS") == "tif tiff"
    end
end


text:
using CircGeometry, Test

p1 = Point(0.0,0.0); q1 = Point(0.0,1.0)
p2 = Point(-0.5,0.5); q2 = Point(0.5,0.5)
@testset "do_intersect tests" begin
    @test CircGeometry.do_intersect(p1,q1,p2,q2) == true
    @test CircGeometry.do_intersect(p1,p2,q1,q2) == false
end

z = Point(0.5+Float64(pi),1.5+Float64(pi))
@testset "get_orientation tests" begin
    @test CircGeometry.get_orientation(p2,q1,z) == 0
    @test CircGeometry.get_orientation(p2,z,q1) == 0
    @test CircGeometry.get_orientation(z,p2,q1) == 0
    @test CircGeometry.get_orientation(p1,q1,Point(1.0,0.0)) == 1
    @test CircGeometry.get_orientation(p1,q1,Point(-1.0,0.0)) == 2
end

a = Point(13.2,12.3)
b = Point(1.0+5*eps(),1.0)
c = Point(-1.0,0.0)
circlewall = CircGeometry.CircleWall(Point(0.0,0.0),1.0)
lwall1 = CircGeometry.LineWall(Point(-1.0,1.0),Point(1.0,1.0))
lwall2 = CircGeometry.LineWall(Point(-1.0,-1.0),Point(-1.0,1.0))
@testset "is_point_in_wall tests" begin
    @test CircGeometry.is_point_in_wall(q1,circlewall) == true
    @test CircGeometry.is_point_in_wall(p2,circlewall) == false
    @test CircGeometry.is_point_in_wall(a,circlewall) == false

    @test CircGeometry.is_point_in_wall(q1,lwall1) == true
    @test CircGeometry.is_point_in_wall(p1,lwall1) == false
    @test CircGeometry.is_point_in_wall(b,lwall1) == false

    @test CircGeometry.is_point_in_wall(c,lwall2) == true
    @test CircGeometry.is_point_in_wall(p1,lwall2) == false
end

# test for on_segment where point is colinear with other two points
polygon = [Point(0.0,0.0),Point(2.0,0.0),Point(2.0,1.0)]
point = Point(1.0,0.0)
@test CircGeometry.is_inside_polygon(polygon, point; extreme = Point(100_000.0, point.y)) == true

# check on_segment
@testset "on_segment tests" begin
    @test CircGeometry.on_segment(Point(0.0,0.0),Point(1.0,0.0),Point(2.0,0.0)) == true
    @test CircGeometry.on_segment(Point(0.0,0.0),Point(1.0,0.1),Point(2.0,0.0)) == false
    @test CircGeometry.on_segment(Point(0.0,0.0),Point(3.0,0.0),Point(2.0,0.0)) == false
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# This file was generated, do not modify it. # hide
df = DataFrame(data)
describe(df)


text:
baremodule ThreadsX

function collect end
function map end
function mapi end

function mapreduce end
function reduce end

# Functions derived from `reduce`
function sum end
function prod end
function count end
function maximum end
function minimum end

# Functions implemented directly in terms of transducers
function any end
function all end
function findfirst end
function findlast end
function findall end
# function findmax end
# function findmin end
# function argmax end
# function argmin end
function extrema end
function issorted end
function unique end
function Set end

function foreach end
function map! end
function copy! end
function copyto! end

function sort end
function sort! end

module Implementations
import SplittablesBase
using ArgCheck: @argcheck, @check
using BangBang: SingletonVector, append!!, push!!, union!!
using Base:
    HasShape, IteratorSize, Ordering, add_sum, mapreduce_empty, mul_prod, reduce_empty
using ConstructionBase: setproperties
using InitialValues: asmonoid
using MicroCollections: EmptyVector
using Referenceables: referenceable
using Setfield: @set
using Transducers:
    Cat,
    Empty,
    Filter,
    Init,
    Map,
    MapSplat,
    NondeterministicThreading,
    OnInit,
    ReduceIf,
    Transducers,
    extract_transducer,
    foldxl,
    opcompose,
    reduced,
    right,
    tcollect,
    tcopy
using ..ThreadsX

@static if VERSION >= v"1.3-alpha"
    using Base.Threads: @spawn
else
    # Mock `@spawn` using `@async`:
    @eval const $(Symbol("@spawn")) = $(Symbol("@async"))
end

if isdefined(Transducers, :foldxt)
    using Transducers: foldxt
else
    foldxt(rf, xf, xs; kw...) = reduce(rf, xf, xs; kw...)
    foldxt(rf, xs; kw...) = reduce(rf, Map(identity), xs; kw...)
end

include("debug.jl")
include("utils.jl")
include("basesizes.jl")
include("reduce.jl")
include("foreach.jl")
include("map.jl")
include("copy.jl")
include("mergesort.jl")
include("quicksort.jl")
include("countingsort.jl")
end  # module Implementations

const MergeSort = Implementations.ParallelMergeSortAlg()
const QuickSort = Implementations.ParallelQuickSortAlg()
const StableQuickSort =
    Implementations.ParallelQuickSortAlg(smallsort = MergeSort.smallsort)

Implementations.define_docstrings()

end # baremodule ThreadX


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rocm_smi_lib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rocm_smi_lib")
JLLWrappers.@generate_main_file("rocm_smi_lib", UUID("920c060a-fbb2-54e8-bd70-fc4b6562f351"))
end  # module rocm_smi_lib_jll


text:
using Test, AbrkgaObop

include("test_dataset_reader.jl")
include("test_utils.jl")
#include("test_parameters.jl")
#include("test_decoder.jl")
include("test_local_search.jl")
#include("test_abrkga.jl")

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
# Adapted from Example 3.99 of [BPT12].
#
# [BPT12] Blekherman, G.; Parrilo, P. & Thomas, R.
# Semidefinite Optimization and Convex Algebraic Geometry
# Society for Industrial and Applied Mathematics, 2012

using Test
using SumOfSquares
using DynamicPolynomials

function BPT12e399_test(optimizer, config::MOIT.TestConfig, remainder::Bool)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x y
    if remainder
        cref = @constraint(model, 10 - (x^2 + α*y) in SOSCone(), newton_of_remainder = true,
                           maxdegree = nothing, domain = @set x^2 + y^2 == 1)
    else
        cref = @constraint(model, 10 - (x^2 + α*y) in SOSCone(), maxdegree = 2, domain = @set x^2 + y^2 == 1)
    end

    @objective(model, Max, α)
    JuMP.optimize!(model)

    α_value = remainder ? 6.0 : 10.0

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ α_value atol=atol rtol=rtol

    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test JuMP.value(α) ≈ α_value atol=atol rtol=rtol

    @test_throws SumOfSquares.ValueNotSupported value(cref)
    p = gram_matrix(cref)
    if remainder
        @test getmat(p) ≈ [1 -3; -3 9] atol=atol rtol=rtol
        @test p.basis.monomials == [y, 1]
    else
        @test getmat(p) ≈ [4.0 0.0 0.0; 0.0 5.0 -5.0; 0.0 -5.0 5.0] atol=atol rtol=rtol
        @test p.basis.monomials == [x, y, 1]
    end

    @test dual_status(model) == MOI.FEASIBLE_POINT
    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ (remainder ? -8/3 : 0.0) atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == x^2
    @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ (remainder ? 1/3 : 1.0) atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1

    μ = moments(cref)
    @test μ isa AbstractMeasure{Float64}
    if remainder
        @test length(moments(μ)) == 3
        @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == y^2
        @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[2]) == y
        @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
        @test monomial(moments(μ)[3]) == 1
    else
        @test length(moments(μ)) == 5
        @test moment_value(moments(μ)[1]) ≈ 0 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == x*y
        @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[2]) == y^2
        @test moment_value(moments(μ)[3]) ≈ 0 atol=atol rtol=rtol
        @test monomial(moments(μ)[3]) == x
        @test moment_value(moments(μ)[4]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[4]) == y
        @test moment_value(moments(μ)[5]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[5]) == 1
    end

    @objective(model, Min, α)
    JuMP.optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ -α_value atol=atol rtol=rtol

    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test JuMP.value(α) ≈ -α_value atol=atol rtol=rtol

    @test_throws SumOfSquares.ValueNotSupported value(cref)
    p = gram_matrix(cref)
    if remainder
        @test getmat(p) ≈ [1 3; 3 9] atol=atol rtol=rtol
        @test p.basis.monomials == [y, 1]
    else
        @test getmat(p) ≈ [4.0 0.0 0.0; 0.0 5.0 5.0; 0.0 5.0 5.0] atol=atol rtol=rtol
        @test p.basis.monomials == [x, y, 1]
    end

    @test dual_status(model) == MOI.FEASIBLE_POINT
    μ = dual(cref)
    @test μ isa AbstractMeasure{Float64}
    @test length(moments(μ)) == 3
    @test moment_value(moments(μ)[1]) ≈ (remainder ? -8/3 : 0.0) atol=atol rtol=rtol
    @test monomial(moments(μ)[1]) == x^2
    @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
    @test monomial(moments(μ)[2]) == y
    @test moment_value(moments(μ)[3]) ≈ (remainder ? 1/3 : 1.0) atol=atol rtol=rtol
    @test monomial(moments(μ)[3]) == 1

    μ = moments(cref)
    @test μ isa AbstractMeasure{Float64}
    if remainder
        @test length(moments(μ)) == 3
        @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == y^2
        @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
        @test monomial(moments(μ)[2]) == y
        @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
        @test monomial(moments(μ)[3]) == 1
    else
        @test length(moments(μ)) == 5
        @test moment_value(moments(μ)[1]) ≈ 0 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == x*y
        @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[2]) == y^2
        @test moment_value(moments(μ)[3]) ≈ 0 atol=atol rtol=rtol
        @test monomial(moments(μ)[3]) == x
        @test moment_value(moments(μ)[4]) ≈ -1 atol=atol rtol=rtol
        @test monomial(moments(μ)[4]) == y
        @test moment_value(moments(μ)[5]) ≈ 1 atol=atol rtol=rtol
        @test monomial(moments(μ)[5]) == 1
    end

end

BPT12e399_rem_test(optimizer, config) = BPT12e399_test(optimizer, config, true)
sd_tests["BPT12e399_rem"] = BPT12e399_rem_test
BPT12e399_maxdegree_test(optimizer, config) = BPT12e399_test(optimizer, config, false)
sd_tests["BPT12e399_maxdegree"] = BPT12e399_maxdegree_test


text:
using QuasiGeostrophy, FFTW, Test
using LinearAlgebra
include(pwd() * "/test/test_utils.jl")
const tol = 1e1

## Define 1D Test
Ω = S¹(0, 2π)
Nx = 2^4;
fourier_grid = FourierGrid(Nx, Ω)
fieldnames = ("ϕ1", "ϕ2", "ϕ3", "ϕ4", "ϕ5")
ϕ1, ϕ2, ϕ3, ϕ4, ϕ5 = create_fields(names = fieldnames, grid = fourier_grid)
∂x = create_operators(fourier_grid)
# initialize FourierField data
x = fourier_grid.grid[1]
ϕ1(sin.(x)); ϕ2(cos.(x));
ϕ3( @. sin(x) * cos(x) );
ϕ4( @. sin(x) + cos(x) );
ϕ5( @. cos(x)^2 - sin(x)^2 );
# Wrap Around Impero Object
@wrapper ϕ1=ϕ1 ϕ2=ϕ2 ϕ3=ϕ3 ϕ4=ϕ4 ϕ5=ϕ5
Δ  = Operator(∂x^2)
∂x = Operator(∂x)

@testset "1D Algebra Test" begin
    println("checking ϕ1 * ϕ2 - ϕ3")
    @test norm(compute(ϕ1 * ϕ2 - ϕ3))/norm(compute(ϕ3)) < eps(tol)
    println("checking ϕ1 + ϕ2 - ϕ4")
    @test norm(compute(ϕ1 + ϕ2 - ϕ4))/norm(compute(ϕ4)) < eps(tol)
    println("checking ϕ1 + ϕ1 - 2 * ϕ1")
    @test norm(compute(ϕ1 + ϕ1- 2*ϕ1))/norm(compute(ϕ1)) < eps(tol)
end

@testset "1D Calculus Test" begin
    println("checking ∂x(ϕ1) - ϕ2")
    @test norm(compute(∂x(ϕ1) - ϕ2))/norm(compute(ϕ2)) < eps(tol)
    println("checking Δ(ϕ1) + ϕ1")
    @test norm(compute(Δ(ϕ1) + ϕ1))/norm(compute(ϕ1)) < eps(tol*10)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy