text:
min_mass = 0.6
max_mass = 1.4
const g = -9.8
#z = [0;0;1]
sigma = 1.0

mass_resolution = 0.2
#mass_resolution = 0.05
masses = collect(min_mass:mass_resolution:max_mass)

np = 20
thetas = 2*pi * (1/np:1/np:1)
circle_ps = map(theta->[cos(theta);sin(theta);0], thetas)
circle_ws = to_wrench_matrix(circle_ps)

# move to separate attachment points
na = 10
thetas = 2*pi * (1/na:1/na:1)
r_attach = 0.8
attachment_ps = map(theta->r_attach*[cos(theta);sin(theta);0], thetas)
attachment_ws = map(point_to_wrench, attachment_ps)

# set up the occupancy prior
resolution = 0.1

# note resolution
interior_q(x) = norm(x) <= 1 - 0.1 * 3/2

normalized_error(err) = sqrt(sum((err.*[0.5;0.5;1/(max_mass - min_mass)]).^2))

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
export is_anagram, anagram_graph

"""
    is_anagram(w1,w2)
Checks if words `w1` and `w2` are anagrams of each other.

All letters are reduced to lower case so `Act` and `cat` are
considered anagrams.
"""
function is_anagram(w1::String, w2::String)::Bool
    c1 = counter(collect(lowercase(w1)))
    c2 = counter(collect(lowercase(w2)))
    return c1 == c2
end

function make_anagrams(word::String)::Vector{String}
    anags = unique(collect(permutations(word)))
    return join.(anags)
end

"""
    anagram_graph(S, len=0, trim=false)
Create a graph whose vertices are words in `S` of length `len` in which two words are
adjacent if they are anagrams of each other.  If `len` is zero, all words are used.
If `trim` is true, vertices of degree zero are removed.

Note that every component of this graph is a clique.
"""
function anagram_graph(
    S::Set{String},
    len::Int = 0,
    trim::Bool = false,
)::SimpleGraph{String}
    G = _bare_graph(S, len)

    # add edges

    VV = deepcopy(G.V)
    n = NV(G)
    PM = Progress(n)

    while length(VV) > 0
        next!(PM)
        v = first(VV)
        anas = [ana for ana ∈ make_anagrams(v) if ana ∈ G.V]
        nana = length(anas)
        for i = 1:nana-1
            for j = i+1:nana
                add!(G, anas[i], anas[j])
            end
            next!(PM)
        end
        for a in anas
            delete!(VV, a)
        end
    end

    # remove vertices of degree zero (if trim is true)
    if trim
        for w in G.V
            if deg(G, w) == 0
                delete!(G, w)
            end
        end
    end

    if len < 1
        name(G, "Anagram Graph on all words")
    else
        name(G, "Anagram Graph on $len-letter words")
    end

    return G

end

ext: jl


text:
"LJH is a module for working with LJH files. The intended interface is to use `LJHGroup` exclusivley, `LJHFile` is for internal use only.
`ljh=LJHGroup(filename)` will open one file. If you instead pass a vector of filenames from the same channel, you will open all the files
and be able to access them as if they were one continuous LJH file.
`ljh[1]` returns the first `LJHRecord`. `LJHRecord` has 3 fields `data`, `rowcount`, `timestamp_usec`. If you want all of the pulse data
but no rowcount or timestamp information do do `[r.data for r in ljh]`. If you want just a few pulse records do `collect(ljh[5:10])`.
Alternativley you can use `get_data_rowcount_timestamp` and `get_data_rowcount_timestamp!`.
Use `record_nsamples`, `pretrig_nsamples`, `frametime`, `filenames`, `lengths`, `column`, `row`, `num_columns`, `num_rows` to access
additional information about the LJH file. If you really need to get access to extra information in the header you can access
`ljh.ljhsfiles[1].headerdict`.
`LJH.writeljhheader` and `LJH.writeljhdata` can be used to write LJH files."
module LJH

export LJHGroup, channel, record_nsamples, pretrig_nsamples, frametime, filenames, lengths, column, row, num_columns, num_rows, get_data_rowcount_timestamp

function ljh_get_header_dict(io::IO)
    headerdict = Dict()
    while true
        line = readline(io)
        if startswith(line, "#End of Header")
            break
        elseif startswith(line, "#")
            continue
        else
            splitline=split(line,":")
            if length(splitline) != 2
              continue
            end
            a,b = splitline
            headerdict[strip(a)]=strip(b)
        end
    end
    headerdict
end

type LJHFile{VersionInt, T<:IO}
    filename             ::AbstractString   # filename
    io               ::IO         # IO to read from LJH file
    headerdict       ::Dict             # LJH file header data
    datastartpos     ::Int
    frametime        ::Float64          # sample spacing (microseconds)
    pretrig_nsamples ::Int64            # nPresample
    record_nsamples  ::Int64            # number of sample per record
    channum          ::Int16            # channel number
    column           ::Int16            # column number
    row              ::Int16            # row number
    num_columns      ::Int16            # number of rows
    num_rows         ::Int16            # number of columns
end
immutable LJHRecord
    data::Vector{UInt16}
    rowcount::Int64
    timestamp_usec::Int64
end

@enum LJHVERSION LJH_21 LJH_22
VERSIONS = Dict(v"2.1.0"=>LJH_21, v"2.2.0"=>LJH_22)

LJHFile(fname::AbstractString) = LJHFile(fname,open(fname,"r"))
function LJHFile(fname::AbstractString,io::IO)
    headerdict = ljh_get_header_dict(seekstart(io))
    datastartpos = position(io)
    ioend = position(seekend(io))
    version = VersionNumber(headerdict["Save File Format Version"])
    version in keys(VERSIONS) || error("$fname has version $version, which is not in VERSIONS $VERSIONS")
    versionint = VERSIONS[version]
    record_nsamples = parse(Int,headerdict["Total Samples"])
    num_rows = parse(Int16, get(headerdict,"Number of rows","0"))
    num_columns = parse(Int16, get(headerdict,"Number of columns","0"))
    LJHFile{versionint, typeof(io)}(
        fname, # filename
        io, # IO
        headerdict, #headerdict
        datastartpos, # datastartpos
        parse(Float64,headerdict["Timebase"]), # frametime
        parse(Int16,headerdict["Presamples"]), # pretrig_nsamples
        parse(Int16,headerdict["Total Samples"]), # record_nsamples
        round(Int16,parse(Float64,headerdict["Channel"])), # channel number
        parse(Int16,get(headerdict,"Column number (from 0-$(num_columns-1) inclusive)","0")),   #column, defaults to zero if not in file
        parse(Int16,get(headerdict,"Row number (from 0-$(num_rows-1) inclusive)","0")), #row, defaults to zero if not in file
        num_columns, # num_columns
        num_rows) # num_rows
end
LJHFile(f::LJHFile) = f

"Return number of bytes per record in LJH file, based on version number"
record_nbytes{T}(f::LJHFile{LJH_21,T}) = 6+2*f.record_nsamples
record_nbytes{T}(f::LJHFile{LJH_22,T}) = 16+2*f.record_nsamples

"ljh_number_of_records(f::LJHFile) Return the number of complete records currently available to read from `f`."
function ljh_number_of_records(f::LJHFile)
    oldpos = position(f.io)
    endpos = position(seekend(f.io))
    nbytes = endpos-f.datastartpos
    seek(f.io,oldpos)
    number_of_records = div(nbytes, record_nbytes(f))
end

# support for ljhfile[1:7] syntax
seekto(f::LJHFile, i::Int) = seek(f.io,f.datastartpos+(i-1)*record_nbytes(f))
Base.getindex(f::LJHFile,indexes::AbstractVector)=LJHGroupSlice(LJHGroup(f), indexes)
function Base.getindex(f::LJHFile,index::Int)
    seekto(f, index)
    pop!(f)
end
function Base.pop!{T}(f::LJHFile{LJH_21,T})
    rowcount, timestamp_usec =  record_row_count_v21(read(f.io, UInt8, 6), f.num_rows, f.row, f.frametime)
    data = read(f.io, UInt16, f.record_nsamples)
    LJHRecord(data, rowcount, timestamp_usec)
end
function Base.pop!{T}(f::LJHFile{LJH_22,T})
    rowcount = read(f.io, Int64)
    timestamp_usec = read(f.io, Int64)
    data = read(f.io, UInt16, f.record_nsamples)
    LJHRecord(data, rowcount, timestamp_usec)
end
Base.size(f::LJHFile) = (ljh_number_of_records(f),)
Base.length(f::LJHFile) = ljh_number_of_records(f)
Base.endof(f::LJHFile) = ljh_number_of_records(f)
# access as iterator
Base.start(f::LJHFile) = (seekto(f,1);1)
Base.next(f::LJHFile,j) = pop!(f),j+1
Base.done(f::LJHFile,j) = j==length(f)+1

# open and close
Base.open(f::LJHFile) = open(f.io)
Base.close(f::LJHFile) = close(f.io)

"""Used only for reading LJH version 2.1.0 files. This parses the ugly
"encoded" version of the frame counter, which is converted into an
approximate time, rounded to 4 microseconds, and divided into the integer
and fractional parts of the millisecond. The latter are stored in bytes
3-6, and the former is in byte 1. Ignore byte 2 (it has some need to be
0 for backward compatibility). Ugly! That's why LJH 2.2.0 does something
totally different.
"""
function record_row_count_v21(header::Vector{UInt8}, num_rows::Integer, row::Integer, frame_time::Float64)
    frac = Int64(header[1])
    ms = UInt64(header[3]) |
         (UInt64(header[4])<<8) |
         (UInt64(header[5])<<16) |
         (UInt64(header[6])<<24)
    count_4usec = Int64(ms*250+frac)
    ns_per_frame = round(Int64,frame_time*1e9)
    ns_per_4usec = Int64(4000)
    count_nsec = count_4usec*ns_per_4usec
    count_frame = cld(count_nsec,ns_per_frame)
    rowcount = count_frame*num_rows+row
    return rowcount, 4*count_4usec
end

"""Represent one or more LJHFiles as a seamless sequence that can be addressed
by record number from 1 to the sum of all records in the group."""
type LJHGroup
    ljhfiles::Vector{LJHFile}
    lengths::Vector{Int}
end
function LJHGroup(x::Vector)
    ljhfiles = LJHFile[LJHFile(f) for f in x]
    LJHGroup(ljhfiles, [length(f) for f in ljhfiles])
end
LJHGroup(x::LJHFile) = LJHGroup(LJHFile[x])
LJHGroup(x::AbstractString) = LJHGroup(LJHFile(x))
Base.length(g::LJHGroup) = sum(g.lengths)
Base.close(g::LJHGroup) = map(close, g.ljhfiles)
Base.open(g::LJHGroup) = map(open, g.ljhfiles)
fieldvalue(g::LJHGroup, s::Symbol) = unique([getfield(f, s) for f in g.ljhfiles])
channel(g::LJHGroup) = (assert(length(fieldvalue(g, :channum))==1);g.ljhfiles[1].channum)
record_nsamples(g::LJHGroup) = (assert(length(fieldvalue(g, :record_nsamples))==1);g.ljhfiles[1].record_nsamples)
pretrig_nsamples(g::LJHGroup) = (assert(length(fieldvalue(g, :pretrig_nsamples))==1);g.ljhfiles[1].pretrig_nsamples)
frametime(g::LJHGroup) = (assert(length(fieldvalue(g, :frametime))==1);g.ljhfiles[1].frametime)
column(g::LJHGroup) = (assert(length(fieldvalue(g, :column))==1);g.ljhfiles[1].column)
row(g::LJHGroup) = (assert(length(fieldvalue(g, :row))==1);g.ljhfiles[1].row)
num_columns(g::LJHGroup) = (assert(length(fieldvalue(g, :num_columns))==1);g.ljhfiles[1].num_columns)
num_rows(g::LJHGroup) = (assert(length(fieldvalue(g, :num_rows))==1);g.ljhfiles[1].num_rows)
filenames(g::LJHGroup) = [f.filename for f in g.ljhfiles]
lengths(g::LJHGroup) = g.lengths
"watch(g::LJHGroup, timeout_s) calls `watch_file` on the last LJH file in `g`, passes `timeout_s` through."
watch(g::LJHGroup, timeout_s) = watch_file(last(g.ljhfiles).filename, timeout_s)
"examine the underlying LJHFiles to determine if any grew. If the last one grew, update `g.lengths`, otherwise throw an error"
function update_num_records(g::LJHGroup)
    old_lengths = copy(g.lengths)
    new_lengths = Int[length(f) for f in g.ljhfiles]
    for i = 1:length(g.lengths)-1
        if old_lengths[i]!=new_lengths[i]
            error("a ljh file other than the last file in grew in length $g it was $(g.ljhfiles[i])")
        end
    end
    g.lengths=new_lengths
end
"filenum_recordnum(g::LJHGroup, j::Int) The record `g[j]` is actually `g.ljhfiles[i][k]`, return `i,k`. "
function filenum_recordnum(g::LJHGroup, j::Int)
    for (i,len) in enumerate(g.lengths)
        j <= len ? (return i,j) : (j-=len)
    end
    1,1 # default return value in case of empty range
end
function Base.getindex(g::LJHGroup, i::Int)
    filenum, recordnum = filenum_recordnum(g,i)
    g.ljhfiles[filenum][recordnum]
end
Base.getindex(g::LJHGroup, slice::AbstractArray) = LJHGroupSlice(g, slice)
Base.endof(g::LJHGroup) = length(g)
function Base.start(g::LJHGroup)
    for f in g.ljhfiles seekto(f,1) end
    filenum, recordnum = filenum_recordnum(g,1)
    donefilenum, donerecordnum = filenum_recordnum(g, length(g))
    (filenum, recordnum, donefilenum, donerecordnum)
