text:
#=
Form factor + edge velocity calculation
=#

##################################################################
"""
$(SIGNATURES)

Find the circumcenter of a triangle.

Created from C source of Jonathan R Shewchuk <jrs@cs.cmu.edu>

Modified to return absolute coordinates.
"""
function tricircumcenter!(circumcenter,a,b,c)

    # Use coordinates relative to point `a' of the triangle.
    xba = b[1] - a[1]
    yba = b[2] - a[2]
    xca = c[1] - a[1]
    yca = c[2] - a[2]

    # Squares of lengths of the edges incident to `a'.
    balength = xba * xba + yba * yba
    calength = xca * xca + yca * yca

    # Calculate the denominator of the formulae.
    # if EXACT
    #    Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html
    #    to ensure a correctly signed (and reasonably accurate) result,
    #    avoiding any possibility of division by zero.
    #  denominator = 0.5 / orient2d((double*) b, (double*) c, (double*) a)

    # Take your chances with floating-point roundoff
    denominator = 0.5 / (xba * yca - yba * xca)

    # Calculate offset (from `a') of circumcenter.
    xcirca = (yca * balength - yba * calength) * denominator
    ycirca = (xba * calength - xca * balength) * denominator

# The result is returned both in terms of x-y coordinates and xi-eta
# coordinates, relative to the triangle's point `a' (that is, `a' is
# the origin of both coordinate systems).	 Hence, the x-y coordinates
# returned are NOT absolute; one must add the coordinates of `a' to
# find the absolute coordinates of the circumcircle.  However, this means
# that the result is frequently more accurate than would be possible if
# absolute coordinates were returned, due to limited floating-point
# precision.  In general, the circumradius can be computed much more
# accurately.

    circumcenter[1] = xcirca+a[1]
    circumcenter[2] = ycirca+a[2]

    return circumcenter
end

################################################
"""
$(SIGNATURES)

Calculate node volume  and voronoi surface contributions for cell.
"""

function cellfactors!(::Type{Edge1D},::Type{Cartesian1D},coord,cellnodes,icell,nodefac,edgefac)
    K=cellnodes[1,icell]
    L=cellnodes[2,icell]
    xK=coord[1,K]
    xL=coord[1,L]
    d=abs(xL[1]-xK[1])
    nodefac[1]=d/2
    nodefac[2]=d/2
    edgefac[1]=1/d
    nothing
end

function cellfactors!(::Type{Edge1D},::Type{<:Polar1D}, coord,cellnodes,icell,nodefac,edgefac)
    K=cellnodes[1,icell]
    L=cellnodes[2,icell]
    xK=coord[1,K]
    xL=coord[1,L]
    r0=xK[1]
    r1=xL[1]
    if r1<r0
        r0=xL[1]
        r1=xK[1]
    end
    rhalf=0.5*(r1+r0);
    # cpar[1]= π*(r1*r1-r0*r0);         # circular volume
    nodefac[1]= π*(rhalf*rhalf-r0*r0);   # circular volume between midline and boundary
    nodefac[2]= π*(r1*r1-rhalf*rhalf);   # circular volume between midline and boundary
    edgefac[1]= 2.0*π*rhalf/(r1-r0);     # circular surface / width
    nothing
end

function cellfactors!(::Type{Triangle2D},::Type{Cartesian2D},coord,cellnodes,icell,npar,epar)
    i1=cellnodes[1,icell]
    i2=cellnodes[2,icell]
    i3=cellnodes[3,icell]

    # Fill matrix of edge vectors
    V11= coord[1,i2]- coord[1,i1]
    V21= coord[2,i2]- coord[2,i1]

    V12= coord[1,i3]- coord[1,i1]
    V22= coord[2,i3]- coord[2,i1]

    V13= coord[1,i3]- coord[1,i2]
    V23= coord[2,i3]- coord[2,i2]

    # Compute determinant
    det=V11*V22 - V12*V21
    vol=0.5*det

    ivol = 1.0/vol

    # squares of edge lengths
    dd1=V13*V13+V23*V23 # l32
    dd2=V12*V12+V22*V22 # l31
    dd3=V11*V11+V21*V21 # l21

    # contributions to \sigma_kl/h_kl
    epar[1]= (dd2+dd3-dd1)*0.125*ivol
    epar[2]= (dd3+dd1-dd2)*0.125*ivol
    epar[3]= (dd1+dd2-dd3)*0.125*ivol

    # contributions to \omega_k
    npar[1]= (epar[3]*dd3+epar[2]*dd2)*0.25
    npar[2]= (epar[1]*dd1+epar[3]*dd3)*0.25
    npar[3]= (epar[2]*dd2+epar[1]*dd1)*0.25
    nothing
end

function cellfactors!(::Type{Tetrahedron3D},::Type{Cartesian3D},coord,cellnodes,icell,npar,epar)
    # Transferred from WIAS/pdelib, (c) J. Fuhrmann, H. Langmach, I. Schmelzer
    es=@MVector zeros(6)
    dd=@MVector zeros(6)
    i1=cellnodes[1,icell]
    i2=cellnodes[2,icell]
    i3=cellnodes[3,icell]
    i4=cellnodes[4,icell]

    #-> celledgenodes
    pp1=(1,1,1,2,2,3)
    pp2=(2,3,4,3,4,4)

    pi1=(5,6,5,1)
    pi2=(6,3,1,4)
    pi3=(4,2,3,2)

    po1=(2,1,2,6)
    po2=(1,4,6,3)
    po3=(3,5,4,5)
    pedge=(6,5,4,3,2,1)

    for i=1:6
        p1=cellnodes[pp1[i],icell]
        p2=cellnodes[pp2[i],icell]
        dx=coord[1,p1]-coord[1,p2]
        dy=coord[2,p1]-coord[2,p2]
        dz=coord[3,p1]-coord[3,p2]
        dd[i]=dx*dx+dy*dy+dz*dz
        es[i]=0.0
        epar[i]=0.0
    end

    x1=coord[1,i2]-coord[1,i1]
    y1=coord[2,i2]-coord[2,i1]
    z1=coord[3,i2]-coord[3,i1]

    x2=coord[1,i3]-coord[1,i1]
    y2=coord[2,i3]-coord[2,i1]
    z2=coord[3,i3]-coord[3,i1]

    x3=coord[1,i4]-coord[1,i1]
    y3=coord[2,i4]-coord[2,i1]
    z3=coord[3,i4]-coord[3,i1]

    det= (x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1))

    if det <0
        det=-det
    end
    vol=det/6
    vv= 96*6*vol

    for  i=1:4
        npar[i]=0.0
        i1=pi1[i]
        i2=pi2[i]
        i3=pi3[i]

        h1=dd[i1]*(dd[i2]+dd[i3]-dd[i1])
        h2=dd[i2]*(dd[i3]+dd[i1]-dd[i2])
        h3=dd[i3]*(dd[i1]+dd[i2]-dd[i3])

        df=h1+h2+h3
        vf=(h1*dd[po1[i]] + h2*dd[po2[i]] + h3*dd[po3[i]]-2*dd[i1]*dd[i2]*dd[i3]) / (vv*df)
        es[i1] += h1*vf
        es[i2] += h2*vf
        es[i3] += h3*vf
    end

    for i=1:6
        epar[pedge[i]] = 6*es[i]/dd[i]
        npar[pp1[i]] += es[i]
        npar[pp2[i]] += es[i]
    end
    nothing
end

function cellfactors!(::Type{Triangle2D},::Type{<:Cylindrical2D},coord,cellnodes,icell,npar,epar)
    function area2d(coord1, coord2, coord3)
        V11= coord2[1]- coord1[1]
        V21= coord2[2]- coord1[2]

        V12= coord3[1]- coord1[1]
        V22= coord3[2]- coord1[2]

        V13= coord3[1]- coord2[1]
        V23= coord3[2]- coord2[2]

        # Compute determinant
        det=V11*V22 - V12*V21
        area=abs(0.5*det)
    end

    i1=cellnodes[1,icell]
    i2=cellnodes[2,icell]
    i3=cellnodes[3,icell]

    # Fill matrix of edge vectors
    V11= coord[1,i2]- coord[1,i1]
    V21= coord[2,i2]- coord[2,i1]

    V12= coord[1,i3]- coord[1,i1]
    V22= coord[2,i3]- coord[2,i1]

    V13= coord[1,i3]- coord[1,i2]
    V23= coord[2,i3]- coord[2,i2]

    # Compute determinant
    det=V11*V22 - V12*V21
    area=0.5*det

    # Integrate R over triangle (via quadrature rule)
    vol=2.0*π*area*(coord[1,i1]+coord[1,i2]+coord[1,i3])/3.0

    # squares of edge lengths
    dd1=V13*V13+V23*V23 # l32
    dd2=V12*V12+V22*V22 # l31
    dd3=V11*V11+V21*V21 # l21

    emid23=[0.5*(coord[1,i3]+coord[1,i2]),
            0.5*(coord[2,i3]+coord[2,i2])]

    emid13=[0.5*(coord[1,i1]+coord[1,i3]),
            0.5*(coord[2,i1]+coord[2,i3])]

    emid12=[0.5*(coord[1,i1]+coord[1,i2]),
            0.5*(coord[2,i1]+coord[2,i2])]

    cc=@MVector zeros(2)
    tricircumcenter!(cc,coord[:,i1],coord[:,i2],coord[:,i3])

    r(p)=p[1]
    z(p)=p[2]
    sq(x)=x*x
    epar[1]= π*(r(cc)+r(emid23))*sqrt(sq(r(cc)-r(emid23))+sq(z(cc)-z(emid23)))/sqrt(dd1);
    epar[2]= π*(r(cc)+r(emid13))*sqrt(sq(r(cc)-r(emid13))+sq(z(cc)-z(emid13)))/sqrt(dd2);
    epar[3]= π*(r(cc)+r(emid12))*sqrt(sq(r(cc)-r(emid12))+sq(z(cc)-z(emid12)))/sqrt(dd3);

    rintegrate(coord1, coord2, coord3)=2.0*π*area2d(coord1,coord2,coord3)*(coord1[1]+coord2[1]+coord3[1])/3.0
    npar[1]=rintegrate(coord[:,i1],cc,emid13)+rintegrate(coord[:,i1],cc,emid12)
    npar[2]=rintegrate(coord[:,i2],cc,emid23)+rintegrate(coord[:,i2],cc,emid12)
    npar[3]=rintegrate(coord[:,i3],cc,emid13)+rintegrate(coord[:,i3],cc,emid23)
    nothing
end

################################################
"""
$(SIGNATURES)

Calculate node volume  contributions for boundary face.
"""

function bfacefactors!(::Type{Vertex0D},::Type{Cartesian1D},coord,bfacenodes,ibface,nodefac)
    nodefac[1]=1.0
    nothing
end

function bfacefactors!(::Type{Vertex0D},::Type{<:Polar1D},coord,bfacenodes,ibface,nodefac)
    inode=bfacenodes[1,ibface]
    r=coord[1,inode]
    nodefac[1]=2*π*r
    nothing
end

function bfacefactors!(::Type{Edge1D},::Type{<:Cartesian2D},coord,bfacenodes,ibface,nodefac)
    i1=bfacenodes[1,ibface]
    i2=bfacenodes[2,ibface]
    dx=coord[1,i1]-coord[1,i2]
    dy=coord[2,i1]-coord[2,i2]
    d=0.5*sqrt(dx*dx+dy*dy)
    nodefac[1]=d
    nodefac[2]=d
    nothing
end

function bfacefactors!(::Type{Edge1D},::Type{<:Cylindrical2D},coord,bfacenodes,ibface,nodefac)
    i1=bfacenodes[1,ibface]
    i2=bfacenodes[2,ibface]
    r1=coord[1,i1]
    r2=coord[1,i2]
    z1=coord[2,i1]
    z2=coord[2,i2]
    dr=r1-r2
    rmid=(r1+r2)/2
    dz=z1-z2
    l=sqrt(dr*dr+dz*dz)
    nodefac[1]=π*(r1+rmid)*l/2
    nodefac[2]=π*(r2+rmid)*l/2
    nothing
end

function bfacefactors!(::Type{Triangle2D},::Type{<:Cartesian3D},coord,bfacenodes,ibface,npar)

    # Transferred from WIAS/pdelib, (c) J. Fuhrmann, H. Langmach, I. Schmelzer
    i1=bfacenodes[1,ibface]
    i2=bfacenodes[2,ibface]
    i3=bfacenodes[3,ibface]
    dd=@MVector zeros(3)
    epar=@MVector zeros(3)
    for j=1:3
        d= coord[j,i1] - coord[j,i3];  dd[2]+= d*d;
        d= coord[j,i2] - coord[j,i1];  dd[3]+= d*d;
        d= coord[j,i3] - coord[j,i2];  dd[1]+= d*d;
    end

    # Kanten-Flaechenanteile (ohne Abschneiden); epar als Hilfsfeld benutzt
    epar[1]= (dd[2]+dd[3]-dd[1])*dd[1];
    epar[2]= (dd[3]+dd[1]-dd[2])*dd[2];
    epar[3]= (dd[1]+dd[2]-dd[3])*dd[3];
    vol= sqrt(epar[1]+epar[2]+epar[3])*0.25;

    d = 1.0/(8*vol);

    # Knoten-Flaechenanteile (ohne Abschneiden)
    npar[1]= (epar[3]+epar[2])*d*0.25;
    npar[2]= (epar[1]+epar[3])*d*0.25;
    npar[3]= (epar[2]+epar[1])*d*0.25;

    # Kantengewichte
    # epar[1]= epar[1]*d/dd[1];
    # epar[2]= epar[2]*d/dd[2];
    # epar[3]= epar[3]*d/dd[3];
    nothing
end

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

#
# TODO: this should be generalized for more quadrules
#
function integrate(coordl,coordr,hnormal,velofunc)
    wl=1.0/6.0
    wm=2.0/3.0
    wr=1.0/6.0
    coordm=0.5*(coordl+coordr)
    (vxl,vyl)=velofunc(coordl[1],coordl[2])
    (vxm,vym)=velofunc(coordm[1],coordm[2])
    (vxr,vyr)=velofunc(coordr[1],coordr[2])
    return (wl*vxl +wm*vxm+wr*vxr)*hnormal[1] + (wl*vyl +wm*vym+wr*vyr)*hnormal[2]
end

