text:
# Methods for porting Indicators.jl functions to TS objects from Temporal.jl package

###### run.jl ######
mode{V}(x::TS{V}) = mode(x.values)
runmean{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runmean(x.values, n, cumulative), x.index, :RunMean)
runsum{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runsum(x.values, n, cumulative), x.index, :RunSum)
runmad{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runmad(x.values, n, cumulative), x.index, :RunMAD)
runvar{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runvar(x.values, n, cumulative), x.index, :RunVar)
runmax{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runmax(x.values, n, cumulative), x.index, :RunMax)
runmin{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runmin(x.values, n, cumulative), x.index, :RunMin)
runsd{V,T}(x::TS{V,T}, n::Int64=10, cumulative::Bool=true) = ts(runsd(x.values, n, cumulative), x.index, :RunSD)
wilder_sum{V,T}(x::TS{V,T}, n::Int64=10) = ts(runsum(x.values, n), x.index, :WilderSum)
function runcov{V,T}(x::TS{V,T}, y::TS{V,T}, n::Int64=10, cumulative::Bool=true)
    @assert size(x,2) == 1 || size(y,2) == 1 "Arguments x and y must both be univariate (have only one column)."
    z = [x y]
    return ts(runcov(z[:,1], z[:,2], n, cumulative), x.index, :RunCov)
end
function runcor{V,T}(x::TS{V,T}, y::TS{V,T}, n::Int64=10, cumulative::Bool=true)
    @assert size(x,2) == 1 || size(y,2) == 1 "Arguments x and y must both be univariate (have only one column)."
    z = [x y]
    ts(runcor(z[:,1], z[:,2], n, cumulative), x.index, :RunCor)
end