end
function Base.next(g::LJHGroup, state)
    filenum, recordnum, donefilenum, donerecordnum = state
    ljhrecord = pop!(g.ljhfiles[filenum])
    recordnum+=1
    recordnum > g.lengths[filenum] && (recordnum-=g.lengths[filenum];filenum+=1)
    ljhrecord, (filenum, recordnum, donefilenum, donerecordnum)
end
function Base.done(g::LJHGroup, state)
    filenum, recordnum, donefilenum, donerecordnum = state
    filenum>donefilenum || filenum==donefilenum && recordnum>donerecordnum
end
function Base.show(io::IO, g::LJHGroup)
    print(io, "LJHGroup with $(length(g.ljhfiles)) files, $(length(g)) records, split as $(lengths(g)),")
    print(io, " record_nsampes $(record_nsamples(g)),\n")
    print(io,"  pretrig_nsamples $(pretrig_nsamples(g)).")
    print(io,"channel $(channel(g)), row $(row(g)), column $(column(g)), frametime $(frametime(g)) s.\n")
    print("  First filename $(g.ljhfiles[1].filename)")
end

"LJHGroupSlice is used to allow acces to ranges of LJH records, eg `[r.data for r in ljh[1:100]]`."
immutable LJHGroupSlice{T<:AbstractArray}
    g::LJHGroup
    slice::T
    function LJHGroupSlice(ljhgroup, slice)
        isempty(slice) || maximum(slice)<=length(ljhgroup) || error("$(maximum(slice)) is greater than nrec=$(length(ljhgroup)) in $ljhgroup")
        new(ljhgroup, slice)
    end
end
Base.length(g::LJHGroupSlice) = length(g.slice)
Base.endof(g::LJHGroupSlice) = length(g.slice)
LJHGroupSlice{T<:AbstractArray}(ljhfile::LJHGroup, slice::T) = LJHGroupSlice{T}(ljhfile, slice)
function Base.start{T<:UnitRange}(g::LJHGroupSlice{T})
    for f in g.g.ljhfiles seekto(f,1) end
    isempty(g.slice) && return (2,2,1,1) # ensure done condition is immediatley met on empty range
    filenum, recordnum = filenum_recordnum(g.g, first(g.slice))
    donefilenum, donerecordnum = filenum_recordnum(g.g, last(g.slice))
    seekto(g.g.ljhfiles[filenum], recordnum)
    (filenum, recordnum, donefilenum, donerecordnum)
end
function Base.next{T<:UnitRange}(g::LJHGroupSlice{T}, state)
    filenum, recordnum, donefilenum, donerecordnum = state
    ljhrecord = pop!(g.g.ljhfiles[filenum])
    recordnum+=1
    recordnum > g.g.lengths[filenum] && (recordnum-=g.g.lengths[filenum];filenum+=1)
    ljhrecord, (filenum, recordnum, donefilenum, donerecordnum)
end
function Base.done{T<:UnitRange}(g::LJHGroupSlice{T}, state)
    filenum, recordnum, donefilenum, donerecordnum = state
    filenum>donefilenum || filenum==donefilenum && recordnum>donerecordnum
end

"Get all data from an `LJHGroupSlice`, returned as a tuple of Vectors `(data, rowcount, timestamp_usec)`."
function get_data_rowcount_timestamp(g::LJHGroupSlice)
    data = Array(Vector{UInt16},length(g))
    rowcount = zeros(Int64, length(g))
    timestamp_usec = zeros(Int64, length(g))
    get_data_rowcount_timestamp!(g,data,rowcount,timestamp_usec)
end
"get_data_rowcount_timestamp!(g,data::Vector{Vector{UInt16}},rowcount::Vector{Int64},timestamp_usec::Vector{Int64})
Get all data from an `LJHGroupSlice`, pass in vectors of length `length(g)` and correct type to be filled with the answers."
function get_data_rowcount_timestamp!(g,data::Vector{Vector{UInt16}},rowcount::Vector{Int64},timestamp_usec::Vector{Int64})
    state = start(g)
    i=0
    @assert length(data)==length(rowcount)==length(timestamp_usec)==length(g) "data, rowcount, timestap_usec: length mismatch: lengths $((length(data), length(rowcount), length(timestamp_usec), length(g)))"
    while !done(g, state)
        i+=1
        record, state = next(g,state)
        data[i] = record.data
        rowcount[i] = record.rowcount
        timestamp_usec[i] = record.timestamp_usec
    end
    @assert i==length(g) "iterated $i times, should have been $(length(g))"
    data,rowcount,timestamp_usec
end
"Get all data from an `LJHGroup`, returned as a tuple of Vectors `(data, rowcount, timestamp_usec)`."
get_data_rowcount_timestamp(g::LJHGroup) = get_data_rowcount_timestamp(g[1:end])
function get_data_rowcount_timestamp!(g::LJHGroup,data::Vector{Vector{UInt16}},rowcount::Vector{Int64},timestamp_usec::Vector{Int64})
    get_data_rowcount_timestamp!(g[1:end],data,rowcount, timestamp_usec)
end

"""Write a header for an LJH file."""
function writeljhheader(filename::AbstractString, dt, npre, nsamp; version="2.2.0")
    open(filename, "w") do f
    writeljhheader(f, dt, npre, nsamp; version=version)
    end #do
end

function writeljhheader(io::IO, dt, npre, nsamp; version="2.2.0")
    write(io,
"#LJH Memorial File Format
Save File Format Version: $(version)
Software Version: Fake LJH file generated by Julia
Software Driver Version: n/a
Date: %(asctime)s GMT
Acquisition Mode: 0
Digitized Word Size in bytes: 2
Location: LANL, presumably
Cryostat: Unknown
Thermometer: Unknown
Temperature (mK): 100.0000
Bridge range: 20000
Magnetic field (mGauss): 100.0000
Detector:
Sample:
Excitation/Source:
Operator: Unknown
SYSTEM DESCRIPTION OF THIS FILE:
USER DESCRIPTION OF THIS FILE:
#End of description
Number of Digitizers: 1
Number of Active Channels: 1
Timestamp offset (s): 1304449182.876200
Digitizer: 1
Description: CS1450-1 1M ver 1.16
Master: Yes
Bits: 16
Effective Bits: 0
Anti-alias low-pass cutoff frequency (Hz): 0.000
Timebase: $(dt)
Number of samples per point: 1
Presamples: $(npre)
Total Samples: $(nsamp)
Trigger (V): 250.000000
Tigger Hysteresis: 0
Trigger Slope: +
Trigger Coupling: DC
Trigger Impedance: 1 MOhm
Trigger Source: CH A
Trigger Mode: 0 Normal
Trigger Time out: 351321
Use discrimination: No
Channel: 1.0
Description: A (Voltage)
Range: 0.500000
Offset: -0.000122
Coupling: DC
Impedance: 1 Ohms
Inverted: No
Preamp gain: 1.000000
Discrimination level (%%): 1.000000
#End of Header\n"
    )
end

"""Write LJH file data to an IO object or a filename given a AbstractString."""
function writeljhdata(filename::AbstractString, a...)
    open(filename, "a") do f
    writeljhdata(f,a...)
    end
end

# Write LJH v2.2+ data, with row # and timestamps
function writeljhdata(io::IO,traces::Array{UInt16,2}, rows::Vector{Int64}, times::Vector{Int64})
    for j = 1:length(times)
        writeljhdata(io, traces[:,j], rows[j], times[j])
    end
end
function writeljhdata(io::IO, trace::Vector{UInt16}, row::Int64, time::Int64)
    write(io, row)
    write(io, time)
    write(io, trace)
end

# Write LJH v2.1 data, with row # but no timestamps
function writeljhdata(io::IO,traces::Array{UInt16,2}, rows::Vector{Int64})
    for j = 1:length(rows)
        writeljhdata(io, traces[:,j], rows[j])
    end
end
function writeljhdata(io::IO, trace::Vector{UInt16}, row::Int64)
    timestamp_us = round(Int32, row*0.32) # made-up line rate of 320 nanoseconds per row.
    timestamp_ms = Int32(div(timestamp_us, 1000))
    subms_part = round(UInt8, mod(div(timestamp_us,4), 250))
    dummy_channum = Int8(0)
    write(io, subms_part)
    write(io, dummy_channum)
    write(io, timestamp_ms)
    write(io, trace)
end

end #module


text:
# Quadratic programming in MATLAB-like style
#
#   minimize x^2 + xy + y^2 + yz + z^2
#
#   s.t.    x + 2 y + 3 z >= 4
#           x +   y       >= 1
#

using Xpress

model = xpress_model(
	name = "qp_02",
	f = [0., 0., 0.],
	H = [2. 1. 0.; 1. 2. 1.; 0. 1. 2.],
	A = -[1. 2. 3.; 1. 1. 0.],
	b = -[4., 1.])

optimize(model)

println(Xpress.get_lp_status(model))
println(Xpress.getq(model))

println("sol = $(get_solution(model))")
println("obj = $(get_objval(model))")


text:
module CUDA

using Logging
@Logging.configure(level=ERROR)

include("errors.jl")

include("base.jl")
include("types.jl")
include("devices.jl")
include("context.jl")
include("module.jl")
include("stream.jl")
include("execution.jl")

include("memory.jl")
include("arrays.jl")

include("native/execution.jl")
include("native/intrinsics.jl")
include("native/arrays.jl")

include("profile.jl")

function initialize()
	initialize_api()
end
initialize()

end


text:
function endowments_test()
    e = init_test_endowments(40);
    @test validate_endowments(e);
end

@testset "Endowments" begin
    endowments_test();
end

# ---------------


text:
abstract type MarkableFloat{N}   <: AbstractFloat    end
abstract type MarkableFloat64{N} <: MarkableFloat{N} end
abstract type MarkableFloat32{N} <: MarkableFloat{N} end
abstract type MarkableFloat16{N} <: MarkableFloat{N} end

primitive type MarkableF64{N} <: MarkableFloat64{N} 64 end
primitive type MarkableF32{N} <: MarkableFloat32{N} 32 end
primitive type MarkableF16{N} <: MarkableFloat16{N} 16 end

const MarkableF640 = MarkableF64{0}
const MarkableF641 = MarkableF64{1}
const MarkableF642 = MarkableF64{2}
const MarkableF643 = MarkableF64{3}

const MarkableF320 = MarkableF32{0}
const MarkableF321 = MarkableF32{1}
const MarkableF322 = MarkableF32{2}
const MarkableF323 = MarkableF32{3}

const MarkableF160 = MarkableF16{0}
const MarkableF161 = MarkableF16{1}
const MarkableF162 = MarkableF16{2}
const MarkableF163 = MarkableF16{3}

MarkableF64(x) = MarkableF640(x)
MarkableF32(x) = MarkableF320(x)
MarkableF16(x) = MarkableF160(x)

MarkableF640(x::Float64) = reinterpret(MarkableF640, x)
MarkableF641(x::Float64) = reinterpret(MarkableF641, x)
MarkableF642(x::Float64) = reinterpret(MarkableF642, x)
MarkableF643(x::Float64) = reinterpret(MarkableF643, x)

MarkableF320(x::Float32) = reinterpret(MarkableF320, x)
MarkableF321(x::Float32) = reinterpret(MarkableF321, x)
MarkableF322(x::Float32) = reinterpret(MarkableF322, x)
MarkableF323(x::Float32) = reinterpret(MarkableF323, x)

MarkableF160(x::Float16) = reinterpret(MarkableF160, x)
MarkableF161(x::Float16) = reinterpret(MarkableF161, x)
MarkableF162(x::Float16) = reinterpret(MarkableF162, x)
MarkableF163(x::Float16) = reinterpret(MarkableF163, x)

const IntFloat = Union{Integer, AbstractFloat}

MarkableF640(x::T) where {T<:IntFloat} = MarkableF640(Float64(x))
MarkableF641(x::T) where {T<:IntFloat} = MarkableF641(Float64(x))
MarkableF642(x::T) where {T<:IntFloat} = MarkableF642(Float64(x))
MarkableF643(x::T) where {T<:IntFloat} = MarkableF643(Float64(x))

MarkableF320(x::T) where {T<:IntFloat} = MarkableF320(Float32(x))
MarkableF321(x::T) where {T<:IntFloat} = MarkableF321(Float32(x))
MarkableF322(x::T) where {T<:IntFloat} = MarkableF322(Float32(x))
MarkableF323(x::T) where {T<:IntFloat} = MarkableF323(Float32(x))

MarkableF160(x::T) where {T<:IntFloat} = MarkableF160(Float16(x))
MarkableF161(x::T) where {T<:IntFloat} = MarkableF161(Float16(x))
MarkableF162(x::T) where {T<:IntFloat} = MarkableF162(Float16(x))
MarkableF163(x::T) where {T<:IntFloat} = MarkableF163(Float16(x))

Float64(x::MarkableF640) = reinterpret(Float64, x)
Float64(x::MarkableF641) = reinterpret(Float64, x)
Float64(x::MarkableF642) = reinterpret(Float64, x)
Float64(x::MarkableF643) = reinterpret(Float64, x)

Float32(x::MarkableF320) = reinterpret(Float32, x)
Float32(x::MarkableF321) = reinterpret(Float32, x)
Float32(x::MarkableF322) = reinterpret(Float32, x)
Float32(x::MarkableF323) = reinterpret(Float32, x)