"""
$(SIGNATURES)

Project velocity onto grid edges,
"""
function edgevelocities(grid,velofunc)

    @assert dim_space(grid)<3

    cn=grid[CellNodes]
    ec=grid[EdgeCells]
    en=grid[EdgeNodes]
    coord=grid[Coordinates]

    velovec=zeros(Float64,num_edges(grid))
    if dim_space(grid)==1
        for iedge=1:num_edges(grid)
            K=en[1,iedge]
            L=en[2,iedge]
            elen=coord[1,L]-coord[1,K]
            vx,vy=velofunc((coord[1,K]+coord[1,L])/2,0)
            velovec[iedge]=-elen*vx
        end
    else
        for iedge=1:num_edges(grid)
            K=en[1,iedge]
            L=en[2,iedge]
            p1=@MVector zeros(2)
            p2=@MVector zeros(2)
            tricircumcenter!(p1,
                             coord[:,cn[1,ec[1,iedge]]],
                             coord[:,cn[2,ec[1,iedge]]],
                             coord[:,cn[3,ec[1,iedge]]])
            if ec[2,iedge]>0
                tricircumcenter!(p2,
                                 coord[:,cn[1,ec[2,iedge]]],
                                 coord[:,cn[2,ec[2,iedge]]],
                                 coord[:,cn[3,ec[2,iedge]]])
            else
                p2.=0.5*(coord[:,K]+coord[:,L])
            end
            hnormal=coord[:,K]-coord[:,L]
            velovec[iedge]=integrate(p1,p2,hnormal,velofunc)
        end
    end
    return velovec
end

function bfacevelocities(grid,velofunc)
    @assert dim_space(grid)<3
    bfacenodes=grid[BFaceNodes]
    coord=grid[Coordinates]
    bfacecells=grid[BFaceCells]
    bfacenormals=grid[BFaceNormals]
    bfr=grid[BFaceRegions]
    velovec=zeros(Float64,2,num_bfaces(grid))
    if dim_space(grid)==1
        for ibface=1:num_bfaces(grid)
            vx,vy=velofunc(coord[1,bfacenodes[1,ibface]])
            velovec[ibface]=vx*bfacenormals[1,ibface]
        end
    else
        for ibface=1:num_bfaces(grid)
            p1=coord[:,bfacenodes[1,ibface]]
            p2=coord[:,bfacenodes[2,ibface]]
            pm=0.5*(p1+p2)
            velovec[1,ibface]=integrate(p1,pm,bfacenormals[:,ibface],velofunc)
            velovec[2,ibface]=integrate(pm,p2,bfacenormals[:,ibface],velofunc)
        end
    end
    return velovec
end


text:
module GnuplotSimple

using Gnuplot
import Gnuplot: PlotElement, DatasetText

export plot

"""
`plot(x,Y; kwargs...)`

A simple plot recipe for `Gnuplot.jl`.

Required: vector-like `x` and vector or matrix-like `Y`.

All keyword arguments are optional:

`title`, `xlab`, `ylab`, `key` should be strings.

`grid`, `scatter`, `xlog`, `ylog`, `mono` are boolean.

`xr` and `yr` are like `[0,1]`. Use `NaN` for unbounded.

`labels`, `linecolor`, `linetype` should be vectors of strings (at least as many as columns in `Y`).

`linewidth` should be a vector of numbers.

```julia
using Gnuplot, GnuplotSimple
x = 0:0.05:2
y = x.^(1:3)'
@gp plot(x,y; title="My Title",xlab="My x",ylab="My y",xr=[0,1],labels=["y1","y2","y3"])
```
"""
function plot(x,Y; grid::Bool=true, scatter::Bool=false, title=nothing, xlab=nothing, ylab=nothing, key=nothing, labels=nothing, xr=[NaN,NaN], yr=[NaN,NaN],
        linetype=nothing, linewidth=nothing, linecolor=nothing, xlog::Bool=false, ylog::Bool=false, mono::Bool=false)
    # assertions and input processing
    @assert !isempty(x)
    @assert !isempty(Y)
    Lx = length(x)
    Ly, Ny = size(Y)
    dimsx = length(size(x))
    dimsY = length(size(Y))
    @assert dimsx == 1 "x must be 1D"
    @assert dimsY <= 2 "Y must be 1D or 2D"
    @assert Lx == Ly "x and Y (or its columns) must have equal number of observations"

    ptype = ifelse(scatter,"p","l")

    # plot setup
    out = Vector{Gnuplot.PlotElement}()
    grid && push!(out,PlotElement(cmds=["set grid"]))
    mono && push!(out,PlotElement(cmds=["set monochrome"]))
    !isnothing(title) && push!(out, PlotElement(;tit=title))
    !isnothing(xlab) && push!(out, PlotElement(;xlabel=xlab))
    !isnothing(ylab) && push!(out, PlotElement(;ylabel=ylab))
    !isnothing(key) && push!(out, PlotElement(;key=key))
    !isnothing(xr) && push!(out, PlotElement(;xrange=xr))
    !isnothing(yr) && push!(out, PlotElement(;yrange=yr))
    push!(out, PlotElement(;xlog=xlog,ylog=ylog))

    # do plotting
    for i in 1:Ny
        plotstr = "w $ptype"
        if !isnothing(labels)
            if length(labels) >= i
                thislab = labels[i]
                plotstr = string(plotstr," t '$thislab'")
            end
        end
        if !isnothing(linetype)
            if length(linetype) >= i
                thislt = linetype[i]
                plotstr = string(plotstr," $thislt")
            end
        end
        if !isnothing(linewidth)
            if length(linewidth) >= i
                thislw = linewidth[i]
                plotstr = string(plotstr," lw $thislw")
            end
        end
        if !isnothing(linecolor)
            if length(linecolor) >= i
                thiscolor = linecolor[i]
                plotstr = string(plotstr," lc rgb '$thiscolor'")
            end
        end

        push!(out, PlotElement(data=DatasetText(x, Y[:,i]), plot=plotstr))
    end
    return out
end

end # module

ext: jl


text:
## ExoplanetsSysSim/src/window_function.jl
## (c) 2018 Darin Ragozzine

# Gather and prepare the window function data
module WindowFunction

import ..cdpp_durations
export setup_window_function, get_window_function_data, get_window_function_id, eval_window_function

#using DataArrays
using DataFrames
#using CSV
#using JLD2
using FileIO
using ExoplanetsSysSim
using ExoplanetsSysSim.SimulationParameters

# Object to hold window function data
struct window_function_data
  window_func_array::Array{Float64,3}      #  Value of window function (window_function_id, duration_id, period_id).  Maybe rename to wf_value or data?
  wf_durations_in_hrs::Array{Float64,1}    # Boundaries for duration bins in window_func_array.  Maybe rename to durations?
  wf_periods_in_days::Array{Float64,1}     # Boundaries for periods bins in window_func_array.  Maybe rename to periods?
  sorted_quarter_strings::Array{Int64,1}   # TODO OPT: Is there a reason to keep this?  Maybe rename to quarter_strings?
  allsortedkepids::Array{Int64,1}          # value is Kepler ID.  Index is same as index to window_function_id_arr
  window_function_id_arr::Array{Int64,1}   # value is index for window_func_array.  Index is same as index to allsortedkepids
  default_wf_id::Int64                     # Index corresponding to the default window function
  already_warned::Array{Bool,1}            # Wether we've already thrown a warning about this kepid
end

function window_function_data()
  window_function_data( Array{Float64,3}(undef,0,0,0), Array{Float64,1}(undef,0),Array{Float64,1}(undef,0), Array{Int64,1}(undef,0),Array{Int64,1}(undef,0),Array{Int64,1}(undef,0), 0, falses(0) )
end

win_func_data = window_function_data()

function setup(sim_param::SimParam; force_reread::Bool = false)
  global win_func_data
  if haskey(sim_param,"read_window_function") && !force_reread
     return win_func_data
  end
  window_function_filename = convert(String,joinpath(dirname(pathof(ExoplanetsSysSim)),
      "..", "data", convert(String,get(sim_param,"window_function","DR25topwinfuncs.jld2")) ) )
  setup(window_function_filename)
  add_param_fixed(sim_param,"read_window_function",true)
  @assert( size(win_func_data.window_func_array,2) == length(win_func_data.wf_durations_in_hrs) )
  @assert( size(win_func_data.window_func_array,3) == length(win_func_data.wf_periods_in_days) )
  @assert( size(win_func_data.window_func_array,1) >= maximum(win_func_data.window_function_id_arr) )
  @assert( size(win_func_data.window_func_array,1) >= win_func_data.default_wf_id )
  return win_func_data
end

function setup(filename::String)
# Reads in the window function data collected from the Kepler Completeness Products
# see Darin Ragozzine's get/cleanDR25winfuncs.jl

  if occursin(r".jld2$",filename)
    try
      wfdata = load(filename)
      window_func_array = wfdata["window_func_array"]
      wf_durations_in_hrs = wfdata["wf_durations_in_hrs"]  # TODO OPT DETAIL: Should we convert units to days here?
      wf_periods_in_days = wfdata["wf_periods_in_days"]
      sorted_quarter_strings = wfdata["sorted_quarter_strings"]
      allsortedkepids = wfdata["allsortedkepids"]
      window_function_id_arr = wfdata["window_function_id_arr"]
      already_warned = falses(length(allsortedkepids))
      global win_func_data = window_function_data(window_func_array, wf_durations_in_hrs, wf_periods_in_days, sorted_quarter_strings,
                                                allsortedkepids, window_function_id_arr, maximum(window_function_id_arr), already_warned )

    catch
      error(string("# Failed to read window function data > ", filename," < in jld2 format."))
    end
  end

  return win_func_data
end

setup_window_function(sim_param::SimParam; force_reread::Bool = false) = setup(sim_param, force_reread=force_reread)
setup_window_function(filename::String; force_reread::Bool = false) = setup(filename, force_reread=force_reread)

function get_window_function_data()::window_function_data
   #global win_func_data
   return win_func_data
end

function get_window_function_id(kepid::Int64; use_default_for_unknown::Bool = true)::Int64
  # takes the quarter string from the stellar catalog and determines the window function id
  # from DR25topwinfuncs.jld2 made by Darin Ragozzine's cleanDR25winfuncs.jl script.
  no_win_func_available::Int64 = -1        # hardcoding this in, should match convention in window function input file

  idx = searchsortedfirst(win_func_data.allsortedkepids,kepid) # all Kepler kepids are in allsortedkepids
  wf_id = win_func_data.window_function_id_arr[idx]

  if wf_id == no_win_func_available && use_default_for_unknown
    # if a target is observed for less than 4 quarters, then it won't have a corresponding
    # window function in this list, so throw a warning and use the last window_function_id
    # which corresponds to an "averaged" window function
    if !win_func_data.already_warned[idx]
       win_func_data.already_warned[idx] = true
       if sum(win_func_data.already_warned) < 20
          @warn "Window function data is not avaialble for kepid $kepid, using default."
       end
    end
    wf_id = win_func_data.default_wf_id
  end
  # TODO SCI DETAIL IMPORTANT? This does not include TPS timeouts or MESthresholds (see DR25 Completeness Products)

  return wf_id
end

function calc_period_idx(P::Float64)::Int64
  @assert(P>zero(P))
  idx = searchsortedlast(win_func_data.wf_periods_in_days,P)
  if idx == 0
     return 1
  elseif idx<length(win_func_data.wf_periods_in_days)
     if P-win_func_data.wf_periods_in_days[idx]>win_func_data.wf_periods_in_days[idx+1]-P
        idx += 1
     end
  end
  return idx
end

function calc_duration_idx(D::Float64)::Int64
  # NOTE: Currently assumes we left wf data in hours, so deal with that conversion here
  @assert(D>=zero(D)) ##### Make sure this function is still doing the right thing if D = 0!
  hours_in_day = 24
  idx = searchsortedlast(win_func_data.wf_durations_in_hrs,D*hours_in_day)
  if idx == 0
     return 1
  elseif idx<length(win_func_data.wf_durations_in_hrs)
     if D*hours_in_day-win_func_data.wf_durations_in_hrs[idx]>win_func_data.wf_durations_in_hrs[idx+1]-D*hours_in_day
        idx += 1
     end
  end
  return idx
end

function eval_window_function(wf_idx::Int64=-1; Duration::Float64=0., Period::Float64=0.)::Float64
  D_idx = calc_duration_idx(Duration)
  P_idx = calc_period_idx(Period)
  wf = eval_window_function(wf_idx,D_idx,P_idx)
  # TODO SCI DETAIL: Improve way deal with missing wf values for some durations. Interpolate?
  while wf<=zero(wf) && D_idx<length(win_func_data.wf_durations_in_hrs)
     D_idx += 1
     wf = eval_window_function(wf_idx,D_idx,P_idx)
  end
  return wf
end

function eval_window_function(wf_idx::Int64, D_idx::Int64, P_idx::Int64)::Float64
   global win_func_data
   #@assert(1<=wf_idx<maximum(win_func_data.window_function_id_arr))
   #@assert(1<=P_idx<=length(win_func_data.wf_periods_in_days))
   #@assert(1<=D_idx<=length(win_func_data.wf_durations_in_hrs))
   return win_func_data.window_func_array[wf_idx,D_idx,P_idx]
end

#Object for storing data necessary for OSD_interpolator
struct OSD_data{T1<:Real, T2<:Real}
    allosds::Array{T1,3}
    kepids::Array{Int64,1}
    #periods_length::Int64
    #durations_length::Int64
    grid::Array{Array{T2,1},1}

    function OSD_data(data::AbstractArray{T1,3}, kepids::AbstractArray{Int64,1}, durations::AbstractArray{T2,1}, periods::AbstractArray{T2,1} ) where {T1<:Real, T2<:Real}
        @assert(size(data,1)==length(kepids))
        @assert(size(data,2)==length(durations))
        @assert(size(data,3)==length(periods))
        @assert issorted(kepids)
        @assert issorted(durations)
        @assert issorted(periods)
        new{T1,T2}(data, kepids, [durations,periods])
    end
end

# Only point of this version is to provide a drop in replacement for Keir's code
function OSD_data(allosds::AbstractArray{T1,3}, kepids::AbstractArray{T3,1}, periods_length::Int64, durations_length::Int64, grid::Array{Array{T2,1},1}) where {T1<:Real, T2<:Real, T3<:Real}
    @assert length(grid) == 2
    @assert durations_length == length(grid[1])
    @assert periods_length == length(grid[2])
    if eltype(kepids) != Int64
        kepids = round.(Int64,kepids)
    end
    @assert grid[1][1] == 1.5  # Checking that durations were passed as hours, since that's what Keir's assumed
    return OSD_data(allosds, kepids, grid[1] ./ 24.0 , grid[2])  # Convert durations to days, since that's units in rest of SysSim
end

num_stars(osd::OSD_data) = size(osd.allosds,1)
num_durations(osd::OSD_data) = size(osd.allosds,2)
num_periods(osd::OSD_data) = size(osd.allosds,3)