##### ma.jl ######
sma{V,T}(x::TS{V,T}, n::Int64=10) = ts(sma(x.values, n), x.index, :SMA)
hma{V,T}(x::TS{V,T}, n::Int64=10) = ts(hma(x.values, n), x.index, :HMA)
mma{V,T}(x::TS{V,T}, n::Int64=10) = ts(mma(x.values, n), x.index, :MMA)
swma{V,T}(x::TS{V,T}, n::Int64=10) = ts(swma(x.values, n), x.index, :SWMA)
kama{V,T}(x::TS{V,T}, n::Int64=10, nfast::Float64=0.6667, nslow::Float64=0.0645) = ts(kama(x.values, n, nfast, nslow), x.index, :KAMA)
alma{V,T}(x::TS{V,T}, n::Int64=9, offset::Float64=0.85, sigma::Float64=6.0) = ts(kama(x.values, n, offset, sigma), x.index, :ALMA)
trima{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(trima(x.values, n; args...), x.index, :TRIMA)
wma{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(wma(x.values, n; args...), x.index, :WMA)
ema{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(ema(x.values, n; args...), x.index, :EMA)
dema{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(dema(x.values, n; args...), x.index, :DEMA)
tema{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(tema(x.values, n; args...), x.index, :TEMA)
mama{V,T}(x::TS{V,T}, fastlimit::Float64=0.5, slowlimit::Float64=0.05) = ts(mama(x.values, fastlimit, slowlimit), x.index, [:MAMA, :FAMA])

##### reg.jl ######
mlr_beta{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr_beta(y.values, n), y.index, :Beta)
mlr_slope{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr_slope(y.values, n), y.index, :Slope)
mlr_intercept{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr_intercept(y.values, n), y.index, :Intercept)
mlr{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr(y.values, n), y.index, :MLR)
mlr_se{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr_se(y.values, n), y.index, :StdErr)
mlr_ub{V,T}(y::TS{V,T}, n::Int64=10, se::Float64=2.0) = ts(mlr_ub(y.values, n, se), y.index, :UB)
mlr_lb{V,T}(y::TS{V,T}, n::Int64=10, se::Float64=2.0) = ts(mlr_lb(y.values, n, se), y.index, :LB)
mlr_bands{V,T}(y::TS{V,T}, n::Int64=10, se::Float64=2.0) = ts(mlr_bands(y.values, n, se), y.index, [:LB, :MLR, :UB])
mlr_rsq{V,T}(y::TS{V,T}, n::Int64=10) = ts(mlr_rsq(y.values, n), y.index, :RSquared)

##### mom.jl ######
momentum{V,T}(x::TS{V,T}, n::Int64=1) = ts(momentum(x.values, n), x.index, :Momentum)
roc{V,T}(x::TS{V,T}, n::Int64=1) = ts(roc(x.values, n), x.index, :ROC)
macd{V,T}(x::TS{V,T}, nfast::Int64=12, nslow::Int64=26, nsig::Int64=9; args...) = ts(macd(x.values, nfast, nslow, nsig; args...), x.index, [:MACD, :Signal, :Histogram])
rsi{V,T}(x::TS{V,T}, n::Int64=10; args...) = ts(rsi(x.values, n; args...), x.index, :RSI)
adx{V,T}(hlc::TS{V,T}, n::Int64=10; args...) = ts(adx(hlc.values, n; args...), hlc.index, [:DiPlus, :DiMinus, :ADX])
psar{V,T}(x::TS{V,T}, af::Float64=0.02, af_max::Float64=0.2, af_min::Float64=0.02) = ts(psar(x.values, af, af_max, af_min), x.index, :SMA)
kst{V,T}(x::TS{V,T}, nroc::Vector{Int64}=[10,15,20,30], navg::Vector{Int64}=[10,10,10,15]; args...) = ts(kst(x.values, nroc, navg; args...), x.index, :KST)
wpr{V,T}(hlc::TS{V,T}, n::Int64=14) = ts(wpr(hlc.values, n), hlc.index, :WPR)
cci{V,T}(hlc::TS{V,T}, n::Int64=20, c::Float64=0.015; args...) = ts(cci(hlc.values, n, c; args...), hlc.index, :WPR)
stoch{V,T}(hlc::TS{V,T}, nK::Int64=14, nD::Int64=3; args...) = ts(stoch(hlc.values, nK, nD; args...), hlc.index, [:Stochastic, :Signal])
smi{V,T}(hlc::TS{V,T}, n::Int64=13, nFast::Int64=2, nSlow::Int64=25, nSig::Int64=9; args...) = ts(smi(hlc.values, n, nFast, nSlow, nSig; args...), hlc.index, [:SMI, :Signal])

##### vol.jl ######
bbands{V,T}(x::TS{V,T}, n::Int64=10, sigma::Float64=2.0; args...) = ts(bbands(x.values, n, sigma; args...), x.index, [:LB, :MA, :UB])
tr{V,T}(x::TS{V,T}, n::Int64=2) = ts(tr(x.values, n), x.index, :TR)
atr{V,T}(x::TS{V,T}, n::Int64=14) = ts(atr(x.values, n), x.index, :ATR)
keltner{V,T}(hlc::TS{V,T}, nema::Int64=20, natr::Int64=10, mult::Int64=2) = ts(keltner(hlc.values, nema, natr, mult), hlc.index, [:LB, :MA, :UB])

##### trendy.jl #####
maxima{V,T}(x::TS{V,T}; args...) = ts(maxima(x.values; args...), x.index, :Maxima)
minima{V,T}(x::TS{V,T}; args...) = ts(minima(x.values; args...), x.index, :Minima)
support{V,T}(x::TS{V,T}; args...) = ts(support(x.values; args...), x.index, :Support)
resistance{V,T}(x::TS{V,T}; args...) = ts(resistance(x.values; args...), x.index, :Resistance)


text:
struct Cursor
    mx::Int
    my::Int

    max::Int
    may::Int

    x::Int
    y::Int

    ax::Int
    ay::Int

    Cursor(mx::Int, my::Int, max::Int, may::Int, x::Int, y::Int, ax::Int, ay::Int) = new(mx, my, max, may, x, y, ax, ay)

    function Cursor(event::eventMouse, camera::Camera, room::Maple.Room)
        mx, my = getMapCoordinates(camera, event.x, event.y)
        max, may = getMapCoordinatesAbs(camera, event.x, event.y)

        x, y = mapToRoomCoordinates(mx, my, room)
        ax, ay = mapToRoomCoordinatesAbs(max, may, room)

        return new(mx, my, max, may, x, y, ax, ay)
    end
end

ext: jl


text:
export imagepaths

"""
    imagepaths(path::AbstractString)
    imagepaths(file::File)
    imagepaths(dir::Dir)

List the image assets found in given markdown `file` or all markdown files in
directory `dir`.
"""
imagepaths(path::AbstractString, dest=homedir()) = imagepaths(Path(path), dest)

function imagepaths(file::File, dest::AbstractString)
    lines = readlines(file.path)
    for line in lines
        if  startswith(line, "![")
            m = match(r"!\[(.*)\]\((.+)\)", line)
            caption, imgpath = m.captures
            println("cp ../", imgpath, " ", dest)
        end
    end
end

imagepaths(dir::Dir, dest::AbstractString) = visitdir(imagepaths, dir.path, dest)

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
"""

Computing GMD contributions for the parts of the solute

"""

using M3GTools
using DelimitedFiles

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Data for my simulation (we only need the topoloy file here, might be PSF or GRO
topology = "./gmdfiles/structure.psf"

# Read GMD file corresponding to solute contributions
println(" Reading GMD solute contribution file... ")
gmd_solute = readdlm("./gmdfiles/gmd-GMD_ATOM_SOLUTE_CONTRIB.dat",comments=true,comment_char='#')

# The distance is the first column of the gmd file, and the total gmd is the second column:
d = gmd_solute[:,1]
gmd_total = gmd_solute[:,2]

# Define the solute selection used for computation of the GMD here:
gmd_backbone = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and (not sidechain)")
gmd_aliphatic = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and aliphatic")
gmd_aromatic = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and aromatic")
gmd_polar = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and polar")

# Plot the results
plot(xlim=[0,8],xlabel="Distance / Angstrom",ylabel="MDDF")
plot!(d,gmd_total,label="Total",linewidth=2)
plot!(d,gmd_backbone,label="Backbone",linewidth=2)
plot!(d,gmd_aliphatic,label="Aliphatic",linewidth=2)
plot!(d,gmd_aromatic,label="Aromatic",linewidth=2)
plot!(d,gmd_polar,label="Polar",linewidth=2)
savefig("./gmd_solute.pdf")

# Of course, you can save the data for further analysis:

output = open("./gmd_solute_contributions.dat","w")
write(output,"# Distance Total Backbone Aliphatic Aromatic Polar\n")
writedlm(output,zip(d,gmd_total,gmd_backbone,gmd_aliphatic,gmd_aromatic,gmd_polar))
close(output)
println(" Wrote file: ./gmd_solute_contribution.dat")


text:
if R123_USE_AESNI

import Random: seed!
using Test: @test

using RandomNumbers
using Random123

import RandomNumbers: split_uint
import Random123: __m128i

x = zero(__m128i)
ctr = __m128i(0x9799b5d54f7b9227, 0xb47607190d0dfefb)
key = __m128i(0x07b8e4b6aa98ec24, 0x5a7da274d3b8146a)
@test rand(ARS1x{1}(x, ctr, key), UInt128) ≡ 0x1a0b14c707b64224e548ef12331396ef
@test rand(ARS1x{2}(x, ctr, key), UInt128) ≡ 0x3ced8e0970690f718336318ba22e8ae1
@test rand(ARS1x{3}(x, ctr, key), UInt128) ≡ UInt128(__m128i(0xb6621a8b006319e8, 0x67c841642c32fc19))
@test rand(ARS1x{10}(x, ctr, key), UInt128) ≡ UInt128(__m128i(0xac35df44f996ed82, 0x4e287697bad2f9a2))

key = rand(UInt128)
r = ARS1x(key)
r1 = ARS4x(split_uint(key, UInt32))
@test seed_type(r) ≡ UInt128
@test seed_type(r1) ≡ NTuple{4, UInt32}
@test copyto!(copy(r), r) == r
@test copyto!(copy(r1), r1) == r1
@test UInt128(r.x) ≡ rand(r1, UInt128)
@test rand(r, UInt128) ≡ rand(r1, UInt128)
set_counter!(r, 0)
set_counter!(r1, 1)
@test rand(r, Tuple{UInt128})[1] ≡ rand(r1, UInt128)

end


text:
function run_weights(returns::Returns, weights::TimeArray{T,1}, transaction_costs::Float64 = 0.005)::Returns  where {T}
    last_weights = nothing
    function mapping(timestamp, value)
        if isnothing(weights[timestamp]) && isnothing(last_weights)
            return (timestamp, NaN)
        elseif isnothing(weights[timestamp])
            return (timestamp, value * last_weights)
        else
            new_weights = values(weights[timestamp])[1]
            cost = abs(new_weights - (isnothing(last_weights) ? 0.0 : last_weights)) * transaction_costs
            last_weights = values(weights[timestamp])[1]
            return (timestamp, value * new_weights - cost)
        end
    end
    return map(mapping, returns) |> dropnan
end

function run_weights(returns::TimeArray, weights::TimeArray{T,2}, transaction_costs::Float64 = 0.005) where {T}
    individual_returns = map(colnames(weights)) do col_name
        return run_weights(returns[col_name], weights[col_name])
    end
    return TimeSeries.merge(individual_returns..., method = :outer) |> sum_rows
end


text:
import PyPlot.plot

################
# Plotting tools
################

# plot one or more time series
function plot(x::TimeSeries; kws...)
  time = get_time(x)
  data = get_data(x)

  n_repeat = vcat(1, collect(size(data))[2:end])
  time_array = repeat(time, outer = n_repeat)

  plot(time_array, data; kws...)
end

function plot_mean(x::TimeSeries; kws...)
  plot(mean(x, dims=2); kws...)
end

# Plot the standard error for a time series as a filled polygon.
function plot_standard_error(x::TimeSeries{<:Any,2}; color="k", alpha=0.2,
                             linewidth=0.0, kws...)
  time = get_time(x)
  data = get_data(x)

  # reverse only handles 1D arrays
  stdes = standard_error(data, 2)[:]
  means = mean(data, dims=2)[:]

  xs = vcat(time, reverse(time))
  ys = vcat(means+stdes, reverse(means-stdes))

  fill(xs, ys; color=color, alpha=alpha, linewidth=linewidth, kws...)
end

# Often, it is to have a single compact plot method for a set of trials
function plot_trials(x::TimeSeries; mean=false, standard_error=false,
                     trials=false, label=Nothing, kws...)
  ret = Dict{Symbol, Any}()

  if trials
    ret[:trials] = plot(x; linestyle=":", kws...)
  end

  if mean
    ret[:mean] = plot_mean(x; label=label, kws...)
  end

  if standard_error
    ret[:standard_error] = plot_standard_error(x; kws...)
  end

  ret
end

###############
# Color schemes
###############

standard_colors = [RGB(0,0,1), RGB(1,0,0), RGB(0,1,0)]

rgb_tuple(color::RGB) = (red(color), green(color), blue(color))

generate_colors(x::AbstractArray) = generate_colors(length(x))
generate_colors(x::Integer) =
  map(rgb_tuple, distinguishable_colors(x, standard_colors))

#############
# File output
#############

# Changes a title to an appropriate file name
to_file_name(s) = replace(lowercase(s), " ", "_")

function save_latex(fig_path, title; kwargs...)
  mkpath(fig_path)

  matplotlib2tikz.save("$(fig_path)/$(title).tex",
                       axis_height="\\figureheight",
                       axis_width="\\figurewidth"; kwargs...)
end

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
# hand-crafted date selector

defaultCalendarHelpText = """
Arrows : move cursor
.      : jump to today
a      : jump to start of cursor's month
e      : jump to end of cursor's month
A      : jump to Jan 1st
E      : jump to Dec 31st
d,D    : add/subtract a day
w,W    : add/subtract a week
( the following would do end-of-month truncation)
m,M    : add/subtract a month
q,Q    : add/subtract a quarter
y,Y    : add/subtract a year
"""
type TwCalendarData
    showHelp::Bool
    helpText::UTF8String
    date::Date
    cursorweekofmonth::Int # cached "nth week" in the current month containing date, 1-based
    geometry::@compat Tuple{Int, Int} # rows x cols in months
    ncalStyle::Bool
    TwCalendarData( dt::Date ) = new( true, defaultCalendarHelpText, dt, 1, (1,1), false )
end

function monthDimension( ncalStyle::Bool )
    ncalStyle ? ( 8, 3*6-1 ): (8, 3*7+1 )
end

function bestfitgeometry( ncalStyle, scr::TwScreen, box::Bool )
    ( parmaxy, parmaxx ) = getwinmaxyx( scr.window )
    leftcols = 0
    monthdim::Tuple{Int,Int} = monthDimension( ncalStyle )
    if ncalStyle
        allowed_geometry = [ (3,4), (2, 3), (1,3), (1,1 ) ]
        leftcols = 2
    else
        allowed_geometry = [ (4,3), (2, 3), (1,3), (1,1 ) ]
    end

    found = false
    finalg = (0,0,0,0)
    for g in allowed_geometry
        h::Int = 1+g[1] * monthdim[1] + (box?2:0) # box borders + contents + year title
        w::Int = leftcols + g[2] * monthdim[2] + (box?2:0)
        if w <= parmaxx && h <= parmaxy
            found = true
            finalg = (h,w, g[1]::Int, g[2]::Int )
            break
        end
    end

    if !found
        throw( "terminal is too small to view even one month")
    end
    return finalg
end

function newTwCalendar( scr::TwScreen, dt::Date;
    posy::Any = :center,posx::Any = :center,
    ncalStyle=true, box=true, showHelp=true, title = "" )
    data = TwCalendarData( dt )
    obj = TwObj( data, Val{:Calendar} )
    registerTwObj( scr, obj )
    obj.box = box
    obj.title = title
    obj.borderSizeV= box ? 1 : 0
    obj.borderSizeH= box ? 1 : 0
    obj.data.showHelp = showHelp
    obj.data.date = dt
    h,w,g1,g2 = bestfitgeometry( ncalStyle, scr, box )
    obj.data.geometry = (g1, g2)
    obj.data.ncalStyle = ncalStyle
    alignxy!( obj, h, w, posx, posy)
    configure_newwinpanel!( obj )
    obj
end

function draw( o::TwObj{TwCalendarData} )
    werase( o.window )
    if o.box
        box( o.window, 0,0 )
    end
    if !isempty( o.title ) && o.box
        mvwprintw( o.window, 0, (@compat round( Int, ( o.width - length(o.title) )/2 )), "%s", o.title )
    end
    starty = o.borderSizeV
    startx = o.borderSizeH
    mvwprintw( o.window, starty, (@compat round( Int, ( o.width - 4 )/2 )), "%s", string(year(o.data.date)))
    starty += 1
    # figure out the start month
    nummonths = o.data.geometry[1] * o.data.geometry[2]
    m = month( o.data.date )
    if nummonths >= 12
        startmonth = 1 # jan
    elseif nummonths >= 6
        startmonth = m >= 7 ? 7 : 1
    elseif nummonths >= 3
        startmonth = m - mod1( m, 3 ) + 1
    else
        startmonth = m
    end

    mth = startmonth
    ncalStyle = o.data.ncalStyle
    monthdim = monthDimension( ncalStyle )
    wkdys = [ "Mo", "Tu", "We", "Th", "Fr", "Sa", "Su" ]
    for i = 1:o.data.geometry[1]
        if ncalStyle # draw the week column on the left
            for (wdidx,wdstr) in enumerate( wkdys )
                mvwprintw( o.window, starty + wdidx, startx, "%s", wkdys[wdidx] )
            end
            startx += 2
        end
        for j = 1:o.data.geometry[2]
            # print the month header
            mvwprintw( o.window, starty, startx + 6, "%s", monthabbr( mth ) )
            mst = Date( year( o.data.date ), mth, 1 )
            men = Date( year( o.data.date ), mth, daysinmonth( year( o.data.date ), mth ) )
            dt = mst
            wkd = dayofweek( dt )
            if ncalStyle
                wcol = 0
                while dt <= men
                    flags = 0
                    if dt == o.data.date
                        flags = COLOR_PAIR( o.hasFocus ? 15 : 30 )
                        o.data.cursorweekofmonth = wcol + 1
                    end
                    if dt == today()
                        flags = flags | A_BOLD
                    end
                    wattron( o.window, flags )
                    if wcol == 0
                        mvwprintw( o.window, starty + wkd, startx, "%2s", string( day( dt ) ) )
                    else
                        mvwprintw( o.window, starty + wkd, startx - 1 + wcol * 3, "%3s", string( day( dt ) ) )
                    end
                    wattroff( o.window, flags )
                    if wkd == 7
                        wkd = 1
                        wcol += 1
                    else
                        wkd += 1
                    end
                    dt = dt + Day(1)
                end
            else
                wrow = 0
                for (wdidx,wdstr) in enumerate( wkdys )
                    mvwprintw( o.window, starty + 1, startx + (wdidx-1)* 3, "%3s", wkdys[wdidx] )
                end
                while dt <= men
                    flags = 0
                    if dt == o.data.date
                        flags = COLOR_PAIR( o.hasFocus ? 15 :30 )
                        o.data.cursorweekofmonth = wrow + 1
                    end
                    if dt == today()
                        flags = flags | A_BOLD
                    end
                    wattron( o.window, flags )
                    mvwprintw( o.window, starty + 2 + wrow, startx + (wkd-1) * 3, "%3s", string( day( dt ) ) )
                    wattroff( o.window, flags )
                    if wkd == 7
                        wkd = 1
                        wrow += 1
                    else
                        wkd += 1
                    end
                    dt = dt + Day(1)
                end
            end
            startx += monthdim[2]
            mth += 1
        end
        starty += monthdim[1]
        startx = o.borderSizeH # reset
    end
end

# returns a range of the dates of that week. n is 1-based
# if n is too large, it'd use the last week of that month
function monthNthWeekRange( y::Int, m::Int, n::Int )
    monthstart = Date( y, m )
    wdaystart = dayofweek( monthstart )
    mdays = daysinmonth( y, m )
    we = 8 - wdaystart + 7 * (n-1)
    ws = we - 6
    while ws > mdays
        n = n-1
        we = 8 - wdaystart + 7 * (n-1)
        ws = we - 6
    end
    return (Date( y, m, max(ws, 1) ),Date( y, m, min( we, mdays )))
end

function inject( o::TwObj{TwCalendarData}, token )
    dorefresh = false
    retcode = :got_it # default behavior is that we know what to do with it

    ncalStyle = o.data.ncalStyle
    if token == :esc
        retcode = :exit_nothing
    elseif token == "."
        o.data.date = today()
        dorefresh = true
    elseif in( token, [ :up, :down, :left, :right ] )
        if ncalStyle
            if token == :left
                o.data.date = o.data.date - Day(7)
            elseif token == :right
                o.data.date = o.data.date + Day(7)
            elseif token == :up
                if dayofweek( o.data.date ) > 1 && day( o.data.date ) > 1
                    o.data.date = o.data.date - Day(1)
                else
                    # last weekday of the same nth week of the previous month
                    # if the nth week of the previous month doesn't exist, it'd
                    # be the weekday of the last week of the previous month
                    prevmonth = o.data.date - Month( ncalStyle ? o.data.geometry[2] : o.data.geometry[1] )
                    o.data.date = monthNthWeekRange( year( prevmonth ), month( prevmonth ), o.data.cursorweekofmonth )[2]
                end
            else # :down
                if dayofweek( o.data.date ) < 7 && day( o.data.date ) < daysinmonth( year( o.data.date ), month( o.data.date ) )
                    o.data.date = o.data.date + Day(1)
                else
                    # last weekday of the same nth week of the "next" month
                    nextmonth = o.data.date + Month( ncalStyle ? o.data.geometry[2] : o.data.geometry[1] )
                    o.data.date = monthNthWeekRange( year( nextmonth ), month( nextmonth ), o.data.cursorweekofmonth )[1]
                end
            end
            dorefresh = true
        else
            if token == :up
                d = day( o.data.date )
                if d >= 7
                    o.data.date = o.data.date - Day(7)
                else
                    currdayofweek = dayofweek( o.data.date )
                    prevmonth = o.data.date - Month( o.data.geometry[2] )
                    (y,m) = (year(prevmonth), month(prevmonth))
                    mds = daysinmonth( y,m )
                    lstdayofwk = dayofweek( Date( y,m,mds ) )
                    o.data.date = Date( y,m, mds - mod(lstdayofwk-currdayofweek, 7))
                end
            elseif token == :down
                (y,m,d) = (year(o.data.date ), month( o.data.date), day( o.data.date ))
                mds = daysinmonth( y,m )
                if d + 7 <= mds
                    o.data.date = o.data.date + Day(7)
                else
                    currdayofweek = dayofweek( o.data.date )
                    nextmonth = o.data.date + Month( o.data.geometry[2] )
                    (y,m) = (year(nextmonth), month(nextmonth))
                    mds = daysinmonth( y,m )
                    fstdayofwk = dayofweek( Date( y,m,1 ) )
                    o.data.date = Date( y,m, 1 + mod(currdayofweek-fstdayofwk, 7))
                end
            elseif token == :left
                if dayofweek( o.data.date ) > 1 && day( o.data.date ) > 1
                    o.data.date = o.data.date - Day(1)
                else
                    # last weekday of the same nth week of the previous month
                    # if the nth week of the previous month doesn't exist, it'd
                    # be the weekday of the last week of the previous month
                    prevmonth = o.data.date - Month( 1 )
                    o.data.date = monthNthWeekRange( year( prevmonth ), month( prevmonth ), o.data.cursorweekofmonth )[2]
                end
            else # :right
                if dayofweek( o.data.date ) < 7 && day( o.data.date ) < daysinmonth( year( o.data.date ), month( o.data.date ) )
                    o.data.date = o.data.date + Day(1)
                else
                    # last weekday of the same nth week of the "next" month
                    nextmonth = o.data.date + Month( 1 )
                    o.data.date = monthNthWeekRange( year( nextmonth ), month( nextmonth ), o.data.cursorweekofmonth )[1]
                end
            end
            dorefresh = true
        end
    elseif token == "a"
        o.data.date = Date( year( o.data.date ), month( o.data.date ) )
        dorefresh = true
    elseif token == "e"
        (y,m) = year(o.data.date), month( o.data.date )
        o.data.date = Date( y,m, daysinmonth( y,m ) )
        dorefresh = true
    elseif token == "A"
        o.data.date = Date( year( o.data.date ), 1 )
        dorefresh = true
    elseif token == "E"
        y = year(o.data.date)
        o.data.date = Date( y, 12, 31 )
        dorefresh = true
    elseif token == "w"
        o.data.date = o.data.date + Day(7)
        dorefresh = true
    elseif token == "W"
        o.data.date = o.data.date - Day(7)
        dorefresh = true
    elseif token == "m"
        o.data.date = o.data.date + Month(1)
        dorefresh = true
    elseif token == "M"
        o.data.date = o.data.date - Month(1)
        dorefresh = true
    elseif token == "d"
        o.data.date = o.data.date + Day(1)
        dorefresh = true
    elseif token == "D"
        o.data.date = o.data.date - Day(1)
        dorefresh = true
    elseif token == "y"
        o.data.date = o.data.date + Year(1)
        dorefresh = true
    elseif token == "Y"
        o.data.date = o.data.date - Year(1)
        dorefresh = true
    elseif token == "q"
        o.data.date = o.data.date + Month(3)
        dorefresh = true
    elseif token == "Q"
        o.data.date = o.data.date - Month(3)
        dorefresh = true
    elseif token == :enter || token == Symbol( "return" )
        o.value = o.data.date
        retcode = :exit_ok
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        refresh(o)
    end

    return retcode
end

function helptext( o::TwObj{TwCalendarData} )
    if o.data.showHelp
        o.data.helpText
    else
        ""
    end
end


text:
"""
    Connectivity{T}

Stores connection data between cells. Connections are stored in a compressed
sparse column (CSC) adjacency matrix: only non-zero values are stored.
Primarily, this consist of two vectors:
* the row value vector holds the cell indices of neighbors
* the column pointers marks the start and end index of the row value vector

This matrix is square: n = m = ncell. nconnection is equal to nnz (the number of
non-zero values).

* ncell: the number of (active) cells in the simulation
* nconnection: the number of connections between cells
* length1: for every connection, the length in the first cell, size nconnection
* length2: for every connection, the length in the second cell, size nconnection
* width: width for every connection, (approx.) perpendicular to length1 and
  length2, size nconnection
* colptr: CSC column pointer (size ncell + 1)
* rowval: CSC row value (size nconnection)
"""
struct Connectivity{T}
    ncell::Int
    nconnection::Int
    length1::Vector{T}
    length2::Vector{T}
    width::Vector{T}
    colptr::Vector{Int}
    rowval::Vector{Int}
end

"""
    connections(C::Connectivity, id::Int)

Returns connections for a single cell, identified by ``id``.
"""
connections(C::Connectivity, id::Int) = C.colptr[id]:(C.colptr[id+1]-1)

"""
    connection_geometry(I, J, Δx, Δy)

Compute geometrical properties of connections for structured input.
"""
function connection_geometry(I, J, Δx, Δy)
    if I[1] != J[1]  # connection in y
        length1 = 0.5 * Δy[I[1]]
        length2 = 0.5 * Δy[J[1]]
        width = Δx[I[2]]
    elseif I[2] != J[2]  # connection in x
        length1 = 0.5 * Δx[I[2]]
        length2 = 0.5 * Δx[J[2]]
        width = Δy[I[1]]
    else
        # TODO: more specific exception? --> Martijn
        error("Inconsistent index")
    end
    return (length1, length2, width)
end

# Define cartesian indices for neighbors
const neighbors = (
    CartesianIndex(0, -1),
    CartesianIndex(-1, 0),
    CartesianIndex(1, 0),
    CartesianIndex(0, 1),
)

# Constructor for the Connectivity structure for structured input
function Connectivity(indices, reverse_indices, Δx::Vector{T}, Δy::Vector{T}) where {T}
    # indices: These map from the 1D internal domain to the 2D external domain.
    # reverse_indices: from the 2D external domain to the 1D internal domain,
    # providing an Int which can be used as a linear index
    nrow, ncol = size(reverse_indices)
    # Pre-allocate output, allocate for full potential number of neighbors (4)
    ncell = length(indices)
    colptr = Vector{Int}(undef, ncell + 1)
    rowval = Vector{Int}(undef, ncell * 4)
    length1 = similar(rowval, T)
    length2 = similar(rowval, T)
    width = similar(rowval, T)

    i = 1  # column index of sparse matrix
    j = 1  # row index of sparse matrix
    for I in indices # loop over active indices
        colptr[j] = i
        # Strictly increasing numbering for any row
        # (Required by a CSCSparseMatrix, if you want to convert)
        for neighbor in neighbors
            J = I + neighbor
            if (1 <= J[1] <= nrow) && (1 <= J[2] <= ncol && reverse_indices[J] != 0) # Check if it's inbounds and neighbor is active
                rowval[i] = reverse_indices[J]
                length1[i], length2[i], width[i] = connection_geometry(I, J, Δx, Δy)
                i += 1
            end
        end
        j += 1
    end
    colptr[j] = i

    nconnection = i - 1
    return Connectivity(
        ncell,
        nconnection,
        length1[1:nconnection],
        length2[1:nconnection],
        width[1:nconnection],
        colptr,
        rowval[1:nconnection],
    )
end


text:
function merge_sort!(
    arr::Vector{T},
    l::Int = 1,
    r::Int = length(arr),
    temp::Vector{T} = Vector{T}(undef, r - l + 1),
) where {T}
    if l >= r
        return
    end
    # split
    mid = (l + r) >> 1
    merge_sort!(arr, l, mid)
    merge_sort!(arr, mid + 1, r)
    # merge
    l_pos = l # pos of the left part
    r_pos = mid + 1 # pos of the right part
    for t_pos in 1:r-l+1
        if l_pos <= mid && (r_pos > r || arr[l_pos] < arr[r_pos])
            temp[t_pos] = arr[l_pos]
            l_pos += 1
        else
            temp[t_pos] = arr[r_pos]
            r_pos += 1
        end
    end
    for i in l:r
        arr[i] = temp[i-l+1]
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
#=
    Goal:
        For the given Hilbert matrix H, and a collection of
        100 random vectors bᵢ,
        compute H⁻¹ bᵢ using
        a) the LU factorization
        b) the automatic factorization
        c) use a SVD factorization

    Hint:
=#

using LinearAlgebra

# construction of a Hilbert matrix
N = 25
H = Float64[1/(i+j-1) for i in 1:N, j in 1:N]

M = 1000
xs = [rand(N) for _ in 1:M]

H_lu = lu(H)
H_fact = factorize(H)

bs = [H*x for x in xs]

ratios = zeros(M)
for i in 1:M
    res_lu   = norm(xs[i] - H_lu\bs[i])
    res_fact = norm(xs[i] - H_fact\bs[i])
    ratios[i] = res_lu/res_fact
end

# optional: visuation of the results
using CairoMakie
begin
    p = hist(ratios, bins=100, axis = (title = "Ratio",))
    vlines!(p.axis, [1.0], color = :red)
    current_figure()
end

norm(xs[1] - lu(H) \ bs[1])
norm(xs[1] - svd(H) \ bs[1])

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
import REPL
using REPL.TerminalMenus

# Make sure JSON is installed
try
    using JSON
catch ex
    # Ask to install JSON
    choice = request(
        "Package JSON not installed - would you like me to install it?",
        RadioMenu(["yes", "no"])
    )
    if choice == 1
        using Pkg
        Pkg.add("JSON")
        using JSON
    else
        println("Canceling!")
        exit()
    end
end

# Constants for paths
AUTOTM_BUILD = joinpath(@__DIR__, "AutoTM", "deps", "build.json")
NGRAPH_BUILD = joinpath(@__DIR__, "deps", "nGraph", "deps", "build.json")

#####
##### Setup for nGraph
#####

options = ["yes", "no"]

# Select NVDIMMs
choice = request(
                 "Would you like to use NVDIMMs? (requires a Cascade Lake system with Optane)",
    RadioMenu(options)
)

# Choice for using NVDIMMs
use_nvdimms = choice == 1

# Select GPU
choice = request(
    "Would you like to build with GPU support?",
    RadioMenu(options)
)

# Choice for using NVDIMMs
use_gpu = choice == 1

ngraph_json = JSON.parsefile(NGRAPH_BUILD)

if use_nvdimms
    ngraph_json["PMDK"] = true
    ngraph_json["NUMA"] = true
else
    ngraph_json["PMDK"] = false
    ngraph_json["NUMA"] = false
end

if use_gpu
    ngraph_json["GPU"] = true
else
    ngraph_json["GPU"] = false
end

open(NGRAPH_BUILD; write = true) do f
    JSON.print(f, ngraph_json)
end

#####
##### Setup AutoTM
#####

# Select GPU
choice = request(
    "Would you like to use Gurobi as the ILP solver?",
    RadioMenu(options)
)

# Choice for using NVDIMMs
use_gurobi = choice == 1

autotm_json = JSON.parsefile(AUTOTM_BUILD)

if use_gurobi
    autotm_json["GUROBI"] = true
else
    autotm_json["GUROBI"] = false
end

open(AUTOTM_BUILD; write = true) do f
    JSON.print(f, autotm_json)
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using RemoteFiles
"""
http://cs.nyu.edu/~silberman/datasets/nyu_depth_v2.html
"""
struct NYUDataset <: RGBDDataset
    path::String
end
function NYUDataset()
    return NYUDataset(joinpath(dirname(pathof(RGBDSegmentation)), "..", "datasets", "NYU"))
end
import Base.length
import Base.iterate
function iterate(s::NYUDataset, state=1)
    if (state > length(s))
        return nothing
    else
        return (get_image(s, state), state + 1)
    end
end
function name(s::NYUDataset)
    return "NYU"
end
function length(s::NYUDataset)
    return 1449
end
function raw_path(ds::NYUDataset)
    return joinpath(ds.path,"raw")
end
function raw_file(ds::NYUDataset)
    return joinpath(raw_path(ds), "nyu_depth_v2_labeled.mat")
end
function transformed_path(ds::NYUDataset)
    return joinpath(ds.path, "transformed")
end
function transformed_file_name(ds::NYUDataset, j::Integer)
    return "image_nyu_v2_$(j).fits"
end
function transformed_file_path(ds::NYUDataset, j::Integer)
    return joinpath(
        transformed_path(ds),
        transformed_file_name(ds, j)
    )
end
function preprocess(ds::NYUDataset, n::Integer=0)
    c_params = camera_params(ds)
    file = matopen(raw_file(ds))
    RGB = read(file,"images")
    depths = read(file,"depths")
    labels = read(file,"labels")
    if n == 0
        n = size(labels,3)
    else
        n = min(n, size(labels, 3))
    end
    mkpath(transformed_path(ds))
    for j=1:n
        @info "Preprocessing image $j of $n"
        t_file_path = transformed_file_path(ds, j)
        if !isfile(t_file_path)
            f = FITS(t_file_path, "w");
            write(f, RGB[:,:,:,j])
            write(f, depth_plane2depth_world(depths[:,:,j],
                                             c_params["c_d"],
                                             c_params["f_d"]))
            write(f, labels[:,:,j])
            close(f)
        end
    end
end
function download(ds::NYUDataset)
    file_path = raw_file(ds)
    if !isfile(file_path)
        @info("Downloading file in $file_path")
        @RemoteFile(dataset_raw_file,
                    "http://horatio.cs.nyu.edu/mit/silberman/nyu_depth_v2/nyu_depth_v2_labeled.mat",
                    dir=raw_path(ds))
        download(dataset_raw_file)

    else
        @info("Raw file already present")
    end
    preprocess(ds)
end
"""
```julia
function get_image(ds::NYUDataset, i::Integer; compute_normals=true)
```
"""
function get_image(ds::NYUDataset, i::Integer; compute_normals=true)
    image_path = joinpath(transformed_path(ds), string("image_nyu_v2_",i,".fits"))
    f = FITS(image_path, "r+")
    N = nothing
    D = read(f[2])
    if (length(f)==3) || compute_normals
        (N1, N, N2) = compute_local_planes(D,window=3,relDepthThresh=0.01)
    else
#        N = read(f[4])
    end
    margin = 15
    RGB = read(f[1]);
    RGB = RGB[margin:end-margin, margin:end-margin,:] ;
    D = D[margin:end-margin, margin:end-margin, :];
    N = N[margin:end-margin, margin:end-margin,:];
    S = read(f[3])
    S = convert.(Integer, S[margin:end-margin, margin:end-margin] .+ 1);
    return LabelledCDNImage(CDNImage(RGB,D,N),S)
end

function camera_params(ds::NYUDataset)
    params = Dict(
        # The maximum depth used, in meters.
        "maxDepth" => 10,
        # RGB Intrinsic Parameters
        "f_rgb" => [5.1885790117450188e+02 5.1946961112127485e+02],
        "c_rgb" => [3.2558244941119034e+02 2.5373616633400465e+02],
        # RGB Distortion Parameters
        "k_rgb" => [ 2.0796615318809061e-01,
                     -5.8613825163911781e-01,
                     4.9856986684705107e-01],
        "p_rgb" => [7.2231363135888329e-04  1.0479627195765181e-03],
        # Depth Intrinsic Parameters
        "f_d" => [5.8262448167737955e+02 5.8269103270988637e+02],
        "c_d" => [3.1304475870804731e+02 2.3844389626620386e+02],

        "k_d" => [-9.9897236553084481e-02,
                  3.9065324602765344e-01,
                  -5.1031725053400578e-01],
        "p_d" => [1.9290592870229277e-03 -1.9422022475975055e-03],
            # 3D Translation
        "t" => [2.5031875059141302e-02,
                -2.9342312935846411e-04,
                6.6238747008330102e-04],
        # Parameters for making depth absolute.
        "depthParam1" => 351.3,
        "depthParam2" => 1092.5

    )
# Rotation
    R = -[ 9.9997798940829263e-01, 5.0518419386157446e-03,
           4.3011152014118693e-03, -5.0359919480810989e-03,
           9.9998051861143999e-01, -3.6879781309514218e-03,
           -4.3196624923060242e-03, 3.6662365748484798e-03,
           9.9998394948385538e-01 ];

    R = reshape(R, (3, 3));
    R = inv(R');
    params["R"]= R
    return params

end


text:
using Images
using DataFrames
using FileIO

## THERE ARE 250 different sketches, with 80 sketches each

dataSetFolder = "/raided/datasets/sketches_png/png/"
df = DataFrame(Any, 0, 1234322)

listOfImgs = open(dataSetFolder * "filelist.txt")
imgPaths = readlines(listOfImgs)
close(listOfImgs)

writetable("/raided/datasets/sketches_png/allSketches.csv", df)

output = open("/raided/datasets/sketches_png/allSketches.csv", "a")

count = 0
for path in imgPaths
    path = path[1:end-1]
    sketchName = split(path, "/")[1]

    println("Loading picture " * dataSetFolder * path * "... ")
    img = load(dataSetFolder * path)
    println("Done loading picture. Now extracting pixel values...\n")

    println("Converting image to array of strings...")
    img = map(x-> round(Int64, x.val), img)
    println("img is now an array of strings\n")

    println("Reshape img to 1d array...")
    flatten = reshape(img, 1234321)
    println("Done reshaping\n")

    println("Appending tranpose of image to data and convert to (1,21234322) array...")
    data = vcat([sketchName], flatten)
    data = reshape(data, (1, 1234322))
    println("Done appending to data\n")

    println("Writing out to csv...")
    writecsv(output, data)
    println("Done writing\n")
    count += 1
    if count == 80
        break
    end
end

close(output)

println("********************************")
println("DONE")
println("********************************")

# 61 839 458


text:
export getchar, set_homedir, get_homedir

const charset = 'a':'z'

getchar(idx) = charset[idx]

const homedir = ["wrong"]

function set_homedir(path)
    global homedir[1] = path
end
get_homedir() = homedir[1]


text:
# using Revise
using ClusterManagers
using Distributed

worker_ids = if "SLURM_JOBID" in keys(ENV)
  wids = ClusterManagers.addprocs_slurm(parse(Int, ENV["SLURM_NTASKS"]))
  wids
else
  addprocs(2)
end

Distributed.@everywhere begin
  using Pkg
  Pkg.activate(".")
  Pkg.instantiate()
end

Distributed.@everywhere begin
  using POMDPs
  using BasicPOMCP
  using POMCPOW
  using VOOTreeSearch
  using ParticleFilters
  using CPUTime
  using DataFrames
  using ARDESPOT
  using Distributions
  using VDPTag2
  using Dates
  using CSV
  using POMDPSimulators
  using VDPTag2Experiments
  using LinearAlgebra: diagm, PosDefException, det, eigvals, eigvecs
  using Random: MersenneTwister, seed!
  using Distances: euclidean

  using MCTS
  using StaticArrays
  using POMDPModelTools: GenerativeBeliefMDP
  include("./scripts/VDPTag2/vdp_next_action.jl")
end

@info "Started $(Distributed.nworkers()) workers..."
@info "Precompiling simulation code..."

max_time = 1.0
max_depth = 10

@show start_mean = [110.0, 30.0, 30.0, 5.0, 100.0, 0.8]
@show start_cov = diagm([100.0^2, 10.0^2, 10.0^2, 10.0^2, 20.0^2, 0.1^2])
d = MvNormal(start_mean, start_cov)
on_cluster = true
save_data = true
parallelize = true

if on_cluster
  max_query = 100_000
  K = 150 # number of parameter samples
  n = 20 # number of evaluation simulations
  m = 30 # number of elite samples
  max_iters = 50
else
  max_query = 1000
  K = 15 # number of parameter samples
  n = 1 # number of evaluation simulations
  m = 10 # number of elite samples
  max_iters = 3
end

# VDPTag2
@show pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 2.8)))

println(stderr, "Max time: ", max_time)
println(stderr, "Max depth: ", max_depth)
println(stderr, "Max query: ", max_query)
println(stderr, "Start mean: ", start_mean)
println(stderr, "Start cov: ", start_cov)
println(stderr, "POMDP: ", pomdp)

function gen_sims(x::Vector{Float64}, n, k, seed)
    c = x[1]
    @assert c >= 0.0
    k_act = x[2]
    @assert k_act >= 1.0
    inv_alpha_act = x[3]
    @assert inv_alpha_act >= 0.1
    k_obs = x[4]
    @assert k_obs >= 1.0
    inv_alpha_obs = x[5]
    @assert inv_alpha_obs >= 0.1
    p_voo = x[6]
    @assert p_voo >= 0.1
    @assert p_voo <= 1.0

    sims = []

    for i in 1:n
        rng = MersenneTwister(7)

        ro = ToNextMLSolver(rng)
        sig = diagm([0.1])
        epsilon = 0.1/10.0
        voo_opt = VOOActionGenerator(p_voo, Distributions.Uniform(), MvNormal(sig), rng, epsilon)

        solver = POMCPOWSolver(tree_queries=max_query,
                               criterion=MaxUCB(c),
                               max_depth=max_depth,
                               max_time=max_time,
                               k_action=k_act,
                               alpha_action=1/inv_alpha_act,
                               k_observation=k_obs,
                               alpha_observation=1/inv_alpha_obs,
                               estimate_value=FORollout(ro),
                               next_action = RootToNextMLFirstVOO(voo_opt),
                               default_action=BasicPOMCP.ReportWhenUsed(TagAction(false, 0.0)),
                               rng=rng
                              )

        planner = solve(solver, pomdp)
        filter = SIRParticleFilter(deepcopy(pomdp),
                                   100_000,
                                   rng=MersenneTwister(i+90_000))

        seed!(planner, i+40000*k)
        sim = Sim(deepcopy(pomdp),
                  planner,
                  filter,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=50,
                  metadata=Dict(:i=>i, :k=>k)
                 )
        push!(sims, sim)
    end

    return sims
end

datestring = Dates.format(now(), "e_d_u_Y_HH_MM")
results_csv = DataFrame(c=Float64[], k_act=Float64[], inv_alpha_act=Float64[], k_obs=Float64[], inv_alpha_obs=Float64[], p=Float64[])
try
  for i in 1:max_iters
      sims = []
      params = Vector{Float64}[]
      println(stderr,"creating $K simulation sets")
      for k in 1:K
          p = rand(d)
          p[1] = max(0.0, p[1])
          p[2] = max(1.0, p[2])
          p[3] = max(0.1, p[3])
          p[4] = max(1.0, p[4])
          p[5] = max(0.1, p[5])
          p[6] = max(0.1, min(1.0, p[6]))
          push!(params, p)
          k_sims = gen_sims(p, n, k, i)
          append!(sims, k_sims)
      end
      @assert length(params) == K
      if parallelize
        results = run_parallel(sims)
      else
        results = run(sims)
      end
      combined = by(results, :k) do df
          DataFrame(mean_reward=mean(df[:reward]))
      end
      order = sortperm(combined[:mean_reward])
      elite = params[combined[:k][order[K-m:end]]]
      elite_matrix = Matrix{Float64}(undef, length(start_mean), m)
      for k in 1:m
          elite_matrix[:,k] = elite[k]
      end
      try
          global d = fit(typeof(d), elite_matrix)
      catch ex
          if ex isa PosDefException
              println(stderr,"pos def exception")
              global d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
          else
              rethrow(ex)
          end
      end
      println(stderr, "Iteration $i")
      println(stderr, "Reward (mean): ", mean(combined[:mean_reward]))
      println(stderr, "Mean: ", mean(d))
      println(stderr, "Cov (det): ", det(cov(d)))
      ev = eigvals(cov(d))
      println(stderr, "Cov (eig): ", ev)
      for j in 1:length(ev)
          println(stderr, "Eigvecs: ", eigvecs(cov(d))[:,j])
      end
      push!(results_csv, mean(d))
      if save_data
        fname_csv = joinpath(@__DIR__, "../data", "VDP_vomcpow_opt"*"_"*datestring*".csv")
        CSV.write(fname_csv, results_csv)
        println(stderr,"Saving parameters to ", fname_csv)
      end
  end
finally
  Distributed.rmprocs(worker_ids)
end


text:
module MITprof_plots

using ArgoData, DataFrames, CSV, Dates, GLMakie, JLD2, Statistics

#trim_cost is used in plot_cost
function trim_cost(cost)
    ii=findall((!ismissing).(cost) .&& (isfinite).(cost) .&& (cost.>0.0))
    med=round(median(cost[ii]), digits=3)
    nb=length(ii)

    (ncost,nii)=(length(cost),length(ii))
    println("keeping $(nii), leaving $(ncost-nii)")
    println("median $(median(cost[ii]))")
    println("mean $(mean(cost[ii]))")
    println("max $(maximum(cost[ii]))")

    cost[ii],med,nb
end

"""
    plot_cost()

```
f=MITprof_plots.plot_cost()
save("cost_pdf.png", f)
```
"""
function plot_cost()
    #
    costT=load("csv/prof_T_stats.jld2","cost")
    costT,medT,nbT=trim_cost(costT)
    #
    costS=load("csv/prof_S_stats.jld2","cost")
    costS,medS,nbS=trim_cost(costS)

    f = Figure()
    hist(f[1, 1], costT, bins = collect(0:0.1:8.0),
        axis = (title = "cost for T (median=$(medT), from $(nbT))",),
        color = :values, strokewidth = 1, strokecolor = :black)

    hist(f[2, 1], costS, bins = collect(0:0.1:8.0),
        axis = (title = "cost for S (median=$(medS), from $(nbS))",),
        color = :values, strokewidth = 1, strokecolor = :black)

    f
end

"""
    plot_stats()

Read file `profile_positions.csv`, pre-process, and then compute
and display basic statistics of the Argo float array.

```
f=MITprof_plots.plot_stats()
save("ArgoDistributions.png", f)
```
"""
function plot_stats()
    df=read_profile_positions()

    #display basic statistics of the Argo float array:
    plot_stats(df)
end

"""
    read_profile_positions(csv_file="csv/profile_positions.csv")

Read "profile_positions.csv" file.

```
df=MITprof_plots.read_profile_positions("csv/profile_positions.csv")
```
"""
function read_profile_positions(csv_file="csv/profile_positions.csv")
df=CSV.read(csv_file,DataFrame)

#time variable = year + (month-1/2)/12 :
df.ym=year.(df.date)+(month.(df.date) .-0.5)/12

#eliminate problematic data points:
#l0=size(df,1)
df=df[df.lat .> -89.99,:]
df=df[df.date .> DateTime(1000,1,1),:]
#l1=size(df,1)
df=df[df.date .< DateTime(2022,4,1),:]
#l2=size(df,1)

df
end

"""
    plot_stats(df::DataFrame; ym1=2004+(12-0.5)/12,ym2=2021+(12-0.5)/12)

Call `plot_stats` with `df` already in memory. Please refer to the `plot_stats()`
source code for documentation on how to load `df`.
"""
function plot_stats(df::DataFrame; ym1=2004+(12-0.5)/12,ym2=2021+(12-0.5)/12)
    txt(ym) = string(Int(floor(ym)))*"/"*string(Int(round((ym-floor(ym))*12+0.5)))

    #ii=findall(peryear[:,:year].>1900)

    gdf=groupby(df,:ym)
    df2=combine(gdf) do df
        (m = length(df.ID), n = length(unique(df.ID)))
    end
    sort!(df2, [:ym])

    ##

    f = Figure()

    ax1 = Axis(f[1,1], title="number of profiles per day")
    lines!(ax1,df2[:,:ym],df2[:,:m]/30)
    lines!(ax1,[ym1,ym1],[0,600],color=:red,linestyle=:dash)
    lines!(ax1,[ym2,ym2],[0,600],color=:red,linestyle=:dash)
    ylims!(ax1, (0, 600))

    ax2 = Axis(f[1,2], title="number of profilers at sea")
    barplot!(ax2,df2.ym,df2.n)
    lines!(ax2,[ym1,ym1],[0,5000],color=:red,linestyle=:dash)
    lines!(ax2,[ym2,ym2],[0,5000],color=:red,linestyle=:dash)
    ylims!(ax2, (0, 5000))

    ax3 = Axis(f[2,1], title="profile positions for "*txt(ym1))
    tmp=df[df.ym.==ym1,:]
    scatter!(ax3,tmp.lon,tmp.lat,markersize=3)
    xlims!(ax3, (-180, 180)); ylims!(ax3, (-90, 90))

    ax4 = Axis(f[2,2], title="profile positions for "*txt(ym2))
    tmp=df[df.ym.==ym2,:]
    scatter!(ax4,tmp.lon,tmp.lat,markersize=3)
    xlims!(ax4, (-180, 180)); ylims!(ax4, (-90, 90))

	f
end

"""
    map_stats(df::DataFrame,G::NamedTuple)

```
G=GriddedFields.load()
MITprof_plots.map_stats(df,G)
```
"""
function map_stats(df::DataFrame,G::NamedTuple)
    gdf=groupby(df,:pos)
    df2=combine(gdf) do df
        (m = log10(length(df.ID)), lon=mean(df.lon), lat=mean(df.lat))
    end

    tmp1=fill(0.0,(90,1170))
    ii=fill(0,2)
    for i in 1:size(df2,1)
        ii.=parse.(Int,split(split(split(df2[i,:pos],"(")[2],")")[1],","))
        tmp1[ii[1],ii[2]]=df2[i,:m]
    end

    #GriddedFields.MeshArrays.read(tmp1,G.msk.grid)
    XC=GriddedFields.MeshArrays.write(G.Γ.XC)
    YC=GriddedFields.MeshArrays.write(G.Γ.YC)
    ii=findall((!).(tmp1.==0))

    f = Figure()

    ax1 = Axis(f[1,1], title="number of profiles per day")
    scatter!(ax1,XC[ii],YC[ii],color=tmp1[ii])
    xlims!(ax1, -180, 180)
    ylims!(ax1, -90, 90)

    #(XC,YC,tmp1,f)
    f
end

end


text:
function lol()
	e3::Int32 = 0
	e3 = 0
end


text:
function solve(case, N, Q, doors, queries)
    answers = Int[]
    for (Si, Ki) in queries
        left = Si - 1
        right = Si
        last = Si

        for _ in 2:Ki
            if right > N - 1 || (left > 0 && doors[left] < doors[right])
                left -= 1
                last = left + 1
            else
                right += 1
                last = right
            end
        end

        push!(answers, last)
    end

    println("Case #$(case): $(join(answers, " "))")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N, Q = parse.(Int, stdin |> readline |> split)
        doors = parse.(Int, stdin |> readline |> split)

        queries = Tuple{Int,Int}[]
        for _ in 1:Q
            Si, Ki = parse.(Int, stdin |> readline |> split)
            push!(queries, (Si, Ki))
        end
        solve(case, N, Q, doors, queries)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
############################################################################ NN
# TODO there is extra hidden layer...

#function getmean(p::NN, features::Vector{Float64})
#   Win, bin, Wh, bh, Wout, bout, _ = split_theta!(p)
#   getmean(Win, bin, Wh, bh, Wout, bout, features)
#end

function getmean!{T<:AbstractFloat}(c::Vector{T}, p::NN,
                                    features::Vector{T})
   Win, bin, Wh, bh, Wout, bout, _ = split_theta!(p)
   #if p.nlayers == 2
   c[:] = getmean(Win, bin, Wh, bh, Wout, bout, features)
end

function getmean(Win, bin, Wh, bh, Wout, bout, features)
   a = tanh.(Win * features + bin)
   nlayers = size(Wh, 3)
   for i=1:nlayers
      a = tanh.(Wh[:,:,i] * a + bh[:,i])
   end
   return Wout * a + bout
end

function getmean2(Win, bin, Wh1, bh1, #Wh2, bh2,
                  Wout, bout, features)
   a = tanh.(Win * features + bin) # m->nhidden
   a[:] = tanh.(Wh1 * a + bh1) # nhidden -> nhidden
   return Wout * a + bout # nhidden -> m
end

# features = n
function getaction!{T<:AbstractFloat}(c::Vector{T}, p::NN,
                                      features::Vector{T})
   Win, bin, Wh, bh, Wout, bout, ls = split_theta!(p)
   #getmean!(c, p, features)
   @. c = c * exp(ls)
   a = tanh.(Win * features + bin) # m->nhidden
   a[:] = tanh.(Wh[:,:,1] * a + bh[:,1]) # nhidden -> nhidden
   c[:] += Wout * a + bout # nhidden -> m
end

function split_theta!(p::NN)
   n = p.n
   m = p.m
   nh = p.nhidden
   nlay = p.nlayers

   # win is matrix from n to nhidden
   Win = reshape( view(p.theta, p.ranges[:Win]), nh, n)
   bin = view( p.theta, p.ranges[:bin])

   # nlayers of nhidden*nhidden
   Wh = reshape( view(p.theta, p.ranges[:Wh]), nh, nh, nlay) # W[:,:,i]
   bh = reshape( view(p.theta, p.ranges[:bh]), nh, nlay) # b[:,i]

   # nhidden -> m actions
   Wout = reshape( view(p.theta, p.ranges[:Wout]), m, nh) # W[:,:]
   bout = reshape( view(p.theta, p.ranges[:bout]), m) # b[:]

   ls = reshape( view(p.theta, p.ranges[:ls]), m) # log_std

   return Win, bin, Wh, bh, Wout, bout, ls
end

function loglikelihood!{T<:AbstractFloat}(ll,
                                          theta::Vector{T}, p::NN,
                                          features::Matrix{T},
                                          actions::Matrix{T})
   const N    = size(features, 2) # N = T*numT
   const n    = p.n
   const m    = p.m
   const nh   = p.nhidden
   const rwin = p.ranges[:Win]
   const rbin = p.ranges[:bin]
   const rwh1 = p.ranges[:Wh1]
   const rbh1 = p.ranges[:bh1]
   const rwou = p.ranges[:Wout]
   const rbou = p.ranges[:bout]
   const rls  = p.ranges[:ls]

   const Win  = reshape( view( theta, rwin), nh, n)
   const bin  = view( theta, rbin)
   const Wh1  = reshape( view( theta, rwh1), nh, nh)
   const bh1  = view( theta, rbh1)
   const Wout = reshape( view( theta, rwou), m, nh)
   const bout = view( theta, rbou)
   const ls   = view( theta, rls)

   term1 = -0.5*m*log(2pi)
   term2 = -sum(ls)
   terms = term1+term2

   if size(p.scratch1) != (nh, N) p.scratch1 = zeros(nh, N) end
   if size(p.scratch2) != (nh, N) p.scratch2 = zeros(nh, N) end
   if size(p.scratch3) != (m, N) p.scratch3 = zeros(m, N) end

   # hidden 1
   A_mul_B!(p.scratch1, Win, convert(Matrix{T}, features))
   broadcast!(+, p.scratch1, p.scratch1, bin)
   p.scratch1 .= tanh.(p.scratch1)

   A_mul_B!(p.scratch2, Wh1, p.scratch1)
   broadcast!(+, p.scratch2, p.scratch2, bh1)
   p.scratch2 .= tanh.(p.scratch2)

   A_mul_B!(p.scratch3, Wout, p.scratch2)
   broadcast!(+, p.scratch3, p.scratch3, bout)

   p.scratch3 .= ((convert(Matrix{T},actions)-p.scratch3) ./ exp.(ls)).^2
   sum!(ll, p.scratch3)
   ll .= ll.*-0.5 .+ terms
end

function loglikelihood(theta, p::NN, features, actions)
   const n    = p.n
   const m    = p.m
   const nh   = p.nhidden
   const rwin = p.ranges[:Win]
   const rbin = p.ranges[:bin]
   const rwh1 = p.ranges[:Wh1]
   const rbh1 = p.ranges[:bh1]
   const rwou = p.ranges[:Wout]
   const rbou = p.ranges[:bout]
   const rls  = p.ranges[:ls]
   Win = reshape( theta[rwin], nh, n)
   bin = theta[rbin]

   Wh1 = reshape( theta[rwh1], nh, nh)
   bh1 = theta[rbh1]

   Wout = reshape( theta[rwou], m, nh)
   bout = theta[rbou]

   ls = theta[rls]

   dt = eltype(theta)
   #terms = dt(-0.5*p.m*log(2pi)) - sum(ls)
   terms = -0.5*p.m*log(2pi) - sum(ls)

   a = tanh.(broadcast(+, Win * features, bin)) # m->nhidden
   a = tanh.(broadcast(+, Wh1 * a, bh1)) # nhidden -> nhidden
   mu = broadcast(+, Wout * a, bout) # nhidden -> m

   zs = ((actions-mu) ./ exp.(ls)).^2
   ll = sum(zs, 1)
   ll .= ll.*-0.5 .+ terms
   return ll
end

function gradLL!{T<:AbstractFloat}(gradll::Matrix{T},
                                   p::NN,
                                   features::Matrix{T},
                                   actions::Matrix{T})
   const n     = p.n
   const m     = p.m
   const nh    = p.nhidden
   const rwin  = p.ranges[:Win]
   const rbin  = p.ranges[:bin]
   const rwh1  = p.ranges[:Wh1]
   const rbh1  = p.ranges[:bh1]
   const rwou  = p.ranges[:Wout]
   const rbou  = p.ranges[:bout]
   const rls   = p.ranges[:ls]
   const N = size(features, 2) # N = T*numT

   function ll_vec(theta, ff, aa)
      loglikelihood(theta, p, ff, aa)
   end

   #println("reverse split")
   BLAS.set_num_threads(1)
   nthread = min(N, Threads.nthreads())
   Threads.@threads for tid=1:nthread
   #@time for tid=1:nthread
      thread_range = Distributed.splitrange(N, nthread)[tid]

      const ff = zeros(T, n)
      const aa = zeros(T, m)

      const ctp = ReverseDiff.GradientTape(ll_vec, (p.theta, ff, aa))

      gr = map(similar, (p.theta, ff, aa))
      for i=thread_range
         ff .= features[:,i] # implicit conversion and memory saving
         aa .= actions[:,i]

         ReverseDiff.gradient!(gr, ctp, (p.theta, ff, aa)) # taped
         gradll[:,i] .= gr[1][:]
      end
   end
   BLAS.set_num_threads(Threads.nthreads())
end


text:
### type definitions ###

"""
A celestial object, such as a planet, star, or moon.
"""
abstract type CelestialObject end

"""
A Keplerian orbit, defined, by semimajor axis (`.a`), eccentricity (`.e`), inclination(`.i`),
right ascension of the ascending node (`.Ω`), argument of the periapsis (`.ω`),
mean anomaly at epoch(`.Mo`), and an epoch (`.epoch`), as well as the celestial body around
which the orbit takes place (.`primary`).
"""
struct Orbit{P}
    a::Float64
    e::Float64
    i::Float64
    Ω::Float64
    ω::Float64
    Mo::Float64
    epoch::Float64
    primary::P
    period::Float64
    basis::RotZXZ{Float64}
end

Base.hash(orb::Orbit, h::UInt) = hash_struct(orb, h, (:a, :e, :i, :Ω, :ω, :Mo, :epoch, :primary))

### helper methods ###
"""
    Δt = period(a, μ)

Compute the period of an orbit with semimajor axis `a` around a body with
gravity parameter `μ`.
"""
period(a, μ) = 2π*√(abs(a)^3 / μ)

## angle-based methods ##
"""
    ϕ = flight_path_angle(θ, e)
    ϕ = flight_path_angle(θ, orbit)

Compute the flight path angle of an orbit with eccentricity `e` at true anomaly `θ`.
"""
flight_path_angle(θ, e) = atan(e*sin(θ)/(1+e*cos(θ)))
flight_path_angle(θ, orb::Orbit) = flight_path_angle(θ, orb.e)

"""
    r = orbital_distance(θ, a, e)
    r = orbital_distance(θ, orbit)

Compute the orbital distance of an orbit with semimajor axis `a` and eccentricity `e`
at true anomaly `θ`, via the [conic-section orbit equation](https://en.wikipedia.org/wiki/Kepler_orbit#Development_of_the_laws).
"""
orbital_distance(θ, a, e) = a*(1-e^2)/(1+e*cos(θ))
orbital_distance(θ, orb::Orbit) = orbital_distance(θ, orb.a, orb.e)

"""
    θ = orbital_angle(r, a, e)
    θ = orbital_angle(r, orbit)

Compute the true anomaly of an orbit with semimajor axis `a` and eccentricity `e`
at orbital distance `r`, via the [conic-section orbit equation](https://en.wikipedia.org/wiki/Kepler_orbit#Development_of_the_laws).
"""
orbital_angle(r, a, e) = wrap_acos((a*(1-e^2)/r - 1)/e)
orbital_angle(r, orb::Orbit) = orbital_angle(r, orb.a, orb.e)

"""
    ψ = orbital_direction(θ, e)
    ψ = orbital_direction(r, a, e)

Compute the direction of velocity of an orbit with semimajor axis `a` at true anomaly `θ`.

Alternatively to true anomaly, orbital distance `r` and eccentricity `e` can be provided.
"""
orbital_direction(θ, e) = θ + π/2 - flight_path_angle(θ, e)
orbital_direction(r, a, e) = orbital_direction(orbital_angle(r,a,e), e)

"""
    r̄ = orbital_position(θ, a, e, basis)
    r̄ = orbital_position(θ, a, e, i, Ω, ω)
    r̄ = orbital_position(θ, orbit)

Compute the position vector at true anomaly `θ` of an orbit with with semimajor axis `a`,
eccentricity `e`, and a rotation `basis` defining the orbit's perifocal plane.

The rotation can be computed from inclinaion `i`, RAAN `Ω`, and the argument of the periapsis `ω`.
"""
function orbital_position(θ, a, e, basis::Rotation{3,Float64})
    r = orbital_distance(θ, a, e)
    ō = r * @SVector([cos(θ),sin(θ),0])
    return perifocal_to_inertial_bases(ō, basis)
end
orbital_position(θ, a, e, i, Ω, ω) = orbital_position(θ, a, e, basis_rotation(Ω, i, ω))
orbital_position(θ, orb::Orbit) = orbital_position(θ, orb.a, orb.e, orb.basis)

"""
    v̄ = orbital_velocity(θ, a, e, μ, basis)
    v̄ = orbital_velocity(θ, a, e, μ, i, Ω, ω)
    v̄ = orbital_velocity(θ, orbit)

Compute the velocity vector at true anomaly `θ` of an orbit with with semimajor axis `a`,
eccentricity `e`, gravity paremter `μ`, and a rotation `basis` defining the orbit's perifocal plane.

The rotation can be computed from inclinaion `i`, RAAN `Ω`, and the argument of the periapsis `ω`.
"""
function orbital_velocity(θ, a, e, μ, basis::Rotation{3,Float64})
    r = orbital_distance(θ, a, e)
    ϕ = orbital_direction(θ, e)
    v = √(μ*(2/r - 1/a))
    dōdt = v * @SVector([cos(ϕ),sin(ϕ),0])
    return perifocal_to_inertial_bases(dōdt, basis)
end
orbital_velocity(θ, a, e, μ, i, Ω, ω) = orbital_velocity(θ, a, e, μ, basis_rotation(Ω, i, ω))
orbital_velocity(θ, orb::Orbit) = orbital_velocity(θ, orb.a, orb.e, orb.primary.μ, orb.basis)

"""
    r̄, v̄ = state_vector(θ, a, e, μ, basis)
    r̄, v̄ = state_vector(θ, a, e, μ, i, Ω, ω)
    r̄, v̄ = state_vector(θ, orbit)

Compute the state vector at true anomaly `θ` of an orbit with with semimajor axis `a`,
eccentricity `e`, gravity paremter `μ`, and a rotation `basis` defining the orbit's perifocal plane.

The rotation can be computed from inclinaion `i`, RAAN `Ω`, and the argument of the periapsis `ω`.
"""
state_vector(θ, a, e, μ, basis) = (orbital_position(θ, a, e, basis), orbital_velocity(θ, a, e, μ, basis))
state_vector(θ, a, e, μ, i, Ω, ω) = state_vector(θ, a, e, μ, basis_rotation(Ω, i, ω))
state_vector(θ, orb::Orbit) = state_vector(θ, orb.a, orb.e, orb.primary.μ, orb.basis)

## time-based methods ##

"""
    r̄ = time_orbital_position(t, a, e, i, Ω, ω, Mo, epoch, μ)
    r̄ = time_orbital_position(t, orbit)

Compute the position vector at time `t` from a full set of orbital paremeters, as well as a gravity parameter `μ`.
"""
time_orbital_position(t, a, e, i, Ω, ω, Mo, epoch, μ
    ) = orbital_position(time_to_true(t, e, period(a,μ), Mo, epoch), a, e, i, Ω, ω)
time_orbital_position(t, orb::Orbit
    ) = orbital_position(time_to_true(t, orb.e, orb.period, orb.Mo, orb.epoch), orb.a, orb.e, orb.basis)

"""
    r̄ = time_orbital_velocity(t, a, e, i, Ω, ω, Mo, epoch, μ)
    r̄ = time_orbital_velocity(t, orbit)

Compute the velocity vector at time `t` from a full set of orbital paremeters, as well as a gravity parameter `μ`.
"""
time_orbital_velocity(t, a, e, i, Ω, ω, Mo, epoch, μ
    ) = orbital_velocity(time_to_true(t, e, period(a,μ), Mo, epoch), a, e, μ, i, Ω, ω)
time_orbital_velocity(t, orb::Orbit
    ) = orbital_velocity(time_to_true(t, orb.e, orb.period, orb.Mo, orb.epoch), orb.a, orb.e, orb.primary.μ, orb.basis)

"""
    r̄ = time_orbital_velocity(t, a, e, i, Ω, ω, Mo, epoch, μ)
    r̄ = time_orbital_velocity(t, orbit)

Compute the state vector at time `t` from a full set of orbital paremeters, as well as a gravity parameter `μ`.
"""
time_state_vector(t, a, e, i, Ω, ω, Mo, epoch, μ
    ) = state_vector(time_to_true(t, e, period(a,μ), Mo, epoch), a, e, μ, i, Ω, ω)
time_state_vector(t, orb::Orbit
    ) = state_vector(time_to_true(t, orb.e, orb.period, orb.Mo, orb.epoch), orb.a, orb.e, orb.primary.μ, orb.basis)

### Alternate Constructors ###

Orbit(a,body::CelestialObject
    ) = Orbit(a,0.,0.,0.,0.,0.,0.,body,period(a,body.μ),basis_rotation(0.,0.,0.))
Orbit(a,e,i,Ω,ω,Mo,body::CelestialObject
    ) = Orbit(a,e,i,Ω,ω,Mo,0.,body,period(a,body.μ),basis_rotation(Ω, i, ω))
Orbit(a,e,i,Ω,ω,Mo,epoch,body::CelestialObject
    ) = Orbit(a,e,i,Ω,ω,Mo,epoch,body,period(a,body.μ),basis_rotation(Ω, i, ω))

"""
    orb = Orbit(t, r̄, v̄, prim, epoch)

Compute an orbit from a time `t`, state vector (`r̄, v̄`), primary body `prim`, and epoch.
"""
function Orbit(t, r̄::AbstractVector{<:Real}, v̄::AbstractVector{<:Real}, prim::CelestialObject, epoch=t)
    r = norm(r̄)
    v = norm(v̄)
    μ = prim.μ
    if r == 0
        throw(ArgumentError("Invalid position"))
    end

    a = 1 / (2/r - (v^2)/μ)
    h̄ = cross(r̄,v̄)
    i = acos(h̄[3]/norm(h̄))

    ē = cross(v̄,h̄)/μ - r̄/r
    e = norm(ē)
    if isapprox(e, 0, atol=1e-15)       # if circular
        e = 0.
        ê = @SVector([1.,0.,0.])        # set periapsis at reference x
    else
        ê = ē/e
    end

    n̄ = cross(@SVector([0.,0.,1.]), h̄)
    n = norm(n̄)
    if isapprox(n, 0, atol=1e-15)       # if uninclined
        n = 0.
        n̂ = ê                           # set Ω at periapsis
    else
        n̂ = n̄/n
    end

    Ω = wrap_angle(copysign(1,n̂[2]) * wrap_acos(n̂[1]))
    ω = wrap_angle(copysign(1,ē[3]) * wrap_acos(dot(n̂,ê)))

    θ = angle_in_plane(r̄, basis_rotation(Ω, i, ω))
    M = true_to_mean(θ, e)
    Δt = t - epoch
    if e < 1
        Mo = wrap_angle(M - time_to_mean(Δt, period(a,μ)))
    else
        Mo = M - time_to_mean(Δt, period(a,μ))
    end

    return Orbit(a,e,i,Ω,ω,Mo,epoch,prim)
end

### comparison and sorting methods ###

function Base.isless(orb1::Orbit, orb2::Orbit)
    fnames = fieldnames(Orbit)
    for fn in fnames
        if orb1.fn != orb2.fn
            return isless(orb1.fn, orb2.fn)
        end
    end
    return false
end

# TO DO: account for possibility of negative angles (i.e. inclination)
function Base.isapprox(orb1::Orbit, orb2::Orbit; atol=1e-6, rtol=1e-6)
    a = isapprox(orb1.a, orb2.a; atol=atol, rtol=rtol)
    e = isapprox(orb1.e, orb2.e; atol=atol, rtol=rtol)
    b = isapprox(orb1.basis, orb2.basis;  atol=atol, rtol=rtol)
    pos = isapprox(time_orbital_position(0,orb1), time_orbital_position(0,orb2); atol=atol, rtol=rtol)
    return a && e && b && pos
end

### descriptive display ###

Base.show(io::IO, orb::Orbit) = println(io,
    "Orbit above $(orb.primary.name):\n",
    "  semi-major axis:\t\t$(orb.a) m\n",
    "  eccentricity:\t\t\t$(orb.e)\n",
    "  inclination:\t\t\t$(rad2deg(orb.i))°\n",
    "  RA of ascending node:\t\t$(rad2deg(orb.Ω))°\n",
    "  argument of periapsis:\t$(rad2deg(orb.ω))°\n",
    "  mean anomaly at epoch:\t$(rad2deg(orb.Mo))°\n",
    "  epoch:\t\t\t$(orb.epoch) s"
)


text:
get_matpower_cols() = Dict(
    "bus" => ["ID", "TYPE", "PD", "QD", "GS", "BS", "BUS_AREA", "VM", "VA", "BASE_KV", "ZONE", "VMAX", "VMIN", "LAM_P",
              "LAM_Q", "MU_VMAX", "MU_VMIN"],
    "gen" => ["ID", "PG", "QG", "QMAX", "QMIN", "VG", "MBASE", "GEN_STATUS", "PMAX", "PMIN", "PC1", "PC2", "QC1MIN",
              "QC1MAX", "QC2MIN", "QC2MAX", "RAMP_AGC", "RAMP_10", "RAMP_30", "RAMP_Q", "APF", "MU_PMAX", "MU_PMIN",
              "MU_QMAX", "MU_QMIN"],
    "branch" => ["SRC", "DST", "BR_R", "BR_X", "BR_B", "RATE_A", "RATE_B", "RATE_C", "TAP", "SHIFT", "BR_STATUS",
                 "ANGMIN", "ANGMAX", "PF", "QF", "PT", "QT", "MU_SF", "MU_ST", "MU_ANGMIN", "MU_ANGMAX"],
    "gencost" => ["MODEL", "STARTUP", "SHUTDOWN", "NCOST", "COST"]
)

function get_matpower_gencost_cols(gencost)
    if size(gencost, 1) == 0
        return ["MODEL", "STARTUP", "SHUTDOWN", "NCOST"]
    end
    nb_costs = size(gencost, 2) - 4
    return vcat(["MODEL", "STARTUP", "SHUTDOWN", "NCOST"], ["COST$i" for i in 1:nb_costs])
end

function _row_to_array_matpower_mat(row)
    row = split(row, ';')[1]
    values = split(row, '\t')[2:end]
    array = map(x -> parse(Float64, x), values)
    return array
end

function _parse_matrix_matpower_mat(lines, start)
    nb_cols = length(split(lines[start], '\t')) - 1
    current_line = start
    matrix = nothing
    while !occursin("];", lines[current_line])
        if isnothing(matrix)
            matrix = transpose(_row_to_array_matpower_mat(lines[current_line]))
        else
            matrix = [matrix; transpose(_row_to_array_matpower_mat(lines[current_line]))]
        end
        current_line += 1
    end
    return matrix
end

function get_data_matpower_m(path::AbstractString)
    file_string = read(open(path, "r"), String)
    lines = split(file_string, '\n')
    lines = [l for l in lines if !startswith(lstrip(l), "%")]
    data = Dict()
    occursin_pattern = ["mpc.bus", "mpc.gen ", "mpc.branch", "mpc.gencost"]
    occursin_key = ["bus", "gen", "branch", "gencost"]
    dims = [17, 25, 21]
    pattern_idx = 1
    for (i, f) in enumerate(lines)
        if occursin("mpc.baseMVA", f)
            data["baseMVA"] = parse(Float64, split(f[15:end], ';')[end - 1])
        elseif occursin(occursin_pattern[pattern_idx], f)
            key = occursin_key[pattern_idx]
            mat = _parse_matrix_matpower_mat(lines, i + 1)
            sizes = size(mat)
            if key == "gencost"
                data[key] = zeros(sizes[1], sizes[2])
            else
                data[key] = zeros(sizes[1], dims[pattern_idx])
            end
            data[key][1:sizes[1], 1:sizes[2]] = mat
            pattern_idx += 1
            if pattern_idx > length(occursin_pattern)
                break
            end
        end
    end
    if !haskey(data, "gencost")
        data["gencost"] = zeros(0, 4)
    end
    return data["bus"], data["gen"], data["branch"], data["gencost"], data["baseMVA"]
end

function get_data_matpower_mat(path::AbstractString)
    error("Not Implemented")
end

nbus_matpower_mat(path::AbstractString) = error("Not Implemented")
nbranch_matpower_mat(path::AbstractString; distinct_pair=false) = error("Not Implemented")
ngen_matpower_mat(path::AbstractString; distinct_pair=false) = error("Not Implemented")
ntransformer_matpower_mat(path::AbstractString; distinct_pair=false) = error("Not Implemented")

function nbus_matpower_m(path::AbstractString)
    file_string = read(open(path, "r"), String)
    buses = split(file_string, "mpc.bus")[2]
    buses = split(buses, "];")[1]
    return length(split(buses, '\n')) - 2
end

function nbranch_matpower_m(path::AbstractString; distinct_pair=false)
    nbranch = nothing
    file_string = read(open(path, "r"), String)
    branches = split(file_string, "mpc.branch")[2]
    branches = split(branches, "];")[1]
    branches = split(branches, '\n')[2:end - 1]
    filter!(b -> !startswith(b, "%"), branches)
    if distinct_pair
        branches = [Set(parse.(Int, split(b, '\t')[2:3])) for b in branches]
        nbranch = length(unique(branches))
    else
        nbranch = length(branches)
    end
    return nbranch
end

function ntransformer_matpower_m(path::AbstractString; distinct_pair=false)
    return 0
end

function ngen_matpower_m(path::AbstractString; distinct_pair=false)
    nbranch = nothing
    file_string = read(open(path, "r"), String)
    branches = split(file_string, "mpc.gen")[2]
    branches = split(branches, "];")[1]
    if distinct_pair
        branches = split(branches, '\n')[2:end - 1]
        branches = [parse(Int, split(b, '\t')[2]) for b in branches]
        nbranch = length(unique(branches))
    else
        nbranch = length(split(branches, '\n')) - 2
    end
    return nbranch
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using OrderedResponse
using Optim
using StatsFuns
using Calculus
using DataFrames
using StatsModels
using StatsBase
using CategoricalArrays
using Test
using CSV
using LinearAlgebra
using StatsFuns
using Distributions: _F1

include("benchmarks.jl")

# read in data
# df = CSV.read(joinpath(dirname(pathof(OrderedResponse) * "/data/testdat.csv")))

@testset "test dlogF for :probit" begin
    for model in (:probit,:logit,)
        for a in [-6.0:5.0..., typemin(1.0),]
            for b in [(max(a,-5)+1):6.0..., typemax(1.0),]
                FF = OrderedResponse.cdf(b, Val{model}) - OrderedResponse.cdf(a, Val{model})
                ff = OrderedResponse.pdf(b, Val{model}) - OrderedResponse.pdf(a, Val{model})
                dlogF  = ff / FF
                dlogF_new = OrderedResponse.dlogcdf_trunc(a,b,Val{model})
                @test dlogF ≈ dlogF_new
                @test dlogF_new != 0.0 || ff == 0.0
                @test FF > 0
            end
        end
    end
end

df = CSV.read("D:/libraries/julia/dev/OrderedResponse" * "/data/testdat.csv")

# test functions in "likelihood" file
include("likelihood.jl")

# test outer wrapper
fm = @formula(y ~ 0 + x1 + x2)   # in StatsModels
mf = ModelFrame(fm, df)
y = OrderedResponse.response_vec(mf)
X = ModelMatrix(mf).m'

orlm(y, X, :logit)
orlm(y, X, :probit)

orlm(fm, df, :logit)
orlm(fm, df, :probit)

#


text:
#=
This file attempts to implement the "imporoved dimension estimation" by Judd:

Judd, K. (1992). An improved estimator of dimension and some comments on providing
confidence intervals. Physica D: Nonlinear Phenomena, 56(2–3), 216–228. https://doi.org/10.1016/0167-2789(92)90025-I

Judd, K. (1994). Estimating dimension from small samples. Physica D: Nonlinear Phenomena,
71(4), 421–429. https://doi.org/10.1016/0167-2789(94)90008-6

So far this attempt has failed.
=#

function interpoint_distances(X, norm = Euclidean())
    N, ds = length(X), eltype(X)[]
    sizehint!(ds, Int(N*(N-1)/2))
    @inbounds for i in 1:N
        @inbounds for j in i+1:N
            push!(ds, evaluate(norm, X[i], X[j]))
        end
    end
    return ds
end

"""
    logspace_histogram(x::Vector, λ = exp(-w0), ϵ₀ = λ^2 * maximum(x)) -> bins, counts
Calculate a logarithmically spaced histogram of the values in `x` (must be positive
definite, i.e. distances) and return the bins and count of values in the bins.

The `i`-th bin is `[ϵ₀*λ^i, ϵ₀*λ^(i+1))`. Elements of `x` greater than `ϵ₀` are disregarded.
"""
function logspace_histogram(x::AbstractVector, λ = exp(-w0(x)), ϵ₀ = λ^2 * maximum(x))
    @assert λ < 1
    bins = [ϵ₀]
    λmin = λ*minimum(x)
    λmin == 0 && error("x cannot contain zero values!")
    for i in 1:length(x)
        push!(bins, bins[end]*λ)
        bins[end] < λmin && break
    end
    reverse!(bins)
    L = length(bins)
    counts = zeros(Int, L)
    # Fill bins with values
    for v in x
        i = searchsortedfirst(bins, v)
        i == L + 1 && continue
        bins[i] += 1
    end
    return reverse!(bins), reverse!(counts)
end
# TODO: test the algorithm with including the bin [∞, ε0)
w0(x) = log(maximum(x)/minimum(x))/√length(x)

# The two papers are really complicated. The way to really do this that I will
# follow is given by the 1992 Judd paper, page 223, second column.
# p is a function ε^d * polynomial of degree `t`, with `t` " = 2 being a good choice."
# Then the formula for L (no equation number, because why write good papers) is
# "minimized" by some routine.
# Notice that pᵢ is just the function of the ϵᵢ , the bin edges.

# To estimate the dimension we have to do an optimization solve of a problem...? wtf!

p(a, bins) = @. (bins^a[1]) * (a[2] + a[3]*bins + a[4]*bins^2)
function judd_estimator(bins, count)
    f = (a) -> -sum( bins .* log.(p(a, bins)) )
    optimize(f, ones(4))
end

# The problem now is OF COURSE, that this optimization problem is ill defined.
# a polynomial can OF COURSE get negative values, and OF COURSE the logarithm
# of a negative number is not real. How do these papers get published my god...


text:
"""
rw_test(p; method = "VR", m = 1)

TEST OF RANDOM WALK

Random walk with uncorrelated increments (RW3) is the weakest form version of the random walk.
The other versions are random walk with iid increments (RW1) and random walk with independent increments (RW2).

Box and Pierce (1970) test of RW(3): the autocorrelation of increments at all leads and lags should be zero under random walk.
By summing the squared autocorrelation, the Box-Pierce Q-statistics is designed to detect departures from zero autocorrelations in either directions and at all lags

Ljung and Box (1978) test of RW(3): the Ljung-Box Q-statistics provides the finite-sample correction which yields a better fit to the chi-squared distribution for small sample size.

Variance ratio test of RW(1): the variance of random walk increments must be a linear function of the time interval under random walk. The variance ratios statistics tests such property.

INPUT
`p': a T x 1 vector of asset prices
`method': method used to test random walk hypthesis with uncorrelated increments:
        - BP: the Box and Pierce (1970) statistics
        - LB: the Ljung and Box (1978) statistics
        - VR: variance ratio test
`m': number of autocorrelation lags used for test

OUTPUT
`test_statistics': test statistics of the random walk test
`pvalue': p-value of the test statistics

"""

function rw_test(p, method = "VR", m = 1)
    r = diff(log.(p), dims = 1)
    T = length(r); # length of the time series

    if method == "BP" # Box and Pierce (1970) test of RW(3)
        Qm = 0.0;
        for k = 1:m
            rho_k = autocor(r, [k]; demean = false);
            Qm = Qm .+ rho_k.^2;
        end
        test_statistics = T * Qm;
        pvalue = 1.0 .- cdf.(Chisq(m), test_statistics);
    elseif method == "LB" # Ljung and Box (1978) test of RW(3)
        Qm = 0.0;
        for k = 1:m
            rho_k = autocor(r, [k]; demean = false);
            Qm = Qm .+ rho_k.^2;
        end
        test_statistics = T * (T + 2) / (T - m) * Qm;
        pvalue = 1.0 .- cdf.(Chisq(m), test_statistics);
    elseif method == "VR" # Variance ratio test of RW(1)
        VR = 1.0;
        for k = 1:m-1
            rho_k = autocor(r, [k]; demean = false);
            VR = VR .+ 2 .* (1 .- k / m) .* rho_k;
        end
        test_statistics = sqrt(T * m) * (VR .- 1.0) ./ sqrt(2.0 .* (m - 1));
        pvalue = 1.0 .- cdf.(Normal(0.0, 1.0), test_statistics);
    end

    return (test_statistics, pvalue)
end


text:
export getparticle,getparticles

function getparticle(Collection,Meta)
	return getparticles(Collection,Meta,1)
end

function getparticles(Collection, Meta,nparticles)
	PArray = Array{Particle,1}(undef, nparticles)
    edgesx,edgesy  = Meta.edgesx, Meta.edgesy
    for particleid=1:nparticles
		foundparticle = false
		r = 0.0
		x,y = 0.0,0.0
		nattempt = 0
		while !foundparticle
			nattempt +=1
			x,y = rand(2)
			i,j  = findindex(x,edgesx), findindex(y,edgesy)
			if length(Collection[i,j].x)==0
				break
			end
			foundparticle = true

			# Make sure the particle is not already inside an obstacle

			dists =  sqrt.((x.-Collection[i,j].x).^2 .+(y.-Collection[i,j].y).^2) .- Collection[i,j].r
			if minimum(dists) < 0.0  foundparticle=false; end
			if nattempt>1000 error("Did not find initial position after $nattempt attempts."); end
		end

		vx,vy = randn(2)
		normv = norm([vx,vy])
		vx /= normv
		vy /= normv
		PArray[particleid] = Particle(x,y,r,vx,vy)
	end
    return PArray

end


text:
mutable struct DoubleMacaulayMatrix
    leftrows::Vector{Vector{Int}}
    rightrows::Vector{Vector{Int}}

    leftcoeffs::Vector{Vector{UInt64}}
    rightcoeffs::Vector{Vector{UInt64}}

    lefthash2col::Vector{Int}
    leftcol2hash::Vector{Int}

    righthash2col::Vector{Int}
    rightcol2hash::Vector{Int}

    leftsize::Int
    rightsize::Int
end

function initialize_double_matrix(ring::PolyRing)
    uprows   = Vector{Vector{Int}}(undef, 0)
    lowrows  = Vector{Vector{Int}}(undef, 0)
    col2hash = Vector{Int}(undef, 0)
    coeffs   = Vector{Vector{UInt64}}(undef, 0)

    size    = 0
    npivots = 0
    nrows   = 0
    ncols   = 0

    nup    = 0
    nlow   = 0
    nleft  = 0
    nright = 0

    up2coef   = Vector{Int}(undef, 0)
    low2coef   = Vector{Int}(undef, 0)

    DoubleMacaulayMatrix(uprows, lowrows, col2hash, coeffs,
            size, npivots, nrows, ncols,
            nup, nlow, nleft, nright,
            up2coef, low2coef)
end

function linear_relation!(
            matrix::DoubleMacaulayMatrix,
            monom::Int, vector::Basis)
    denserow = convert_vector_to_dense_row(matrix, monom, vector)
    reduce_double_dense_row_by_known_pivots_sparse!(matrix, denserow)

    if there is a relation
        return true
    else
        false
end


text:
"An isomorphism is a bijection between types that preserves norms."
abstract type Isomorphism{T,U} <: TypedMap{T,U} end

show(io::IO, m::Isomorphism{T,U}) where {T,U} = print(io, "x : $(T) -> x : $(U)")
Display.object_parentheses(m::Isomorphism) = true

"Map a length 1 vector `x` to `x[1]`."
struct VectorToNumber{T} <: Isomorphism{SVector{1,T},T}
end
"Map a number `x` to a length 1 vector `[x]`."
struct NumberToVector{T} <: Isomorphism{T,SVector{1,T}}
end

mapsize(::VectorToNumber) = (1,1)
mapsize(::NumberToVector) = (1,)

inverse(::VectorToNumber{T}) where {T} = NumberToVector{T}()
inverse(::NumberToVector{T}) where {T} = VectorToNumber{T}()
inverse(m::VectorToNumber, x) = inverse(m)(x)
inverse(m::NumberToVector, x) = inverse(m)(x)

applymap(::VectorToNumber, x) = x[1]
applymap(::NumberToVector, x) = SVector(x)

jacobian(::VectorToNumber{T}, x) where {T} = transpose(SVector(one(T)))
jacobian(::VectorToNumber{T}) where {T} = ConstantMap{SVector{1,T}}(transpose(SVector(one(T))))

jacobian(::NumberToVector{T}, x) where {T} = SVector(one(T))
jacobian(::NumberToVector{T}) where {T} = ConstantMap{T}(SVector(one(T)))

"Map a length 2 vector `x` to `x[1] + im*x[2]`."
struct VectorToComplex{T} <: Isomorphism{SVector{2,T},Complex{T}}
end
"Map a complex number `x` to the length 2 vector `[real(x); imag(x)]`."
struct ComplexToVector{T} <: Isomorphism{Complex{T},SVector{2,T}}
end

mapsize(::VectorToComplex) = (1,2)
mapsize(::ComplexToVector) = (2,)

applymap(::VectorToComplex, x) = x[1] + im*x[2]
applymap(::ComplexToVector, x) = SVector(real(x), imag(x))

inverse(::VectorToComplex{T}) where {T} = ComplexToVector{T}()
inverse(::ComplexToVector{T}) where {T} = VectorToComplex{T}()
inverse(m::VectorToComplex, x) = inverse(m)(x)
inverse(m::ComplexToVector, x) = inverse(m)(x)

jacobian(::VectorToComplex{T}, x) where {T} = transpose(SVector(one(T),one(T)*im))
jacobian(::VectorToComplex{T}) where {T} = ConstantMap{SVector{2,T}}(transpose(SVector(one(T),one(T)*im)))

"Map a static vector to a tuple."
struct VectorToTuple{N,T} <: Isomorphism{SVector{N,T},NTuple{N,T}}
end
"Map a tuple to a static vector."
struct TupleToVector{N,T} <: Isomorphism{NTuple{N,T},SVector{N,T}}
end

inverse(::VectorToTuple{N,T}) where {N,T} = TupleToVector{N,T}()
inverse(::TupleToVector{N,T}) where {N,T} = VectorToTuple{N,T}()
inverse(m::VectorToTuple, x) = inverse(m)(x)
inverse(m::TupleToVector, x) = inverse(m)(x)

applymap(::VectorToTuple, x) = tuple(x...)
applymap(::TupleToVector, x) = SVector(x)

"Map a nested vector or tuple to a flat vector."
struct NestedToFlat{N,T,U,DIM} <: Isomorphism{U,SVector{N,T}}
end
"Map a flattened vector to a nested one."
struct FlatToNested{N,T,U,DIM} <: Isomorphism{SVector{N,T},U}
end

inverse(::NestedToFlat{N,T,U,DIM}) where {N,T,U,DIM} = FlatToNested{N,T,U,DIM}()
inverse(::FlatToNested{N,T,U,DIM}) where {N,T,U,DIM} = NestedToFlat{N,T,U,DIM}()
inverse(m::NestedToFlat, x) = inverse(m)(x)
inverse(m::FlatToNested, x) = inverse(m)(x)

applymap(::NestedToFlat{N,T,U,DIM}, x) where {N,T,U,DIM} = convert_tocartesian(x, Val{DIM}())
applymap(::FlatToNested{N,T,U,DIM}, x) where {N,T,U,DIM} = convert_fromcartesian(x, Val{DIM}())


text:
using Colors: RGB, N0f8, Gray
using FastAI
using FastAI: ParamGroups, IndexGrouper, getgroup, DiscriminativeLRs, decay_optim
import FastAI: Image, Keypoints, Mask, testencoding, Label, OneHot, ProjectiveTransforms,
    encodedblock, decodedblock, encode, decode, mockblock, checkblock, Block, Encoding,
    imagedatasetstats
import Makie
using FilePathsBase
using FastAI.Datasets
using FastAI.Models
using DLPipelines
import DataAugmentation
import DataAugmentation: getbounds, NormalizeRow
using Flux
using Flux.Optimise: Optimiser, apply!
using StaticArrays
using Test
using TestSetExtensions
using DataFrames
using Tables
using CSV
import ReTest

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
include("testdata.jl")


text:
import GBIF
using DataFrames
import CSV

sp_name = "Amphiprion"
sp = GBIF.taxon(sp_name, strict=false)

# Get 900 occurrences
occ = GBIF.occurrences(sp, "hasCoordinate" => "true", "limit" => 300)
GBIF.occurrences!(occ)
GBIF.occurrences!(occ)
GBIF.occurrences!(occ)
GBIF.occurrences!(occ)
GBIF.occurrences!(occ)
GBIF.occurrences!(occ)

raw_data = DataFrame(occ)

# JUDGE NOT LEST YE BE JUDGED
raw_data = raw_data[raw_data.rank .== "SPECIES", :]
# Seriously
ok_names = "Amphiprion ".*["melanopus", "ocellaris", "clarkii"]
# It's 10pm
ok_index = map(f -> in(f, ok_names), raw_data.name)
# I'm so tired
raw_data = raw_data[ok_index, :]
# Can 2020 be over yet, please?
raw_data = raw_data[raw_data.longitude.>=0.0, :]
# I really should have used Query.jl or DataFramesMeta.jl is what I'm trying to say

# Save only the columns we care about
data = select(raw_data, [:name, :longitude, :latitude])
CSV.write("occurrences.csv", data; writeheader=true)


text:
"Abstract type defining backends which determine the behavior of functions that
act on TensorNetworkCircuits"
abstract type AbstractBackend end

"""
A struct to record various costs of storing and contracting a tensor network.
"""
mutable struct Metrics
    # Maximum number of elements in a single tensor
    max_tensor_size::Int

    # Total number of elements required to store the tensor network and
    # intermidiate tensors during contraction.
    total_space_allocated::Int

    # Number of multiplications required to contract a network
    flops::Int

    Metrics() = new(0, 0, 0)
end

"""
    record_compute_costs!(backend::AbstractBackend,
                                   open_dims::Array{Int, 1},
                                   contracted_dims::Array{Int, 1})

Function to calculate the space and time costs of contracting two tensors
and add them to the metrics contained in the given backend.
"""
function record_compute_costs!(backend::AbstractBackend,
                               open_dims::Array{Int, 1},
                               contracted_dims::Array{Int, 1})

    space_cost = prod(open_dims)
    record_memory_costs!(backend, space_cost)

    flops_cost = space_cost * prod(contracted_dims)
    backend.metrics.flops += flops_cost

    flops_cost, space_cost
end

"""
    record_memory_costs!(backend::AbstractBackend,
                                  tensor_memory::Int)

Function to add the space cost of a tensor to the metrics contained in
the given backend.
"""
function record_memory_costs!(backend::AbstractBackend,
                              tensor_memory::Int)

    backend.metrics.total_space_allocated += tensor_memory
    if tensor_memory > backend.metrics.max_tensor_size
       backend.metrics.max_tensor_size = tensor_memory
    end
end

include("backends/interactive.jl")
include("backends/dsl.jl")

export AbstractBackend
export save_tensor_data, load_tensor_data, save_output
export contract_tensors, record_costs!
export reshape_tensor, permute_tensor
export decompose_tensor!, view_tensor!, delete_tensor!


text:
# ----- Computational metrics -----
# Relative norm between 2 vectors
"""
    relativenorm(x, x₀[, p]) where T<:Number

Returns the relative norm of base p between original signal x₀ and noisy signal
x.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Signal with noise.
- `x₀::AbstractArray{T} where T<:Number`: Reference signal.
- `p::Real`: (Default: 2) ``p``-norm to be computed.

# Returns
`::AbstractFloat`: The relative norm between x and x₀.

# Examples
```@repl
using WaveletsExt

x = randn(8)
y = randn(8)

relativenorm(x, y)
```

**See also:** [`psnr`](@ref), [`snr`](@ref), [`ssim`](@ref)
"""
function relativenorm(x::AbstractArray{T},
                      x₀::AbstractArray{T},
                      p::Real = 2) where T<:Number
    @assert length(x) == length(x₀)             # ensure same lengths
    return norm(x-x₀,p)/norm(x₀,p)
end

# PSNR between 2 vectors
@doc raw"""
    psnr(x, x₀)

Returns the peak signal to noise ratio (PSNR) between original signal x₀ and noisy signal x.

PSNR definition: ``10 \log_{10} \frac{\max{(x_0)}^2}{MSE(x, x_0)}``

# Arguments
- `x::AbstractVector{T} where T<:Number`: Signal with noise.
- `x₀::AbstractVector{T} where T<:Number`: Reference signal.

# Returns
`::AbstractFloat`: The PSNR between x and x₀.

# Examples
```@repl
using WaveletsExt

x = randn(8)
y = randn(8)

psnr(x, y)
```

**See also:** [`relativenorm`](@ref), [`snr`](@ref), [`ssim`](@ref)
"""
function psnr(x::AbstractVector{T}, x₀::AbstractVector{T}) where T<:Number
    @assert length(x) == length(x₀)              # ensure same lengths
    sse = zero(T)
    for i in eachindex(x)
        @inbounds sse += (x[i] - x₀[i])^2
    end
    mse = sse/length(x)
    return 20 * log(10, maximum(x₀)) - 10 * log(10, mse)
end

# SNR between 2 vectors
@doc raw"""
    snr(x, x₀)

Returns the signal to noise ratio (SNR) between original signal x₀ and noisy signal x.

SNR definition: ``20 \log_{10} \frac{||x_0||_2}{||x-x_0||_2}``

# Arguments
- `x::AbstractVector{T} where T<:Number`: Signal with noise.
- `x₀::AbstractVector{T} where T<:Number`: Reference signal.

# Returns
`::AbstractFloat`: The SNR between x and x₀.

# Examples
```@repl
using WaveletsExt

x = randn(8)
y = randn(8)

snr(x, y)
```

**See also:** [`relativenorm`](@ref), [`psnr`](@ref), [`ssim`](@ref)
"""
function snr(x::AbstractVector{T}, x₀::AbstractVector{T}) where T<:Number
    @assert length(x) == length(x₀)             # ensure same lengths
    return 20 * log(10, norm(x₀,2)/norm(x-x₀,2))
end

# SSIM between 2 arrays
"""
    ssim(x, x₀)

Wrapper for `assess_ssim` function from ImageQualityIndex.jl.

Returns the Structural Similarity Index Measure (SSIM) between the original signal/image x₀
and noisy signal/image x.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Signal with noise.
- `x₀::AbstractArray{T} where T<:Number`: Reference signal.

# Returns
`::AbstractFloat`: The SNR between x and x₀.

# Examples
```@repl
using WaveletsExt

x = randn(8)
y = randn(8)

ssim(x, y)
```

**See also:** [`relativenorm`](@ref), [`psnr`](@ref), [`snr`](@ref)
"""
function ssim(x::AbstractArray{T}, x₀::AbstractArray{T}) where T<:Number
    return assess_ssim(x, x₀)
end


text:
module Lucas

include("string.jl")
export strcut

if Pkg.installed("DICOM") != nothing
    import Base.show
    include("dicom.jl")
end

end # module


text:
C=[1 1 1 1; 0 1 0 1; 0 0 1 1; 0 0 0 1]
P=[1/3 1/3 1/3 0; 1/3 1/3 0 1/3; 1/3 0 1/3 1/3; 0 1/3 1/3 1/3 ]

C1=[1 1 1 1 1 1; 0 1 0 1 1 1; 0 0 1 1 1 1; 0 0 0 1 0 1; 0 0 0 0 1 1; 0 0 0 0 0 1]
P3=[1/2 1/6 0 1/3 0 0; 1/3 1/6 1/6 1/3 0 0; 1/3 1/6 0 1/3 1/6 0; 1/6 1/6 1/6 0 1/6 1/3; 1/6 0 1/6 1/6 1/3 1/6; 0 1/6 1/3 0 1/6 1/3]
P6=[17/24 0 0 1/8 1/8 1/24; 1/8 5/16 5/16 1/12 1/12 1/12; 1/8 5/16 5/16 1/12 1/12 1/12; 1/12 1/12 1/12 5/16 5/16 1/8; 1/12 1/12 1/12 5/16 5/16 1/8; 1/24 1/8 1/8 0 0 17/24]
P7=[1/3 1/3 1/3 0 0 0; 1/3 1/3 0 1/3 0 0; 1/3 0 1/3 1/3 0 0; 0 1/3 1/3 1/3 0 0; 0 1/6 1/6 1/6 1/6 1/3; 0 1/6 1/6 1/6 1/6 1/3]
P11=[1/3 1/6 1/6 1/6 1/6 0; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 0 1/6 1/6 1/6 1/6 1/3]
P12=[1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 0 1/3 1/6 1/6 1/6 1/6]
P13=[1/6 1/6 1/6 1/6 1/3 0; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6]
P14=[1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6; 1/6 1/6 1/6 1/6 1/6 1/6]
P15=[17/24 0 0 1/8 1/8 1/24; 1/8 5/16 5/16 1/12 1/12 1/12; 1/8 5/16 5/16 1/12 1/12 1/12; 1/12 1/12 1/12 5/16 5/16 1/8; 1/12 1/12 1/12 5/16 5/16 1/8; 1/24 1/16 1/16 1/16 1/16 17/24]
P16=[17/24 1/16 1/16 1/16 1/16 1/24; 1/8 5/16 5/16 1/12 1/12 1/12; 1/8 5/16 5/16 1/12 1/12 1/12; 1/12 1/12 1/12 5/16 5/16 1/8; 1/12 1/12 1/12 5/16 5/16 1/8; 1/24 1/8 1/8 0 0 17/24]

C2=[1 1 1 1 1 1 1 1; 0 1 0 0 1 1 0 1; 0 0 1 0 1 0 1 1; 0 0 0 1 0 1 1 1; 0 0 0 0 1 0 0 1; 0 0 0 0 0 1 0 1; 0 0 0 0 0 0 1 1; 0 0 0 0 0 0 0 1]
P1=[3/8 1/8 1/8 1/8 0 1/8 1/8 0; 1/4 1/4 0 1/8 1/8 0 1/8 1/8; 1/4 0 3/8 0 0 1/4 0 1/8; 1/4 0 1/8 0 1/8 1/4 1/8 1/8; 0 1/8 1/8 0 0 1/4 1/4 1/4; 0 0 1/4 1/4 0 1/8 0 3/8; 0 0 1/8 1/8 1/4 1/4 1/8 1/8; 0 0 1/8 0 1/8 1/8 0 5/8]
P2=[5/8 0 1/8 1/8 0 1/8 0 0; 1/8 1/8 1/4 1/4 1/8 1/8 0 0; 3/8 0 1/8 0 1/4 1/4 0 0; 1/4 1/4 1/4 0 0 1/8 1/8 0; 1/8 1/8 1/4 1/8 0 1/8 0 1/4; 1/8 0 1/4 0 0 3/8 0 1/4; 1/8 1/8 0 1/8 1/8 0 1/4 1/4; 0 1/8 1/8 0 1/8 1/8 1/8 3/8]

C3=[1 1 1 1 1; 0 1 0 0 1; 0 0 1 0 1; 0 0 0 1 1; 0 0 0 0 1]
P4=[2/5 0 2/5 1/5 0; 0 2/5 1/5 2/5 0; 0 1/5 2/5 2/5 0; 0 1/5 1/5 1/5 2/5; 0 1/5 1/5 1/5 2/5]
P5=[2/5 1/5 1/5 1/5 0; 2/5 1/5 1/5 1/5 0; 0 2/5 2/5 1/5 0; 0 2/5 1/5 2/5 0; 0 1/5 2/5 0 2/5]
P8=[2/5 1/5 1/5 1/5 0; 2/5 1/5 1/5 1/5 0; 2/5 0 1/5 1/5 1/5; 1/5 1/5 2/5 1/5 0; 0 2/5 1/5 0 2/5]
P9=[2/5 0 1/5 2/5 0; 0 1/5 2/5 1/5 1/5; 1/5 1/5 1/5 0 2/5; 0 1/5 1/5 1/5 2/5; 0 1/5 1/5 1/5 2/5]

is_mobius_down(C,P)
is_mobius_up(C,P)
is_weak_down(C,P)
is_weak_up(C,P)
is_stochastic(C,P)
is_realizable(C,P)

which_is(C,P)

which_is(C1,P3)
which_is_not_real(C1,P6)
which_is(C1,P7)
which_is(C1,P11)
which_is(C1,P12)
which_is(C1,P13)
which_is(C1,P14)
which_is_not_real(C1,P15)
which_is_not_real(C1,P16)

which_is(C2,P1)
which_is(C2,P2)

which_is(C3,P4)
which_is(C3,P5)
which_is(C3,P8)
which_is(C3,P9)

searching_for_stoch_mob_not_real(100)

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
struct LossLog{time,time_titer,data,M,control,titer,shift} <: DiffEqBase.DECostFunction
  t::time #Time points for data
  t_titer::time_titer #Time points for virus count
  data::data # LFC of data
  measured::M #Indices in model with data
  mock::control #raw data for control
  vtiter::titer #Viral titer estimates
  shift::shift #Shift in data from noise
end

mutable struct NN
  sizeLayers
  numLayers
  biases
  weights
end

function (f::LossLog)(sol::DESolution)

t = f.t
data = f.data
measured = f.measured
mock = f.mock
vtiter = f.vtiter
shift=f.shift

  if sol.retcode != :Success
      return Inf
  end

  subSol = sol(t)

  sumsq = 0.0
  maxlog2(x) = log2(max(0,x)+1)
  @inbounds for i = 1:length(subSol) #Loop through each data time point
    for (j,m) in enumerate(measured) #loop though states
      sumsq +=(maxlog2(subSol[m,i]+shift[m]) - maxlog2(mock[i,j]) - data[i,j])^2
    end
  end

subSol = sol(t_titer)

@inbounds for i = 1:length(subSol) #Loop through each virus time point
  sumsq += 1000000*(titer[i,1]-subSol[7,i])^2 #Virus ssqe, weighted
end

return sumsq
end

max_stars_repo_path: Julia Monte Carlo Fitting/RIGI/LossFunction.jl


text:
# generate 2D boxy canonical graph pattern

export generateCanonicalFG_Boxes2D!

_eps_w_int(::Type{T}) where {T <: Real} = eps(T)
_eps_w_int(::Type{T}) where {T <: Integer} = 1e-17

function _calcOdoBox( t4::T;
                      slew_x::Real=2/3,
                      length_x::Real = 15,
                      length_y::Real = length_x ) where {T <: Real}
  #

  # select direction
  t_ = t4 % 4
  if isapprox(t_, 0, atol=_eps_w_int(T))
    # drive along x
    return [length_x;0.0], :POSITIVE_X
  elseif isapprox(t_, 1, atol=_eps_w_int(T))
    # drive along y
    return [0.0;length_y], :POSITIVE_Y
  elseif isapprox(t_, 2, atol=_eps_w_int(T))
    # drive along -x
    return [-slew_x*length_x;0.0], :NEGATIVE_X
  elseif isapprox(t_, 3, atol=_eps_w_int(T))
    # drive along -y
    return [0.0;-length_y], :NEGATIVE_Y
  else
    error("non integer values of ")
  end
end

function driveLeg!( fg,
                    lastPose,
                    odo::AbstractVector{<:Real},
                    direction::Symbol;
                    graphinit::Bool=true,
                    Qd::AbstractArray{<:Real}=[1.0; 1.0],
                    factor::AbstractRelative=Point2Point2(MvNormal(odo, Qd)),
                    overridePPE=nothing,
                    postpose_cb::Function=(fg_,latestpose)->() )
  #
  # fg - factor graph object
  # lastpose
  # start -
  # v - displacement vector (e.g. [length_x; 0.0] for direction= :NORTH) ASSUMES ONE ELEMENT ALWAYS ZERO
  # direction - direction symbol (:NORTH,:SOUTH,:EAST,:WEST)

  # add end pose of the leg
  newPose = incrSuffix(lastPose)

  v_n = RoME._addPoseCanonical!(fg, lastPose, -1, factor, genLabel=newPose, srcType=RoME.Point2,
                                graphinit=false, variableTags=[:POSE;direction],
                                factorTags=[:ODOMETRY; direction],
                                overridePPE=overridePPE,
                                postpose_cb=postpose_cb  )
  #

  #
  # addFactor!(fg, [newPose;], PartialPrior(Normal(trueY, 0.1),(2,)), graphinit=graphinit)

  # return new pose name
  return v_n
end

# drive clockwise, x is North, y is East (NED convention).
# boxes start bottom left, spine of boxy helix is on x-axis
function driveOneBox!(fg;
                      lastPose=sortDFG(ls(fg, tags=[:POSE]))[end],
                      slew_x = 2/3,
                      length_x = 15,
                      length_y = length_x,
                      postpose_cb::Function=(fg_,latestpose)->()  )
  #
  # (g,lp) -> callbackSequenceScalar_ex4(g, lp, start_x, v )

  for leg in 0:3
    odo, direction = _calcOdoBox(leg, slew_x=slew_x, length_x=length_x, length_y=length_y)
    lastPose = driveLeg!(fg, lastPose, odo, direction,
                          postpose_cb=postpose_cb  ) |> getLabel
    #
  end

  #
  nothing
end

"""
    $SIGNATURES

Canonical graph generator for box shapes along x and y axes, with default behaviour slewing each new box along x-axes with 2/3 overlap on edges.

Notes
- Function is still experimental
- TODO currently `numposes` works in multiples of 4, for each corner of a new box -- i.e. 8 --> 2 boxes.
- Use `postpose_cb(fg, lastpose)` for additional behaviour after each pose variable has been added to the graph.

Related

[`generateCanonicalFG_ZeroPose`](@ref), [`generateCanonicalFG_Helix2DSlew!`](@ref), [`generateCanonicalFG_Hexagonal`](@ref)
"""
function generateCanonicalFG_Boxes2D!(numposes::Integer=16;
                                      graphinit::Bool=false,
                                      dfg::AbstractDFG = LightDFG{SolverParams}(solverParams=SolverParams(graphinit=graphinit)),
                                      useMsgLikelihoods::Bool=getSolverParams(dfg).useMsgLikelihoods,
                                      length_x::Real=15,
                                      length_y::Real=length_x,
                                      slew_x::Real=2/3,
                                      poseRegex::Regex=r"x\d+",
                                      refKey::Symbol=:simulated,
                                      Qd::Matrix{<:Real}=diagm( [0.1;0.1;0.05].^2 ),
                                      postpose_cb::Function=(fg_,latestpose)->()   )
  #

  # actually start adding nodes
  generateCanonicalFG_ZeroPose(dfg=dfg, varType=RoME.Point2, variableTags=[:POSE;], postpose_cb=postpose_cb)

  numboxes = ceil(Int, numposes/4)

  for tr in 0:(numboxes-1)
    # drive boxes
    driveOneBox!(dfg, slew_x=slew_x, length_x=length_x, postpose_cb=postpose_cb )
  end

  return dfg
end

#