Float16(x::MarkableF160) = reinterpret(Float16, x)
Float16(x::MarkableF161) = reinterpret(Float16, x)
Float16(x::MarkableF162) = reinterpret(Float16, x)
Float16(x::MarkableF163) = reinterpret(Float16, x)

Float64(x::MarkableF320) = Float64(Float32(x))
Float64(x::MarkableF321) = Float64(Float32(x))
Float64(x::MarkableF322) = Float64(Float32(x))
Float64(x::MarkableF323) = Float64(Float32(x))

Float64(x::MarkableF160) = Float64(Float16(x))
Float64(x::MarkableF161) = Float64(Float16(x))
Float64(x::MarkableF162) = Float64(Float16(x))
Float64(x::MarkableF163) = Float64(Float16(x))

Float32(x::MarkableF640) = Float32(Float64(x))
Float32(x::MarkableF641) = Float32(Float64(x))
Float32(x::MarkableF642) = Float32(Float64(x))
Float32(x::MarkableF643) = Float32(Float64(x))

Float32(x::MarkableF160) = Float32(Float16(x))
Float32(x::MarkableF161) = Float32(Float16(x))
Float32(x::MarkableF162) = Float32(Float16(x))
Float32(x::MarkableF163) = Float32(Float16(x))

Float16(x::MarkableF640) = Float16(Float64(x))
Float16(x::MarkableF641) = Float16(Float64(x))
Float16(x::MarkableF642) = Float16(Float64(x))
Float16(x::MarkableF643) = Float16(Float64(x))

Float16(x::MarkableF320) = Float16(Float32(x))
Float16(x::MarkableF321) = Float16(Float32(x))
Float16(x::MarkableF322) = Float16(Float32(x))
Float16(x::MarkableF323) = Float16(Float32(x))

MarkableF641(x::MarkableF640) = reinterpret(MarkableF641, x)
MarkableF641(x::MarkableF641) = x
MarkableF641(x::MarkableF642) = reinterpret(MarkableF641, x)
MarkableF641(x::MarkableF643) = reinterpret(MarkableF641, x)

MarkableF321(x::MarkableF320) = reinterpret(MarkableF321, x)
MarkableF321(x::MarkableF321) = x
MarkableF321(x::MarkableF322) = reinterpret(MarkableF321, x)
MarkableF321(x::MarkableF323) = reinterpret(MarkableF321, x)

MarkableF161(x::MarkableF160) = reinterpret(MarkableF161, x)
MarkableF161(x::MarkableF161) = x
MarkableF161(x::MarkableF162) = reinterpret(MarkableF161, x)
MarkableF161(x::MarkableF163) = reinterpret(MarkableF161, x)

MarkableF641(x::MarkableF320) = MarkableF641(Float64(Float32(x)))
MarkableF641(x::MarkableF321) = MarkableF641(Float64(Float32(x)))
MarkableF641(x::MarkableF322) = MarkableF641(Float64(Float32(x)))
MarkableF641(x::MarkableF323) = MarkableF641(Float64(Float32(x)))

MarkableF641(x::MarkableF160) = MarkableF641(Float64(Float16(x)))
MarkableF641(x::MarkableF161) = MarkableF641(Float64(Float16(x)))
MarkableF641(x::MarkableF162) = MarkableF641(Float64(Float16(x)))
MarkableF641(x::MarkableF163) = MarkableF641(Float64(Float16(x)))

MarkableF321(x::MarkableF640) = MarkableF321(Float32(Float64(x)))
MarkableF321(x::MarkableF641) = MarkableF321(Float32(Float64(x)))
MarkableF321(x::MarkableF642) = MarkableF321(Float32(Float64(x)))
MarkableF321(x::MarkableF643) = MarkableF321(Float32(Float64(x)))

MarkableF321(x::MarkableF160) = MarkableF321(Float32(Float16(x)))
MarkableF321(x::MarkableF161) = MarkableF321(Float32(Float16(x)))
MarkableF321(x::MarkableF162) = MarkableF321(Float32(Float16(x)))
MarkableF321(x::MarkableF163) = MarkableF321(Float32(Float16(x)))

MarkableF161(x::MarkableF640) = MarkableF161(Float16(Float64(x)))
MarkableF161(x::MarkableF641) = MarkableF161(Float16(Float64(x)))
MarkableF161(x::MarkableF642) = MarkableF161(Float16(Float64(x)))
MarkableF161(x::MarkableF643) = MarkableF161(Float16(Float64(x)))

MarkableF161(x::MarkableF320) = MarkableF161(Float16(Float32(x)))
MarkableF161(x::MarkableF321) = MarkableF161(Float16(Float32(x)))
MarkableF161(x::MarkableF322) = MarkableF161(Float16(Float32(x)))
MarkableF161(x::MarkableF323) = MarkableF161(Float16(Float32(x)))

MarkableF642(x::MarkableF640) = reinterpret(MarkableF642, x)
MarkableF642(x::MarkableF641) = reinterpret(MarkableF642, x)
MarkableF642(x::MarkableF642) = x
MarkableF642(x::MarkableF643) = reinterpret(MarkableF642, x)

MarkableF322(x::MarkableF320) = reinterpret(MarkableF322, x)
MarkableF322(x::MarkableF321) = reinterpret(MarkableF322, x)
MarkableF322(x::MarkableF322) = x
MarkableF322(x::MarkableF323) = reinterpret(MarkableF322, x)

MarkableF162(x::MarkableF160) = reinterpret(MarkableF162, x)
MarkableF162(x::MarkableF161) = reinterpret(MarkableF162, x)
MarkableF162(x::MarkableF162) = x
MarkableF162(x::MarkableF163) = reinterpret(MarkableF162, x)

MarkableF642(x::MarkableF320) = MarkableF642(Float64(Float32(x)))
MarkableF642(x::MarkableF321) = MarkableF642(Float64(Float32(x)))
MarkableF642(x::MarkableF322) = MarkableF642(Float64(Float32(x)))
MarkableF642(x::MarkableF323) = MarkableF642(Float64(Float32(x)))

MarkableF642(x::MarkableF160) = MarkableF642(Float64(Float16(x)))
MarkableF642(x::MarkableF161) = MarkableF642(Float64(Float16(x)))
MarkableF642(x::MarkableF162) = MarkableF642(Float64(Float16(x)))
MarkableF642(x::MarkableF163) = MarkableF642(Float64(Float16(x)))

MarkableF322(x::MarkableF640) = MarkableF322(Float32(Float64(x)))
MarkableF322(x::MarkableF641) = MarkableF322(Float32(Float64(x)))
MarkableF322(x::MarkableF642) = MarkableF322(Float32(Float64(x)))
MarkableF322(x::MarkableF643) = MarkableF322(Float32(Float64(x)))

MarkableF322(x::MarkableF160) = MarkableF322(Float32(Float16(x)))
MarkableF322(x::MarkableF161) = MarkableF322(Float32(Float16(x)))
MarkableF322(x::MarkableF162) = MarkableF322(Float32(Float16(x)))
MarkableF322(x::MarkableF163) = MarkableF322(Float32(Float16(x)))

MarkableF162(x::MarkableF640) = MarkableF162(Float16(Float64(x)))
MarkableF162(x::MarkableF641) = MarkableF162(Float16(Float64(x)))
MarkableF162(x::MarkableF642) = MarkableF162(Float16(Float64(x)))
MarkableF162(x::MarkableF643) = MarkableF162(Float16(Float64(x)))

MarkableF162(x::MarkableF320) = MarkableF162(Float16(Float32(x)))
MarkableF162(x::MarkableF321) = MarkableF162(Float16(Float32(x)))
MarkableF162(x::MarkableF322) = MarkableF162(Float16(Float32(x)))
MarkableF162(x::MarkableF323) = MarkableF162(Float16(Float32(x)))

MarkableF643(x::MarkableF640) = reinterpret(MarkableF643, x)
MarkableF643(x::MarkableF641) = reinterpret(MarkableF643, x)
MarkableF643(x::MarkableF642) = reinterpret(MarkableF643, x)
MarkableF643(x::MarkableF643) = x

MarkableF323(x::MarkableF320) = reinterpret(MarkableF323, x)
MarkableF323(x::MarkableF321) = reinterpret(MarkableF323, x)
MarkableF323(x::MarkableF322) = reinterpret(MarkableF323, x)
MarkableF323(x::MarkableF323) = x

MarkableF163(x::MarkableF160) = reinterpret(MarkableF163, x)
MarkableF163(x::MarkableF161) = reinterpret(MarkableF163, x)
MarkableF163(x::MarkableF162) = reinterpret(MarkableF163, x)
MarkableF163(x::MarkableF163) = x

MarkableF643(x::MarkableF320) = MarkableF643(Float64(Float32(x)))
MarkableF643(x::MarkableF321) = MarkableF643(Float64(Float32(x)))
MarkableF643(x::MarkableF322) = MarkableF643(Float64(Float32(x)))
MarkableF643(x::MarkableF323) = MarkableF643(Float64(Float32(x)))

MarkableF643(x::MarkableF160) = MarkableF643(Float64(Float16(x)))
MarkableF643(x::MarkableF161) = MarkableF643(Float64(Float16(x)))
MarkableF643(x::MarkableF162) = MarkableF643(Float64(Float16(x)))
MarkableF643(x::MarkableF163) = MarkableF643(Float64(Float16(x)))

MarkableF323(x::MarkableF640) = MarkableF323(Float32(Float64(x)))
MarkableF323(x::MarkableF641) = MarkableF323(Float32(Float64(x)))
MarkableF323(x::MarkableF642) = MarkableF323(Float32(Float64(x)))
MarkableF323(x::MarkableF643) = MarkableF323(Float32(Float64(x)))

MarkableF323(x::MarkableF160) = MarkableF323(Float32(Float16(x)))
MarkableF323(x::MarkableF161) = MarkableF323(Float32(Float16(x)))
MarkableF323(x::MarkableF162) = MarkableF323(Float32(Float16(x)))
MarkableF323(x::MarkableF163) = MarkableF323(Float32(Float16(x)))

MarkableF163(x::MarkableF640) = MarkableF163(Float16(Float64(x)))
MarkableF163(x::MarkableF641) = MarkableF163(Float16(Float64(x)))
MarkableF163(x::MarkableF642) = MarkableF163(Float16(Float64(x)))
MarkableF163(x::MarkableF643) = MarkableF163(Float16(Float64(x)))

MarkableF163(x::MarkableF320) = MarkableF163(Float16(Float32(x)))
MarkableF163(x::MarkableF321) = MarkableF163(Float16(Float32(x)))
MarkableF163(x::MarkableF322) = MarkableF163(Float16(Float32(x)))
MarkableF163(x::MarkableF323) = MarkableF163(Float16(Float32(x)))


text:
using Statistics, SpecialFunctions

"""
    tnmom2(a, b)

Second moment of the truncated standard normal distribution.
"""
function tnmom2(a::Real, b::Real)
    #return tnmom2c(0, a, b)

    if !(a ≤ b)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)^2
    elseif abs(a) > abs(b)
        return tnmom2(-b, -a)
    elseif isinf(a) && isinf(b)
        return one(middle(a, b))
    elseif isinf(b)
        return 1 + √(2 / π) * a / erfcx(a / √2)
    end

    @assert a < b < Inf && abs(a) ≤ abs(b)
    @assert a ≤ 0 ≤ b || 0 ≤ a ≤ b

    if a ≤ 0 ≤ b
        ea = √(π/2) * erf(a / √2)
        eb = √(π/2) * erf(b / √2)
        fa = ea - a * exp(-a^2 / 2)
        fb = eb - b * exp(-b^2 / 2)
        m2 = (fb - fa) / (eb - ea)
        fb ≥ fa && eb ≥ ea || error("error: a=$a, b=$b")
        0 ≤ m2 ≤ 1 || error("error: a=$a, b=$b")
        return m2
    else # 0 ≤ a ≤ b
        exΔ = exp((a - b)middle(a, b))
        ea = √(π/2) * erfcx(a / √2)
        eb = √(π/2) * erfcx(b / √2)
        fa = ea + a
        fb = eb + b
        m2 = (fa - fb * exΔ) / (ea - eb * exΔ)
        a^2 ≤ m2 ≤ b^2 || error("error: a=$a, b=$b")
        return m2
    end
end

"""
    tnmom2(a, b, μ, σ)

Second moment of the truncated normal distribution, where μ, σ are the mean and
standard deviation of the untruncated distribution.
"""
function tnmom2(a, b, μ, σ)
    if σ > 0
        α = (a - μ) / σ
        β = (b - μ) / σ
        #return σ^2 * tnmom2c(-μ / σ, α, β)
        return μ^2 + σ^2 * tnmom2(α, β) + 2μ * σ * tnmean(α, β)
    elseif iszero(σ) && a ≤ b
        return clamp(μ^2 / one(μ), a, b)
    else
        return oftype(middle(a, b), NaN)
    end
end

"""
    tnmom2i(a, b)

Second moment of the normal distribution with variance -1 and mean 0,
truncated to [a,b].
"""
function tnmom2i(a::Real, b::Real)
    if !(-Inf < a ≤ b < Inf)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)^2
    elseif abs(a) > abs(b)
        return tnmom2i(-b, -a)
    end

    @assert -Inf < a < b < Inf && abs(a) ≤ abs(b)
    @assert a ≤ 0 ≤ b || 0 < a < b

    Δ = (b - a) * middle(a, b)
    exΔ = exp(-Δ)
    da = dawson(a/√2)
    db = dawson(b/√2)
    m = ((da - a/√2)exΔ - (db - b/√2)) / (db - da * exΔ)
    return m