function setup_OSD(sim_param::SimParam; force_reread::Bool = false)			#reads in 3D table of OSD values and sets up global variables to be used in interpolation
  global OSD_setup
  if haskey(sim_param,"read_OSD_function") && !force_reread
     return OSD_setup
  end
  #OSD_file = load(joinpath(Pkg.dir(), "ExoplanetsSysSim", "data", convert(String,get(sim_param,"osd_file","allosds.jld"))))
  #OSD_file = load(joinpath(Pkg.dir(), "ExoplanetsSysSim", "data", convert(String,get(sim_param,"osd_file","allosds.jld"))))
  #OSD_file = load(joinpath(dirname(pathof(ExoplanetsSysSim)),"data",convert(String,get(sim_param,"osd_file","allosds.jld"))))
  #OSD_file = load(joinpath(dirname(pathof(ExoplanetsSysSim)),"..","data",convert(String,get(sim_param,"osd_file","dr25fgk_relaxcut_osds.jld"))))
  OSD_file = load(joinpath(dirname(pathof(ExoplanetsSysSim)),"..","data",convert(String,get(sim_param,"osd_file","dr25fgk_small_osds.jld2"))))
  allosds = OSD_file["allosds"]			#table of OSDs with dimensions: kepids,durations,periods
  periods = OSD_file["periods"][1,:]		#1000 period values corresponding to OSD values in the third dimension of the allosds table
  kepids = OSD_file["kepids"]			#kepids corresponding to OSD values in the first dimension of the allosds table
  OSD_file = 0 # unload OSD file to save memory
  #durations = [1.5,2.,2.5,3.,3.5,4.5,5.,6.,7.5,9.,10.5,12.,12.5,15.] #14 durations corresponding to OSD values in the first dimension of theh allosds table
  periods_length = length(allosds[1,1,:])
  durations_length = length(allosds[1,:,1])
  @assert length(cdpp_durations) == durations_length
  grid = Array{Float64,1}[]			#grid used in OSD_interpolator
  push!(grid, cdpp_durations)
  push!(grid, periods)
  #global compareNoise = Float64[]		#testing variable used to make sure OSD_interpolator is producing reasonable snrs
  OSD_setup = OSD_data(allosds, kepids, periods_length, durations_length, grid)
  allosds = 0 # unload OSD table to save memory
  add_param_fixed(sim_param,"read_OSD_function",true)
  return OSD_setup
end

setup_OSD_interp(sim_param::SimParam; force_reread::Bool = false) = setup_OSD(sim_param, force_reread=force_reread)

function find_index_lower_bounding_point(grid::AbstractArray{T1,1}, x::T2; verbose::Bool = false) where {T1<:Real, T2<:Real}
    if verbose
        @assert issorted(grid)
        @assert length(grid)>=2
        println("# ", grid[1], " <= ", x, " <= ", grid[end])
        @assert grid[1] <= x <= grid[end]
    end
    idx = searchsortedlast(grid,x)
    if idx == 0
        idx = 1
    #elseif idx >= length(grid) # should never happen
    #    idx = length(grid)-1
    end
    return idx
end

function interp_OSD_from_table(kepid::Int64, period::T2, duration::T3; verbose::Bool = false) where {T2<:Real, T3<:Real}
  @assert eltype(OSD_setup.kepids) == Int64      # otherwise would need something like kepid = convert(eltype(OSD_setup.kepids),kepid)
  kepid_idx = searchsortedfirst(OSD_setup.kepids, kepid)
  if (kepid_idx > num_stars(OSD_setup)) || (OSD_setup.kepids[kepid_idx] != kepid) # if we don't find the kepid in allosds.jld, then we make a random one
     kepid_idx = rand(1:size(OSD_setup.allosds,1))
     if verbose
            println("# picked random kepid = ", OSD_setup.kepids[kepid_idx])
     end
  end
  idx_duration = find_index_lower_bounding_point(OSD_setup.grid[1], duration)
  idx_period   = find_index_lower_bounding_point(OSD_setup.grid[2], period)
  #z = view(OSD_setup.allosds,kepid_idx,idx_duration:(idx_duration+1),idx_period:(idx_period+1))     # use correct kepid index to extract 2D table from 3D OSD table

  #= value = z[1,1] * w_dur * w_per +
          z[2,1] * (1-w_dur) * w_per +
          z[1,2] * w_dur * (1-w_per) +
          z[2,2] * (1-w_dur) * (1-w_per)  =#

  if idx_duration < length(OSD_setup.grid[1]) && idx_period < length(OSD_setup.grid[2])
      w_dur = (duration-OSD_setup.grid[1][idx_duration]) / (OSD_setup.grid[1][idx_duration+1]-OSD_setup.grid[1][idx_duration])
      w_per = (period  -OSD_setup.grid[2][idx_period])   / (OSD_setup.grid[2][idx_period+1]  -OSD_setup.grid[2][idx_period])

      value = OSD_setup.allosds[kepid_idx,idx_duration,  idx_period  ] * w_dur * w_per +
          OSD_setup.allosds[kepid_idx,idx_duration+1,idx_period  ] * (1-w_dur) * w_per +
          OSD_setup.allosds[kepid_idx,idx_duration,  idx_period+1] * w_dur * (1-w_per) +
          OSD_setup.allosds[kepid_idx,idx_duration+1,idx_period+1] * (1-w_dur) * (1-w_per)
  elseif idx_period < length(OSD_setup.grid[2])
      w_per = (period  -OSD_setup.grid[2][idx_period])   / (OSD_setup.grid[2][idx_period+1]  -OSD_setup.grid[2][idx_period])

      value = (OSD_setup.allosds[kepid_idx,idx_duration,  idx_period+1] - OSD_setup.allosds[kepid_idx,idx_duration,  idx_period  ]) * w_per +
           OSD_setup.allosds[kepid_idx,idx_duration,  idx_period  ]
  elseif idx_duration < length(OSD_setup.grid[1])
      w_dur = (duration-OSD_setup.grid[1][idx_duration]) / (OSD_setup.grid[1][idx_duration+1]-OSD_setup.grid[1][idx_duration])

      value = (OSD_setup.allosds[kepid_idx,idx_duration+1,idx_period  ] - OSD_setup.allosds[kepid_idx,idx_duration,  idx_period  ]) * w_dur +
          OSD_setup.allosds[kepid_idx,idx_duration,idx_period  ]
  else
      value = OSD_setup.allosds[kepid_idx,idx_duration,  idx_period  ]
  end
end

# function interp_OSD_from_table(kepid::Int64, period::Real, duration::Real)
#   kepid = convert(Float64,kepid)
#   meskep = OSD_setup.kepids			#we need to find the index that this planet's kepid corresponds to in allosds.jld
#   kepid_index = findfirst(meskep, kepid)
#   if kepid_index == 0
#      kepid_index = rand(1:88807)		#if we don't find the kepid in allosds.jld, then we make a random one
#   end
#   olOSD = OSD_setup.allosds[kepid_index,:,:]    #use correct kepid index to extract 2D table from 3D OSD table
#   # olOSD = convert(Array{Float64,2},olOSD)
#   @time lint = Lininterp(olOSD, OSD_setup.grid)	#sets up linear interpolation object
#   osd = ApproXD.eval2D(lint, [duration*24,period])[1]	#interpolates osd
#   return osd
# end

# function cdpp_vs_osd(ratio::Float64, cuantos::Int64)
# #testing function that takes ratios of cdpp_snr/osd_snr and plots a histogram to make sure the results are reasonable.
#   global compareNoise
#   push!(compareNoise,ratio)
#   if length(compareNoise) == cuantos
#     PyPlot.plt[:hist](compareNoise,100)
#     println("MES median: ",median(compareNoise)," MES mean: ",mean(compareNoise), " Standard deviation: ",std(compareNoise))
#     cuantos = 100000000
#   end
#   return cuantos
# end

end  # module WindowFunction


text:
## Model Generators for 2D Gaussian PSFs

function genmodel!(s::StateFlatBg, sz::Int, psf::MicroscopePSFs.PSF,
    model::Array{Float32,2})

    for ii=1:sz^2
        model[ii]=s.bg;
    end

    for jj = 1:sz, ii = 1:sz,nn = 1:s.n
                model[ii+sz*(jj-1)] += s.photons[nn] *
                MicroscopePSFs.pdf(psf,(jj,ii),(s.x[nn],s.y[nn]))
    end
end


text:
#=
    Unit tests for EodData external utility functions
=#
using Base.Test
using EodData
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

dt = DateTime("20140101", "yyyymmdd")
d = Date("20140101", "yyyymmdd")

# set_date_string() with a ::DateTime
@test set_date_string(dt) == "20140101"
@test set_date_string(dt) != 20140101
@test set_date_string(dt) != dt
@test set_date_string(dt) != d
@test typeof(set_date_string(dt)) == typeof("20140101")

# set_date_string() with a ::Date
@test set_date_string(d) == "20140101"
@test set_date_string(d) != 20140101
@test set_date_string(d) != dt
@test set_date_string(d) != d
@test typeof(set_date_string(d)) == typeof("20140101")


text:
"""
    modularity(g, c)

Return a value representing Newman's modularity `Q` for the undirected graph
`g` given the partitioning vector `c`.
"""
function modularity end
@traitfn function modularity(g::::(!IsDirected), c::Vector)
    m = 2 * ne(g)
    m == 0 && return 0.
    nc = maximum(c)
    a = zeros(Int, nc)
    Q = 0
    for u in vertices(g)
        for v in neighbors(g, u)
            if u <= v
                c1 = c[u]
                c2 = c[v]
                if c1 == c2
                    Q += 2
                end
                a[c1] += 1
                a[c2] += 1
            end
        end
    end
    Q = Q * m
    @inbounds for i = 1:nc
        Q -= a[i] * a[i]
    end
    return Q / m / m
end


text:
# Hash table acceleration index
struct UniqueHashIndex{D <: Dict} <: AbstractUniqueIndex
    dict::D
end

function UniqueHashIndex(a::AbstractArray)
	dict = Dict{eltype(a), SingleVector{eltype(keys(a))}}()

    @inbounds for i in keys(a)
        value = a[i]
        index = Base.ht_keyindex2!(dict, value)
	    if index > 0 # `value` found in `dict`
	        error("Input not unique") # TODO Use appropriate Exception
	    else # `value` is ready to be inserted into `dict` at slot `-index`
	        @inbounds Base._setindex!(dict, SingleVector(i), value, -index)
	    end
    end
    return UniqueHashIndex(dict)
end

Base.summary(::UniqueHashIndex) = "UniqueHashIndex"

# Accelerations
Base.in(x, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}) = haskey(a.index.dict, x)