end

"""
    tnmom2i(a, b, μ, σ)

Second moment of the normal distribution with variance -σ^2 and mean μ,
truncated to [a,b].
"""
function tnmom2i(a, b, μ, σ)
    α = (a - μ) / σ
    β = (b - μ) / σ
    return μ^2 + σ^2 * tnmom2i(α, β) + 2μ * σ * tnmom1i(α, β)
end


text:
using Zygote
using LichessBot
using LichessBot.Chess
using DelimitedFiles
using StatsBase
tab = copy(LichessBot.psttable)
filename = "newspeed.fen"
games = @time readdlm(filename, ',')
boards = @time fromfen.(String.(games[:, 1]))
scores = games[:, 2]
data = @time collect(zip(boards, scores))
K = 100

sigmoid(s, K = K) = 1 / (1 + 10^(-K * s / 400))

function err(tup::Tuple, table = tab, K = K)
    tup[2] - sigmoid(LichessBot.pstFull(tup[1], 1, table), K)
end

function err(t::Real, table = tab, K=K)
    -sigmoid(t, K)
end
#=function error(v, table=tab, K=K)
    length(v) <= 2 && return error(tuple(v..., table, K))
    mean([error(t, table, K)^2 for t in v])
end=#

#error(K::Real) = error((startboard(), 0.5), tab, K)

#d = 10
#lastErr = error(K)
atexit(()->println(round.(Int, tab)))
atexit(()->writedlm("src/pstNeu.txt", round.(Int, tab)))
accu = ones(size(tab))
evals = @time [LichessBot.reallySimpleEval(b, 2) for b in boards]
println("length of date: $(length(data))")
for e in length(evals):1
    if abs(evals[e]) >= 10000
        deleteat!(data, e)
    end
end
println("length of data after: $(length(data))")
mut = ReentrantLock()
@time gradient(err, data[1], tab)[2]
while sum(abs.(accu)) > 0.001
    global accu, tab
    global accu = zeros(size(tab))
    @time#= Threads.@threads=# for #=(i, =#t in #=enumerate(=#rand(data, 100000)
        #print("\r$i")
        #print("\r$(err(t, tab))")
        locAcc = gradient(err, t, tab)[2]
        replace!(x->isnan(x) ? zero(x) : x, locAcc)
        lock(mut)
        accu += locAcc
        unlock(mut)
        #println(maximum(accu))
        #sleep(0.1)
        #isnan(accu) && break
    end
    println("\n Summe der Werte: $(sum(accu)), Summe der Beträge: $(sum(abs.(accu)))")
    nancount = count(isnan, accu)
    zerocount = count(x->!iszero(x), accu)
    #@time replace!(accu, NaN=>0)
    @time replace!(x->isnan(x) ? zero(x) : x, accu)
    println("$nancount NaNs replaced, $zerocount neither zero nor NaN\n")
    if !any(isnan, accu)
        tab += accu
    else
        accu = ones(size(tab))
    end
    accu = 0
    #println(K)
end
#println(tab)


text:
#!/usr/bin/env julia

## =========================== ## ## =========================== ##
## 	                         2D  L-system                      	 ##
##	                Ollin Demian Langle Chimal					 ##
##	                         						             ##
##                         09 / 23 / 2018	    		         ##
## =========================== ## ## =========================== ##

"""
    evolve(s,d)

Changes a string s to the value of a dictionary counterpart if it is on it
"""
function evolve(s::String,d::Dict)
    if s ∈ keys(d) |> collect
        s = d[s]
    end
    s
end

"""
    paste_string(a)

Makes a single string from an array of strings
"""
function paste_string(a::Array)
    whole = a[1]
    if length(a) > 1
        for i ∈ 2:length(a)
            whole = whole * a[i]
        end
    end
    whole
end

"""
    lsystem(axiom,rules)

Evolves a L-system from the axiom to the rules given as a dictionary
"""
function lsystem(axiom::String,rules::Dict)
    splitted = split(axiom,"")
    to_string = map(String,splitted)
    evolved_seed = map(x -> evolve(x,rules),to_string) |> paste_string
end

"""
    mult_evol_lsystem(axiom,rules,times)

Evolves a L-system from the axiom to the rules given as a dictionary
and does the same the times given
"""
function mult_evol_lsystem(axiom::String,rules::Dict,times::Int)
    ls = lsystem(axiom,rules)
    if times > 1
        for i ∈ 1:times-1
            ls = lsystem(ls,rules)
        end
    end
    ls
end

"""
    turtle_step(x,y,typ,angle,size)

Makes a L-system turtle step given the current position, an angle, the step size
and an typ which is the axiom which can currently only be F or G and for coloring
"""
function turtle_step(x,y,typ,angle,size=1)
    if typ=="F"
        col = :brown
    elseif typ=="G"
        col = :green
    end
    nx = x + size*cos(angle)
    ny = y + size*sin(angle)
    line = ([y,ny],[x,nx])
    plot!(line,color=col,w=2,leg=false,aspect_ratio=:equal,ticks=nothing,border=false),nx,ny
end

"""
    update_angle(angle,amount,times,direction)

Updates the angle of the turtle step by the times of an amount with a positive or negative
change.
"""
function update_angle(angle,amount,times,direction)
    ex = Expr(:call, Meta.parse(direction), angle, times*amount)
    eval(ex)
end


text:
abstract OrderingStrategy
type LogisticOrdering <: OrderingStrategy end
type PartialCorrOrdering <: OrderingStrategy end
type HDPSOrdering <: OrderingStrategy end
type ManualOrdering <: OrderingStrategy
  indexes::Vector{Int}
  ManualOrdering(indexes) = new(copy(indexes))
end

"""
Defines how the next covariate added to g is chosen

* `ForwardStepwise` - default, tries each unused covariate to find the one that helps the most. Slow, O(p^2) time for p covariates.
* `PreOrdered` - Orders covariates once then adds them in that order. Takes O(p) time for p covariates. Specify with `PreOrdered(ordering)` where `ordering` is
    * `LogisticOrdering()` - ranks covariates based on reduction in log likelihood with intital Q as an offset if added one at a time. Basically the first step of `ForwardStepwise`
    * `PartialCorrOrdering()` - ranks covariates based on partial correlation of each covariate and the residual (y - predicted y from initial Qbar) given a
    * `HDPSOrdering()` - ranks covariates based on bias reduction potential as in the hdps procedure. This ignores intial Qbar.
    * `ManualOrdering(ord)` - orders covars based on specified `indexes`, a `Vector{Int}`, of column indexes of W, not including 1, the intercept.
"""
abstract SearchStrategy
type ForwardStepwise <: SearchStrategy end
type PreOrdered <: SearchStrategy
    ordering::OrderingStrategy
    k_at_once::Int
    covar_order::Vector{Int}
    PreOrdered(ordering::OrderingStrategy, k_at_once=1) = new(ordering, k_at_once, Array(Int, 0))
end
Base.show(io::IO, strat::PreOrdered) = print(io, "Preordered($(strat.ordering))")

#strategies may store some state and we need to clear it out sometimes
init!(strat::SearchStrategy) = strat
function init!(strat::PreOrdered)
    strat.covar_order = Int[]
    strat
end

# added covar(s) are added to used_covars and deleted from unused_covars
# q is fluctuated
#return best_risk, gfit, new_fluc
function add_covars!{T<:AbstractFloat}(::ForwardStepwise,
                                       q::Qmodel{T},
                                       param::Parameter{T},
                                       W::Matrix{T},
                                       A::Vector{T},
                                       Y::Vector{T},
                                       used_covars::IntSet,
                                       unused_covars::IntSet,
                                       prev_risk::T,
                                       gbounds::Vector{T},
                                       penalize::Bool)
    @debug("unused_covars: $(unused_covars)")
    @debug("used_covars: $(used_covars)")

    next_covar_risk = fill(Inf, maximum(unused_covars))
    g_fluc_dict = Dict{IntSet, Tuple{LR, Fluctuation}}()

    #Remove most recent fluctuation (and save for later use)
    #So that we can add one covar to the current fluctuation to see if it works
    fluc_now = defluctuate!(q)

    #try adding each covariate to the current fluctuation
    for j in unused_covars
        #estimate g with additional covariate j
        #println("Vars: $(union(used_covars, IntSet(j)))")
        current_covars = union(used_covars, IntSet(j))
        gfit = lreg(W, A, subset=collect(current_covars)) #sslreg(w, a, current_covars)
        #fluctuate and get risk
        fluc = computefluc(q, param, bound!(predict(gfit, W), gbounds), A, Y)
        fluctuate!(q, fluc)
        next_covar_risk[j] = risk(q, A, Y, param, penalize)
        if isnan(next_covar_risk[j])
            @warn("Risk is NaN. used covars: $used_covars")
        end
        #remove newest fluctuation
        defluctuate!(q)
        g_fluc_dict[current_covars] = (gfit, fluc)
    end
    best_risk, best_j = findmin(next_covar_risk)

    if best_risk < prev_risk
        @debug("Adding to old fluc, best_j:$best_j")
        new_fluc = false
        #adding one covariate to current fluc does better, then use that fluctuation
        push!(used_covars, best_j)
        @debug("used_covars after adding best_j: $used_covars")
        delete!(unused_covars, best_j)
        gfit, fluc = g_fluc_dict[used_covars]
        fluctuate!(q, fluc)
    else
        new_fluc = true
        #otherwise, keep current fluctuation find best new fluctuation adding a
        #single covariate
        fluctuate!(q, fluc_now) #reuse previous best fluctuation
        next_covar_risk = fill(Inf, maximum(unused_covars))
        for j in unused_covars
            current_covars = union(used_covars, IntSet(j))
            gfit = g_fluc_dict[current_covars][1]
            fluc = computefluc(q, param, bound!(predict(gfit, W), gbounds), A, Y)
            fluctuate!(q, fluc)
            next_covar_risk[j] = risk(q, A, Y, param, penalize)
            if isnan(next_covar_risk[j])
                @warn("Risk is NaN. used covars: $used_covars")
            end
            defluctuate!(q)
            g_fluc_dict[current_covars] = (gfit, fluc)
        end
        @debug("next_covar_risk: $next_covar_risk")
        best_risk, best_j = findmin(next_covar_risk)
        if !penalize && best_risk > prev_risk
            @warn("Best risk is worse than previous best risk. This doesn't make sense. Continuing anyway...")
        end
        @debug("best_j: $best_j, unused_covars: $unused_covars")
        push!(used_covars, best_j)
        delete!(unused_covars, best_j)
        gfit, fluc = g_fluc_dict[used_covars]
        fluctuate!(q, fluc)
    end
    best_risk, gfit, new_fluc
end

# added covar(s) are added to used_covars and deleted from unused_covars
# q is fluctuated
#return best_risk, gfit, new_fluc
function add_covars!{T<:AbstractFloat}(strategy::PreOrdered,
                                       q::Qmodel{T},
                                       param::Parameter{T},
                                       W::Matrix{T},
                                       A::Vector{T},
                                       Y::Vector{T},
                                       used_covars::IntSet,
                                       unused_covars::IntSet,
                                       prev_risk::T,
                                       gbounds::Vector{T},
                                       penalize::Bool)

    if isempty(strategy.covar_order)
        append!(strategy.covar_order, order_covars(strategy.ordering, q, param, W, A, Y, unused_covars, gbounds, penalize))
    end

    ordered_unused_covars = filter(x -> x ∉ used_covars, strategy.covar_order)
    next_covars = length(ordered_unused_covars) >= strategy.k_at_once ?
        ordered_unused_covars[1:strategy.k_at_once] :
        ordered_unused_covars

    fluc_now = defluctuate!(q)

    union!(used_covars, next_covars)
    setdiff!(unused_covars, next_covars)

    g_fit = lreg(W, A, subset=collect(used_covars))
    gn1 = bound!(predict(g_fit, W), gbounds)

    fluc = computefluc(q, param, gn1, A, Y)
    fluctuate!(q, fluc)
    new_risk = risk(q, A, Y, param, penalize)

    if new_risk < prev_risk
        return new_risk, g_fit, false
    end

    defluctuate!(q)
    fluctuate!(q, fluc_now)
    fluctuate!(q, param, gn1, A, Y)
    return risk(q, A, Y, param, penalize), g_fit, true
end

function order_covars(ordering::LogisticOrdering, q, param, W, A, Y, available_covars, gbounds, penalize)
    logitQnAA = linpred(q, A)
    scores = Dict{Int, Float64}()
    for i in available_covars
        g_fit = lreg(W, A, subset=[i])
        gn1 = bound!(predict(g_fit, W), gbounds)
        fluc = computefluc(q, param, gn1, A, Y)
        fluctuate!(q, fluc)
        scores[i] = risk(q, A, Y, param, penalize)
        defluctuate!(q)
    end
    sort!(collect(keys(scores)), by = x -> scores[x])
end

function order_covars(ordering::PartialCorrOrdering, q, param, W, A, Y, available_covars, gbounds, penalize)
    resid = Y .- predict(q, A)
    W_available = W[:, collect(available_covars)]
    scores = Dict(collect(zip(available_covars, abs(pcor(resid, W_available, A)))))
    return sort!(collect(keys(scores)), by = x -> scores[x], rev=true)
end

function order_covars(ordering::HDPSOrdering, q, param, W, A, Y, available_covars, gbounds, penalize)
    error()
end

function order_covars(ordering::ManualOrdering, q, param, W, A, Y, available_covars, gbounds, penalize)
    collect(available_covars) == sort(ordering.indexes) ||
      error("ManualOrdering indexes do not match available covariates")
    ordering.indexes
end


text:
############################################################################################
############################################################################################
############################################################################################
# Tempering

@testset "Tempering - Iteration Tempering: " begin
    #Define Parameter
    itertemper_type = Float16
    itertemper_adaption = UpdateTrue()
    itertemper_val = Float32(0.5)
    itertemper_iter = 1000
    ## Adaptive Case
    itertemper_tune = IterationTempering(
        itertemper_type,
        itertemper_adaption,
        itertemper_val,
        itertemper_iter
    )
    temps = BaytesCore.checktemperature(itertemper_tune.parameter, itertemper_iter)
    @test all( temps[iter] >= temps[iter-1] for iter in 2:itertemper_iter)
    @test itertemper_tune.initial.current isa itertemper_type
    @test itertemper_tune.parameter.L isa itertemper_type
    # Update temperature
    @test initial(itertemper_tune) == itertemper_tune.initial.current
    itertemper_tmp = update!(itertemper_tune, 2)
    @test itertemper_tmp > itertemper_tune.initial.current
    ## Fixed stepsize
    itertemper_tune = IterationTempering(
        itertemper_type,
        UpdateFalse(),
        itertemper_val,
        itertemper_iter
    )
    @test itertemper_tune.initial.current isa itertemper_type
    @test itertemper_tune.initial.current == itertemper_val
    @test itertemper_tune.parameter.L isa itertemper_type
    # Update temperature
    @test initial(itertemper_tune) == itertemper_tune.initial.current
    itertemper_tmp = update!(itertemper_tune, 2)
    @test itertemper_tmp == itertemper_tune.initial.current
end

@testset "Tempering - Joint Tempering: " begin
    # Initial Parameter
    jointtemper_type = Float16
    jointtemper_adaption = UpdateTrue()
    jointtemper_unnormalizedweights = [.1, .2, .3, .4, .5]
    jointtemper_weights = jointtemper_unnormalizedweights ./ sum(jointtemper_unnormalizedweights)
    jointtemper_λ₀ = 0.5
    jointtemper_ESS = Float64(length(jointtemper_weights))
    # Test if weights are incrementally increasing
    @test BaytesCore.weightedESS(jointtemper_weights, 1., 1.) == length(jointtemper_weights)
    jointtemper_essdiff = BaytesCore.get_ESSDifference(jointtemper_weights, jointtemper_λ₀, jointtemper_ESS)
    @test jointtemper_essdiff(jointtemper_λ₀ + 0.01) >= 0.0
    jointtemper_newval = update(jointtemper_weights, jointtemper_λ₀, jointtemper_ESS)
    @test jointtemper_newval >= jointtemper_λ₀
    @test update(jointtemper_weights, 1.0, jointtemper_ESS) ≈ 1.0 atol = _TOL
    ## Adaptive joint tempering
    jointtemper_tune = JointTempering(
        jointtemper_type,
        jointtemper_adaption,
        jointtemper_λ₀,
        jointtemper_ESS,
        length(jointtemper_weights)
    )
    initial(jointtemper_tune)
    @test typeof(jointtemper_tune.val.current) == typeof(jointtemper_tune.ESSTarget.current) == jointtemper_type
    @test eltype(jointtemper_tune.weights) == jointtemper_type
    # Check for Updates
    update!(jointtemper_tune, log.(jointtemper_unnormalizedweights))
    @test jointtemper_tune.val.current > jointtemper_λ₀
    @test sum(jointtemper_tune.weights) ≈ 1.0

    ## Fixed joint tempering
    jointtemper_tune = JointTempering(
        jointtemper_type,
        UpdateFalse(),
        jointtemper_λ₀,
        jointtemper_ESS,
        length(jointtemper_weights)
    )
    @test typeof(jointtemper_tune.val.current) == typeof(jointtemper_tune.ESSTarget.current) == jointtemper_type
    @test eltype(jointtemper_tune.weights) == jointtemper_type
    # Check for Updates
    update!(jointtemper_tune, log.(jointtemper_unnormalizedweights))
    @test jointtemper_tune.val.current == jointtemper_λ₀
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1Fields <: SwaggerModel

    function IoK8sApimachineryPkgApisMetaV1Fields(;)
        o = new()
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1Fields

const _property_map_IoK8sApimachineryPkgApisMetaV1Fields = Dict{Symbol,Symbol}()
const _property_types_IoK8sApimachineryPkgApisMetaV1Fields = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1Fields }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1Fields))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1Fields }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1Fields[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1Fields }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1Fields[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1Fields)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1Fields }, name::Symbol, val)
end


text:
"""Mutual information from the joint probability mass function p
of two variables X,Y. p is an NxM array which elements represent
the probabilities of observing x~X and y~Y."""
function mutual_information(p::AbstractArray{T,2},base::Real=2) where T
    @assert sum(p) ≈ one(T)  "Entries in p have to sum to 1"
    @assert all(p .>= zero(T))  "Entries in p have to be >= 1"

    nx,ny = size(p)

    mpx = sum(p,dims=2)     # marginal probabilities of x
    mpy = sum(p,dims=1)

    M = zero(T)
    for j in 1:ny
        for i in 1:nx
            if p[i,j] > zero(T)
                M += p[i,j]*log(p[i,j]/mpx[i]/mpy[j])
            end
        end
    end

    M /= log(base)
end

"""Mutual bitwise information of the elements in input arrays A,B.
A and B have to be of same size and eltype."""
function bitinformation(A::AbstractArray{T},
                        B::AbstractArray{T};
                        set_zero_insignificant::Bool=true,
                        confidence::Real=0.99
                        ) where {T<:Union{Integer,AbstractFloat}}

    @assert size(A) == size(B)
    nelements = length(A)
    nbits = 8*sizeof(T)

    # array of counters
    C = zeros(Int,nbits,2,2)

    for (a,b) in zip(A,B)   # run through all elements in A,B pairwise
        bitcount!(C,a,b)    # count the bits and update counter array C
    end

    # P is the join probabilities mass function
    P = [C[i,:,:]/nelements for i in 1:nbits]
    M = [mutual_information(p) for p in P]

    # remove information that is insignificantly different from a random 50/50
    if set_zero_insignificant
        p = binom_confidence(nelements,confidence)  # get chance p for 1 (or 0) from binom distr
        M₀ = 1 - entropy([p,1-p],2)                 # free entropy of random 50/50 at trial size
        M[M .<= M₀] .= 0                            # set zero insignificant
    end

    return M
end

"""Mutual bitwise information of the elements in input arrays A,B.
A and B have to be of same size and eltype."""
function bitinformation(A::AbstractArray{T},
                        B::AbstractArray{T},
                        n::Int) where {T<:Union{Integer,AbstractFloat}}

    @assert size(A) == size(B)
    nelements = length(A)
    nbits = 8*sizeof(T)

    # array of counters, first dim is from last bit to first
    C = [zeros(Int,2^min(i,n+1),2) for i in nbits:-1:1]

    for (a,b) in zip(A,B)   # run through all elements in A,B pairwise
        bitcount!(C,a,b,n)  # count the bits and update counter array C
    end

    # P is the join probabilities mass function
    P = [C[i]/nelements for i in 1:nbits]
    M = [mutual_information(p) for p in P]

    # filter out rounding errors
    M[isapprox.(M,0,atol=10eps(Float64))] .= 0
    return M
end

"""Update counter array of size nbits x 2 x 2 for every
00|01|10|11-bitpairing in a,b."""
function bitcount!(C::Array{Int,3},a::T,b::T) where {T<:AbstractFloat}
    uia = reinterpret(Base.uinttype(T),a)
    uib = reinterpret(Base.uinttype(T),b)
    bitcount!(C,uia,uib)
end

"""Update counter array of size nbits x 2 x 2 for every
00|01|10|11-bitpairing in a,b."""
function bitcount!(C::Vector{Matrix{Int}},a::T,b::T,n::Int) where {T<:AbstractFloat}
    uia = reinterpret(Base.uinttype(T),a)
    uib = reinterpret(Base.uinttype(T),b)
    bitcount!(C,uia,uib,n)
end

"""Update counter array of size nbits x 2 x 2 for every
00|01|10|11-bitpairing in a,b."""
function bitcount!(C::Array{Int,3},a::T,b::T) where {T<:Integer}
    nbits = 8*sizeof(T)
    mask = one(T) << (nbits-1)
    for i in 1:nbits
        j = 1+((a & mask) >>> (nbits-i))
        k = 1+((b & mask) >>> (nbits-i))
        C[i,j,k] += 1
        mask >>>= 1
    end
end

"""Update counter array of size nbits x 2 x 2 for every
00|01|10|11-bitpairing in a,b."""
function bitcount!(C::Vector{Matrix{Int}},a::T,b::T,n::Int) where {T<:Integer}

    nbits = 8*sizeof(T)

    for i = nbits:-1:1
        @boundscheck size(C[end-i+1],1) == 2^min(i,n+1) || throw(BoundsError())
    end

    maskb = one(T) << (nbits-1)
    maska = reinterpret(T,signed(maskb) >> n)

    for i in 1:nbits
        j = 1+((a & maska) >>> max(nbits-n-i,0))
        k = 1+((b & maskb) >>> (nbits-i))
        C[i][j,k] += 1

        maska >>>= 1    # shift through nbits
        maskb >>>= 1
    end
end

"""Calculate the bitwise redundancy of two arrays A,B. Redundancy
is a normalised measure of the mutual information 1 for always
identical/opposite bits, 0 for no mutual information."""
function redundancy(A::AbstractArray{T},
                    B::AbstractArray{T}) where {T<:Union{Integer,AbstractFloat}}
    mutinf = bitinformation(A,B)    # mutual information
    HA = bitcount_entropy(A)        # entropy of A
    HB = bitcount_entropy(B)        # entropy of B
    R = 2mutinf./(HA+HB)            # redundancy (symmetric)

    R[iszero.(HA+HB)] .= 0.0         # HA+HB = 0 creates NaN
    return R
end

"""Calculate the bitwise redundancy of two arrays A,B.
Multi-bit predictor version which includes n lesser significant bit
as additional predictors in A for the mutual information to account
for round-to-nearest-induced carry bits. Redundancy is normalised
by the entropy of A. To be used for A being the reference array
and B an approximation of it."""
function redundancy(A::AbstractArray{T},
                    B::AbstractArray{T},
                    n::Int) where {T<:Union{Integer,AbstractFloat}}
    mutinf = bitinformation(A,B,n)  # mutual information
    HA = bitcount_entropy(A)        # entropy of A
    R = mutinf./HA                  # redundancy (asymmetric)

    R[iszero.(HA)] .= 0.0           # HA = 0 creates NaN
    return R
end


text:
# Solution to the first problem 1.
import JLD2
import Random
import CSV
import DataFrames
import FreqTables
using Statistics

dfs = DataFrames

"""
Runs the first part of the PS1.
"""
function q1()
    # 1.a
    Random.seed!(1234);

    A = Random.rand(10, 7) .* (10 - (-5)) .+ (-5);
    B = Random.randn(10, 7) .* (-15) .+ (-2);
    C = zeros(5, 7);
    C[1:5, 1:5] = A[1:5, 1:5];
    C[1:5, end - 1:end] = B[1:5, end - 1:end];
    D = A .* (A .≤ 0);

    # Printing out the arrays.
    display(A)
    display(B)
    display(C)
    display(D)

    # 1.b
    display(length(A))

    # 1.c
    display(length(unique(D)))

    # 1.d
    E = reshape(B, :, 1);
    display(E[1:10])

    # 1.e
    F = cat(A, B, dims=3)
    display(F)

    # 1.f
    F = permutedims(F, [3, 2, 1]);
    display(F)

    # 1.g
    G = kron(B, C);
    display(G)
    #H = kron(C, F)  # Doesn't work because they are of different
                     # dimensions

    # 1.h
    @JLD2.save "matrixpractice.jld" A B C D E F G

    # 1.i
    @JLD2.save "firstmatrix.jld" A B C D

    # 1.j
    C_df = DataFrames.DataFrame(C);
    CSV.write("Cmatrix.csv", C_df)

    # 1.k
    D_df = DataFrames.DataFrame(D)
    CSV.write("Dmatrix.dat", D_df; delim="|")

    return A, B, C, D
end

"""
Runs the second part of the PS1.
"""
function q2(A, B, C, D)
    # Part 2.
    # 2.a
    AB = zeros(size(A)...);

    for i = 1:size(A)[1], j = 1:size(A)[2]
        AB[i, j] = A[i, j] * B[i, j];
    end

    AB2 = A .* B;

    display(AB)
    display(AB2)

    # 2.b
    Cprime = [];

    for c in C
        if -5 ≤ c ≤ 5
            push!(Cprime, c);
        end
    end

    Cprime2 = C[-5 .≤ C .≤ 5];

    display(Cprime)
    display(Cprime2)

    # 2.c
    N = 15_169;
    K = 6;
    T = 5;

    X = zeros(N, K, T);

    # Creating the binornd since it is not in Julia base.
    binord(;n=20, p=0.6) = sum(Random.rand(n) .≤ p);

    for t in 1:T
        X[:, 1, t] .= 1;
        X[:, 2, t] = Random.rand(N, 1) .≤ 0.75 * (6 - t) / 5;
        X[:, 3, t] = Random.randn(N, 1) .* 5(t - 1) .+ (15 + t - 1);
        X[:, 4, t] = Random.randn(N, 1) .* (1 / ℯ) .+ π * (6 - t) / 3;
        X[:, 5, t] = [binord() for _ in 1:N];
        X[:, 6, t] = [binord(p=0.5) for _ in 1:N];
    end

    # 2.d
    β = zeros(K, T);
    β[1, :] = [1 + 0.25(t - 1) for t in 1:T];
    β[2, :] = [log(t) for t in 1:T];
    β[3, :] = [-sqrt(t) for t in 1:T];
    β[4, :] = [exp(t) - exp(t + 1) for t in 1:T];
    β[5, :] = [t for t in 1:T];
    β[6, :] = [t / 3 for t in 1:T];

    # 2.e
    ε = Random.randn(N, T) .* 0.36;
    Y = [sum(X[n, :, t] .* β[:, t]) + ε[n, t] for n in 1:N, t in 1:T];

    return nothing