function Base.count(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
    index = Base.ht_keyindex(a.index.dict, f.x)
    if index < 0
        return 0
    else
        return 1
    end
end

function Base.findall(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
	index = Base.ht_keyindex(a.index.dict, f.x)
	if index < 0
		return MaybeVector{eltype(keys(a))}()
	else
		return MaybeVector(@inbounds a.index.dict.vals[index][])
	end
end

function Base.filter(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
	index = Base.ht_keyindex(a.index.dict, f.x)
	if index < 0
		return MaybeVector{eltype(a)}()
	else
		return MaybeVector{eltype(a)}(f.x)
	end
end

function SplitApplyCombine.group2(a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}, b::AbstractArray)
    return Dict((key, SingleVector(@inbounds b[inds[]])) for (key, inds) in a.index.dict)
end

function SplitApplyCombine.groupreduce(::typeof(identity), f, op, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}; kw...)
	return Dict((k, mapreduce(i -> f(@inbounds a[i]), op, v; kw...)) for (k, v) in a.index.dict)
end

function SplitApplyCombine._groupinds(a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
	return a.index.dict
end

function SplitApplyCombine._innerjoin!(out, left::AbstractArray, right::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}, v::AbstractArray, ::typeof(isequal))
    @boundscheck if (axes(l)..., axes(r)...) != axes(v)
        throw(DimensionMismatch("innerjoin arrays do not have matching dimensions"))
    end

    dict = right.index.dict

    @inbounds for i ∈ keys(left)
        dict_indUniqueex = Base.ht_keyindex(dict, left(i_l))
        if dict_index > 0 # -1 if key not found
            i_r = dict.vals[dict_index][]
            push!(out, v[Tuple(i_l)..., Tuple(i_r)...])
        end
    end

    return out
end


text:
using YaoExtensions, Yao, YaoSym

function sym_variational_circuit(nbit, nlayer; pairs = pair_ring(nbit), entangler = cnot)
    circuit = chain(nbit)
    k = 0
    genθ() = (k += 1; Basic(Symbol(:θ, k)))

    ent = chain(nbit, entangler(nbit, i, j) for (i, j) in pairs)
    has_param = nparameters(ent) != 0
    for i in 1:(nlayer+1)
        i != 1 && push!(circuit, has_param ? deepcopy(ent) : ent)
        r = chain(nbit)
        for j in 1:nbit
            ci =
                i == 1 ? (i == nlayer + 1 ? Rx(genθ()) : chain(Rx(genθ()), Rz(genθ()))) :
                (
                    i == nlayer + 1 ? chain(Rz(genθ()), Rx(genθ())) :
                    chain(Rz(genθ()), Rx(genθ()), Rz(genθ()))
                )
            push!(r, put(nbit, j => ci))
        end
        push!(circuit, r)
    end
    circuit
end

circ = sym_variational_circuit(2, 1)

reg = ArrayReg(Basic, bit"00") => circ
op = put(2, 2 => Z)
ex = expect'(op, reg)[2]
params = randn(nparameters(circ))
assign = Dict(zip(parameters(circ), params))
ex2 = map(x -> ComplexF64(subs(x, assign...)), ex)
circ2 = variational_circuit(2, 1, pair_ring(2); mode = :Merged)
dispatch!(circ2, params)
real(ex2) ≈ real(expect'(op, ArrayReg(bit"00") => circ2)[2])


text:
#=
   IMPORTANT

   Julia v0.4.0 does not provide elementary functions that always round to nearest.
   For example:  exp( 0.05489354701373524 ) evaluates 1.0564281487582248
                                            actually  1.0564281487582246

   Consequently, the bounds provided below must be expanded to allow for this.

   This file should be used with functions that evaluate to nearest.
   It is part of the release so it is available for use with a library that rounds to nearest.
=#

# functions not exported by CRlibm or available with extraCRlibm.jl
for (fn) in (:acsc, :asec, :acot,
             :csch, :sech, :coth,
             :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
             :erf, :erfinv)
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            lo = ($fn)(lo)
            hi = ($fn)(hi)
            if lo > hi
               lo,hi = hi,lo
            end
            lo = prevFloat(lo)
            hi = nextFloat(hi)

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x) such that f(x+dx) > f(x) for dx>0 and f(x) defined for -Inf..x..Inf
=#

for (fn) in (:exp, :expm1, :atan, :sinh, :tanh)
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            loIsOpen, hiIsOpen = boundries(S)
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(x.lo, RoundDown)
            else
                lo = ($fn)(x.lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(x.hi, RoundUp)
            else
                hi = ($fn)(x.hi, RoundNearest)
            end
            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x), f(x) defined for -Inf..x..Inf (f(x) may be +/-Inf)
=#

for (fn) in (:sin, :cos, :tan, :csc, :sec, :cot, :cosh)
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            loIsOpen, hiIsOpen = boundries(S)
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(lo, RoundDown)
            else
                lo = ($fn)(lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(hi, RoundUp)
            else
                hi = ($fn)(hi, RoundNearest)
            end
            # given the e.g. periodicity of (fn), it is possible that fn(x.lo) > fn(x.hi)
            if lo>hi
                hiIsOpen, loIsOpen = boundries(S)
                hi, lo = value(x)
                if loIsOpen
                    lo = ($fn)(lo, RoundDown)
                else
                    lo = ($fn)(lo, RoundNearest)
                end
                if hiIsOpen
                    hi = ($fn)(hi, RoundUp)
                else
                    hi = ($fn)(hi, RoundNearest)
                end
            end

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x) such that f(x+dx) > f(x) for dx>0 and f(x) defined for domainMin..x..Inf
=#

for (fn, domainMin) in ((:log, 0.0), (:log1p, -1.0))
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            if x.lo < ($domainMin)
                throw(ErrorException("DomainError: $($fn) expected x.lo>=($($domainMin)), got $(x)."))
            end
            loIsOpen, hiIsOpen = boundries(S)
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(lo, RoundDown)
            else
                lo = ($fn)(lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(hi, RoundUp)
            else
                hi = ($fn)(hi, RoundNearest)
            end

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x) such that f(x+dx) > f(x) for dx>0 and f(x) defined for domainMin..x..domainMax
=#

for (fn, domainMin, domainMax) in ((:asin, -1.0, 1.0), (:atanh, -1.0, 1.0), (:erfinv, -1.0, 1.0))
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            if (x.lo < ($domainMin)) | (x.hi > ($domainMax))
                throw(ErrorException("DomainError: $($fn) expected x.lo,x.hi within [($($domainMin),$($domainMax))], got $(x)."))
            end
            loIsOpen, hiIsOpen = boundries(S)
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(lo, RoundDown)
            else
                lo = ($fn)(lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(hi, RoundUp)
            else
                hi = ($fn)(hi, RoundNearest)
            end

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x) such that f(x+dx) < f(x) for dx>0 and f(x) defined for domainMin..x..domainMax
=#

for (fn, domainMin, domainMax) in ((:acos, -1.0, 1.0),)
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            if (x.lo < ($domainMin)) | (x.hi > ($domainMax))
                throw(ErrorException("DomainError: $($fn) expected x.lo,x.hi within [($($domainMin),$($domainMax))], got $(x)."))
            end
            # function is decreasing, so hi,lo sense is reversed
            hiIsOpen, loIsOpen = boundries(S)
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(lo, RoundDown)
            else
                lo = ($fn)(lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(hi, RoundUp)
            else
                hi = ($fn)(hi, RoundNearest)
            end

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#=
   functions y=f(x) and f(x) defined for -Inf..x..domainLoMax, domainHiMin..x..+Inf
   (there is a region of x where f(x) is undefined between domainLoMax and domainHiMin)
=#
#=
for (fn, domainLoMax, domainHiMin) in ((:acsc, -1.0, 1.0), (:asec, -1.0, 1.0), (:acoth, -1.0, 1.0))
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            if (($domainLoMax) < x.lo < ($domainHiMin)) | (($domainLoMax) < x.lo < ($domainHiMin))
                throw(ErrorException("DomainError: $($fn) expected x.lo,x.hi outside [($($domainLoMax),$($domainHiMin))], got $(x)."))
            end
            # function is decreasing, so hi,lo sense is reversed
            loIsOpen, hiIsOpen = boundries(S)
            lo, hi = value(x)
            if loIsOpen
                lo = ($fn)(lo, RoundDown)
            else
                lo = ($fn)(lo, RoundNearest)
            end
            if hiIsOpen
                hi = ($fn)(hi, RoundUp)
            else
                hi = ($fn)(hi, RoundNearest)
            end
            # given the nature of (fn), it is quite possible that fn(x.lo) > fn(x.hi)
            if lo>hi
                hiIsOpen,loIsOpen = boundries(S)
                hi,lo = value(x)
                if loIsOpen
                    lo = ($fn)(lo, RoundDown)
                else
                    lo = ($fn)(lo, RoundNearest)
                end
                if hiIsOpen
                    hi = ($fn)(hi, RoundUp)
                else
                    hi = ($fn)(hi, RoundNearest)
                end
            end
            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end
=#


text:
@testset "read" begin
    dir = @__DIR__
    file = joinpath(dir,  "..", "example", "data", "Malitson.csv")

    d0 = SellmeierFit.read(file)
    d1 = SellmeierFit.read(file, unit_prefix=1.0)
    @test (d0.λ == micro .* d1.λ) && d0.ε==d1.ε

    d2 = SellmeierFit.read(file, unit_prefix=1.0, s_col=2, n_col=1)
    @test (d2.λ .^2 == sort(d1.ε)) && (sort(d2.ε) == d1.λ .^ 2)

    d3 = SellmeierFit.read(file, s_var=FREQUENCY)
    @test (d3.λ == reverse(c₀ ./ d0.λ)) && (d3.ε == reverse(d0.ε))

    d4 = SellmeierFit.read(file, s_var=ENERGY)
    @test (d4.λ == reverse(ℎc₀ ./ d0.λ)) && (d4.ε == reverse(d0.ε))

    url = "https://refractiveindex.info/tmp/data/main/SiO2/Malitson.csv"
    try
        d5 = SellmeierFit.read(url)
        @test d5 == d0
    catch e
        @test e isa Downloads.RequestError
    end
end  # @testset "read"


text:
abstract AdevectionsSolver

type ResolveLL <: AdevectionsSolver
  n_off :: Int64
end

type ResolveLS <: AdevectionsSolver
  n_off :: Int64
end

type ResolveSL  <: AdevectionsSolver
  n_off :: Int64
end

type ResolveSS <: AdevectionsSolver
  n_off :: Int64
end

export ResolveSL, ResolveSS, ResolveLL , ResolveLS, AdevectionsSolver
export ExternalAdectionResolver
abstract ExternalAdectionResolver <: AdevectionsSolver
function advection(  m:: NonlinearSpectralModel, x :: SpectralDataField ,TR :: AdevectionsSolver)
  error("No function for abstract types")
end
export advection


text:
try
@testset "Lint Self" begin
    msgs = lintpkg("Lint")
    if !isempty(msgs)
        display(msgs)
    end
    # TODO: reenable when #200 fixed
    @test_broken isempty(msgs)
    @test_broken length(msgs) === 0
    @test_broken size(msgs) === (0,)
    # @test_throws BoundsError msgs[1]
end
end


text:
# manipulate keyframe images

using Caesar
using CloudGraphs
using Images, ImageMagick
using JSON
using Colors
using Gadfly
using IncrementalInference, KernelDensityEstimate

addrdict = nothing
include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
# prepare the factor graph with just one node
# (will prompt on stdin for db credentials)
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

session = "SESSTURT38"
frd = fetchrobotdatafirstpose(cloudGraph, session)
camf = frd["CAMK"][1,1]
camc = frd["CAMK"][1,3]
tagwidth = 0.172

dcamjl = DepthCamera(frd["CAMK"])
buildmesh!(dcamjl)

# convenience quick peek function
@async cloudimshow(cloudGraph, session, :x3)

# prematurely get subgraph around all landmarks
ids = getExVertexNeoIDs(cloudGraph.neo4j.connection, label="LANDMARK", session=session)

sfg = initfg(cloudgraph=cloudGraph, sessionname=session)
fetchsubgraph!(sfg, Int[ids[i][2] for i in 1:length(ids)], numneighbors=2 )

xx,ll = ls(sfg)

vsym = xx[5]  #:x19

# fetch and show images
cv = getCloudVert(cloudGraph, session, sym=vsym, bigdata=true);
sleep(0.5)
imdata = Caesar.getBigDataElement(cv, "keyframe_rgb").data;
img = ImageMagick.readblob(imdata);

@async imshowhackpng(imdata)
@async imshowhack(img)  # prefered

neis = CloudGraphs.get_neighbors(cloudGraph, cv)

neis[2].properties["label"]
jd = JSON.parse(neis[2].properties["frtend"])
meas = parse.(split(jd["meas"]))

jd

atbeari = meas[1]
atdepth = meas[2]#*520/570
deptstd = meas[5]
measimgcol = round(Int,camc-camf*tan(atbeari)) # measurement center
atvertsides = tagwidth/2/atdepth*camf

r, c = size(img)
for i in 1:r
  img[i, measimgcol] = RGB(1.0,0.0,0.0)
  side1 = round(Int,measimgcol-atvertsides)
  side2 = round(Int,measimgcol+atvertsides)
  if 0 < side1 <= c
    img[i, side1] = RGB(1.0,1.0,0.0)
  end
  if 0 < side2 <= c
    img[i, side2] = RGB(1.0,1.0,0.0)
  end
end

filename = joinpath("/tmp","testimg.png")
ImageMagick.save_(filename, img)

@async run(`eog $(filename)`)

# Look at xtion depth information

depthdata = Caesar.getBigDataElement(cv, "depthframe_image").data;
arr = bin2arr(depthdata, dtype=Float32) # should also store dtype for arr in Mongo
depthimg = Array{Float64,2}(reshape(arr, c, r)')

depthcloud = reconstruct(dcamjl, depthimg)

# get slam vertex elements
lms = landmarks(sfg, vsym)
@show lsym = lms[2]
# SLAM estimates of distance
slamdist = getRangeKDEMax2D(cloudGraph, session, vsym, lsym)

deptdist = Base.mean(depthcloud[260:270,measimgcol,3])

plot(
  layer(y=depthcloud[:,measimgcol,3],Geom.line),
  # layer(y=depthcloud[:,measimgcol,3]+deptstd,Geom.line),
  # layer(y=depthcloud[:,measimgcol,3]-deptstd,Geom.line),
  Guide.title("AprilTag says: $(round(atdepth,3))\n
  SLAM says is: $(round(slamdist,3))\n
  Depth cam is: $(round(deptdist,3))\n
  AprilTag id:$(lsym)")
)

# investigate belief on vertices
ls(sfg, lsym)
vsym
plotKDEresiduals(sfg, :x105l200545, api=localapi)
getVal(sfg, lsym, api=localapi)
plotKDE(sfg, lsym, dims=[1;2], api=localapi)

pro, = localProduct(sfg, lsym)
plotKDE(pro)
lvert = getVert(sfg, lsym, api=localapi)
setVal!(lvert, getPoints(pro), getBW(pro)[:,1])
plotKDE(getBelief(lvert))

plotLocalProduct(sfg, lsym, api=localapi)

# via full graph copy
# fg = initfg(cloudgraph=cloudGraph, sessionname=session)
# fullLocalGraphCopy!(fg)

# for l in ll
fsym = ls(sfg, ll[1])
fs = fsym[1]
if length(lsf(sfg, fs)) > 1
  @show fs
  ty = getfnctype(sfg, sfg.fIDs[fs])
  zij = ty.range.μ
  pij = getRangeKDEMax2D(sfg, ll[1])
end

# end
#


text:
#finitemps.jl
struct FiniteMPS{S<:EuclideanSpace} <: AbstractMPS{S}
    Ts::Vector{<:AbstractTensorMap{S,2,1}}
    Ss::Vector{<:AbstractTensorMap{S,1,1}}
    canonical_form::Vector{Symbol}
    function FiniteMPS{S}(Ts::Vector{<:AbstractTensorMap{S,2,1}},
                            Ss::Vector{<:AbstractTensorMap{S,1,1}},
                            canonical_form::Vector{Symbol}) where {S<:EuclideanSpace}
        (dim(space(Ts[1],1)) == 1 && dim(space(Ts[end],3)) == 1) ||
            error("The left bond of the first tensor and the right bonnd of the
                    last tensor should have dimension 1 for finite MPS!")
        canonical_form[1] ∈ [:left_canonical,:right_canonical,:site_canonical,
                                :bond_canonical,:GammaLambda_canonical,:non_canonical] ||
            error("There is no such canonical form!")
        for i in 1:length(Ts)-1
            space(Ts[i],3) == dual(space(Ts[i+1],1)) ||
                error("The contracted spaces of tensors do not match!")
        end
        new(Ts,Ss,canonical_form)
    end
end

FiniteMPS(Ts,Ss,canonical_form) = FiniteMPS{spacetype(Ts[1])}(Ts, Ss, canonical_form)

function make_left_canonical(mps::FiniteMPS)
    for i in 1:length(mps.Ts)-1
        Q, R = leftorth(mps.Ts[i])
        mps.Ts[i] = Q
        A_tmp = TensorMap(undef, codomain(R)*codomain(mps.Ts[i+1],2), domain(mps.Ts[i+1]))
        @tensor A_tmp[a b;c] = R[a d]*mps.Ts[i+1][d b;c]
        mps.Ts[i+1] = A_tmp
    end
    mps_norm = tr(mps.Ts[L]'*mps.Ts[L])
    mps.Ts[L] = mps.Ts[L]/sqrt(mps_norm)
    mps.canonical_form[1] = :left_canonical
    return
end

function make_right_canonical(mps::FiniteMPS)
    for i in length(mps.Ts):-1:2
        L, Q = rightorth(mps.Ts[i], (1,), (2,3))
        mps.Ts[i] = permute(Q,(1,2),(3,))
        A_tmp = TensorMap(undef, codomain(mps.Ts[i-1]), domain(L))
        @tensor A_tmp[a b;d] = mps.Ts[i-1][a b;c]*L[c;d]
        mps.Ts[i-1] = A_tmp
    end
    T_tmp = TensorMap(undef, codomain(mps.Ts[1],1),codomain(mps.Ts[1],1))
    @tensor T_tmp[a;d] := mps.Ts[1][a b;c]*conj(mps.Ts[1][d b;c])
    mps_norm = tr(T_tmp)
    mps.Ts[1] = mps.Ts[1]/sqrt(mps_norm)
    mps.canonical_form[1] = :right_canonical
    return
end

function LinearAlgebra.norm(mps::FiniteMPS)
    A = mps.Ts[1]'*mps.Ts[1]
    for i in 2:length(mps.Ts)
        A_tmp1 = TensorMap(undef,codomain(A)*codomain(mps.Ts[i],2),domain(mps.Ts[i]))
        @tensor A_tmp1[d b;c] = A[d;a]*mps.Ts[i][a b;c]
        A_tmp2 = TensorMap(undef,domain(mps.Ts[i]),domain(mps.Ts[i]))
        @tensor A_tmp2[d;c] =  conj(mps.Ts[i][a b;d])*A_tmp1[a b;c]
        A = A_tmp2
    end
    return tr(A)
end

function make_site_canonical(psi::FiniteMPS, n::Int)
end

function bond_canonical(psi::FiniteMPS, b::Int)
end


text:
include("silicon_runners.jl")

@testset "Silicon LDA (small, Float64)" begin
    run_silicon_lda(Float64, Ecut=7, test_tol=0.03, n_ignored=0, grid_size=17, scf_tol=1e-5)
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Silicon LDA (large, Float64)" begin
        run_silicon_lda(Float64, Ecut=25, test_tol=5e-6, n_ignored=0,
                        grid_size=33, scf_tol=1e-7)
    end
end

@testset "Silicon LDA (small, Float32)" begin
    run_silicon_lda(Float32, Ecut=7, test_tol=0.03, n_ignored=1, grid_size=19, scf_tol=1e-4,
                    n_noconv_check=1)
end


text:
using SparseSensors
using Test

@testset "SparseSensors.jl" begin
    # Write your tests here.
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ de098278-8e74-11eb-13e3-49c1b6e06e7d
using Pkg; Pkg.activate("MLJ_env", shared=true)

# ╔═╡ ac463e7a-8b59-11eb-229e-db560e17c5f5
begin
	using Test
	using Random
	using PlutoUI
	using LinearAlgebra
	using CSV
	using StringEncodings
	using Printf
	using Plots

	push!(LOAD_PATH, "./src")
	using YaWorkingData: DT, pca, transform
	using YaDistances: cosine_similarity
end

# ╔═╡ 8c80e072-8b59-11eb-3c21-a18fe43c4536
md"""
## Recommender Systems

ref. from book **"Data Science from Scratch"**, Chap 23

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ e7373726-8b59-11eb-2a2b-b5138e4f5268
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f5ee64b2-8b59-11eb-2751-0778efd589cd
html"""

"""

# ╔═╡ 81290d1c-8ce2-11eb-3340-337957fd81b7
html"""

"""

# ╔═╡ 8ff1bb20-8ce2-11eb-1de6-fd84daec8930
md"""
#### Recommending What is Popular
"""

# ╔═╡ d3ee2138-8ce2-11eb-0b29-659a3be01512
const Users_Interests = [
    ["Hadoop", "Big Data", "HBase", "Java", "Spark", "Storm", "Cassandra"],
    ["NoSQL", "MongoDB", "Cassandra", "HBase", "Postgres"],
    ["Python", "scikit-learn", "scipy", "numpy", "statsmodels", "pandas"],
    ["R", "Python", "statistics", "regression", "probability"],
    ["machine learning", "regression", "decision trees", "libsvm"],
    ["Python", "R", "Java", "C++", "Haskell", "programming languages"],
    ["statistics", "probability", "mathematics", "theory"],
    ["machine learning", "scikit-learn", "Mahout", "neural networks"],
    ["neural networks", "deep learning", "Big Data", "artificial intelligence"],
    ["Hadoop", "Java", "MapReduce", "Big Data"],
    ["statistics", "R", "statsmodels"],
    ["C++", "deep learning", "artificial intelligence", "probability"],
    ["pandas", "R", "Python"],
    ["databases", "HBase", "Postgres", "MySQL", "MongoDB"],
    ["libsvm", "regression", "support vector machines"]
]

# ╔═╡ 0fcb268a-8ea8-11eb-309c-ef92de740db2
length(Users_Interests)

# ╔═╡ 51cc0044-8e95-11eb-1d0b-e7e3e537710c
begin
	const VVS = Vector{Vector{S}} where S <: AbstractString
	const VVF = Vector{Vector{Float64}}
	const VS = Vector{String}
	const VP = Vector{Pair{String, Integer}}
	const VT_IF = Vector{Tuple{Integer, Float64}}
	const DSI = Dict{String, Integer}
	const AVT = AbstractArray{T, N} where {T <: Any, N}
end

# ╔═╡ 13a0940e-8ce4-11eb-231c-f331a607203c
function gen_popular_interests(ds::VVS)::DSI
  count_word = Dict{String, Integer}()
  for ui_s ∈ ds, w ∈ ui_s
    count_word[w] = get(count_word, w, 0) + 1
  end
  count_word
end

# ╔═╡ 6b5de662-8e95-11eb-1054-5144e8e0d5db
function most_popular_new_interests(ui::VS, pi::DSI; max_res=5)::VP
	suggestions = filter(((i, _f)=t) -> i ∉ ui,
		sort(collect(pi), by=t -> t[2], rev=true)
	)[1:max_res]
end

# ╔═╡ 19a89dd2-8eac-11eb-2faa-375f52037a94
pop_interests = gen_popular_interests(Users_Interests)

# ╔═╡ f9716ebe-8e96-11eb-1f7f-9f61e16a374b
typeof(Users_Interests)

# ╔═╡ 12d52658-8e9b-11eb-1825-870032e6ad57
md"""
If you are user 2 (`Users_Interests[2]`) we recommend:
`[("R" => 4), ("Python" => 4), ("Java" => 3), ("Big Data" => 3), ("statistics" => 3)]`:
"""

# ╔═╡ 930441a6-8e96-11eb-1374-7b8b4b5199af
begin
	@test most_popular_new_interests(Users_Interests[2], pop_interests) == [
		("R" => 4), ("Python" => 4), ("Java" => 3), ("Big Data" => 3),
		("statistics" => 3)
	]
end

# ╔═╡ ee423fb0-8e9a-11eb-161a-3f3ba95287de
md"""
If you are user 4, who is already interested in many of those things, you would get instead `[("HBase" => 3), ("Java" => 3), ("Big Data" => 3), ("Hadoop" => 2),  ("statsmodels" => 2)]`:
"""

# ╔═╡ 537f09c6-8e9b-11eb-342c-fdcf36135acf
@test most_popular_new_interests(Users_Interests[4], pop_interests) == [
	("HBase" => 3), ("Java" => 3), ("Big Data" => 3), ("Hadoop" => 2),
	("statsmodels" => 2)
]

# ╔═╡ 80ddf0f4-8ce2-11eb-3046-331119a0dc9b
html"""

"""

# ╔═╡ 80a71408-8ce2-11eb-3978-75a4a2df9116
md"""
#### User-Based Collaborative Filtering

One way of taking a user’s interests into account is to look for users who
are somehow similar to him/her, and then suggest the things that those users are
interested in.

In order to do that we need to measure how similar two users are and for this we need a "distance". A typical one for this case is the cosine distance.
To apply it, we will binary encode the user preference's vector. With this encoding, similar user (which means users with vectors pointing nearly in the same direction of the feature space) will have a similarity close to 1. Conversely, very dissimilar users will a similarity  of 0. Otherwise the similarity will fall between 0 and 1.

To achieve this, let us start collecting the known (unique) interests.
"""

# ╔═╡ d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
const Uniq_Interests = Iterators.flatten(Users_Interests) |>
	collect |>
	sort |>
	unique;

# ╔═╡ 1ede2ea4-8ea0-11eb-3cb2-6f9f009cc47c
@test Uniq_Interests[1:6] == String["Big Data", "C++", "Cassandra", "HBase", "Hadoop", "Haskell"]

# ╔═╡ 2b256fe4-8ead-11eb-094e-319b84541a47
length(Uniq_Interests), length(Users_Interests)

# ╔═╡ d2b198b6-8ce2-11eb-170e-0f17904c9f2c
begin

function gen_binary_vect(user_int::VS, uniq_int::VS)::BitArray
	v = BitArray(undef, length(uniq_int))
	for (ix, ui) ∈ enumerate(uniq_int)
		v[ix] = ui ∈ user_int ? 1 : 0
	end
	v
end

gen_binary_vect(user_int::VVS, uniq_int::VS) = gen_binary_vect.(user_int, uniq_int)

end

# ╔═╡ 5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
## NOTE: we need to prevent broadcasting on the second array => Ref()
const User_Interest_Vect = gen_binary_vect.(Users_Interests, Ref(Uniq_Interests))

# ╔═╡ 1df61abc-8ead-11eb-0613-73865afa45f4
length(User_Interest_Vect[1])

# ╔═╡ 394093fa-8cea-11eb-0071-eff3045a012b
md"""
Now `user_interest_vect[i][j]` is 1 if user i has interest in j and 0 otherwise.

Check cosine similarity definition in `YaDistances` module.
"""

# ╔═╡ 526acb22-8ea2-11eb-158e-63c3ec6df336
#function cosine_similarity(v₁::AVT, v₂::AVT)::Float64 # where {T <: Any}
#	dot(v₁, v₂) / (norm(v₁) * norm(v₂))
#end

# ╔═╡ 66a3df24-8ea5-11eb-309c-f549d788b65d
# begin
# 	@test cosine_similarity([1., 1, 1], [2., 2, 2]) ≈ 1. ## "same direction"
# 	@test cosine_similarity([-1., -1], [2., 2]) ≈ -1.    ## "opposite direction"
# 	@test cosine_similarity([1., 0], [0., 1]) ≈ 0.       ## "orthogonal"
# end

# ╔═╡ 1a305bc2-8ea3-11eb-128a-27139da33ae3
md"""
Because we have a small dataset, it is not a problem to compute the pairwise
similarities between all of our users (*it is symetric*) :
"""

# ╔═╡ 435c4444-8dfd-11eb-24e8-5543f905f199
# user_similarities = [
# 	[cosine_similarity(u_i, u_j) for (i, u_i) ∈ enumerate(user_interest_vect)
# 			for (j, u_j) ∈ enumerate(user_interest_vect) if j > i]
# ]

User_Similarities = [
	[cosine_similarity(u_i, u_j) for u_i ∈ User_Interest_Vect]
			for u_j ∈ User_Interest_Vect
]

# ╔═╡ 1ab295c8-8ea4-11eb-166c-8de4f02f0fce
md"""
`user_similarities[i][j]` gives us the similarity between users i and j (i < j)
"""

# ╔═╡ 7365ac88-8eaa-11eb-04da-0ff2fe7e8361
typeof(User_Similarities)

# ╔═╡ e95983ba-8ceb-11eb-38fd-ed92cdcf754c
## Users 1 and 10 share interests in Hadoop, Java, and Big Data
@test 0.56 ≤ User_Similarities[1][10] ≤ 0.58  ## several shared interests

## @test 0.56 ≤ user_similarities[1][9] ≤ 0.58  ## several shared interests

# ╔═╡ 88c98664-8ea6-11eb-1f63-abf82c75db63
## Users 1 and 9 share only one interest: Big Data
@test 0.18 ≤ User_Similarities[1][9] ≤ 0.20  ## "only one shared interest"

## @test 0.18 ≤ user_similarities[1][8] ≤ 0.20  ## "only one shared interest"

# ╔═╡ d3a749a2-8cec-11eb-1f06-b568f244b576
md"""
In particular, `user_similarities[i]` is the vector of user i’s similarities
to every other user. We can use this to write a function that finds the most
similar users to a given user (avoiding including the user himself/herself,
nor any users with zero similarity). Then we will sort the results from most
similar to least similar:
"""

# ╔═╡ 31d3226a-8cf4-11eb-0897-39989ba76b58
function most_similar_users(user_simil::VVF, user_id::Integer)::VT_IF
	type_ = eltype(user_simil[user_id])
	filter(
		((uid, sim)=t) -> uid ≠ user_id && sim > zero(type_),
		collect(enumerate(user_simil[user_id]))
	) |> a -> sort(a, by=t -> t[end], rev=true)
end

# ╔═╡ c99f7c3c-8eaa-11eb-0d3c-8def70eef326
most_similar_users(User_Similarities, 1)

# ╔═╡ 683fb19e-8eaf-11eb-1bed-67aa24fcd830
typeof(most_similar_users(User_Similarities, 1))

# ╔═╡ 6e7e7896-8cf8-11eb-062e-99492ec8cff8
md"""
How do we use this to suggest new interests to a user?

For each interest we can add up the user similarities of the other iuser interested in it:
"""

# ╔═╡ 70d84470-8cf6-11eb-23b6-c7ba506d8552
function user_based_suggestions(user_id::Integer, user_simil, user_inter::VVS;
		incl_curr_interests=false)
	#
	suggs = Dict{String, Float64}()
	for (o_uid, simil) ∈ most_similar_users(user_simil, user_id),
		inter ∈ user_inter[o_uid]
		suggs[inter] = get(suggs, inter, 0.) + simil
	end

	suggs = sort(collect(suggs), by=t -> t[end], rev=true)  ## Vector of pairs

	incl_curr_interests && (return suggs)

	[(sugg, w) for (sugg, w) ∈ suggs if sugg ∉ user_inter[user_id]]
end

# ╔═╡ c451451a-8cf7-11eb-183b-e358c9d618e0
user_based_suggestions(1, User_Similarities, Users_Interests)

# ╔═╡ ebd531e6-8cf7-11eb-228f-d530e7a48151
md"""
This seems to make sense for someone whose stated interests are "Big Data" and database related.

Note however that this approach does not scale well and in large dimensional vector spaces, vectors are far apart.
"""

# ╔═╡ 6f3fc602-8ed5-11eb-050b-977276cabfee
html"""

"""

# ╔═╡ 6f006b24-8ed5-11eb-11d5-c32ce953dcaf
md"""
#### Item-Based Collaborative Filtering
An alternative approach is to compute similarities between interests directly. We can then generate suggestions for each user by aggregating interests that are similar to his/her current interests.

First let us transpose our user_interest_vect (a list of column vector) into a matrix of row vector.
"""

# ╔═╡ 30bdc788-8d01-11eb-2f8f-9fd79593ddb8
size(User_Interest_Vect), typeof(User_Interest_Vect)

# ╔═╡ 6f36902c-8ed6-11eb-08c4-d545f9adab62
const Interest_User_Matrix = hcat.(User_Interest_Vect...)

# ╔═╡ df8f698e-8ed6-11eb-12c9-3b495c1356f3
size(Interest_User_Matrix)

# ╔═╡ eb391ec2-8ed6-11eb-306b-0d5c56ef4063
@test Interest_User_Matrix[1] == [1 0 0 0 0 0 0 0 1 1 0 0 0 0 0]

## user 1, 9 and 10 indicated their interest in Big Data

# ╔═╡ 84946000-8d02-11eb-3160-571efee8fb0b
md"""
With this we can now use the cosine similarity again. If same users are interested in two topics, their similarity will be 1 and conversely when two users are interested in different topices tehir similarity will be 0.
"""

# ╔═╡ bdf358d8-8ed8-11eb-263e-d9adffd4a1db
const Interest_Similarities = [
	[cosine_similarity(uvectᵢ, uvectⱼ) for uvectⱼ ∈ Interest_User_Matrix]
	for uvectᵢ ∈ Interest_User_Matrix
]

# ╔═╡ d05f1c40-8d17-11eb-2724-4d989c4c6b92
function most_similar_interests(inter_id::Integer, inter_similarities::VVF,
		uniq_inters)
	simils = inter_similarities[inter_id]

	[(uniq_inters[o_inter_id], sim) for (o_inter_id, sim) ∈ enumerate(simils)
	 		if inter_id ≠ o_inter_id && sim > 0.] |>
	pairs -> sort(pairs, by=pair -> pair[end], rev=true)
end

# ╔═╡ 70b220c4-8edb-11eb-0db8-97a36ff4675b
md"""
We can, for example, find the interests most similar to "big Data" (interest 1) using  the following, which suggest these similar interests:
"""

# ╔═╡ 848fd812-8d02-11eb-01d6-75965b08bcc5
msit1 = most_similar_interests(1, Interest_Similarities, Uniq_Interests)

# ╔═╡ 7d35f85e-8e07-11eb-228c-81d28a444b52
begin
	@test msit1[1][1] == "Hadoop"
	@test 0.815 < msit1[1][2] < 0.817
	@test msit1[2][1] == "Java"
	@test 0.666 < msit1[2][2] < 0.667
end

# ╔═╡ 6e7c18b4-8e09-11eb-0238-2515d39a89cd
md"""
Thus now we can create recommendations for a user by summing up the similarities of the interests similar to his/her.
"""

# ╔═╡ 80dff352-8d02-11eb-06b1-5f5e325046f5
function item_based_suggestions(user_id::Integer, user_inter_vect, users_inter,
		inter_similarities::VVF, uniq_inters;
		incl_curr_interests=false)
	#
	suggs = Dict{String, Float64}()
	user_inter_v = user_inter_vect[user_id]
	for (inter_id, is_inter) ∈ enumerate(user_inter_v)
		if is_inter == 1
			similar_inters = most_similar_interests(inter_id, inter_similarities,
				uniq_inters)
			for (inter, sim) ∈ similar_inters
				suggs[inter] = get(suggs, inter, 0.) + sim
			end
		end
	end
	# sort...
	suggs = collect(suggs) |>
		p -> sort(p, by=pair -> pair[end], rev=true)
	#
	incl_curr_interests && (return suggs)
	[(sugg, w) for (sugg, w) ∈ suggs if sugg ∉ users_inter[user_id]]
end

# ╔═╡ 34811c32-8edd-11eb-0018-e73666ea7667
ibs1 = item_based_suggestions(1, User_Interest_Vect, Users_Interests,
	Interest_Similarities, Uniq_Interests)

# ╔═╡ d878e9b6-8edd-11eb-03ab-f3a908ee47fe
begin
	@test ibs1[1][1] == "MapReduce"
	@test 1.86 < ibs1[1][2] < 1.87

	@test ibs1[2][1] ∈ ("Postgres", "MongoDB")  # A tie
	@test 1.31 < ibs1[2][2] < 1.32
end

# ╔═╡ d9e37240-8f31-11eb-1caf-dbbf84b424f0
html"""

"""

# ╔═╡ d962551e-8f31-11eb-3ddf-8bc8b04b3086
md"""
#### Matrix Factorization
"""

# ╔═╡ 1a9f18a4-8f32-11eb-056a-b3ee5e093957
begin
	const MOVIES = "./recommender_data/ml-100k/u.item"
	const RATINGS = "./recommender_data/ml-100k/u.data"
end

# ╔═╡ 1a84db5e-8f32-11eb-2724-23abff284462
struct Rating
	user_id:: String
	movie_id::String
	rating::Float32
end

# ╔═╡ 1a69eff8-8f32-11eb-1513-63f5a01b3ec0
begin
	movies = Dict{String, String}()

	for row ∈ CSV.File(open(read, MOVIES, enc"ISO-8859-1"); delim="|")
    	movies[string(row[1])] = row[2]
	end
end

# ╔═╡ d3183c3a-8f38-11eb-21ac-bdb7bb7d21c7
@test length(movies) == 1681

# ╔═╡ 40705800-8f39-11eb-1f91-b3a013911ef8
collect(keys(movies))[1:10]

# ╔═╡ 5f543412-8f39-11eb-3a1f-4d7b9960e772
movies["519"]

# ╔═╡ 1a4d76a0-8f32-11eb-2293-51b004ec835f
begin
	ratings= Rating[]

	for row ∈ CSV.File(open(read, RATINGS, enc"ISO-8859-1"); delim="\t")
    	push!(ratings, Rating(string(row[1]), string(row[2]), Float32(row[3])))
	end
end

# ╔═╡ db2dbb10-8f31-11eb-259b-2f0f3ca5a29d
@test length(rating.user_id for rating ∈ ratings) == 99999

# ╔═╡ 7237b428-8f39-11eb-2b91-effb877198e5
##
## Create a dictionary for ratings by movie ids
##
# starwars_ratings = Dict{String, Vector{Float32}}(
# 	movie_id => Float32[] for (movie_id, title) ∈ movies if occursin(r"Star Wars|Empire Strikes|Jedi", title)
# )
starwars_ratings = reduce(
	(hsh, (movie_id, title)=r) -> (
		occursin(r"Star Wars|Empire Strikes|Jedi", title) &&
			(hsh[movie_id] = Float32[]); hsh
		),
	movies;
	init=Dict{String, Vector{Float32}}()
)

# ╔═╡ 71c9a9fe-8f39-11eb-1d2e-5994b912b366
begin
	## Iterate over ratings, accumulating the Star Wars ones
	for rating ∈ ratings
		rating.movie_id ∈ keys(starwars_ratings) &&
			(push!(starwars_ratings[rating.movie_id], rating.rating))
	end

	## Compute avg rating for each movie
	avg_ratings = [
		(sum(title_ratings) / length(title_ratings), movie_id)
			for (movie_id, title_ratings) ∈ starwars_ratings
	]

	## then print them in order
	with_terminal() do
		for (avg_rating, movie_id) ∈ sort(avg_ratings, by=t -> t[1], rev=true)
			@printf("%.2f %5s\n", avg_rating, movies[movie_id])
		end
	end
end

# ╔═╡ fb52634a-8f3c-11eb-173f-fdd2dfe5d9d3
md"""
Ok, let us try to come up with a model to predict these ratings. Our first step is to split the ratings data into 3 subsets: train, validation and test
"""

# ╔═╡ f093174a-8f43-11eb-154c-e54e019c2a93
typeof(ratings)

# ╔═╡ 61fb5714-8f3d-11eb-265d-97ac6dc8f150
begin
	Random.seed!(42)
	Random.shuffle!(ratings)
	s₁, s₂ = round(Int, length(ratings) * .7), round(Int, length(ratings) * .05)

	train, valid = ratings[1:s₁], ratings[s₁+1:s₁+1+s₂]
	test = ratings[s₁+2+s₂:end]

	length(train), length(valid), length(test)
end

# ╔═╡ 61df4254-8f3d-11eb-1b73-d5045174f781
@test length(train) + length(valid) + length(test) == length(ratings)

# ╔═╡ 3adaca50-8f3f-11eb-2881-a1bc37f221b0
train[1], typeof(train)

# ╔═╡ 61c6f168-8f3d-11eb-311c-47ad4dc9bb5e
md"""
Let us define a baseline model which will predict the average rating. We will use MSE (Mean Squared Error) as our metric and check how the baseline does on our test set.
"""

# ╔═╡ 71ac1ecc-8f39-11eb-0489-7bda4977893d
begin
	bl_avg_ratings = map(row -> row.rating, train) |>
		sum |> s -> s/length(train)

	bl_error = map(row -> (row.rating - bl_avg_ratings) ^ 2, test) |>
		sum |> s -> s/length(test);
end

# ╔═╡ ea84aeb2-8f3f-11eb-09e2-c7cb101a48d1
@test 1.26 < bl_error < 1.28   ## this is what we hope to beat...

# ╔═╡ ea67ca68-8f3f-11eb-28ce-5b807376fa85
md"""
Given our embeddings, the predicted ratings are given by the matrix product of the user embeddings and the movie embeddings.
For a given user and movie, that value is just the dot product of the corresponding
embeddings.

Let us start by creating the embeddings. We will represent them as dictionaries
where the keys are IDs and the values are vectors, which will allow us to easily retrieve the embedding for a given ID:
"""

# ╔═╡ 2580fa72-8f40-11eb-1894-8dc53902cff1
begin
	const EMB_DIM = 2

	user_ids = map(r -> r.user_id, ratings) |> unique
	movie_ids = map(r -> r.movie_id, ratings) |> unique

	user_vects = Dict{String, Vector{Float32}}(
		user_id => rand(Float32, EMB_DIM) for user_id ∈ user_ids
	)
	movie_vects = Dict{String, Vector{Float32}}(
		movie_id => rand(Float32, EMB_DIM) for movie_id ∈ movie_ids
	)
end

# ╔═╡ 25666c1c-8f40-11eb-265b-070dc6e71052
md"""
Now, it is time to write our taining loop.
"""

# ╔═╡ 4d52b86a-8f41-11eb-07de-632ac5b2c513
function train_loop(ds::Vector{Rating}, ds_name::Symbol;
		η::Union{Float32, Nothing}=nothing)
	loss = 0.
	for (ix, rating) ∈ enumerate(ds)
		movie_v = movie_vects[rating.movie_id]
		user_v = user_vects[rating.user_id]
		ŷ = dot(user_v, movie_v)
		err = ŷ - rating.rating
		loss += err ^ 2
		if !isnothing(η)
			## we have ŷ ≡ m₀ × u₀	+ m₁ × u₁ + ... +mₖ × uₖ
			## thus each uⱼ contributes to output with coefficient mⱼ
			## and conv. each mⱼ contributes to output with coefficient uⱼ
			∇user = movie_v * err # [err * mⱼ for mⱼ ∈ movie_v]
			∇movie = user_v * err
			#
			user_v .-= η * ∇user
			movie_v .-= η * ∇movie
		end
		ix % 10_000 == 0 && @printf("\t%5d: avg_loss: %3.6f\n", ix, loss / ix)
	end
	avg_loss = loss / length(ds)
	@printf("\tFinal avg_loss(%12s): %3.6f\n", ds_name, avg_loss)
	avg_loss
end

# ╔═╡ 25489752-8f40-11eb-129b-e52fde513eed
begin
	η = Float32(0.05)
	avg_train_losses, avg_valid_losses = [], []
	avg_test_loss = nothing

	with_terminal() do
		for epoch ∈ 1:20
			global η *= Float32(0.9)
			@printf("%2d => η: %2.5f", epoch, η)
			#
			avg_train_loss = train_loop(train, :Training; η)
			avg_valid_loss = train_loop(valid, :Validation)
			push!(avg_train_losses, avg_train_loss)
			push!(avg_valid_losses, avg_valid_loss)
		end
		#
		global avg_test_loss = train_loop(test, :Test)
	end
end

# ╔═╡ a1bc5930-8f55-11eb-0fdb-af0f194e92b3
(avg_test_loss=string(round(avg_test_loss * 100., digits=2), "%"),)

# ╔═╡ 27dd9746-8f55-11eb-28a7-d3037bd284e8
begin
	plot(avg_train_losses, title="avg loss/epoch", label="train loss")
	plot!(avg_valid_losses, label="valid loss")
end

# ╔═╡ 46a227c8-8f55-11eb-07f9-f7eaff5b38de
md"""
Now, we will inspect  the learned vectors. Because there is no reason to expect the two components to be particularly meaningful, we are going to use PCA.
"""

# ╔═╡ f82b1e48-8f71-11eb-1988-c970a4a0a5a3
begin
	orig_vects = values(movie_vects) |> collect
	comps = pca(orig_vects, 2)
end

# ╔═╡ 143bd95e-8f74-11eb-10af-49df13561522
md"""
Let us transform our vectors to represent the principal components and join
in the movie IDs and average ratings:
"""

# ╔═╡ 21dcc152-8f76-11eb-1511-674e702e4072
begin
	ratings_by_movie = Dict{String, DT}()
	for rating in ratings
		id = rating.movie_id
		ary = get(ratings_by_movie, id, DT{Float32}())
		push!(ary, rating.rating)
		ratings_by_movie[id] = ary
	end
	vs = [
		(
			movie_id,
			sum(ratings_by_movie[movie_id]) / length(ratings_by_movie[movie_id]),
			get(movies, movie_id, "_"),
			vect
		)
		for (movie_id, vect) ∈ zip(keys(movie_vects), transform(orig_vects, comps))
	]
end

# ╔═╡ f1c75d7c-8f7c-11eb-1e75-dbc0b7d8225d
# top 25 by first principal component
sort(vs, by=t -> t[end][1], rev=true)[1:25]

# ╔═╡ 51e6fadc-8f7d-11eb-1f1f-93756139d30b
# bottom 25 by first principal component
sort(vs, by=t -> t[end][1], rev=false)[1:25]

# ╔═╡ 81cb1cec-8f7d-11eb-0b7d-4f4c7dd01a2c
md"""
The top 25 are all highly rated, while the bottom 25 are *mostly* low-rated (or unrated in the training data), which suggests that the first principal component is mostly capturing “how good is this movie?”

It is hard to make much sense of the second component; and, indeed the 2-dimensional embeddings performed only slightly better than the one-dimensional embeddings, suggesting that whatever the second component captured is possibly very subtle...
"""

# ╔═╡ Cell order:
# ╟─8c80e072-8b59-11eb-3c21-a18fe43c4536
# ╠═de098278-8e74-11eb-13e3-49c1b6e06e7d
# ╠═ac463e7a-8b59-11eb-229e-db560e17c5f5
# ╟─e7373726-8b59-11eb-2a2b-b5138e4f5268
# ╟─f5ee64b2-8b59-11eb-2751-0778efd589cd
# ╟─81290d1c-8ce2-11eb-3340-337957fd81b7
# ╟─8ff1bb20-8ce2-11eb-1de6-fd84daec8930
# ╠═d3ee2138-8ce2-11eb-0b29-659a3be01512
# ╠═0fcb268a-8ea8-11eb-309c-ef92de740db2
# ╠═51cc0044-8e95-11eb-1d0b-e7e3e537710c
# ╠═13a0940e-8ce4-11eb-231c-f331a607203c
# ╠═6b5de662-8e95-11eb-1054-5144e8e0d5db
# ╠═19a89dd2-8eac-11eb-2faa-375f52037a94
# ╠═f9716ebe-8e96-11eb-1f7f-9f61e16a374b
# ╟─12d52658-8e9b-11eb-1825-870032e6ad57
# ╠═930441a6-8e96-11eb-1374-7b8b4b5199af
# ╟─ee423fb0-8e9a-11eb-161a-3f3ba95287de
# ╠═537f09c6-8e9b-11eb-342c-fdcf36135acf
# ╟─80ddf0f4-8ce2-11eb-3046-331119a0dc9b
# ╟─80a71408-8ce2-11eb-3978-75a4a2df9116
# ╠═d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
# ╠═1ede2ea4-8ea0-11eb-3cb2-6f9f009cc47c
# ╠═2b256fe4-8ead-11eb-094e-319b84541a47
# ╠═d2b198b6-8ce2-11eb-170e-0f17904c9f2c
# ╠═5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
# ╠═1df61abc-8ead-11eb-0613-73865afa45f4
# ╟─394093fa-8cea-11eb-0071-eff3045a012b
# ╠═526acb22-8ea2-11eb-158e-63c3ec6df336
# ╠═66a3df24-8ea5-11eb-309c-f549d788b65d
# ╟─1a305bc2-8ea3-11eb-128a-27139da33ae3
# ╠═435c4444-8dfd-11eb-24e8-5543f905f199
# ╟─1ab295c8-8ea4-11eb-166c-8de4f02f0fce
# ╠═7365ac88-8eaa-11eb-04da-0ff2fe7e8361
# ╠═e95983ba-8ceb-11eb-38fd-ed92cdcf754c
# ╠═88c98664-8ea6-11eb-1f63-abf82c75db63
# ╟─d3a749a2-8cec-11eb-1f06-b568f244b576
# ╠═31d3226a-8cf4-11eb-0897-39989ba76b58
# ╠═c99f7c3c-8eaa-11eb-0d3c-8def70eef326
# ╠═683fb19e-8eaf-11eb-1bed-67aa24fcd830
# ╟─6e7e7896-8cf8-11eb-062e-99492ec8cff8
# ╠═70d84470-8cf6-11eb-23b6-c7ba506d8552
# ╠═c451451a-8cf7-11eb-183b-e358c9d618e0
# ╟─ebd531e6-8cf7-11eb-228f-d530e7a48151
# ╟─6f3fc602-8ed5-11eb-050b-977276cabfee
# ╟─6f006b24-8ed5-11eb-11d5-c32ce953dcaf
# ╠═30bdc788-8d01-11eb-2f8f-9fd79593ddb8
# ╠═6f36902c-8ed6-11eb-08c4-d545f9adab62
# ╠═df8f698e-8ed6-11eb-12c9-3b495c1356f3
# ╠═eb391ec2-8ed6-11eb-306b-0d5c56ef4063
# ╟─84946000-8d02-11eb-3160-571efee8fb0b
# ╠═bdf358d8-8ed8-11eb-263e-d9adffd4a1db
# ╠═d05f1c40-8d17-11eb-2724-4d989c4c6b92
# ╟─70b220c4-8edb-11eb-0db8-97a36ff4675b
# ╠═848fd812-8d02-11eb-01d6-75965b08bcc5
# ╠═7d35f85e-8e07-11eb-228c-81d28a444b52
# ╟─6e7c18b4-8e09-11eb-0238-2515d39a89cd
# ╠═80dff352-8d02-11eb-06b1-5f5e325046f5
# ╠═34811c32-8edd-11eb-0018-e73666ea7667
# ╠═d878e9b6-8edd-11eb-03ab-f3a908ee47fe
# ╟─d9e37240-8f31-11eb-1caf-dbbf84b424f0
# ╟─d962551e-8f31-11eb-3ddf-8bc8b04b3086
# ╠═1a9f18a4-8f32-11eb-056a-b3ee5e093957
# ╠═1a84db5e-8f32-11eb-2724-23abff284462
# ╠═1a69eff8-8f32-11eb-1513-63f5a01b3ec0
# ╠═d3183c3a-8f38-11eb-21ac-bdb7bb7d21c7
# ╠═40705800-8f39-11eb-1f91-b3a013911ef8
# ╠═5f543412-8f39-11eb-3a1f-4d7b9960e772
# ╠═1a4d76a0-8f32-11eb-2293-51b004ec835f
# ╠═db2dbb10-8f31-11eb-259b-2f0f3ca5a29d
# ╠═7237b428-8f39-11eb-2b91-effb877198e5
# ╠═71c9a9fe-8f39-11eb-1d2e-5994b912b366
# ╟─fb52634a-8f3c-11eb-173f-fdd2dfe5d9d3
# ╠═f093174a-8f43-11eb-154c-e54e019c2a93
# ╠═61fb5714-8f3d-11eb-265d-97ac6dc8f150
# ╠═61df4254-8f3d-11eb-1b73-d5045174f781
# ╠═3adaca50-8f3f-11eb-2881-a1bc37f221b0
# ╟─61c6f168-8f3d-11eb-311c-47ad4dc9bb5e
# ╠═71ac1ecc-8f39-11eb-0489-7bda4977893d
# ╠═ea84aeb2-8f3f-11eb-09e2-c7cb101a48d1
# ╟─ea67ca68-8f3f-11eb-28ce-5b807376fa85
# ╠═2580fa72-8f40-11eb-1894-8dc53902cff1
# ╟─25666c1c-8f40-11eb-265b-070dc6e71052
# ╠═4d52b86a-8f41-11eb-07de-632ac5b2c513
# ╠═25489752-8f40-11eb-129b-e52fde513eed
# ╠═a1bc5930-8f55-11eb-0fdb-af0f194e92b3
# ╠═27dd9746-8f55-11eb-28a7-d3037bd284e8
# ╟─46a227c8-8f55-11eb-07f9-f7eaff5b38de
# ╠═f82b1e48-8f71-11eb-1988-c970a4a0a5a3
# ╟─143bd95e-8f74-11eb-10af-49df13561522
# ╠═21dcc152-8f76-11eb-1511-674e702e4072
# ╠═f1c75d7c-8f7c-11eb-1e75-dbc0b7d8225d
# ╠═51e6fadc-8f7d-11eb-1f1f-93756139d30b
# ╟─81cb1cec-8f7d-11eb-0b7d-4f4c7dd01a2c

max_stars_repo_path: Data Science From Scratch/23-RecomSys.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
const J2000_TO_JULIAN = 2.451545e6days
const J2000_TO_MJD = 51544.5days

"""
    Epoch{S}(jd1::T, jd2::T=zero(T); origin=:j2000) where {S, T<:AstroPeriod}

Construct an `Epoch` with time scale `S` from a Julian date
(optionally split into `jd1` and `jd2`). `origin` determines the
variant of Julian date that is used. Possible values are:

- `:j2000`: J2000 Julian date, starts at `2000-01-01T12:00`
- `:julian`: Julian date, starts at `-4712-01-01T12:00`
- `:modified_julian`: Modified Julian date, starts at `1858-11-17T00:00`

### Examples ###

```jldoctest; setup = :(using AstroTime)
julia> Epoch{InternationalAtomicTime}(0.0days, 0.5days)
2000-01-02T00:00:00.000 TAI

julia> Epoch{InternationalAtomicTime}(2.451545e6days, origin=:julian)
2000-01-01T12:00:00.000 TAI
```
"""
function Epoch{S}(jd1::T, jd2::T=zero(T), args...; origin=:j2000) where {S, T<:AstroPeriod}
    if jd2 > jd1
        jd1, jd2 = jd2, jd1
    end

    u = unit(jd1)

    if origin == :j2000
        # pass
    elseif origin == :julian
        jd1 -= u(J2000_TO_JULIAN)
    elseif origin == :modified_julian
        jd1 -= u(J2000_TO_MJD)
    else
        throw(ArgumentError("Unknown Julian epoch: $origin"))
    end

    jd1v = jd1 |> seconds |> value
    jd2v = jd2 |> seconds |> value

    sum, residual = two_sum(jd1v, jd2v)
    epoch = floor(Int64, sum)
    offset = (sum - epoch) + residual
    return Epoch{S}(epoch, offset)
end

"""
    julian_period([T,] ep::Epoch; origin=:j2000, scale=timescale(ep), unit=days)

Return the period since Julian Epoch `origin` within the time scale `scale` expressed in
`unit` for a given epoch `ep`. The result is an [`AstroPeriod`](@ref) object by default.
If the type argument `T` is present, the result is converted to `T` instead.

### Example ###

```jldoctest; setup = :(using AstroTime)
julia> ep = TAIEpoch(2018, 2, 6, 20, 45, 0.0)
2018-02-06T20:45:00.000 TAI

julia> julian_period(ep; scale=TT)
6611.364955833334 days

julia> julian_period(ep; unit=years)
18.100929728496464 years

julia> julian_period(Float64, ep)
6611.364583333333
```
"""
function julian_period(ep::Epoch; origin=:j2000, scale=timescale(ep), unit=days)
    ep1 = Epoch(ep, scale)
    jd1 = unit(ep1.second * seconds)
    jd2 = unit(ep1.fraction * seconds)

    if origin == :j2000
        # pass
    elseif origin == :julian
        jd1 += unit(J2000_TO_JULIAN)
    elseif origin == :modified_julian
        jd1 += unit(J2000_TO_MJD)
    else
        throw(ArgumentError("Unknown Julian epoch: $origin"))
    end

    return jd2 + jd1
end

function julian_period(::Type{T}, ep::Epoch; kwargs...) where T
    jd = julian_period(ep; kwargs...)
    return T(value(jd))
end

"""
    j2000(ep)

Return the J2000 Julian Date for epoch `ep`.

### Example ###

```jldoctest; setup = :(using AstroTime)
julia> j2000(TAIEpoch(2000, 1, 1, 12))
0.0 days
```
"""
j2000(ep::Epoch) = julian_period(ep)

"""
    julian(ep)

Return the Julian Date for epoch `ep`.

### Example ###

```jldoctest; setup = :(using AstroTime)
julia> julian(TAIEpoch(2000, 1, 1, 12))
2.451545e6 days
```
"""
julian(ep::Epoch) = julian_period(ep; origin=:julian)

"""
    modified_julian(ep)

Return the Modified Julian Date for epoch `ep`.

### Example ###

```jldoctest; setup = :(using AstroTime)
julia> modified_julian(TAIEpoch(2000, 1, 1, 12))
51544.5 days
```
"""
modified_julian(ep::Epoch) = julian_period(ep; origin=:modified_julian)

"""
    julian_twopart(ep)

Return the two-part Julian Date for epoch `ep`, which is a tuple consisting
of the Julian day number and the fraction of the day.

### Example ###

```jldoctest; setup = :(using AstroTime)
julia> julian_twopart(TAIEpoch(2000, 1, 2))
(2.451545e6 days, 0.5 days)
```
"""
function julian_twopart(ep::Epoch)
    sec_in_days = ep.second / SECONDS_PER_DAY
    frac_in_days = ep.fraction / SECONDS_PER_DAY
    j2k1, j2k2 = divrem(sec_in_days, 1)
    jd1 = j2k1 * days + J2000_TO_JULIAN
    jd2 = (frac_in_days + j2k2) * days
    return jd1, jd2
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
module Signals

using SignalAnalysis: signal
using Statistics: mean
using DataFrames: DataFrame, DataFrameRow
using TimeZones: ZonedDateTime, localzone, @tz_str, astimezone
using Dates: unix2datetime
using Base64: base64decode
import PooledArrays: PooledArray

function read(filename; tz=localzone())
  df = DataFrame(
    time = ZonedDateTime[],
    rxtime = Int64[],
    rssi = Union{Float64,Missing}[],
    preamble = Int64[],
    channels = Int64[],
    fc = Float64[],
    fs = Float64[],
    len = Int64[],
    lno = Int64[],
    filename = PooledArray(String[]),
    dtype = DataType[])
  for (lno, line) ∈ enumerate(eachline(filename))
    if contains(line, "|RxBasebandSignalNtf:INFORM[")
      m = match(r"^(\d+)\|RxBasebandSignalNtf:INFORM\[(.*) \((\d+) .*samples\)\]", line)
      t = astimezone(ZonedDateTime(unix2datetime(parse(Int64, m[1]) / 1000), tz"UTC"), tz)
      len = parse(Int64, m[3])
      ps = Dict(=>(split(p, ':')...) for p ∈ split(m[2], ' '; keepempty=false))
      rxtime = parse(Int64, ps["rxTime"])
      rssi = "rssi" ∈ keys(ps) ? parse(Float64, ps["rssi"]) : missing
      pre = "preamble" ∈ keys(ps) ? parse(Int64, ps["preamble"]) : 0
      ch = "channels" ∈ keys(ps) ? parse(Int64, ps["channels"]) : 1
      fc = parse(Float64, ps["fc"])
      fs = parse(Float64, ps["fs"])
      dt = fc == 0 ? Float32 : ComplexF32
      push!(df, (t, rxtime, rssi, pre, ch, fc, fs, len, lno+1, filename, dt))
    end
  end
  df
end

function read(filenames::AbstractVector; tz=localzone())
  df = vcat(read.(filenames; tz)...)
  df.filename = PooledArray(df.filename)
  sort!(df, :time)
end

function read(row::DataFrameRow)
  for (n, line) ∈ enumerate(eachline(row.filename))
    if n == row.lno
      raw = ntoh.(Array{row.dtype}(reinterpret(row.dtype, base64decode(line))))
      x = signal(collect(transpose(reshape(raw, row.channels, :))), row.fs)
      return x .- mean(1.0x; dims=1)
    end
  end
  throw(ErrorException("Signal not found"))
end

read(df::DataFrame, i) = read(df[i,:])

end # module


text:
include("reader.jl")

struct Treebank{T<:Dependency}
    files::Vector{String}
    add_id::Bool
    allow_nonprojective::Bool
    allow_multiheaded::Bool
end

"""
    Treebank{T}(treebank; pattern=r".", add_ad=false, allow_nonprojective=false,kwargs...) whereT

Create a treebank from a `treebank` directory from files matching `pattern`.
"""
function Treebank{T}(treebank::String; pattern = r".", add_id=false,
                     allow_nonprojective=true, allow_multiheaded=true) where T
    if isfile(treebank)
        Treebank{T}([treebank], add_id, allow_nonprojective, allow_multiheaded)
    elseif isdir(treebank)
        treebank_files = String[]
        for (root, dirs, files) in walkdir(treebank), file in files
            if occursin(pattern, file)
                push!(treebank_files, joinpath(root, file))
            end
        end
        Treebank{T}(treebank_files, add_id, allow_nonprojective, allow_multiheaded)
    else
        error("Couldn't read '$treebank'")
    end
end

"""
    Treebank{T}(files; add_id=false, allow_nonprojective=false, kwargs...)

Create a treebank from `files`.
"""
function Treebank{T}(files::Vector{String}; add_id=false, allow_nonprojective=true, allow_multiheaded=true, kwargs...) where T
    Treebank{T}(files, add_id, allow_nonprojective, allow_multiheaded)
end

function Treebank(treebank; kwargs...)
    endswith(treebank, ".conllu") ? Treebank{CoNLLU}(treebank; kwargs...) :
        endswith(treebank, "conll") ? Treebank{CoNLLU}(treebank; kwargs...) :
        error("error reading '$treebank'")
end

deptype(::Type{<:Treebank{T}}) where T = T
deptype(treebank::Treebank) = deptype(typeof(treebank))

trees(treebank::Treebank) = TreebankIterator(treebank)

function Base.iterate(treebank::Treebank)
    iter = TreebankIterator(treebank)
    graph, state = iterate(iter)
    return (graph, (iter, state))
end
function Base.iterate(treebank::Treebank, state)
    (iter, st) = state
    next = iterate(iter, st)
    if next == nothing
        return nothing
    else
        graph, next_state = next
        return (graph, (iter, next_state))
    end
end

Base.IteratorSize(treebank::Treebank) = Base.SizeUnknown()

Base.length(treebank::Treebank) = length(collect(treebank))

function Base.show(io::IO, treebank::Treebank)
    T = deptype(treebank)
    len = length(treebank.files)
    print(io, "Treebank{$T} of $len file(s)")
end

mutable struct TreebankIterator{T}
    t::Treebank{T}
    i::Int
    reader::TreebankReader
end

function TreebankIterator(tb::Treebank)
    r = TreebankReader{deptype(tb)}(first(tb.files), add_id=tb.add_id)
    TreebankIterator(tb, 1, r)
end

function Base.iterate(t::TreebankIterator)
    graph, state = iterate(t.reader)
    if !t.t.allow_nonprojective && !isprojective(graph)
        return iterate(t, (state, 1))
    end
    return (graph, (state, 1))
end
function Base.iterate(t::TreebankIterator, state)
    next = iterate(t.reader, state)
    if next != nothing
        (graph, (st, i)) = next
        if !t.t.allow_nonprojective && !isprojective(graph)
            return iterate(t, (st, i+1))
        end
        return (graph, (st, i+1))
    else
        if t.i < length(t.t.files)
            t.i += 1
            t.reader = TreebankReader{deptype(t.t)}(t.t.files[t.i]; allow_nonprojective=t.t.allow_nonprojective, allow_multiheaded=t.t.allow_multiheaded)
            graph, state = iterate(t.reader)
            return (graph, (state, 1))
        else
            return nothing
        end
    end
end

Base.IteratorSize(t::TreebankIterator) = Base.SizeUnknown()


text:
# module Bukdu

function config{AE<:ApplicationEndpoint}(app::Symbol, ::Type{AE}; kw...)
    # app, endpoint
    merge!(Configuration.env, Dict(map(kw) do kv
        (k,v) = kv
        (k,Dict(v))
    end))
end

function getindex{AE<:ApplicationEndpoint}(::Type{AE}, sym::Symbol)
    Configuration.env[sym]
end

module Configuration
env = Dict{Symbol, Any}()
end # module Bukdu.Configuration


text:
# Julia wrapper for header: /usr/local/include/libavutil/mem.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_malloc,
    av_realloc,
    av_free,
    av_mallocz,
    av_strdup,
    av_freep

function av_malloc(size::Csize_t)
    ccall((:av_malloc,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_realloc(ptr,size::Csize_t)
    ccall((:av_realloc,libavutil),Ptr{Void},(Ptr{Void},Csize_t),ptr,size)
end

function av_free(ptr)
    ccall((:av_free,libavutil),Void,(Ptr{Void},),ptr)
end

function av_mallocz(size::Csize_t)
    ccall((:av_mallocz,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_strdup(s)
    ccall((:av_strdup,libavutil),Ptr{Uint8},(Ptr{Uint8},),s)
end

function av_freep(ptr)
    ccall((:av_freep,libavutil),Void,(Ptr{Void},),ptr)
end


text:
@otestset "Allocations" begin
    # ensure that there are 0 allocations for all the benchmarks
    for b in Benchmarks.all_benchmarks()
        @test Benchmarks.runbenchmark(b, samples = 1, evals = 1).allocs == 0
    end
end # testset Allocations


text:
# kuramoto_7-26-19.jl

# Continuing to read through the Julia documentation.

using BenchmarkTools

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

str = "∀ x ∈ S : x ∈ T";

println("************")
print("Proof this ran again: \t\t\t\t\t")
println(convert(Int64, trunc(rand()*1000)))
println("String to print: ", str)

println("Idiomatic method:")

for s in str
    print(s)
end
# According to @btime: 606.084 μs (63 allocations: 1.22 KiB)

println("\nInefficient but working method:")

for i = firstindex(str):lastindex(str)
    try
        print(str[i])
    catch
        # Do nothing
    end
end
# According to @btime: 691.052 μs (76 allocations: 1.63 KiB).
# Not a huge difference, but we did like a thousand trials of it,
# so it matters.

println("\nDone")

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# One more @benchmark test, because I love that software.

@benchmark begin
    str = "∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀"
    for s in str
        s
    end
end

@benchmark begin
    str = "∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀∀"
    i = 1
    for i = firstindex(str):lastindex(str)
        try
            str[i]
        catch
            # Do nothing
        end
    end
end

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# Fun with numpads.

println("Strings to Arrays, 1")

x = "123"
y = "456"
z = "789"

A = [x y z]

display(A)

A_ = [collect(x) collect(y) collect(z)]

display(A_)
A__ = collect.(A)

display(A__)

A___ = collect.(collect.(A))

display(A___)

# Interestingly, this displays very different behavior than when I
# initialize A using punctuation. What gives?

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

println("Strings to Arrays, 2")

x = "123"
y = "456"
z = "789"

A = [x; y; z]

display(A)

A_ = [collect(x) collect(y) collect(z)]

display(A_)
A__ = collect.(A)

display(A__)

A___ = collect.(collect.(A))

display(A___)

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

println("Strings to Arrays, 3")

x = "123"
y = "456"
z = "789"

A = [x, y, z]

display(A)

A_ = [collect(x) collect(y) collect(z)]

display(A_)
A__ = collect.(A)

display(A__)

A___ = collect.(collect.(A))

display(A___)

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

println(x * y * z)
println(x, y, z)
println("$x$y$z")
println("$x\t$y\t$z")
println("$x\n$y\n$z")

ffp = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

for x in ffp
    for y in ffp
        println("$x\t⋅\t$y\t=\t$(x*y)")
    end
end

print("""

There are $(lastindex(ffp)) prime numbers under 100.
So I wish to memorize $(lastindex(ffp)^2) semiprimes.

At a rate of 10 per day,
this will take me about

$(convert(Int64,round(lastindex(ffp)^2 / 10)))

days.
""")

open("diary\\semiprimes-from-primes-under-100.txt", "w") do file
    ffp = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

    for x in ffp
        for y in ffp
            if x*y < 10000
                write(file, "$x ⋅ $y = $(x*y)\n")
            end
        end
    end
end

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# I think that's it for today.


text:
using Plots

global NPLOTS = 0
pgfplots_flag = true

# function to plot
function plot_bound(f, a, b, xp, yp; existing_plot=nothing, saveflag=false)
	"""
	This function plots f and its overapproximator
		defined by poins xp and yp
	"""
	x = range(a, stop=b, length=1000)
	y = f.(x)
	if pgfplots_flag
		if isnothing(existing_plot)
			p = PGFPlots.Axis([PGFPlots.Linear(x,  y, style="black, thick", mark="none", legendentry=L"f(x)"),
			PGFPlots.Linear(xp, yp, style="blue, thick", mark=:o, legendentry=L"g(x)")],
			legendPos="east")
			xlabel!(p, L"x")
			#display(p)
			global NPLOTS
			NPLOTS += 1
			# PGFPlots.save("plots/bound_"*string(NPLOTS)*".pdf")
			# PGFPlots.save("plots/bound_"*string(NPLOTS)*".tex") # relative to run dir of top level file...dumb...
			#return p
		else
			# push!(existingAxis.plots, new_plot)
			plot!(existing_plot, x,  y, color="red", linewidth=2, label=L"f(x)");
			plot!(existing_plot, xp, yp, color="blue", marker=:o, linewidth=2, label=L"g(x)", legend=:right);
			xlabel!(existing_plot, L"x")
			display(existing_plot);
			global NPLOTS
			NPLOTS += 1
			# savefig("plots/bound_"*string(NPLOTS)*".pdf")
			#return existing_plot
		end
	else
		if isnothing(existing_plot)
			p = plot(x,  y, color="red", linewidth=2, label=L"f(x)");
			p = plot!(p, xp, yp, color="blue", marker=:o, linewidth=2, label=L"g(x)", legend=:right);
			xlabel!(p, L"x")
			display(p)
			global NPLOTS
			NPLOTS += 1
			#savefig("plots/bound_"*string(NPLOTS)*".pdf") # relative to run dir of top level file...dumb...
			#return p
		else
			plot!(existing_plot, x,  y, color="red", linewidth=2, label=L"f(x)");
			plot!(existing_plot, xp, yp, color="blue", marker=:o, linewidth=2, label=L"g(x)", legend=:right);
			xlabel!(existing_plot, L"x")
			display(existing_plot);
			global NPLOTS
			NPLOTS += 1
			if saveflag
				savefig(existing_plot, "plots/bound_"*string(NPLOTS)*".html")
			end
			#return existing_plot
		end
	end

	print("\u1b[1F")
end


text:
using PathwiseSampling
using AbstractGPs
using SparseGPs
using Distributions
using LinearAlgebra

k = 3 * (SqExponentialKernel() ∘ ScaleTransform(3.0))
gp = GP(k)
x = rand(9)
fx = gp(x, 0.01)
y = rand(fx)
C = rand(length(x), length(x))
ap = approx_posterior(SVGP(), fx, MvNormal(C'C))

n_feats = Int(1e6)
n_samples = 100

ϕ = prior_basis(k, n_feats)
k_approx(ϕ, x, y) = ϕ(x)'ϕ(y)
k_matrix_approx(ϕ, X, Y) = hcat(ϕ.(RowVecs(X))...)'hcat(ϕ.(RowVecs(Y))...)
# k_matrix_approx(ϕ, X, Y) = reshape(hcat(k_approx.(ϕ, X, Y')...), (length(X), length(Y)))
k_matrix_approx(ϕ, X) = k_matrix_approx(ϕ, X, X)

km_true = kernelmatrix(k, x)
km_approx = k_matrix_approx(ϕ, x) # TODO: too slow?

using Test
@test km_true ≈ km_approx rtol=1e-3 atol=1e-3

n_feats = 10
ϕ = prior_basis(k, n_feats)

f, w = sample_prior_functions(ϕ, n_feats, n_samples)

f_post = sample_posterior_functions(ap, n_feats, n_samples)

using Plots
x_plot = sort(rand(100))
plot(x_plot, hcat(f.(x_plot)...)', color=:blue, linealpha=0.1, label="")
plot!(x_plot, gp, label="True prior")

plot(x_plot, hcat(f_post.(x_plot)...)', label="", color=:red, linealpha=0.2)
plot!(x_plot, ap, color=:green, label="True posterior")

# # 2 dimensional input
x2 = RowVecs(rand(9, 2))
ψ = prior_basis(k, 10, 2)

p = hcat(ψ.(x2)...)

p*p'


text:
using Submodular
using FactCheck

Tol = 1e-3

facts("Combinatorial Sets") do

  context("intersection") do
    A = SetVariable(2)
    B = SetVariable([1, 4, 5])
    C = [1, 10, 11]
    fix!(A, [1, 2])
    fix!(B, [1, 4])
    D = intersect(A, B, C)
    @fact Set(get_elements(D)) == Set([1]) --> true
    @fact Set(D.baseset) == Set([1, 2, 4, 5]) --> true
  end

  context("set difference") do
    A = SetVariable(2)
    B = SetVariable([1, 4, 5])
    fix!(A, [1, 2])
    fix!(B, [1, 4])
    C = setdiff(A, B)
    @fact Set(get_elements(C)) == Set([2]) --> true
    @fact Set(C.baseset) == Set([1, 2, 4, 5]) --> true

    B₁ = setdiff(A, [1, 3])
    @fact Set(get_elements(B₁)) == Set([2]) --> true
    @fact Set(B₁.baseset) == Set([1, 2]) --> true

    B₂ = setdiff([1, 3], A)
    @fact Set(get_elements(B₂)) == Set([3]) --> true
    @fact Set(B₂.baseset) == Set([1, 2]) --> true
  end

  context("union") do
    A = SetVariable(3)
    fix!(A, [1, 2])
    B = union(A, [5])
    @fact Set(get_elements(B)) == Set([1, 2, 5]) --> true
    @fact Set(B.baseset) == Set([1, 2, 3, 5]) --> true

    A₁ = SetVariable(5)
    fix!(A₁, [5])
    B₁ = union(A, A₁)
    @fact Set(get_elements(B₁)) == Set([1, 2, 5]) --> true
    @fact Set(B₁.baseset) == Set([1, 2, 3, 4, 5]) --> true
  end

end


text:
module SACFiles

export HeaderEnum, Header, readsachdr
export AbstractSACData, AbstractSpectrum, AbstractTimeSeries
export EvenTimeSeries, UnevenTimeSeries, AmplitudeSpectrum, ComplexSpectrum, GeneralXY
export readsac

include("header.jl")
include("data_types.jl")
include("io.jl")

end # module


text:
module TestRegressionLGBM

using MLJBase
using Test
using Random: seed!

import LightGBM

## Regression -- shamelessly adapted from the other tests

Nsamples = 3000
seed!(0)
calc_rmse(p, t) = sqrt(sum((p - t) .^ 2) / length(p))

model = LightGBM.MLJInterface.LGBMRegressor(num_iterations=100)

X       = rand(Nsamples, 5)
y       = sqrt.(sum(X .^ 2, dims=2)) # make the targets the L2 norm of the vectors
weights = rand(Nsamples)

# fit once, without weights
train, test              = MLJBase.partition(MLJBase.eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train];)
yhat                     = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

# fit again with weights
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train], weights[train])
yhat_with_weights        = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

rmse                     = calc_rmse(y[test], yhat)
rmse_weights             = calc_rmse(y[test], yhat_with_weights)

@test rmse < 0.05
@test !isapprox(rmse, rmse_weights, rtol=0.1) # check that they differ by at least around 10% with/without weights

# Cache contains iterations counts history
@test cache isa NamedTuple
@test cache.num_boostings_done == [100]

@test isa(report, NamedTuple)

expected_return_type = Tuple{
    LightGBM.LGBMRegression,
    Vector{Any}, # blep
    LightGBM.MLJInterface.LGBMRegressor,
}

@test isa(fitresult, expected_return_type)

# Provided by Anthony Blaom as a simple integration test
X, y = @load_boston;
model = LightGBM.MLJInterface.LGBMRegressor()
yhat = fit!(machine(model, X, y); verbosity=0) |> predict;
scitype(yhat) == AbstractVector{Continuous}
@test rms(yhat, y) < 6  # or do proper out-of-sample test

end # module
true


text:
function binary(n::Int)
	ans=Int[]
	while n>=1

		push!(ans,n%2)
		n=floor(Int,n/2)

	end

	return reverse(ans)
end

function is_decimal_binary_palindrome(n::Int)
	s=string(n)
	if s!=reverse(s)
		return false
	end

	bin=binary(n)

	return bin==reverse(bin)
end

function compute()

	ans= sum(i for i=1:1000000 if is_decimal_binary_palindrome(i))

	return ans
end

println(compute())