end

"""
Runs the third part of the PS1.
"""
function q3()
    # Part 3.
    # 3.a
    data = CSV.read("nlsw88.csv");
    @JLD2.save "nlsw88.jld" data

    # 3.b
    perc_no_marriage = mean(data["married"]) * 100.0;
    println("Percentage that hasn't been married: $perc_no_marriage %")

    # 3.c
    # Alternative way.
    nobs = dfs.nrow(data);
    races = dfs.combine(dfs.groupby(data, :race), dfs.nrow => :count);
    races = dfs.transform(races, :count => x -> x * 100 / nobs);
    races = dfs.rename(races, :count_function => "percentage (%)");
    display(races)

    # Expected way (?).
    races = (FreqTables.prop ∘ FreqTables.freqtable)(data["race"]);
    display(races)

    # 3.d
    summarystats = dfs.describe(data);
    missing_grade = summarystats[summarystats["variable"] .== :grade, :]["nmissing"][1];
    println("Missing observations in the grade columns: $missing_grade")

    # 3.e
    industry_occupation = FreqTables.freqtable(data, :industry, :occupation);
    display(industry_occupation)

    # 3.f
    tab_iow = dfs.combine(dfs.groupby(data[["industry", "occupation", "wage"]],
                                      [:industry, :occupation]),
                          :wage => mean);
    display(tab_iow)
end

"""
Runs the fourth part of the PS1.
"""
function q4()
    data = CSV.read("nlsw88.csv");
    # Part 4
    # 4.a
    @JLD2.load "firstmatrix.jld" A B C D

    # 4.b, 4.c and 4.e
    """
    Takes as input a matrix A and B and returns three outputs

    i) the elemnt-by-element product of the inputs
    ii) the product A'B
    iii) the sum of all the elements of A + B
    """
    function matrixops(A, B)
        if size(A) != size(B)
            throw(DimensionMismatch("size of A not equal to size of B"))
        end
        return A .* B, A' * B, sum(A + B)
    end

    # 4.d
    display(matrixops(A, B))

    # 4.f
    try
        display(matrixops(C, D))
    catch e
        println("size of A not equal to size of B")
    end

    # 4.g
    # No need to convert it into an array.
    display(matrixops(data["ttl_exp"], data["wage"]))
end

A, B, C, D = q1();
q2(A, B, C, D);
q3();
q4();


text:
using Documenter
using LaplacianOpt
using DocumenterTools: Themes

for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "src/assets/themes/style.scss"), String)
    theme = read(joinpath(@__DIR__, "src/assets/themes/$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "src/assets/themes/$(w).scss"), header*"\n"*theme)
end

Themes.compile(joinpath(@__DIR__, "src/assets/themes/light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "src/assets/themes/dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

makedocs(
    modules = [LaplacianOpt],
    sitename = "LaplacianOpt",
    format = Documenter.HTML(mathengine = Documenter.MathJax(),
                             prettyurls = get(ENV, "CI", nothing) == "true"),
    strict = true,
    authors = "Harsha Nagarajan",
    pages = [
        "Introduction" => "index.md",
        "Quick Start guide" => "quickguide.md",
        "Function References" => "function_references.md",
    ],
)

deploydocs(
    repo = "github.com/harshangrjn/LaplacianOpt.jl.git",
    push_preview = true
)


text:
#=
Settings: tests for `CommandLineEquationsOfState.Settings` module
- Julia version: 1.0
- Author: singularitti
- Date: 2019-08-27
=#
using Test

@testset "Test 1 dataset with multiple equations of state" begin

end # testset


text:
using Zygote

@doc raw"
    einsum_grad(::EinCode{ixs, iy}, xs, size_dict, cdy, i)

return the gradient of the result of evaluating the `EinCode` w.r.t
the `i`th tensor in `xs`. `cdy` is the result of applying the `EinCode`
to the `xs`.

# example
```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: einsum_grad, get_size_dict

julia> a, b = rand(2,2), rand(2,2);

julia> c = einsum(EinCode((('i','j'),('j','k')), ('i','k')), (a,b));

julia> sd = get_size_dict((('i','j'),('j','k')), (a,b));

julia> einsum_grad(EinCode((('i','j'),('j','k')), ('i','k')), (a,b), sd, c, 1) ≈ c * transpose(b)
true
```
"
function einsum_grad(::EinCode{ixs, iy}, xs, size_dict, cdy, i) where {ixs, iy}
    nixs = TupleTools.insertat(ixs, i, (iy,))
    nxs  = TupleTools.insertat( xs, i, (cdy,))
    niy = ixs[i]
    y = conj!(einsum(EinCode(nixs, niy), nxs, size_dict))
    typeof(y) == typeof(xs[i]) && return y
    xs[i] isa Array{<:Real} && return convert(typeof(xs[i]), real(y))
    convert(typeof(xs[i]), y)
end

@Zygote.adjoint function einsum(code::EinCode{ixs, iy}, xs::NTuple{N,T where T}, size_dict::IndexSize) where {N, ixs, iy}
    y = einsum(code, xs, size_dict)
    return y, dy -> let cdy = map(conj,dy)
                (
                    nothing,
                    ntuple(i -> einsum_grad(code, xs, size_dict, cdy, i), N),
                    nothing
                )
            end
end

@doc raw"
    bpcheck(f, args...; η = 1e-5, verbose=false)
returns a `Bool` indicating whether Zygote calculates the gradient of `f(args...) -> scalar`
correctly using the relation `f(x - ηg) ≈ f(x) - η|g|²` with a relative tolerance
of 1e-2 and an absolute tolerance of 1e-8.
If `verbose=true`, print `f(x) - f(x - ηg)`and `η|g|²`.

# example

```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: bpcheck

julia> a, b = rand(2,2), rand(2,2);

julia> bpcheck(sum ∘ ein\"ij,jk -> ik\", a, b)
true
```
"
function bpcheck(f, args...; η = 1e-5, verbose = false)
    g = gradient(f, args...)
    all(==(nothing), g) && error()
    dy_ref = 0
    for x in g
        x === nothing && continue
        x isa Tuple && (dy_ref += η * mapreduce(y -> y == nothing ? 0 : sum(abs2,y), +, x))
        x isa AbstractArray && (dy_ref += η * sum(abs2,x))
    end
    dy = f(args...) - f([gi == nothing ? arg : arg .- η .* gi for (arg, gi) in zip(args,g)]...)

    verbose && @show dy
    verbose && @show dy_ref

    isapprox(dy, dy_ref, rtol=1e-2, atol=1e-8)
end

@Zygote.nograd get_size_dict


text:
#!/usr/bin/env julia

# include("common.jl")

function analyze(host=gethostname(), dst=nothing, suite="julia_cuda")
    measurements = CSV.read("measurements_$host.dat")
    grouped = summarize(measurements)

    # add device totals for each benchmark
    grouped_kernels = filter(entry->!startswith(entry[:target], '#'), grouped)
    by(grouped_kernels, [:benchmark, :target]) do df
        if df[df[:suite] .== baseline, :target_iterations] != df[df[:suite] .== suite, :target_iterations]
            @error "Kernel $(first(df[:target])) of $(first(df[:benchmark])) not executed identically across suites"
        end
    end
    grouped = vcat(grouped,
                   by(grouped_kernels, [:suite, :benchmark],
                      dt->DataFrame(target = "#device",
                                    time = sum(dt[:time]),
                                    target_iterations = missing,
                                    benchmark_iterations = allequal(dt[:benchmark_iterations]))))

    # XXX: add compilation time to host measurements
    # grouped[(grouped[:suite].==suite) .& (grouped[:target].=="#host"), :time] +=
    #     grouped[(grouped[:suite].==suite) .& (grouped[:target].=="#compilation"), :time]

    # select values for the requested suite, adding ratios against the baseline
    analysis = by(grouped, [:benchmark, :target]) do df
        baseline_data = df[df[:suite] .== baseline, :]
        suite_data = df[df[:suite] .== suite, :]

        entry_or_missing(data, col) = size(data, 1) == 1 ? data[1, col] : missing
        get_baseline(col) = entry_or_missing(baseline_data, col)
        get_suite(col) = entry_or_missing(suite_data, col)

        DataFrame(reference = get_baseline(:time), time = get_suite(:time),
                  ratio = get_suite(:time) / get_baseline(:time),
                  target_iterations = get_suite(:target_iterations),
                  benchmark_iterations = get_suite(:benchmark_iterations))
    end

    # calculate ratio grand totals
    geomean(x) = prod(x)^(1/length(x))  # ratios are normalized, so use geomean
    totals = by(filter(row->startswith(row[:target], '#'), analysis), [:target]) do df
        DataFrame(reference = sum(df[:reference]), time = sum(df[:time]),
                  ratio = geomean(df[:ratio]),
                  target_iterations = missing,
                  benchmark_iterations = missing)
    end
    totals = hcat(DataFrame(benchmark=fill("#all", size(totals,1))), totals)
    analysis = vcat(analysis, totals)

    @info "Analysis complete"
    println(filter(row->startswith(row[:target], '#'), analysis))

    # prepare measurements for PGF
    analysis[:speedup] = 1 - analysis[:ratio]
    function decompose(df)
        df = copy(df)
        for column in names(df)
            if eltype(df[column]) <: Union{Measurement, Union{Missing,<:Measurement}}
                df[Symbol("$(column)_val")] = map(datum->ismissing(datum)?missing:datum.val,
                                                        df[column])
                df[Symbol("$(column)_err")] = map(datum->ismissing(datum)?missing:datum.err,
                                                        df[column])
                delete!(df, column)
            end
        end
        return df
    end

    # per-benchmark totals
    let analysis = filter(row->!startswith(row[:benchmark], '#') && startswith(row[:target], '#'), analysis)
        # move from time and reference time for every benchmark x target,
        # to different time columns for every benchmark
        # (easier for PGF to plot)
        analysis = DataFrame(
            benchmark            = unique(analysis[:benchmark]),
            kernels              = analysis[analysis[:target] .== "#cudanative", :target_iterations],
            cuda_host            = analysis[analysis[:target] .== "#host",       :reference],
            cuda_device          = analysis[analysis[:target] .== "#device",     :reference],
            julia_host           = analysis[analysis[:target] .== "#host",       :time],
            julia_device         = analysis[analysis[:target] .== "#device",     :time],
            jit_julia            = analysis[analysis[:target] .== "#julia",      :time],
            jit_cudanative       = analysis[analysis[:target] .== "#cudanative", :time],
        )
        sort!(analysis, cols=:julia_host)
        if dst != nothing
            CSV.write(joinpath(dst, "perf.csv"), decompose(analysis); header=true)
        end

        analysis[:jit_julia] = analysis[:jit_julia] ./ 1000
        analysis[:jit_cudanative] = analysis[:jit_cudanative] ./ 1000
        @info("JIT compilation",
              data=analysis[[:benchmark, :kernels, :jit_julia, :jit_cudanative]],
              julia=mean(analysis[:jit_julia]), cudanative=mean(analysis[:jit_cudanative]))
    end

    # per-benchmark kernel timings
    let analysis = filter(row->!startswith(row[:benchmark], '#') && row[:target] == "#device", analysis)
        analysis = analysis[[:benchmark, :speedup]]
        sort!(analysis, cols=:speedup; rev=true)
        if dst != nothing
            CSV.write(joinpath(dst, "perf_device.csv"), decompose(analysis); header=true)
        end
    end

    # total kernel speedup
    let analysis = analysis[(analysis[:benchmark] .== "#all") .& (analysis[:target] .== "#device"), :speedup]
        if dst != nothing
            open(joinpath(dst, "perf_device_total.csv"), "w") do io
                println(io, analysis[1].val)
            end
        end
    end

    return
end


text:
"""
Node is the basic structural datatype of a Binary tree
"""
struct Node{T}
    payload::T
    left::Union{Node, Nothing}
    right::Union{Node, Nothing}

    function Node(payload)
        new{typeof(payload)}(payload, nothing, nothing)
    end

    function Node(payload, l::Node, r::Node)
        new{typeof(payload)}(payload, l, r)
    end
end

"""
Encodes the traversal and visitation of a Binayr tree.

Note that the actual visitor is the function visitor_fn. This is one
of the nice things about Julia's runtime multiple dispatch.
"""
function traverse(n::Node, visitor_fn)
    visitor_fn(n.payload)
    if !isnothing(n.left)
        traverse(n.left, visitor_fn)
    end
    if !isnothing(n.right)
        traverse(n.right, visitor_fn)
    end
end

function myprint(value::String)
    println("I am a Sting $value")
end

function myprint(value::Int)
    println("I am an int $value")
end

l = Node(34)
r = Node("Hi there")
r = Node(1, l, r)

traverse(r, myprint)

max_stars_repo_name: julitopower/JuliaLangAdventures

max_issues_repo_name: julitopower/JuliaLangAdventures

max_forks_repo_name: julitopower/JuliaLangAdventures


text:
import SemiparametricAUC: simulate_one_predictor

ds_betas = @time simulate_one_predictor(iter = 500, m = 120, p = 100)
@time simulate_one_predictor(iter = 1000, m = 220, p = 200)
@time simulate_one_predictor(iter = 1000, m = 400, p = 300)
@time simulate_one_predictor(iter = 100000, m = 120, p = 100)

meanbeta = DataFrames.colwise(mean, DataFrames.DataFrame(ds_betas[1])) # mean betas
meanvar = DataFrames.colwise(var, DataFrames.DataFrame(ds_betas[2]))  # mean variances
meansd = DataFrames.colwise(sqrt, DataFrames.DataFrame(ds_betas[2]))

#Calculating 95% CI coverage
low_bo  = 0.15 .>= ds_betas[3][:,1]
high_bo = 0.15 .<= ds_betas[3][:,2]
# sum(low_bo & high_bo)/ds_betas[4]

low_b1  = 0.50 .>= ds_betas[3][:,3]
high_b1 = 0.50 .<= ds_betas[3][:,4]
# sum(low_b1 & high_b1)/ds_betas[4]

low_b2  = 1.00 .>= ds_betas[3][:,5]
high_b2 = 1.00 .<= ds_betas[3][:,6]
# sum(low_b2 & high_b2)/ds_betas[4]

all_coverage = vcat(sum(low_bo & high_bo)/ds_betas[4],sum(low_b1 & high_b1)/ds_betas[4],sum(low_b2 & high_b2)/ds_betas[4])


text:
function exchange_accept_reject_entry!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    γ_curr::Float64,
    i::Int,
    mcmc::T,
    posterior::MultigraphSnfPosterior,
    aux_data::Vector{Vector{Int}},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    ) where {T<:Union{SnfExchangeRandGibbs}}

    aux_mcmc = mcmc.aux_mcmc
    data = posterior.data
    dist = posterior.dist

    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:mcmc.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability

    # Auxiliary sampling
    G_prop = vec_to_adj_mat(
        x_prop,
        directed=posterior.directed, self_loops=posterior.self_loops
    )
    aux_model = SNF(G_prop, γ_curr)

    if aux_init_at_prev
        tmp = copy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = γ_curr * (
        sum(x->dist(x,x_curr), aux_data) - sum(x->dist(x,x_prop), aux_data)
    )

    suff_stat_prop = sum(x->dist(x,x_prop), data)

    log_lik_ratio = γ_curr * (suff_stat_curr - suff_stat_prop)

    log_prior = G_prior.γ * (dist(x_curr, G_prior.mode) - dist(x_prop, G_prior.mode))

    log_α = log_lik_ratio + log_prior + aux_log_lik_ratio

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        return true, suff_stat_prop
    else
        @inbounds x_prop[i] = w_curr
        return false, suff_stat_curr
    end
end

# Optimised version for the hamming distances between mutlisets
function accept_reject_entry_hamming!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    i::Int,
    mcmc::T,
    model::MultigraphSNF,
    x_mode::Vector{Int} # Must also pass vectorised mode
    ) where {T<:Union{SnfMcmcRandGibbs,SnfMcmcSysGibbs}}

    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:mcmc.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability
    γ, d = (model.γ, model.d)
    @inbounds w_mode = x_mode[i]
    log_α = γ * (abs(w_curr - w_mode) - abs(w_prop - w_mode))

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        return true
    else
        @inbounds x_prop[i] = w_curr
        return false
    end
end


text:
include("roblib.jl")

function main(r::Robot)

    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    side = North

    while side != "border"
        side = sneak!(r, side)
    end

    for i in (North, East, South, West)
        while isborder(r, clockwise_side(i) )
            putmarker!(r)
            move!(r, i)
        end
        putmarker!(r)
        move!(r, clockwise_side(i) )
    end

    go_to_west_south_corner_and_return_path!(r)
    go_by_path!(r, path)
end

function sneak!(r, side)

    while !isborder(r, side)
        (isborder(r, East)) ? break : move!(r,side)
    end
    if isborder(r, East)
        my_ans = "border"
    else
        move!(r, East)
        my_ans = inverse_side(side)
    end

    return my_ans
end

r=Robot( "6.sit" ,animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
using Pkg
using Pkg.TOML
Pkg.activate(@__DIR__)

# ---------------------------------------------------------------------------------------------
# Parsing ARGS

try; @eval import ArgParse
catch
    Pkg.add("ArgParse")
    @eval import ArgParse
end
const AP = ArgParse

argset = AP.ArgParseSettings()
AP.@add_arg_table! argset begin
    "--deep", "-d"
        help = "define how many extra version of each pkg to install, " *
                "from newer to older."
        arg_type = Int
        default = 0
    "--root", "-r"
        help = "The root folder for searching projects."
        arg_type = String
        default = "."
    "--dry-run"
        help = "Run without consequences."
        action = :store_true
end

if isinteractive()
    # Dev values
    deep = 0
    dry_run = true
    root = abspath(".")
else
    parsed_args = AP.parse_args(ARGS, argset)
    deep = parsed_args["deep"]
    dry_run = parsed_args["dry-run"]
    root = abspath(parsed_args["root"])
end

## ---------------------------------------------------------------------------------------------
# Tools
function load_deps(proj_file)
    raw_deps = get(TOML.parsefile(proj_file), "deps", Dict())
    deps = Dict()
    for (name, uuid) in raw_deps
        deps[uuid] = name
    end
    deps
end

function find_projects(root = ".")
    projs = String[]
    c = 1
    root = abspath(root)
    for (dir, subdirs, _) in walkdir(root)
        for subdir in subdirs, projname in Base.project_names
            try
                projfile = joinpath(dir, subdir, projname)
                !isfile(projfile) && continue

                # ignores
                occursin("/.git/", projfile) && continue

                push!(projs, projfile)
                println(c, ": ", relpath(projfile, root))
                c += 1
            catch err; end
        end
    end
    projs
end

reg_dir() = joinpath(first(DEPOT_PATH),"registries")

function findin_regs(name::String, reg_root = reg_dir())
    firstletter = string(first(name))
    founds = []
    for dir in readdir(reg_root)
        reg_dir = joinpath(reg_root, dir, firstletter, name)
        !isdir(reg_dir) && continue
        push!(founds, reg_dir)
    end
    founds
end

function extract_versions(regpath)
    ctx = Pkg.Types.Context()
    di=Pkg.Operations.load_versions(ctx, regpath)
    keys(di) |> collect
end

function extract_uuid(regpath)
    pkg_file = joinpath(regpath, "Package.toml")
    dat = TOML.parsefile(pkg_file)
    return dat["uuid"]
end

function create_temp_proj(proj_dir)
    temp_proj_dir = tempdir()
    mkpath(temp_proj_dir)
    for file_names in [Base.project_names, Base.manifest_names]
        for file_name in file_names
            file = joinpath(file_name, proj_dir)
            !isfile(file) && continue
            cp(file, temp_proj_dir; force = true, follow_symlinks = true)
        end
    end
    return temp_proj_dir
end

## ---------------------------------------------------------------------------------------------
# Looking for projects
println("\n", "-" ^ 45)
@info("Looking for projects", root)
const projs = find_projects(root);

## ---------------------------------------------------------------------------------------------
# Install source projects
const pkgs_pool = Dict()
let
    println("\n", "-" ^ 45)
    @info("Installing source projects", root)

    for proj_file in projs
        try
            # Install and update
            println("\n", "-" ^ 45)
            proj_file = relpath(proj_file, root)
            @info("Found on disk", proj_file)
            proj_file = abspath(proj_file)

            # Collect devs
            pkgs = load_deps(proj_file)
            merge!(pkgs_pool, pkgs)

            # Installing
            if !dry_run

                # Temp env
                proj_dir = create_temp_proj(dirname(proj_file))
                try
                    proj_file = joinpath(proj_dir, basename(proj_file))
                    Pkg.activate(proj_file)

                    for manf_name in Base.manifest_names
                        manf_file = joinpath(proj_dir, manf_name)
                        isfile(manf_file) && Pkg.instantiate()
                    end
                    Pkg.resolve()
                    Pkg.update()
                    Pkg.build()

                finally
                    rm(proj_dir; force = true, recursive = true)
                end
            end

        catch err
            @warn("ERROR", proj_file, err)
        end
    end
end

## ---------------------------------------------------------------------------------------------
# Install version range
let
    pkgs_count = length(pkgs_pool)
    println("\n", "-" ^ 45)
    @info("Installing pkg versions", pkgs_count, deep)

    deep <= 0 && return
    for (uuidpkg, name) in pkgs_pool
        try
            versions = VersionNumber[]
            for path in findin_regs(name)
                uuidfile = extract_uuid(path)
                uuidpkg != uuidfile && continue
                push!(versions, extract_versions.(path)...)
            end
            sort!(unique!(versions), rev = true)

            println("\n" ^ 3, "-" ^ 45)
            @info("Found in register", name, uuidpkg, length(versions))

            # Installing version range
            c = 0
            for version in versions

                c >= deep && break
                c += 1

                println("\n", "-" ^ 45)
                @info("Installing", name, uuidpkg, version)

                if !dry_run
                    tempenv = tempdir()
                    mkpath(tempenv)
                    try
                        Pkg.activate(tempenv)
                        # pkg = PackageSpec(name, Base.UUID(uuidpkg), version)
                        # Pkg.add(pkg)
                        # Pkg.build()

                    catch err
                        @warn("ERROR", name, uuidpkg, version, err)
                    finally
                        rm(tempenv; force = true, recursive = true)
                    end
                end
            end
        catch err
            @warn("ERROR", name, uuidpkg, err)
        end
    end
end


text:
# (letter10 kind letter color gray shape q screen-x 338 screen-y 515 value q :attended nil bottom-up activation 0.6 top-down activation 1 activation 1.26)
# (letter11 kind letter color gray shape q screen-x 351 screen-y 193 value q :attended nil bottom-up activation 1.1 top-down activation 1 activation 1.81)
# (letter12 kind letter color gray shape q screen-x 511 screen-y 462 value q :attended nil bottom-up activation 0.6 top-down activation 1 activation 1.26)
# (letter13 kind letter color gray shape q screen-x 233 screen-y 303 value q :attended nil bottom-up activation 0.85 top-down activation 1 activation 1.54)
# (letter14 kind letter color gray shape q screen-x 305 screen-y 350 value q :attended nil bottom-up activation 0.81 top-down activation 1 activation 1.49)
# (letter15 kind letter color gray shape p screen-x 301 screen-y 193 value p :attended t bottom-up activation 0.92 top-down activation 2 activation 2.22)
# (letter16 kind letter color black shape p screen-x 296 screen-y 134 value p :attended nil bottom-up activation 0.84 top-down activation 1 activation 1.52)
# (letter17 kind letter color black shape p screen-x 135 screen-y 136 value p :attended nil bottom-up activation 0.66 top-down activation 1 activation 1.33)
# (letter18 kind letter color black shape p screen-x 236 screen-y 186 value p :attended nil bottom-up activation 0.85 top-down activation 1 activation 1.54)
# (letter19 kind letter color black shape p screen-x 398 screen-y 300 value p :attended nil bottom-up activation 0.89 top-down activation 1 activation 1.58)
# (letter20 kind letter color black shape p screen-x 400 screen-y 133 value p :attended nil bottom-up activation 0.79 top-down activation 1 activation 1.47)
# time: 2.6009866538579667 location: 324 324
# Iconic Memory Size: 3

vals = [(color = :gray, shape = :q, x = 338, y =515, value = :q),
(color = :gray, shape = :q, x = 351 , y = 193, value = :q),
(color = :gray, shape = :q, x = 511 , y = 462, value = :q),
(color = :gray, shape = :q, x = 233 , y = 303, value = :q),
(color = :gray, shape = :q, x = 305 , y = 350, value = :q),
(color = :gray, shape = :p,  x = 301 , y = 193, value = :p),
(color = :black, shape = :p,  x = 296 , y = 134, value = :p),
(color = :black, shape = :p,  x = 135 , y = 136, value = :p),
(color = :black, shape = :p,  x = 236 , y = 186, value = :p),
(color = :black, shape = :p,  x = 398 , y = 300, value = :p),
(color = :black, shape = :p,  x = 400 , y = 133, value = :p)]
cd(@__DIR__)
using Pkg
Pkg.activate("..")
using Revise, PAAV, Plots, DataFrames, Statistics, StatsPlots

experiment = Experiment(n_trials=10^1,
    trace=true, visible=true, speed=.2)

target = (color=:gray,shape=:p)
visicon = [VisualObject(features=(color=Feature(;value=x.color), shape=Feature(;value=x.shape)), width=32.0,
location=[x.x-0.,x.y-0.]) for x in vals ]

model = Model(;target=target, iconic_memory=visicon, noise=0.0)
PAAV.compute_angular_size!(model)
model.focus = [324.0,324.0]
#PAAV.orient!(model, experiment)
PAAV.update_decay!(model)
PAAV.update_finst!(model)
PAAV.update_visibility!(model)
PAAV.compute_activations!(model)

map(x->(bottomup=x.bottomup_activation,topdown=x.topdown_activation), visicon)
#run_trial!(experiment)


text:
# Test residual block
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using LinearAlgebra, InvertibleNetworks, Test, Flux
import Flux.Optimise.update!

# Input
nx = 28
ny = 28
n_in = 4
n_hidden = 8
batchsize = 2

# Flux networks
model = Chain(
    Conv((3,3), n_in => n_hidden; pad=1),
    Conv((3,3), n_hidden => n_in; pad=1)
)

model0 = Chain(
    Conv((3,3), n_in => n_hidden; pad=1),
    Conv((3,3), n_hidden => n_in; pad=1)
)

# Input
X = glorot_uniform(nx, ny, n_in, batchsize)
X0 = glorot_uniform(nx, ny, n_in, batchsize)
dX = X - X0

# Flux blocks
RB = FluxBlock(model)
RB0 = FluxBlock(model0)

# Observed data
Y = RB.forward(X)

function loss(RB, X, Y)
    Y_ = RB.forward(X)
    ΔY = Y_ - Y
    f = .5f0*norm(ΔY)^2
    ΔX = RB.backward(ΔY, X)
    return f, ΔX, RB.params[1].grad, RB.params[2].grad
end

###################################################################################################
# Gradient tests

# Gradient test w.r.t. input
f0, ΔX = loss(RB, X0, Y)[1:2]
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = loss(RB, X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
RB_ini = deepcopy(RB0)
dW1 = RB.params[1].data - RB0.params[1].data
dW2 = RB.params[2].data - RB0.params[2].data

f0, ΔX, ΔW1, ΔW2 = loss(RB0, X, Y)
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    RB0.params[1].data[:] = RB_ini.params[1].data + h*dW1
    RB0.params[2].data[:] = RB_ini.params[2].data + h*dW2
    f = loss(RB0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
using HDF5,JLD,KUparser,KUnet,Base.Test
@date @load "acl11.trn.jld4"
@show net=newnet("bp10acl11eager13.jld")
ft = Flist.acl11eager
pt = ArcEager13
ncpu = 20
ndeps = length(deprel)
nbatch = 2000
@date pxy1=gparse(pt, corpus, ndeps, ft, net, nbatch, ncpu; xy=true)
@date pxy2=gparse(pt, corpus, ndeps, ft, net, nbatch; xy=true)
@test @show pxyequal(pxy1, pxy2)

# The 12-cpu version takes 150 seconds on ilac.
# The RAM usage does not exceed 12GB according to free.
# The 1-cpu version takes 212 seconds on ilac.
# The RAM usage goes up to 22GB according to free.

# GPU mem:

# For 12-cpu nbatch=500: main=673 children=320 (peak at 380MB)
# (could unload main net for multi-cpu work.)
# For 1-cpu nbatch=500: main<=717MB.

# On iui:
# For 20-cpu nbatch=2000: main=681 children=460 time=70s (20s for init/term)
# For 1-cpu nbatch=2000: main<=900 time=147s
# Same with julia4, except init takes 10s longer.
# pxyequal() runs out of RAM and gets killed.


text:
#reml.jl
 Base.@propagate_inbounds function subutri!(a, b)
    s = size(a,1)
    if s == 1 return a[1,1] -= b[1,1] end
    @simd for m = 1:s
        @simd for n = m:s
            a[m,n] -= b[m,n]
        end
    end
    a
end

function fillzeroutri!(a::AbstractArray{T})  where T
    tr = UpperTriangular(a)
    fill!(tr, zero(T))
end

################################################################################
#                     REML without provided β
################################################################################
function reml_sweep_β(lmm, θ::Vector{T}; syrkblas::Bool = false) where T <: Number
    data = LMMDataViews(lmm)
    reml_sweep_β(lmm, data, θ; syrkblas = syrkblas)
end
function reml_sweep_β(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}; syrkblas::Bool = false) where T <: Number
    n             = length(lmm.covstr.vcovblock)
    N             = length(lmm.data.yv)
    c             = (N - lmm.rankx)*log(2π)
    #---------------------------------------------------------------------------
    V⁻¹           = Vector{AbstractArray{T}}(undef, n)
    θ₁            = zero(T)
    θ₂            = zeros(T, lmm.rankx, lmm.rankx)
    θ₃            = zero(T)
    βm            = zeros(T, lmm.rankx)
    β             = Vector{T}(undef, lmm.rankx)
    #---------------------------------------------------------------------------
    logdetθ₂      = zero(T)
    θs₂           = Symmetric(θ₂)
    #θ₂ut          = UpperTriangular(θ₂)
    #akk = Vector{T}(undef, lmm.covstr.maxn + lmm.rankx) #temp for sweep
    noerror       = true

        ncore     = min(Polyester.num_cores(), n)
        accθ₁     = zeros(T, ncore)
        accθ₂     = Vector{Matrix{T}}(undef, ncore)
        accβm     = Vector{Vector{T}}(undef, ncore)
        erroracc  = trues(ncore)
        d, r = divrem(n, Polyester.num_cores())
        Base.Threads.@threads for t = 1:ncore

        #@batch per=core for i = 1:n #@fastmath
            # Vp - matrix for sweep operation
            # [V  X
            #  X' 0]
            offset = min(t-1, r) + (t-1)*d
            accθ₂[t] = zeros(T, lmm.rankx, lmm.rankx)
            accβm[t] = zeros(T, lmm.rankx)

            for j ∈ 1:d+(t ≤ r)
                i =  offset + j

                q    = length(lmm.covstr.vcovblock[i])
                qswm = q + lmm.rankx
                Vp   = Matrix{T}(undef, qswm, qswm)
                V    = view(Vp, 1:q, 1:q)
                Vx   = view(Vp, 1:q, q+1:qswm)
                Vc   = view(Vp, q+1:qswm, q+1:qswm)
                fillzeroutri!(V)
                copyto!(Vx, data.xv[i])
                fillzeroutri!(Vc)
            #-------------------------------------------------------------------
            # Make V matrix
                vmatrix!(V, θ, lmm, i)
            #-----------------------------------------------------------------------
            #swr  = sweepb!(view(akk, 1:qswm), Vp, 1:q; logdet = true, syrkblas = syrkblas)
                swm, swr, ne  = sweepb!(Vector{T}(undef, qswm), Vp, 1:q; logdet = true, syrkblas = syrkblas)
                V⁻¹[i] = V
            #-----------------------------------------------------------------------
                if ne == false erroracc[t] = false end
                accθ₁[t] += swr
                subutri!(accθ₂[t], Vc)
                mulαtβinc!(accβm[t], Vx, data.yv[i])
            end
            #-----------------------------------------------------------------------
        end
        θ₁ = sum(accθ₁)
        map(x->θ₂ .+= x, accθ₂) # make upper triangular use only!
        map(x->βm .+= x, accβm)
        noerror = all(erroracc)
        # Cholesky decomposition for matrix inverse θs₂ - Summetric(θ₂); C = θ₂⁻¹
        cθs₂ = cholesky(θs₂)
        # β calculation
        mul!(β, inv(cθs₂), βm)
        # θ₃
        @inbounds @simd for i = 1:n
            θ₃ += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
        end
        # final θ₂
        logdetθ₂ = logdet(cθs₂)

    return   θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerror #REML, β, iC, θ₃, errors
end
################################################################################
#                     REML with provided β
################################################################################
function reml_sweep_β(lmm, θ::Vector{T}, β::Vector) where T <: Number
    data = LMMDataViews(lmm)
    reml_sweep_β(lmm, data, θ, β)
end
function reml_sweep_β(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}, β::Vector) where T <: Number
    n             = length(lmm.covstr.vcovblock)
    N             = length(lmm.data.yv)
    c             = (N - lmm.rankx)*log(2π)
    #---------------------------------------------------------------------------
    θ₁            = zero(T)
    θ₂            = zeros(T, lmm.rankx, lmm.rankx)
    θ₃            = zero(T)
    #---------------------------------------------------------------------------
    logdetθ₂      = zero(T)
    akk           = zeros(T, lmm.covstr.maxn + lmm.rankx) #temp for sweep
    #Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
    θs₂           = Symmetric(θ₂)
    noerror       = true

        l = Base.Threads.SpinLock()
        Base.Threads.@threads for i = 1:n
            q    = length(lmm.covstr.vcovblock[i])
            qswm = q + lmm.rankx
            Vp   = Matrix{T}(undef, qswm, qswm)
            V    = view(Vp, 1:q, 1:q)
            Vx   = view(Vp, 1:q, q+1:qswm)
            Vc   = view(Vp, q+1:qswm, q+1:qswm)
            #Vp   = view(Vm, 1:q + lmm.rankx, 1:q + lmm.rankx)
            #V    = view(Vm, 1:q, 1:q)
            #Vx   = view(Vm, 1:q, q+1:q+lmm.rankx)
            #Vc   = view(Vm, q + 1:qswm, q + 1:qswm)
            fillzeroutri!(V)
            copyto!(Vx, data.xv[i])
            fillzeroutri!(Vc)
            vmatrix!(V, θ, lmm, i)
            #-----------------------------------------------------------------------
            swm, swr, ne = sweepb!(Vector{T}(undef, qswm), Vp, 1:q; logdet = true)
            #-----------------------------------------------------------------------
            θ₃t = mulθ₃(data.yv[i], data.xv[i], β, V)
            lock(l) do
                if ne == false noerror = false end
                θ₁  += swr
                subutri!(θ₂, view(Vp, q + 1:qswm, q + 1:qswm))
                θ₃  += θ₃t
            end
        end
        logdetθ₂ = logdet(θs₂)

    return   θ₁ + logdetθ₂ + θ₃ + c, θs₂, θ₃, noerror #REML, iC, θ₃
end
################################################################################
#                     REML AI-like / scoring part
################################################################################
function sweep_ai(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}, β::Vector) where T <: Number
    n             = length(lmm.covstr.vcovblock)
    θ₃            = zero(T)
    #akk           = zeros(T, lmm.covstr.maxn + lmm.rankx) #temp for sweep
    #Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
    l = Base.Threads.SpinLock()
    @inbounds Base.Threads.@threads for i = 1:n
        q    = length(lmm.covstr.vcovblock[i])
        qswm = q + lmm.rankx
        Vp   = Matrix{T}(undef, qswm, qswm)
        V    = view(Vp, 1:q, 1:q)
        Vx   = view(Vp, 1:q, q+1:qswm)
        Vc   = view(Vp, q+1:qswm, q+1:qswm)
        #Vp   = view(Vm, 1:q + lmm.rankx, 1:q + lmm.rankx)
        #V    = view(Vm, 1:q, 1:q)
        #Vx   = view(Vm, 1:q, q+1:q+lmm.rankx)
        #Vc   = view(Vm, q + 1:qswm, q + 1:qswm)
        fillzeroutri!(V)
        copyto!(Vx, data.xv[i])
        fillzeroutri!(Vc)
        vmatrix!(V, θ, lmm, i)
        sweepb!(Vector{T}(undef, qswm), Vp, 1:q)
        θ₃t = mulθ₃(data.yv[i], data.xv[i], β, V)
        lock(l) do
            θ₃  += θ₃t
        end
    end
    return  θ₃
end
function sweep_score(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}, β::Vector) where T <: Number
    n             = length(lmm.covstr.vcovblock)
    θ₁            = zero(T)
    θ₃            = zero(T)
    #---------------------------------------------------------------------------
    #akk           = zeros(T, lmm.covstr.maxn + lmm.rankx) #temp for sweep
    #Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
    l = Base.Threads.SpinLock()
    @inbounds Base.Threads.@threads for i = 1:n
        q    = length(lmm.covstr.vcovblock[i])
        qswm = q + lmm.rankx
        Vp   = Matrix{T}(undef, qswm, qswm)
        V    = view(Vp, 1:q, 1:q)
        Vx   = view(Vp, 1:q, q+1:qswm)
        Vc   = view(Vp, q+1:qswm, q+1:qswm)

        #Vp   = view(Vm, 1:q + lmm.rankx, 1:q + lmm.rankx)
        #V    = view(Vm, 1:q, 1:q)
        #Vx   = view(Vm, 1:q, q+1:q+lmm.rankx)
        #Vc   = view(Vm, q + 1:qswm, q + 1:qswm)

        fillzeroutri!(V)
        copyto!(Vx, data.xv[i])
        fillzeroutri!(Vc)
        vmatrix!(V, θ, lmm, i)
        #-----------------------------------------------------------------------
        swr  = sweepb!(Vector{T}(undef, qswm), Vp, 1:q; logdet = true)

        θ₃t  = mulθ₃(data.yv[i], data.xv[i], β, V)
        lock(l) do
            θ₁  += swr[2]
            θ₃  += θ₃t
        end
    end

    return   -θ₁ + θ₃
end

################################################################################
#                     variance-covariance matrix of β
################################################################################
function sweep_β_cov(lmm, θ::Vector{T}, β::Vector) where T <: Number
    n             = length(lmm.covstr.vcovblock)
    θ₂            = zeros(T, lmm.rankx, lmm.rankx)
    akk           = zeros(T, lmm.covstr.maxn + lmm.rankx) #temp for sweep
    Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
    @inbounds @simd for i = 1:n
        q    = length(lmm.covstr.vcovblock[i])
        qswm = q + lmm.rankx
        Vp   = view(Vm, 1:q + lmm.rankx, 1:q + lmm.rankx)
        V    = view(Vm, 1:q, 1:q)
        fillzeroutri!(V)
        Vx   = view(Vm, 1:q, q+1:q+lmm.rankx)
        copyto!(Vx,  view(lmm.data.xv,  lmm.covstr.vcovblock[i],:))
        Vc   = view(Vm, q + 1:qswm, q + 1:qswm)
        fillzeroutri!(Vc)
        vmatrix!(V, θ, lmm, i)
        #-----------------------------------------------------------------------
        sweepb!(fill!(view(akk, 1:qswm), zero(T)), Vp, 1:q)
        subutri!(θ₂, view(Vp, q + 1:qswm, q + 1:qswm))
    end
    return Symmetric(θ₂)
end