text:
# Rounding of vectors to a given number of significant digits
vcode(v::Verts,prec=10^5) = zeros!(map(round, v*prec)/prec)
vcode(v::Vector,prec=10^5) = zeros!(map(round, v*prec)/prec)

# random 2D point with given number of digits.
# A single 2D random point, codified in floating point format, and with a fixed (quite small)
# number of digits, is returned by the \texttt{rand_point2d()} function, with no input
# parameters.
rand_point2d(prec=10^15) = vcode(rand(2,1))

# Generation of a random line segment.
# A single random segment, scaled about its centroid by the \texttt{scaling} parameter, is
# returned by the \texttt{redge()} function, as a tuple ot two random points in the unit
# square.
function rand_edge(scaling=1)
	v1,v2 = (rand_point2d(), rand_point2d())
    c = (v1+v2)/2
    pos = rand_point2d()
    v1 = (v1-c)*scaling + pos
    v2 = (v2-c)*scaling + pos
    (vcode(v1), vcode(v2))
end

# Generation of random lines
# Return a LAR pair (V,EV)
function randomLines(lineNumber,scaling=1)
    L = hcat([hcat(rand_edge(scaling)...)  for k=1:lineNumber]...)
    xmin, ymin = min(L[1,:]...),min(L[2,:]...)
    v = [-xmin,-ymin]
    lines = [L[:,k]+v  for k=1:size(L,2)]
	V = hcat(lines...)
	EV = reshape(collect(1:2*lineNumber),(2,lineNumber))
	V,EV
end

# Generation of containment boxes of LAR cells
# Input: a LAR pair V,EV (or V,FV etc)
# Output:
# Containment boxes
function containment2DBoxes(V,EV)
	boxes = Vector{Float64}[]
	for k=1:size(EV,2)
		(v1,v2) = EV[:,k]
		(x1,y1),(x2,y2) = V[:,v1],V[:,v2]
		box = [min(x1,x2),min(y1,y2),max(x1,x2),max(y1,y2)]
		boxes = push!(boxes,vcode(box))
	end
	boxes
end

function containment2DBoxes(V,EV)
	boxes = Vector{Float64}[]
	for k=1:size(EV,2)
		(v1,v2) = EV[:,k]
		(x1,y1),(x2,y2) = V[:,v1],V[:,v2]
		box = [min(x1,x2),min(y1,y2),max(x1,x2),max(y1,y2)]
		boxes = push!(boxes,vcode(box))
	end
	boxes
end

# generate an HPC value from data of `Verts` type
polyline(verts::Verts) = p.POLYLINE(PyObject(map(PyObject,verts')))

# produce the 5 vertices of a closed rectangular polyline
# input is a box, i.e. a quadruple xmin,ymin,xmax,ymax()
function box2rect(box)
    x1,y1,x2,y2 = box
    verts = [[x1,y1] [x2,y1] [x2,y2] [x1,y2] [x1,y1]]
end

# linesFromLineArray(V,EV) = Any[[V[:,EV[:,k][1]] V[:,EV[:,k][2]]  ] for k=1:size(EV,2)]
# Implemented indirect indexing via multidimensional array
function linesFromLineArray(V,EV)
	result = Array(Float64, size(V,1), 2, size(EV,2))
	for k = 1 : size(EV,2)
		for j = 1 : 2
			for i = 1 : size(V,1)
				result[i,j,k] = V[i,EV[j,k]]
			end
		end
	end
	result
end

# Computation of the 1D centroid of a list of 2D boxes. The direction
# is chosen depending on the value of the `xy in Set([1,2])` parameter.
function centroid(boxes::Array{Float64,2},xyz)
	m = size(boxes,1)
	n = Int(m/2)
	average = zeros(m)
	for h=1:m
		average[h] = mean(boxes[h,:])
	end
	median = (average[xyz] + average[xyz+n])/2
end

function centroid(boxes::Array{Array{Float64,1},1},xyz)
	myboxes = hcat(boxes...)
	return centroid(myboxes::Array{Float64,2},xyz)
end

# Splitting the input above and below a median threshold
function splitOnThreshold(boxes,subset,coord)
	if subset==Set{Int64}()
		return Set{Int64}(), Set{Int64}()
	end
    theBoxes = [boxes[:,k] for k in subset]
    threshold = centroid(theBoxes,coord)
    ncoords = Int(floor(length(boxes[:,1])/2))
    a = coord % ncoords +1
    b = Int(a + ncoords)
    below,above = Int[],Int[]
    for k in subset
        if boxes[:,k][a] <= threshold
        	push!(below, k) end
    end
    for k in subset
        if boxes[:,k][b] >= threshold
        	push!(above, k) end
    end
    Set{Int64}(below), Set{Int64}(above)
end

function geomPartitionate(boxes,buckets)
    geomInters = [Set{Int64}() for h=1:length(boxes)]
    for bucket in buckets
        for k in bucket
            geomInters[k] = union(geomInters[k],bucket)
        end
    end
    for (h,inters) in enumerate(geomInters)
        geomInters[h] = setdiff(geomInters[h],Set([h]))
    end
    geomInters
end

# local function
function splitting(bucket,below,above, finalBuckets,splittingStack)
	a = (length(below)<4) & (length(above)<4)
	b = length(setdiff(bucket,below)) < 7
	c = length(setdiff(bucket,above)) < 7
	if ( a | b | c )
		push!(finalBuckets, below)
		push!(finalBuckets, above)
		println("pushed finalBuckets")
	else
		push!(splittingStack, below)
		push!(splittingStack, above)
		println("pushed splittingStack")
	end
end

function lar2boxes(V,CV::Array{Int64,2})
	m = size(V,1)*2
	boxes = zeros(Float64,(m,size(CV,2)))
	for k=1:size(CV,2)
	cell = CV[:,k]
		verts = hcat([V[:,v] for v in cell]...)
		pmin,pmax = Float64[],Float64[],Float64[]
		for h=1:size(verts,1)
			coords = verts[h,:]
			push!(pmin,minimum(coords))
			push!(pmax,maximum(coords))
		end
		box = vcat(pmin,pmax)
		boxes[:,k] = box
	end
	return boxes
end

function lar2boxes(V,CV::Array{Array{Int64,1},1})
	m = size(V,1)*2
	n = length(CV)
	boxes = zeros(Float64,(m,n))
	for k=1:n
		cell = CV[k]
		verts = hcat([V[:,v] for v in cell]...)
		pmin,pmax = Float64[],Float64[]
		for h=1:size(verts,1)
			coords = verts[h,:]
			push!(pmin,minimum(coords))
			push!(pmax,maximum(coords))
		end
		box = vcat(pmin,pmax)
		boxes[:,k] = box
	end
	return boxes
end

# Intersection of two line segments
function segmentIntersect(boxes,lineArray,lineStorage)
	V,EV = lineArray
    function segmentIntersect0(h)
    	h1,h2 = EV[:,h]
        p1,p2 = V[:,h1],V[:,h2]
        line1 = vcode([p1 p2])
        #line1 = [p1 p2]
        (x1,y1),(x2,y2) = p1,p2
        B1,B2,B3,B4 = boxes[:,h]
        function segmentIntersect1(k)
			k1,k2 = EV[:,k]
			p3,p4 = V[:,k1],V[:,k2]
			line2 = vcode([p3 p4])
			(x3,y3),(x4,y4) = p3,p4
			b1,b2,b3,b4 = boxes[:,k]
            if ! |(b3<B1 , B3<b1 , b4<B2 , B4<b2)
                m23 = [p2 p3]
                m14 = [p1 p4]
                m = m23 - m14
                v = p3-p1
                a=m[1,1]; b=m[1,2]; c=m[2,1]; d=m[2,2]
                det = a*d-b*c
                if abs(det) > 10^-6.
                    m_inv = inv(m)
                    alpha, beta = m_inv*v
                    if (-0.0<=alpha<=1.0) & (-0.0<=beta<=1.0)
                    	if 0.0 != vcode([alpha])[1] != 1.0
                        	push!(lineStorage[line1], alpha)
                        end
                    	#if 0.0 != vcode([alpha])[1] != 1.0
                    	if 0.0 != vcode([beta])[1] != 1.0
	                        push!(lineStorage[line2], beta)
	                    end
                        return p1+alpha*(p2-p1)
                    end
                end
            end
            #return nothing
        end
    end
end

# Brute force intersection within the buckets
function lineBucketIntersect(boxes,lineArray, h,bucket, lineStorage)
    intersect0 = segmentIntersect(boxes,lineArray,lineStorage)
    intersectionPoints = Array{Float64,1}[]
    intersect1 = intersect0(h)
    for line in bucket
    	#@show line
        point = intersect1(line)
        if point != nothing
            push!(intersectionPoints, vcode(point))
        end
    end
    intersectionPoints
end

function boxBucketing(boxes::Array{Float64,2})
	nboxes = size(boxes,2)
	dim = Int(size(boxes,1)/2)
	trees,boxes1D = Any[],Any[]
	# preparation of d interval-trees
	for d=1:dim
		intervals = [ (boxes[d,k], boxes[d+dim,k], k) for k=1:nboxes ];
		push!(boxes1D, intervals)
		tree1D = IntervalTree{Float64, IntervalValue{Float64, Int64}}();
		for triple in intervals
			push!(tree1D, IntervalValue{Float64, Int64}(triple...));
		end
		push!(trees, tree1D)
	end
	# execution of spatial queries
	buckets,queries = Any[],Any[]
	for k=1:nboxes
		sets = Any[]
		for d=1:dim
			push!(sets, Set(Int64[]))
			query = boxes1D[d][k]
			for item in intersect(trees[d],query[1:2])
				push!(sets[d], item.value)
			end
		end
		push!(buckets,intersect(sets...))
	end
	buckets = [sort(collect(bucket)) for bucket in buckets]
	return buckets
end

# Accelerate intersection of lines
function lineIntersection(lineArray)
	V,EV = lineArray
	lineStorage = Dict{Array{Float64,2},Array{Float64,1}}()
	for k = 1:size(EV,2)
		v1,v2 = EV[:,k]
		p1,p2 = V[:,v1],V[:,v2]
		key = vcode([p1 p2])
		lineStorage[key] = Int[]
	end
	boxes = lar2boxes(lineArray...)
	buckets = boxBucketing(boxes)
    for (h,bucket) in enumerate(buckets)
    	#@show h,bucket
        pointBucket = lineBucketIntersect(boxes,lineArray, h,bucket, lineStorage)
    end
    frags = keys(lineStorage)
    params = [sort([x for x in Set(vcode(group))]) for group in values(lineStorage)]
    lineFrags = Dict(zip(frags,params))
end

function check4DoubleEdges(EV)
	for k=1:(size(EV,2)-1)
		if EV[:,k]==EV[:,k+1]

		end
	end
end

function larModelCheck(V,EV)
	vertexIndex,k = Dict{Int64,Int64}(),0
	W = sort(collect(Set(EV)))
	for (k,v) in enumerate(W)
		vertexIndex[v] = k
	end
	EW = deepcopy(EV)
	for k=1:size(EV,2)
		EW[1,k] = vertexIndex[EV[1,k]]
		EW[2,k] = vertexIndex[EV[2,k]]
	end
	vkeys = p.TRANS(sort(collect(vertexIndex)))'[:,1]
	V = hcat([V[:,k] for k in vkeys]...)
	EV = hcat([sort(EW[:,k]) for k=1:size(EW,2)]...)
	EV = sortcols(EV, by=x->(x[1],x[2]))
	EV = sortrows(EV)
	#EV = check4Edges(EV)
	V,EV
end

# Transform a lineArray (Array of pairs of 2D points) into a 1-complex (V,EV)
function lines2lar(lineArray,prec=10^15)
	lineFrags = lineIntersection(lineArray)

	verts = zeros(2,1)
	index = 0
	EV = Array{Any,1}()
	for k in 1:length(lineFrags)
		ps = collect(keys(lineFrags))[k]
		alphas = vcat([[0.],collect(values(lineFrags))[k],[1.]]...)
		pts = [ (ps[:,1]*(1.-alphas[i])+ps[:,2]*alphas[i]) for i=1:length(alphas) ]
		verts = hcat(verts,hcat(pts...))
		for h=(index+1):(index+length(pts)-1)
			push!(EV, [h,h+1])
		end
		index += length(pts)
	end
	V,EV = verts[:,2:end],hcat([EV[k] for k=1:length(EV)]...)

	W,close,clusters,vmap = p.pruneVertices(V')
	if typeof(W) == Array{Any,1}
		Z = Array{Float64,2}(hcat(W...))
	else
		Z = Array{Float64,2}(W)
	end
	vmap = vmap + 1
	EW = [[vmap[EV[1,k]],vmap[EV[2,k]]] for k in 1:size(EV,2)]
	EZ = hcat(EW...)
	V,EV = biconnectedComponents(Z,EZ)
	if size(V,2)==2 V=V' end
	V,EV = larModelCheck(V,EV)
end

# Adjacency lists of 1-complex vertices
function vertices2vertices(V,EV)
	lar = Lar()
	lar.EV = [EV[:,k] for k=1:size(EV,2)]
	ev = cellComplex(lar.EV)
	vv = ev'*ev
	triples = hcat(findnz(vv)...)'
	VV = [Int64[] for k=1:size(vv,1)]
	for k=1:size(triples,2)
		row,col,datum = triples[:,k]
		if row != col
			push!(VV[col],row)
		end
	end
	return VV
end

# Main procedure for biconnected components
function biconnectedComponents(W,EV)
    V = 1:length(W)
    count = 0
    stack,out = [],[]
    visited = [false for v in V]
    parent = [-1 for v in V]
    d = [-1 for v in V]
    low = [-1 for v in V]
    VV = vertices2vertices(W,EV)
    for u in V
        if !visited[u]
            DFV_visit( VV,out,count,visited,parent,d,low,stack, u )
        end
    end
    out = [component for component in out if length(component) > 1]
    EV = [[u,v] for (u,v) in vcat(out...)]
    W,hcat(EV...)
end

# Hopcroft-Tarjan algorithm
function DFV_visit( VV,out,count,visited,parent,d,low,stack,u )
    visited[u] = true
    count += 1
    d[u] = count
    low[u] = d[u]
    if u > length(VV) return() end
    for v in VV[u]
        if !visited[v]
            push!(stack, (u,v))   # possible bug: []
            #println("stack = $stack")
            parent[v] = u
            DFV_visit( VV,out,count,visited,parent,d,low,stack, v )
            if low[v] >= d[u]
                push!(out,outputComp(stack,u,v))
                #println("\tout = $out")
                #println("\tv = $v")
                if v == length(VV) return() end
            end
            low[u] = min( low[u], low[v] )
        else
            if ! (parent[u]==v) & (d[v] < d[u])
                push!(stack, (u,v))   # possible bug: []
                low[u] = min( low[u], d[v] )
            end
        end
    end
end

# Output of biconnected components
function outputComp(stack,u,v)
    out = []
    while true
        elem = pop!(stack)
        push!(out, elem)
        if elem == (u,v)
        	break
        end
    end
    out
end

function larFromLines(datafile) # lineArray'
	V = reshape(datafile',(size(datafile',1)÷2,size(datafile',2)*2))
	len = length(datafile)
	EV = collect(reshape(1:(len÷2), 2,(len÷4)))
	W,EW = lines2lar((V,EV))
	chains = boundary(W,EW)
	operator = boundaryOp(EW,chains)
	FW = [sort(collect(Set(vcat([EW[:,abs(e)] for e in face]...)))) for face in chains]
	W,FW,EW
end


text:
__precompile__(true)

module ObjectiveFunctions

using Compat: @compat
using Reexport
@reexport using LearnBase
@reexport using LossFunctions
@reexport using Transformations
@reexport using PenaltyFunctions

import LearnBase: transform!, grad, grad!, params, update!
import Transformations: input_node, output_node, totalcost, InputNode, SumNode, OutputNode

export
    AbstractLossTransform,
    NoLoss,
    LossTransform,
    CrossEntropy,
    RegularizedObjective,
    objective

@compat abstract type AbstractLossTransform{T} <: Transformation end

immutable NoLoss <: AbstractLossTransform{Void} end

# this assumes that input_value(lt) is pre-populated... copy the target in, then pass to no-arg version
function transform!(lt::AbstractLossTransform, target::AbstractVector)
    copy!(value(lt.target), target)
    transform!(input_node(lt))
    transform!(lt)
end

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

# this is a transformation which stores the calculated loss, as well as the gradients of the inputs
# it can be linked to a transformation to allow easier forward/backward calcs
immutable LossTransform{T,L<:Loss} <: AbstractLossTransform{T}
    loss::L
    nin::Int
    input::SumNode{T,1}
    target::SumNode{T,1}
    output::OutputNode{T,1}

    function (::Type{LossTransform{T, L}}){T, L <: Loss}(loss::Loss, nin::Int)
        input = InputNode(T, nin)
        target = InputNode(T, nin)
        output = OutputNode(T, 1)
        grad(output)[1] = one(T)  # ∂L/∂L == 1
        new{T, L}(loss, nin, input, target, output)
    end
end

LossTransform{L<:Loss}(loss::L, nin::Int) = LossTransform{Float64, L}(loss, nin)

# input and target are pre-populated... compute the output value as: ∑ loss(targetᵢ, inputᵢ)
function transform!(lt::LossTransform)
    lt.output.val[1] = value(lt.loss, value(lt.target), input_value(lt), AvgMode.Sum())
    lt
end

# TODO: is this right??  what is LossFunctions.deriv?
function grad!(lt::LossTransform)
    # update the input gradient using the values of target and input
    deriv!(grad(lt.input), lt.loss, value(lt.target), value(lt.input))
end

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

# NOTE: special handling for the Softmax/CrossEntropyLoss layer.
#   Instead of a proper activation function layer, we do both the
#   softmax activation and cross entropy loss together.

immutable CrossEntropy{T} <: AbstractLossTransform{T}
    n::Int
    input::SumNode{T,1}
    target::SumNode{T,1}
    output::OutputNode{T,1}

    function (::Type{CrossEntropy{T}}){T}(n::Int)
        input = InputNode(T, n)
        target = InputNode(T, n)
        output = OutputNode(T, 1)
        grad(output)[:] = one(T)  # ∂L/∂L == 1
        new{T}(n, input, target, output)
    end
end
CrossEntropy(n::Int) = CrossEntropy{Float64}(n)

function transform!{T}(ce::CrossEntropy{T})
    ce.output.val[1] = -sum(ce.target.val[i] * log(ce.input.val[i]) for i=1:ce.n)
end

# ** NOTE **: This is really solving for the gradient with respect to the softmax activation inputs,
#             but we diverge from correctness to allow for the simplified gradient calc.
function grad!{T}(ce::CrossEntropy{T})
    for i=1:ce.n
        ce.input.∇[i] = ce.input.val[i] - ce.target.val[i]
    end
end

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

# TODO: should we support penalties on outputs as well as parameters?

# convenience type for Empirical Risk Minimization and similar
immutable RegularizedObjective{T<:Transformation, L<:AbstractLossTransform, P<:Penalty} <: Minimizable
    transformation::T
    loss::L
    penalty::P
end

# first link the transformation to the loss, then construct
function objective(transformation::Transformation, lt::AbstractLossTransform, penalty::Penalty = NoPenalty())
    # link_nodes!(output_node(transformation), input_node(lt))
    link_nodes!(transformation, lt)
    RegularizedObjective(transformation, lt, penalty)
end

# first create the LossTransform to wrap the Loss, then construct
function objective(transformation::Transformation, loss::Loss, penalty::Penalty = NoPenalty())
    lt = LossTransform{Float64,typeof(loss)}(loss, output_length(transformation))
    objective(transformation, lt, penalty)
end

# convenience when creating chains... auto-choose loss given final layer
function objective(chain::Chain, penalty::Penalty = NoPenalty())
    T = typeof(chain[end])
    loss = if T <: Affine || T <: Activation{:identity}
        L2DistLoss()
    elseif T <: Activation{:logistic}
        CrossentropyLoss()
    elseif T <: Activation{:softmax}
        CrossEntropy(output_length(chain[end]))
    else
        error("Can't pick a default loss for $T... choose it explicitly.")
    end
    objective(chain, loss, penalty)
end

function objective(t::Transformation, nl::NoLoss, penalty::Penalty = NoPenalty())
    RegularizedObjective(t,nl,penalty)
end

input_node(obj::RegularizedObjective) = input_node(obj.transformation)
output_node(obj::RegularizedObjective) = output_node(isa(obj.loss, NoLoss) ? obj.transformation : obj.loss)

params(obj::RegularizedObjective) = params(obj.transformation)
grad(obj::RegularizedObjective) = grad(obj.transformation)
totalcost(obj::RegularizedObjective) = norm(value(output_node(obj)))

function transform!(obj::RegularizedObjective, target::AbstractArray, input::AbstractArray)
    # forward pass through the transformation... assuming output has been linked to loss already
    transform!(obj.transformation, input)

    # loss input should be populated already... forward pass with target
    transform!(obj.loss, target)

    # add the penalty
    total_penalty = value(obj.penalty, params(obj.transformation))
    output_value(obj.loss)[1] += total_penalty
end

# if our target is a single value, lets convert it to a length-1 vector
function transform!{T}(obj::RegularizedObjective, target::T, input::AbstractArray{T})
    transform!(obj, [target], input)
end

# function apply_penalty(penalty::Penalty, θ, ∇)
#     for (i,j) in zip(eachindex(θ), eachindex(∇))
#         ∇[j] += deriv(penalty, θ[i])
#     end
#     ∇
# end

# we don't need data args because they were given or computed in transform (check this fact!)
function grad!(obj::RegularizedObjective)
    grad!(obj.loss)
    grad!(output_node(obj.transformation))
    grad!(obj.transformation)
    addgrad!(grad(obj), obj.penalty, params(obj))
end

# # handle the no-data case... probably just minimizing a function
# function update!(obj::RegularizedObjective, ::Void)
#     transform!(obj.transformation)
#
#     # output grad is ones, backwards pass
#     output_grad(obj.transformation)[:] = 1
#     grad!(obj.transformation)
#
#     apply_penalty(obj.penalty, params(obj), grad(obj))
# end

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

end # module

ext: jl


text:
using OrdinaryDiffEq
using Trixi
using Plots

########## functions that need to be loaded ##########################
function WR2_initial_condition_constant(x, t, equations::CompressibleEulerEquations2D)
    rho = 1.0
    rho_v1 = 1.5
    rho_v2 = 0
    rho_e = 25.0
    return SVector(rho, rho_v1, rho_v2, rho_e)
end

function boundary_condition_noslip_isothermal(u_inner, q1_inner, q2_inner, normal::AbstractVector,
        direction, x, t, surface_flux_function, equations::Trixi.AbstractEquations)
        # Freeslip wall Conditions
    rho, rho_v1, rho_v2, rho_e = u_inner
    u_boundary = SVector(rho, -rho_v1, -rho_v2, rho_e)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, normal, equations)
    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, normal, equations)
    end

    if typeof(equations) == Trixi.CompressibleEulerEquations2D{Float64} && equations.viscous

        p_0 = Trixi.cons2prim([1.0, 1.5, 0.0, 25], equations)[4] # initial pressure with constant initial condition, v1 !
        ρ_0 = 1.0 # initial density
        T_0 = p_0/ρ_0 # initial temperature
        γ = equations.gamma
        # computation of the exact pressure solution of the riemann problem
        rho, v1, v2, p = cons2prim([rho, rho_v1, rho_v2, rho_e], equations) # inner pressure
        c = sqrt(γ * (p/rho)) # sound speed
        v_norm = sqrt((normal[1]* v1)^2 + (normal[2] * v2)^2) # norm of velocity in normal direction
        Mach_n = v_norm / c # Mach number in normal direction
        hmpf = Mach_n*(γ+1)/4 # hmpf
        if v_norm > 0
            pstar_p = 1 + γ* Mach_n* (hmpf + sqrt(hmpf^2 + 1)) # exact pressure solution of the riemann problem
        else
            pstar_p = (1 + 0.5* (γ - 1)* Mach_n)^(2* γ/(γ - 1)) # exact pressure solution of the riemann problem
        end

        rho_viscous = pstar_p / T_0 # 𝑝 = 𝜌* 𝑅* T perfect gas model temperature relation with R = 1 universal gas constant
        vector_dummy = SVector(rho_viscous, 0, 0, pstar_p)
        u_viscous = prim2cons(vector_dummy, equations) # u_boundary for the viscous computation
        viscous = Trixi.viscous_flux(u_viscous, q1_inner, q2_inner, normal, equations)
        flux += viscous # add viscous part to boundary flux computation
    end
    return flux
end

function boundary_condition_stream(u_inner, q1, q2, orientation, direction, x, t,
        surface_flux_function,
        equations::Trixi.AbstractEquations)

    u_boundary = SVector(1, 1.5, 0, 25)
        # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end

    return flux
end
function mapping1zu1(xi_, eta_)
    x = xi_
    y = eta_
    return SVector(x, y)
end
    # mapping O-mesh
function mapping(xi_, eta_)
    ξ = xi_
    η = eta_
    size = 2    # gebietsgröße * 3
    d = 1.5       # relativ zum durchmesser
    x = size * (d + ξ) * cos(π * (η + 1))
    y = size * (d + ξ) * sin(π * (η + 1))
    return SVector(x, y)
end

#########################   Einstellungen  ############################################
"Note : if you change the velocity, it needs to be changed in the boundaries as well !
∣∣v0∣∣ D ρ0/µ = RE
turbulenter:        mu = 0.0008, v1 = 0.1
weniger turbulent:  mu = 0.0008, v1 = 1.5  "
CFL = 0.8
tspan = (0.0, 500)
N = 3
cells_per_dimension = (32, 64)
mu = 0.0008
#visualization = VisualizationCallback(interval = 500, plot_creator=Trixi.save_plot)# clims=(-0.5,0.5),

initial_condition = WR2_initial_condition_constant
boundary_conditions = (x_neg = boundary_condition_noslip_isothermal,# boundary_condition_freeslip,
                       x_pos = boundary_condition_stream,
                       y_neg = boundary_condition_periodic,
                       y_pos = boundary_condition_periodic)

mesh = CurvedMesh(cells_per_dimension, mapping, periodicity=(false, true))
# , periodicity=(false, true))#(-1.0, -1.0), (1.0, 1.0))#, mapping)
# mesh = CurvedMesh(cells_per_dimension, (-1.0, -1.0), (1.0, 1.0))

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

equations = CompressibleEulerEquations2D(1.4, viscous=true, mu=mu)

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive))
basis = LobattoLegendreBasis(N)
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
solver = DGSEM(basis, surface_flux, volume_integral)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                      # extra_analysis_errors=(:conservation_error,),
                                      # extra_analysis_integrals=(entropy, energy_total,
                                      # energy_kinetic, energy_internal)
                                      )

visualization = VisualizationCallback(interval=150, variable_names=["v1", "v2"], plot_creator=Trixi.save_plot)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval = 500,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        # alive_callback,
                        # save_solution,
                        visualization,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

pd = PlotData2D(sol)
b = plot(pd["rho"])
plot!(getmesh(pd))

max_stars_repo_path: src/WR2/MS5 ZylinderFlow.jl

max_issues_repo_path: src/WR2/MS5 ZylinderFlow.jl

max_forks_repo_path: src/WR2/MS5 ZylinderFlow.jl


text:
## Start with `julia --project=.. -i -- example01_sim.jl`

# Make sure necessary packages (as defined in ../Project.toml) are installed
import Pkg; Pkg.instantiate()

using Distributed
## If you want to run on multiple cores, uncomment the following line:
#addprocs()

@everywhere begin
    using DataFrames, CSV
    using DataFramesMeta
    using ProgressMeter
    import Base.Iterators: product
    include("../src/mcmc.jl")
end
# Temperature range
Ts = range(2.0, 2.5, length=100)
# System sizes
Ls = 20:10:60

it = product(Ls, Ts)
lit = length(it)

## This gives you a dataframe filled with the full time series
raw_rows = @showprogress pmap(it) do (L,T)
    begin
        τ = 10    # time between samples in MCS
        N = 50000 # No. of samples

        (L=L, T=T, m=metropolis_timeseries(L, 1/T, N*τ*L^2; sample_interval=τ*L^2)) # take samples
    end
end
d = DataFrame(reshape(raw_rows, lit));

## Compute averages
## m, |m|, m^2, m^4, e, e^2
e = @linq d |>
    select(:L,:T, e=mean.(:e), e2=mean.(map(x->x.^2, :e)),
         m=mean.(:m), m2=mean.(map(x->x.^2, :m)), m4=mean.(map(x->x.^4,:m)), absem=mean.(map(x->abs.(x),:m)));

CSV.write("mydata.csv", e);

max_stars_repo_name: skleinbo/Ising

max_issues_repo_name: Yaren-Github/Ising

max_forks_repo_name: Yaren-Github/Ising


text:
using QPSReader, QuadraticModels, QuadraticModelsXpress
using Test

@testset "QuadraticModelsXpress.jl" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = xpress(QuadraticModel(qps1), presolve=0, crossover=0)
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = xpress(QuadraticModel(qps2))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = xpress(QuadraticModel(qps3))
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)

    qps4 = readqps("AFIRO.SIF") # test lp
    stats4 = xpress(QuadraticModel(qps4))
    @test isapprox(stats4.objective, -4.6475314286e2, atol=1e-2)
end

max_stars_repo_name: geoffroyleconte/QuadraticModelsXpress

max_issues_repo_name: geoffroyleconte/QuadraticModelsXpress

max_forks_repo_name: geoffroyleconte/QuadraticModelsXpress


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0, 1f0,360))
end
# last argument can be used to control the granularity of the resulting mesh
sphere = GLNormalMesh(Sphere(Point3f0(0.5), 0.5f0), 24)
c = collect(linspace(0.1f0,1.0f0,10f0))
rotation = map(rotationmatrix_z, const_lift(*, timesignal, 2f0*pi))
# create a visualisation
vis = visualize((sphere, c), model=rotation, scale_y=0.1f0)
_view(vis, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
module SpSkipGramB
# share sparse representations

using CUDArt: HostArray, CudaArray, synchronize, null_stream
using JLD: jldopen
using CycUtils: invertable_unique, get, set!
using Veclib: gamma_kernel!, gamma_kernel_loss!, normalize_columns!

export SparseSkipGramB, train, init, close_device

include("config.jl")
include("array_handler.jl")
include("samples.jl")
include("utils.jl")

type SparseSkipGramB{Tensor} # Tensor could be HostArray{Float32, 2}
    α::Tensor
    β::Tensor
    B::Tensor
    C::Tensor

    # the following fields are training parameters.
    m::Int
    k::Int
    lr::Float32
    lr0::Float32
    λ::Float32

    # maps
    idx2word::Array{String, 1}
    word2idx::Dict{String, UInt32}

    # the following fields are named memory for training
    ∇αβ::Tensor
    ᾱ::Tensor
    β̄::Tensor
    a⃗::Tensor
    b⃗::Tensor
    γ::Tensor
    σ::Tensor
    ∇α::Tensor
    N::Tensor
    M::Tensor
    ∇β::Tensor
    Q::Tensor
    ∇B::Tensor
    ∇C::Tensor
    idxes::Tensor
    sign::Tensor

    function SparseSkipGramB(d::Integer, # vector dimension
        num_base::Integer,     # number of atom words
        num_vocab::Integer,    # number of vocabulary
        mini_batch::Integer,     # mini-batch size
        negative::Integer,        # number of negative samples per positive sample
        lr::AbstractFloat=0.01,       # learning rate
        λ::AbstractFloat=0.1,         # sparseness weight
        init_model::String=""   # dumped init_model.
    )
        @assert num_base ≤ num_vocab "The number of atom words must be less than the vocabulary size"

        α_ = randn(Float32, num_base, num_vocab) ./ sqrt(num_base)
        B_ = randn(Float32, d, num_base) ./ sqrt(d)
        C_ = randn(Float32, d, num_base) ./ sqrt(d)

        i2w = Array{String, 1}()
        w2i = Dict{String, UInt32}()

        if init_model != ""
            e, c, i2w = jldopen(init_model, "r") do f
                (read(f, "embeddings"), read(f, "contexts"),
                 read(f, "idx2word"))
            end
            @assert num_vocab ≤ length(i2w) "Initial model must contains as many vocabularies as the model to be trained."
            @assert d == size(e, 1) "Initial vector dim. and specified dim. must agree."

            i2w = i2w[1:num_vocab]
            for (i, w) in enumerate(i2w)
                w2i[w] = i
            end

            α_[1:num_base, 1:num_base] = eye(num_base)
            copy!(B_, @view e[:, 1:num_base])
            copy!(C_, @view c[:, 1:num_base])
        end

        α = Tensor(α_)
        β = α
        B = Tensor(B_)
        C = Tensor(C_)

        ∇αβ = Tensor(Float32, num_base*(negative+2)*mini_batch)
        ∇α, ∇β = splitarray(∇αβ, (num_base, mini_batch),
                                   (num_base, negative+1, mini_batch))

        ᾱ = Tensor(Float32, num_base, mini_batch)
        β̄ = Tensor(Float32, num_base, negative+1, mini_batch)
        a⃗ = Tensor(Float32, d, mini_batch)
        b⃗ = Tensor(Float32, d, negative+1, mini_batch)
        γ = Tensor(Float32, negative+1, mini_batch)
        σ = Tensor(Float32, negative+1, mini_batch)
        # ∇α = Tensor(Float32, num_base, mini_batch)
        N = Tensor(Float32, d, mini_batch)
        M = Tensor(Float32, num_base, mini_batch)
        # ∇β = Tensor(Float32, num_base, negative+1, mini_batch)
        Q = Tensor(Float32, num_base, mini_batch)
        ∇B = Tensor(Float32, d, num_base)
        ∇C = Tensor(Float32, d, num_base)
        idxes = Tensor(UInt32, (negative+2)*mini_batch)

        fill!(∇αβ, 0.0)
        fill!(ᾱ, 0.0)
        fill!(β̄, 0.0)
        fill!(a⃗, 0.0)
        fill!(b⃗, 0.0)
        fill!(γ, 0.0)
        fill!(σ, 0.0)
        # fill!(∇α, 0.0)
        fill!(N, 0.0)
        fill!(M, 0.0)
        # fill!(∇β, 0.0)
        fill!(Q, 0.0)
        fill!(∇B, 0.0)
        fill!(∇C, 0.0)
        fill!(idxes, 0)

        tmp = -ones(Float32, (negative+1, mini_batch))
        tmp[1, :] = 1.0f0
        s = Tensor(tmp)

        new{Tensor}(α, β, B, C,                               # model parameters
            mini_batch, negative, lr, lr, λ,              # training parameters
            i2w, w2i,                                 # maps
            ∇αβ, ᾱ, β̄, a⃗, b⃗, γ, σ, ∇α, N, M, ∇β, Q, ∇B, ∇C, idxes, s
        )
    end
end

function save(model::SparseSkipGramB{Array}, fname::String; save_maps=false)
    jldopen(fname, "w") do f
        write(f, "α", model.α)
        write(f, "B", model.B)
        write(f, "C", model.C)

        # the following fields are training parameters.
        write(f, "m", model.m)
        write(f, "k", model.k)
        write(f, "lr", model.lr0)
        write(f, "λ", model.λ)

        # maps
        if save_maps == true
            write(f, "idx2word", model.idx2word)
            write(f, "word2idx", model.word2idx)
        end
    end
end

num_vocab(model::SparseSkipGramB) = size(model.α, 2)
num_base(model::SparseSkipGramB) = size(model.α, 1)
num_vecdim(model::SparseSkipGramB) = size(model.B, 1)
num_minibatch(model::SparseSkipGramB) = size(model.a⃗, 2)
num_negative(model::SparseSkipGramB) = size(model.b⃗, 2) - 1

function to_device{M <: Union{SparseSkipGramB{Array},
                    SparseSkipGramB{HostArray}}}(model::M)
    d = num_vecdim(model)
    num_base_ = num_base(model)
    num_vocab_ = num_vocab(model)
    mini_batch = num_minibatch(model)
    negative = num_negative(model)

    d_model = SparseSkipGramB{CudaArray}(d, num_base_, num_vocab_, mini_batch,
                                        negative, model.lr, model.λ)
    for fn in fieldnames(typeof(d_model))
        if fn == :β
            continue
        end
        dest = getfield(d_model, fn)
        if isa(dest, CudaArray)
             src = getfield(model, fn)
             copy!(dest, src)
        end
    end
    return d_model
end

@generated function get_training_model(model::SparseSkipGramB)
    if BACKEND == "GPU"
        return :(to_device(model))
    else
        return :(model)
    end
end

function compute_grad_sparse(model, batch)
    const d = num_vecdim(model)
    const m = model.m
    const k = model.k
    const B = model.B
    const C = model.C
    const N = model.N
    const M = model.M
    w, c = batch
    const nb = num_base(model)
    const ᾱ = takelast!(model.ᾱ, model.α, w)
    const β̄ = takelast!(model.β̄, model.β, c)
    const a⃗ = model.a⃗
    const b⃗ = model.b⃗
    const γ = model.γ
    const ∇α = model.∇α
    const ∇β = model.∇β

    flat_b⃗ = reshape(b⃗, (d, (k+1)*m))
    gemm!('N', 'N', 1.0f0, B, ᾱ, 0.0f0, a⃗)
    gemm!('N', 'N', 1.0f0, C, reshape(β̄, (nb, (k+1)*m)), 0.0f0, flat_b⃗)

    a⃗_ = reshape(a⃗, (d, 1, m))
    γ_ = reshape(γ, (k+1, 1, m))
    gemm_strided_batched!('T', 'N', 1.0f0, b⃗, a⃗_, 0.0f0, γ_)
    gamma_kernel_loss!(γ, model.σ, model.sign)

    #∇α
    N_ = reshape(N, (d, 1, m))
    gemm_strided_batched!('N', 'N', 1.0f0, b⃗, γ_, 0.0f0, N_)
    gemm!('T', 'N', 1.0f0, B, N, 0.0f0, ∇α)

    # ∇β
    gemm!('T', 'N', 1.0f0, C, a⃗, 0.0f0, M)
    M_ = reshape(M, (num_base(model), 1, m))
    gemm_strided_batched!('N', 'T', 1.0f0, M_, γ_, 0.0f0, ∇β)
    nothing
end

function compute_grad_dict(model)
    const m = model.m
    const k = model.k

    gemm!('N', 'T', 1.0f0, model.N, model.ᾱ, 0.0f0, model.∇B)
    β̄ = model.β̄                                   # |B|×(k+1)×m
    γ = reshape(model.γ, (k+1, 1, m))             # (k+1) × 1 × m
    Q = reshape(model.Q, (num_base(model), 1, m)) # |B| × 1 × m
    gemm_strided_batched!('N', 'N', 1.0f0, β̄, γ, 0.0f0, Q)
    gemm!('N', 'T', 1.0f0, model.a⃗, model.Q, 0.0f0, model.∇C)
    nothing
end

function update_dict_u1f1!(model)
    const m = model.m
    const k = model.k
    const lr_ = model.lr

    gemm!('N', 'T', -lr_, model.N, model.ᾱ, 1.0f0, model.B)
    β̄ = model.β̄  # |B|×(k+1)×m
    γ = reshape(model.γ, (k+1, 1, m))  # (k+1) × 1 × m
    Q = reshape(model.Q, (num_base(model), 1, m)) # |B| × 1 × m
    gemm_strided_batched!('N', 'N', 1.0f0, β̄, γ, 0.0f0, Q)
    gemm!('N', 'T', -lr_, model.a⃗, model.Q, 1.0f0, model.C)
    nothing
end

function update_dict_u2f1!(model)
    compute_grad_dict(model)
    update_dict_u2!(model)
end

function accumulate_grad_dict(model, t)
    const m = model.m
    const k = model.k
    w1 = (Float32(t)/(t+1)) :: Float32
    w2 = (1.0f0/(t+1)) :: Float32

    gemm!('N', 'T', w2, model.N, model.ᾱ, w1, model.∇B)
    β̄ = model.β̄                                   # |B|×(k+1)×m
    γ = reshape(model.γ, (k+1, 1, m))             # (k+1) × 1 × m
    Q = reshape(model.Q, (num_base(model), 1, m)) # |B| × 1 × m
    gemm_strided_batched!('N', 'N', 1.0f0, β̄, γ, 0.0f0, Q)
    gemm!('N', 'T', w2, model.a⃗, model.Q, w1, model.∇C)
    nothing
end

function update_dict_u1!(model)
    const lr_::Float32 = model.lr
    axpy!(-lr_, vec(model.∇B), vec(model.B))
    axpy!(-lr_, vec(model.∇C), vec(model.C))
    nothing
end

const ∇B_cache = Ref{HostArray}()
const ∇C_cache = Ref{HostArray}()
function update_dict_u2!(model)
    global ∇B_cache
    global ∇C_cache
    ∇B = ∇B_cache[]
    ∇C = ∇C_cache[]
    lr_::Float32 = model.lr

    copy!(∇B, model.∇B)
    CUDArt.synchronize(CUDArt.null_stream)
    normalize_columns!(∇B)
    copy!(model.∇B, ∇B)
    copy!(∇C, model.∇C)
    CUDArt.synchronize(CUDArt.null_stream)
    axpy!(-lr_, vec(model.∇B), vec(model.B))
    normalize_columns!(∇C)
    copy!(model.∇C, ∇C)
    axpy!(-lr_, vec(model.∇C), vec(model.C))
    nothing
end

function fetch_grad_sparse!(grads, model)
    copy!(grads, model.∇αβ)
    return grads
end

# function cpu_sparse_update!(arr, grad, idxes, lr, lambda)
#     for (j, idx) in zip(1:size(grad,2), idxes)
#         @inbounds for i in 1:size(arr, 1)
#             x = arr[i, idx] -  lr * grad[i, j]
#             s = (x > 0.0f0 ? 1.0f0 : -1.0f0)
#             x *= s
#             x -= lambda*lr
#             x = x > 0.0f0 ? x : 0.0f0
#             x = x < 1.0f0 ? x : 1.0f0
#             arr[i, idx] = x * s
#         end
#     end
# end

function gpu_sparse_update!(arr, grad, idxes, lr, lambda)
    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)
    ccall((:sparse_update, vec_lib), Void, (Ptr{Float32}, Ptr{Float32},
        Cint, Cint, Ptr{Cuint}, Float32, Float32, Ptr{Void}, Cint),
        arr.ptr, grad.ptr, size(grad, 1), size(grad, 2),
        idxes.ptr, lr, lambda, null_stream.inner.handle, 256)
end

@generated function proximal_update!(model, ∇αβ, idxes)
    if BACKEND == "GPU"
        quote
            lr_::Float32 = model.lr
            dev = Int(DEVICE)
            siz_αβ = (num_base(model), length(idxes))
            d_∇αβ = CudaArray(model.∇αβ.ptr, siz_αβ, dev)
            copy!(d_∇αβ, ∇αβ)

            ptr_idxes = model.idxes.ptr
            d_idxes = CudaArray(ptr_idxes, (length(idxes),), dev)
            copy!(d_idxes, idxes)
            gpu_sparse_update!(model.α, d_∇αβ, d_idxes, lr_, model.λ)
        end
    else
        quote
            error("Not implemented")
            # lr_::Float32 = model.lr
            # cpu_sparse_update!(model.α, ∇αβ, idxes, lr_, model.λ)
        end
    end
end

function update_sparse_u1!(model, ∇αβ, grad_αβ, uniq_batch)
    idx, inv_idxes = uniq_batch

    ∇αβ = @view ∇αβ[:, 1:length(idx)]
    fill!(∇αβ, 0.0f0)

    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)

    for (i, inv_idx) in enumerate(inv_idxes)
        ∇αβ[:, inv_idx] += grad_αβ[:, i]
    end

    proximal_update!(model, ∇αβ, idx)
    nothing
end

function update_sparse_u2!(model, ∇αβ, grad_αβ, uniq_batch)
    idx, inv_idxes = uniq_batch

    ∇αβ = @view ∇αβ[:, 1:length(idx)]
    fill!(∇αβ, 0.0f0)

    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)

    for (i, inv_idx) in enumerate(inv_idxes)
        ∇αβ[:, inv_idx] += grad_αβ[:, i]
    end
    normalize_columns!(∇αβ)

    proximal_update!(model, ∇αβ, idx)
    nothing
end

# function update_all!(model, grad_cache, grad, uniq_batch)
#     lr_::Float32 = model.lr
#     update_αβ!(model, grad_cache, grad, uniq_batch)
#     axpy!(-lr_, vec(model.∇B), vec(model.B))
#     axpy!(-lr_, vec(model.∇C), vec(model.C))
#     nothing
# end

@generated function sync_params!(cpu_model, model)
    if BACKEND == "CPU"
        return :(cpu_model)
    else
        quote
            for fn in (:α, :B, :C)  #fieldnames(typeof(model))
                src = getfield(model, fn)
                dest = getfield(cpu_model, fn)
                copy!(dest, src)
            end
            CUDArt.synchronize(CUDArt.null_stream)
            cpu_model
        end
    end
end

const _start_time_loginfo = Ref{Float64}(0.0)
const _next_report_time = Ref{Float64}(0.0)
const _sigmas_loginfo = Ref{HostArray{Float32, 2}}()

function log_progress(model, current, progress=progress)
    global _start_time_loginfo
    global _next_report_time
    global _sigmas_loginfo
    global sampleFactory
    factory = get(sampleFactory)

    if current >= _next_report_time[]
        σ = _sigmas_loginfo[]
        copy!(σ, model.σ)
        CUDArt.synchronize(CUDArt.null_stream)
        loss = 0.0f0
        for x in σ
            loss += -log(x+eps(typeof(x)))
        end
        ppl = exp(loss/(model.m*(model.k+1)))

        _next_report_time[] += 3.0
        st = _start_time_loginfo[]
        factory = get(sampleFactory)
        prog = progress(factory) * 100.0
        elapsed = current - st
        etc = (100.0 - prog) / prog * elapsed
        hours = div(etc, 3600)
        rem = etc-hours*3600
        mins = div(rem, 60)
        secs = rem - mins * 60

        mean_speed = factory.consumed[]/elapsed
        info(@sprintf("ETA: %02d:%02d:%02d - loss: %.4f - ppl: %.3f - progress: %5.2f%% - speed: %.1f %s/s - lr: %g", hours, mins, secs, loss, ppl, prog, mean_speed, factory.corpus_type == String ? "lines" : "bytes", model.lr))
    end
    nothing
end

function init_logger(model, start_time)
    global _start_time_loginfo
    global _next_report_time
    global _sigmas_loginfo
    _start_time_loginfo[] = start_time
    _next_report_time[] = start_time
    _sigmas_loginfo[] = HostArray(Float32, (model.k+1), model.m)
end

function train(model; min_lr=1f-4, normalize_sp=false,
        normalize_dict=false, every=500, preratio=0.1,
        saveratio=0.05, part=1.0, save_basename="B")
    if every == 1
        train_f1(model, min_lr=min_lr, normalize_sp=normalize_sp,
        normalize_dict=normalize_dict, preratio=preratio,
        saveratio=saveratio, part=part, save_basename=save_basename)
        return nothing
    end
    global sampleFactory

    if normalize_dict
        global ∇B_cache
        global ∇C_cache
        ∇B_cache[] = HostArray(Float32, (num_vecdim(model), num_base(model)))
        ∇C_cache[] = HostArray(Float32, (num_vecdim(model), num_base(model)))
    end

    cpu_model = model
    model = to_device(cpu_model)

    beta = model.lr0 / min_lr - 1.0f0

    info(@sprintf("Train with lr=%g, normalize_sp=%s, normalize_dict=%s. Update dictionary parameters every %d updates of sparse parameters. Pretrain sparse parameters with %.2f%% amount of training data. Save after every %.2f%% of training. Use %.2f%% of the total training data. Save the trained model to \"%s*.jld\"", model.lr, string(normalize_sp), string(normalize_dict), every, preratio*100, saveratio*100, part*100, save_basename))

    st = time()
    init_logger(model, st)

    grads = HostArray(Float32, (num_base(model), (model.k+2)*model.m))
    grad_cache = HostArray(Float32, (num_base(model), (model.k+2)*model.m))

    factory = get(sampleFactory)
    const producer = factory.producer
    batch_set = Task(() -> producer(factory, model.m, model.k))

    const update_sparse! = normalize_sp ? update_sparse_u2! : update_sparse_u1!
    const update_dict! = normalize_dict ? update_dict_u2! : update_dict_u1!
    const progress_ = factory -> progress(factory)/part

    batch = consume(batch_set)
    old_batch = batch
    compute_grad_sparse(model, old_batch) # this method immediately returns when
                                  # the computations are taking place on a GPU.
                                  # this makes the gpu and cpu runs concurrently
    for batch in batch_set
        uniq_batch = invertable_unique(old_batch[1], old_batch[2])
        fetch_grad_sparse!(grads, model)
        update_sparse!(model, grad_cache, grads, uniq_batch)

        old_batch = batch
        prog = progress_(factory)
        if prog >= preratio
            if preratio > 0.0f0
                info("Pre-train finished.")
            end
            break
        end
        compute_grad_sparse(model, old_batch)
        log_progress(model, time(), progress_)
        model.lr = model.lr0 / (1.0f0 + beta*prog)
    end

    next_save_point = preratio + saveratio
    compute_grad_sparse(model, old_batch)
    nb_batch = 0
    for batch in batch_set
        uniq_batch = invertable_unique(old_batch[1], old_batch[2])
        # synchronize(CUDArt.null_stream)
        # TODO: maybe another stream to handle the data transferring
        fetch_grad_sparse!(grads, model)

        accumulate_grad_dict(model, nb_batch)
        if nb_batch == every - 1
            update_dict!(model)
        end
        update_sparse!(model, grad_cache, grads, uniq_batch)

        old_batch = batch
        compute_grad_sparse(model, old_batch)
        prog = progress_(factory)
        if prog >= next_save_point
            if prog >= 1.0f0
                break
            end
            # synchronize
            sync_params!(cpu_model, model)
            sn = joinpath(project_root,
                @sprintf("%s-checkpoint-%03.0f.jld", save_basename, prog*100.0f0))
            info("Save point reached. saveing to $sn")
            save(cpu_model, sn)
            next_save_point += saveratio
        end
        log_progress(model, time(), progress_)
        model.lr = model.lr0 / (1.0f0 + beta*prog)
        nb_batch = (nb_batch+1) % every
    end

    sync_params!(cpu_model, model)
    sn = joinpath(project_root,
        @sprintf("%s-checkpoint-100.jld", save_basename))
    info("Training finished. Saving to $sn")
    save(cpu_model, sn, save_maps=true)
    nothing
end

function train_f1(model; min_lr=1f-4, normalize_sp=false,
        normalize_dict=false, preratio=0.1,
        saveratio=0.05, part=1.0, save_basename="B")
    global sampleFactory

    if normalize_dict
        global ∇B_cache
        global ∇C_cache
        ∇B_cache[] = HostArray(Float32, (num_vecdim(model), num_base(model)))
        ∇C_cache[] = HostArray(Float32, (num_vecdim(model), num_base(model)))
    end

    cpu_model = model
    model = to_device(cpu_model)

    beta = model.lr0 / min_lr - 1.0f0

    info(@sprintf("Train with lr=%g, normalize_sp=%s, normalize_dict=%s. Update dictionary parameters every 1 updates of sparse parameters. Pretrain sparse parameters with %.2f%% amount of training data. Save after every %.2f%% of training. Use %.2f%% of the total training data. Save the trained model to \"%s*.jld\"", model.lr, string(normalize_sp), string(normalize_dict), preratio*100, saveratio*100, part*100, save_basename))

    st = time()
    init_logger(model, st)

    grads = HostArray(Float32, (num_base(model), (model.k+2)*model.m))
    grad_cache = HostArray(Float32, (num_base(model), (model.k+2)*model.m))

    factory = get(sampleFactory)
    const producer = factory.producer
    batch_set = Task(() -> producer(factory, model.m, model.k))

    const update_sparse! = normalize_sp ? update_sparse_u2! : update_sparse_u1!
    const update_dict! = normalize_dict ? update_dict_u2f1! : update_dict_u1f1!
    const progress_ = factory -> progress(factory)/part

    batch = consume(batch_set)
    old_batch = batch
    compute_grad_sparse(model, old_batch) # this method immediately returns when
                                  # the computations are taking place on a GPU.
                                  # this makes the gpu and cpu runs concurrently
    for batch in batch_set
        uniq_batch = invertable_unique(old_batch[1], old_batch[2])
        fetch_grad_sparse!(grads, model)
        update_sparse!(model, grad_cache, grads, uniq_batch)

        old_batch = batch
        prog = progress_(factory)
        if prog >= preratio
            if preratio > 0.0f0
                info("Pre-train finished.")
            end
            break
        end
        compute_grad_sparse(model, old_batch)
        log_progress(model, time(), progress_)
        model.lr = model.lr0 / (1.0f0 + beta*prog)
    end

    next_save_point = preratio + saveratio
    compute_grad_sparse(model, old_batch)
    for batch in batch_set
        uniq_batch = invertable_unique(old_batch[1], old_batch[2])
        # synchronize(CUDArt.null_stream)
        # TODO: maybe another stream to handle the data transferring
        fetch_grad_sparse!(grads, model)
        update_dict!(model)
        update_sparse!(model, grad_cache, grads, uniq_batch)

        old_batch = batch
        compute_grad_sparse(model, old_batch)
        prog = progress_(factory)
        if prog >= next_save_point
            if prog >= 1.0f0
                break
            end
            # synchronize
            sync_params!(cpu_model, model)
            sn = joinpath(project_root,
                @sprintf("%s-checkpoint-%03.0f.jld", save_basename, prog*100.0f0))
            info("Save point reached. saveing to $sn")
            save(cpu_model, sn)
            next_save_point += saveratio
        end
        log_progress(model, time(), progress_)
        model.lr = model.lr0 / (1.0f0 + beta*prog)
    end

    sync_params!(cpu_model, model)
    sn = joinpath(project_root,
        @sprintf("%s-checkpoint-100.jld", save_basename))
    info("Training finished. Saving to $sn")
    save(cpu_model, sn, save_maps=true)
    nothing
end

end # end of module


text:
@testset "Deprecation" begin
Lint.addDummyDeprecates()

s = """
function testDep1(x)
    x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E211
@test occursin("generic deprecate message", msgs[1].message)

# THIS DOESN'T TRIGGER LINT WARNING SINCE THE SIGNATURE DOESN'T MATCH
s = """
function testDep2(x::Int, y)
    x + y
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function testDep2(x::Complex{Int})
    x
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function testDep2(x::Int...)
    x
end
"""
msgs = lintstr(s)
@test isempty(msgs)
# END OF NON-MATCHING SIGNATURES

s = """
function testDep2(x::Int)
    x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E211
@test occursin("generic deprecate message", msgs[1].message)

s = """
function testDep3{T <: Real}(x::Array{T,1})
    x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E211
@test occursin("generic deprecate message", msgs[1].message)

s = """
function testDep4(x::Int, y::Int...)
    x + length(y)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E211
@test occursin("generic deprecate message", msgs[1].message)

s = """
function testDep4(x::Int, y::Int)
    x + length(y)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E211
@test occursin("generic deprecate message", msgs[1].message)

s = """
function testDep5{T <: AbstractString}(x::Array{T,1})
    x
end
"""
msgs = lintstr(s)
@test isempty(msgs)
end


text:
module GBPirate

import SparseArrays
using SuiteSparseGraphBLAS
import LinearAlgebra: mul!

export gbmul!
function gbmul!(C::StridedVecOrMat{T}, A::SparseArrays.SparseMatrixCSC{T, Int64}, B::StridedVecOrMat{T}, α::Number, β::Number) where {T}

    Agb = GBMatrix{T}(size(A))
    SuiteSparseGraphBLAS._packcscmatrix!(Agb, A.colptr, A.rowval, A.nzval)
    SuiteSparseGraphBLAS._makeshallow!(Agb)

    Bgb = GBMatrix{T}(size(B))
    SuiteSparseGraphBLAS._packdensematrix!(Bgb, B)
    SuiteSparseGraphBLAS._makeshallow!(Bgb)

    Cgb = GBMatrix(size(C), convert(T, α))
    gbset(Cgb, :format, :bycol)
    gbset(Cgb, :sparsity_control, :full)
    mul!(Cgb, Agb, Bgb; accum=*)

    # Done with Agb and Bgb now.
    # They do not need to be unpacked since they're shallow.
    # But packing modified the index vectors of A.
    A.colptr .+= 1
    A.rowval .+= 1
    return Cgb
    Ctemp = SuiteSparseGraphBLAS._unpackdensematrix!(Cgb)

    if iszero(β)
        copyto!(C, Ctemp)
    else
        C .*= Ctemp
    end
    ccall(:jl_free, Cvoid, (Ptr{T},), pointer(Ctemp))
    return C
end
end


text:
module CDLagrangianRefFEsTests

using Gridap
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Test

using Gridap.ReferenceFEs: _CDLagrangianRefFE

T = Float64
reffe = _CDLagrangianRefFE(T,SEGMENT,(2,),(DISC,))
test_lagrangian_reference_fe(reffe)

T = Float64
reffe = _CDLagrangianRefFE(T,QUAD,1,(CONT,DISC))
test_lagrangian_reference_fe(reffe)

conf = CDConformity( (CONT,DISC) )
@test conf == Conformity(reffe)
@test get_face_own_dofs(QUAD4,conf) == get_face_own_dofs(reffe)
@test get_face_own_dofs_permutations(QUAD4,conf) == get_face_own_dofs_permutations(reffe)

T = VectorValue{2,Float64}
reffe = _CDLagrangianRefFE(T,QUAD,(2,0),(CONT,DISC))
test_lagrangian_reference_fe(reffe)

@test get_face_own_dofs(reffe,L2Conformity()) == [[], [], [], [], [], [], [], [], [1, 2, 3, 4, 5, 6]]
@test get_face_own_nodes(reffe,L2Conformity()) == [[], [], [], [], [], [], [], [], [1, 2, 3]]

T = VectorValue{3,Float64}
reffe = _CDLagrangianRefFE(T,HEX,(2,2,2),(CONT,CONT,DISC))
test_lagrangian_reference_fe(reffe)

T = VectorValue{3,Float64}
reffe = _CDLagrangianRefFE(T,HEX,2,(CONT,CONT,DISC))
test_lagrangian_reference_fe(reffe)

T = VectorValue{2,Float64}
reffe = LagrangianRefFE(T,QUAD,(2,0))
@test Conformity(reffe) == CDConformity((CONT,DISC))
test_lagrangian_reference_fe(reffe)

T = VectorValue{3,Float64}
reffe = LagrangianRefFE(T,HEX,(0,2,0))
@test Conformity(reffe) == CDConformity((DISC,CONT,DISC))
test_lagrangian_reference_fe(reffe)

reffe = ReferenceFE(HEX,lagrangian,T,(0,2,0))
@test Conformity(reffe) == CDConformity((DISC,CONT,DISC))
test_lagrangian_reference_fe(reffe)

reffe = LagrangianRefFE(T,QUAD,(2,2))

end # module


text:
using Statistics
using StatsBase: sample
using Base.Threads: @threads
using BenchmarkTools
using Revise
using EvoTrees

n_obs = Int(1e6)
n_vars = 100
n_bins = 255
𝑖 = collect(1:n_obs)
𝑗 = collect(1:n_vars)
δ = rand(n_obs)
δ² = rand(n_obs)

hist_δ = zeros(n_bins, n_vars)
hist_δ² = zeros(n_bins, n_vars)
X_bin = rand(UInt8, n_obs, n_vars)

function iter_1(X_bin, hist_δ, δ, 𝑖, 𝑗)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            hist_δ[X_bin[i,j], j] += δ[i]
        end
    end
end

@time iter_1(X_bin, hist_δ, δ, 𝑖, 𝑗)

# takeaway : significant speedup from depth 3 if building all hit simultaneously
𝑖_4 = sample(𝑖, Int(n_obs/4), ordered=true)
@btime iter_1($X_bin, $hist_δ, $δ, $𝑖, $𝑗)
@btime iter_1($X_bin, $hist_δ, $δ, $𝑖_4, $𝑗)

# try adding all info on single array rather than seperate vectors
function iter_1B(X_bin, hist_δ, hist_δ², δ, δ², 𝑖, 𝑗)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            @inbounds hist_δ[X_bin[i,j], j] += δ[i]
            @inbounds hist_δ²[X_bin[i,j], j] += δ²[i]
        end
    end
end

@btime iter_1B($X_bin, $hist_δ, $hist_δ², $δ, $δ², $𝑖, $𝑗)

# try adding all info on single array rather than seperate vectors
δ2 = rand(2, n_obs)
hist_δ2 = zeros(n_bins, 2, n_vars)
function iter_2(X_bin, hist_δ2, δ2, 𝑖, 𝑗)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            # view(hist_δ2, X_bin[i,j], j, :) .+= view(δ2, i, :)
            @inbounds for k in 1:2
                hist_δ2[X_bin[i,j], k, j] += δ2[k, i]
                # @inbounds hist_δ2[X_bin[i,j], 1, j] += δ2[i, 1]
                # @inbounds hist_δ2[X_bin[i,j], 2, j] += δ2[i, 2]
            end
        end
    end
end
@time iter_2(X_bin, hist_δ2, δ2, 𝑖_4, 𝑗)
@btime iter_2($X_bin, $hist_δ2, $δ2, $𝑖, $𝑗)

# integrate a leaf id
hist_δ2 = zeros(n_bins, 2, n_vars, 8);
@time hist_δ2 .= 0;
@time hist_δ2 .* 0;
function iter_3(X_bin, hist, δ2, 𝑖, 𝑗, leaf)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            @inbounds for k in 1:2
                # view(hist_δ2, X_bin[i,j], j, :) .+= view(δ2, i, :)
                @inbounds hist[X_bin[i,j], k, j, leaf[i]] += δ2[k, i]
            end
        end
    end
end

leaf_vec = ones(UInt8, n_obs)
@time iter_3(X_bin, hist_δ2, δ2, 𝑖_4, 𝑗, leaf_vec);
@btime iter_3($X_bin, $hist_δ2, $δ2, $𝑖, $𝑗, $leaf_vec);


text:
"""
Package: Forecast

    forecast(xar, n; levels = (0.8,.95))

Forecast a univariate or multivariate autoregressive model.

The forecasting follows the model:

```math
Xt = \\Phi_0 + \\sum_{i=1}^p \\Phi_i \\cdot X_{t-i} + E
```

# Arguments
`xar`           AR struct coming from the `ar` function.
`n`             Number of time periods to be forecasted.
`alpha`         Prediction intervals levels; its default value is (0.8, 0.95)
`fixMean`       Fixes the mean in the forecast with a DataFrame which first
                column is a timestamp type and missing values indicate values
                to be estimated. Default value is `nothing`.
`fixΣ2`         fixes Σ2 values in the forecast

# Returns
A FORECAST struct
"""
function forecast(xar, n::Integer;
                  alpha::Tuple{Real,Real} = (0.8,.95),
                  fixMean::Union{Nothing,DataFrame} = nothing,
                  fixΣ2::Union{T,AbstractMatrix{T}} = xar.Σ2) where T<:Real

    @assert n > 0 "n must be greater than 0"

    m = xar.ndims
    np = xar.order

    Φ,Φ0 = xar.Φ, xar.Φ0

    dfts = xar.x = tots(xar.x)
    names_x = names(dfts)[2:end]
    x = Array(dfts[:,2:end])
    ts = dfts[:,1]
    name_ts = names(dfts)[1]

    x0 = compact(x[end:-1:end-np+1,:]')
    PT = promote_type(eltype(Φ),eltype(Φ0))
    x0 = PT.(x0)

    Σ20 = compact(zeros(m,m))
    !isnothing(fixMean) && (fixMean = Array(fixMean[:,2:end]))
    mu = arsim(Φ,Φ0,x0,n; Σ2=Σ20, fix=fixMean)
    se = sqrt.(fvar(Φ,compact(fixΣ2),max(n,np)))[1:n,:]

    # Prediction Intervals
    a1 = alpha[1]
    a2 = alpha[2]

    z1 = mapslices(x -> Distributions.quantile.(Normal.(0,x),
                                               a1 + (1-a1)/2),
                                               reshape(se,n,m), dims = 2)
    z1 = compact(z1)
    lower1 = mu .- z1
    upper1 = mu .+ z1
    z2 = mapslices(x -> Distributions.quantile.(Normal.(0,x),
                                               a2 + (1-a2)/2),
                                               reshape(se,n,m), dims = 2)
    z2 = compact(z2)
    lower2 = mu .- z2
    upper2 = mu .+ z2

    # mu
    mu_names = ["mean_" * n  for n in names_x]
    mu_df = hcat(DataFrame(nΔt(ts,n),[name_ts]),
                 DataFrame(reshape(mu,:,size(mu,2)),mu_names))

    # upper
    u1_names = ["upper1_" * n  for n in names_x]
    upper1_df = DataFrame(reshape(upper1,:,size(upper1,2)),u1_names)
    u2_names = ["upper2_" * n  for n in names_x]
    upper2_df = DataFrame(reshape(upper2,:,size(upper2,2)),u2_names)
    upper_df = hcat(mu_df[:,1:1], upper1_df, upper2_df)

    # lower
    l1_names = ["lower1_" * n  for n in names_x]
    lower1_df = DataFrame(reshape(lower1,:,size(lower1,2)),l1_names)
    l2_names = ["lower2_" * n  for n in names_x]
    lower2_df = DataFrame(reshape(lower2,:,size(lower2,2)),l2_names)
    lower_df = hcat(mu_df[:,1:1],lower1_df, lower2_df)

    # se
    se_names = ["se_" * n  for n in names_x]
    se_df = DataFrame(reshape(se,:,size(se,2)),se_names)
    se_df = hcat(mu_df[:,1:1],se_df)

    return FORECAST(xar, alpha,
                    mu_df,
                    upper_df,
                    lower_df,
                    se_df,
                    "Forecasting " * xar.call)
end

"""
Forecast recursive variance/covariance
"""
function fvar(Φ::AbstractArray{T},
              Σ2::AbstractMatrix{T},
              n::Integer) where T<:Real

    m = size(Φ,1)
    np = size(Φ,3)

    Φ = reshape(Φ,m,m,np)

    Σ2i = zeros(T,m,m)
    v = Array{T}(undef,n,m)

    Φ = reshape(Φ,m,m*np)
    for i in 1:n
        #Σ2i = Φ * Σ2i * Φ' + Σ2
        Σ2i = Φ * repeat(repeat(Σ2i,1,np) * Φ',np,1) + Σ2
        v[i,:] = diag(Σ2i)
    end

    v
end

function fvar(Φ::AbstractVector{T},
              Σ2::T,
              n::Integer) where T<:Real

    np = length(Φ)
    Σ2i = T(0)
    v = Vector{T}(undef,n)
    Φ2 = Φ.^2

    for i in 1:n
        v[i] = Σ2i = sum(Φ2 .* Σ2i) + Σ2
    end

    v
end

function fvar(Φ::T,
              Σ2::T,
              n::Integer) where T<:Real

    Σ2i = T(0)
    v = Vector{T}(undef,n)
    Φ2 = Φ^2

    for i in 1:n
        v[i] = Σ2i = Φ2 * Σ2i + Σ2
    end

    v
end


text:
include("settings.jl")

makedocs(; sitename=sitename, settings...)

deploydocs(; repo="github.com/ITensor/ITensors.jl.git")


text:
# # Machine Learning in Julia (continued)

# An introduction to the
# [MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
# toolbox.

# ### Set-up

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate("env")
Pkg.instantiate()

# ## General resources

# - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
# - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
# - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
# - [Data Science Tutorials in Julia](https://juliaai.github.io/DataScienceTutorials.jl/)

# ## Part 2 - Selecting, Training and Evaluating Models

# > **Goals:**
# > 1. Search MLJ's database of model metadata to identify model candidates for a supervised learning task.
# > 2. Evaluate the performance of a model on a holdout set using basic `fit!`/`predict` work-flow.
# > 3. Inspect the outcomes of training and save these to a file.
# > 3. Evaluate performance using other resampling strategies, such as cross-validation, in one line, using `evaluate!`
# > 4. Plot a "learning curve", to inspect performance as a function of some model hyper-parameter, such as an iteration parameter

# The "Hello World!" of machine learning is to classify Fisher's
# famous iris data set. This time, we'll grab the data from
# [OpenML](https://www.openml.org):

using MLJ
OpenML.describe_dataset(61)

#-

iris = OpenML.load(61); # a column dictionary table

import DataFrames
iris = DataFrames.DataFrame(iris);
first(iris, 4)

# **Main goal.** To build and evaluate models for predicting the
# `:class` variable, given the four remaining measurement variables.

# ### Step 1. Inspect and fix scientific types

schema(iris)

# These look fine.

# ### Step 2. Split data into input and target parts

# Here's how we split the data into target and input features, which
# is needed for MLJ supervised models. We can randomize the data at
# the same time:

y, X = unpack(iris, ==(:class), rng=123);
scitype(y)

# This puts the `:class` column into a vector `y`, and all remaining
# columns into a table `X`.

# Here's one way to access the documentation (at the REPL, `?unpack`
# also works):

@doc unpack #!md

# <display omitted, as not markdown renderable> #md

# ### On searching for a model

# Here's how to see *all* models (not immediately useful):

all_models = models()

# If you already have an idea about the name of the model, you could
# search by string or regex:

some_models = models("LinearRegressor")

# Each entry contains metadata for a model whose defining code is not
# yet loaded:

meta = some_models[1]

#-

targetscitype = meta.target_scitype

#-

scitype(y) <: targetscitype

# So this model won't do. Let's find all pure julia classifiers:

filter_julia_classifiers(meta) =
    AbstractVector{Finite} <: meta.target_scitype &&
    meta.is_pure_julia

models(filter_julia_classifiers)

# Find all (supervised) models that match my data!

models(matching(X, y))

# ### Step 3. Select and instantiate a model

# To load the code defining a new model type we use the `@load` macro:

NeuralNetworkClassifier = @load NeuralNetworkClassifier

# Other ways to load model code are described
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/loading_model_code/#Loading-Model-Code).

# We'll instantiate this type with default values for the
# hyperparameters:

model = NeuralNetworkClassifier()

#-

info(model)

# In MLJ a *model* is just a struct containing hyper-parameters, and
# that's all. A model does not store *learned* parameters. Models are
# mutable:

model.epochs = 12

# And all models have a key-word constructor that works once `@load`
# has been performed:

NeuralNetworkClassifier(epochs=12) == model

# ### On fitting, predicting, and inspecting models

# In MLJ a model and training/validation data are typically bound
# together in a machine:

mach = machine(model, X, y)

# A machine stores *learned* parameters, among other things. We'll
# train this machine on 70% of the data and evaluate on a 30% holdout
# set. Let's start by dividing all row indices into `train` and `test`
# subsets:

train, test = partition(1:length(y), 0.7)

# Now we can `fit!`...

fit!(mach, rows=train, verbosity=2)

# ... and `predict`:

yhat = predict(mach, rows=test);  # or `predict(mach, Xnew)`
yhat[1:3]

# We'll have more to say on the form of this prediction shortly.

# After training, one can inspect the learned parameters:

fitted_params(mach)

#-

# Everything else the user might be interested in is accessed from the
# training *report*:

report(mach)

# You save a machine like this:

MLJ.save("neural_net.jlso", mach)

# And retrieve it like this:

mach2 = machine("neural_net.jlso")
yhat = predict(mach2, X);
yhat[1:3]

# If you want to fit a retrieved model, you will need to bind some data to it:

mach3 = machine("neural_net.jlso", X, y)
fit!(mach3)

# Machines remember the last set of hyper-parameters used during fit,
# which, in the case of iterative models, allows for a warm restart of
# computations in the case that only the iteration parameter is
# increased:

model.epochs = model.epochs + 4
fit!(mach, rows=train, verbosity=2)

# For this particular model we can also increase `:learning_rate`
# without triggering a cold restart:

model.epochs = model.epochs + 4
model.optimiser.eta = 10*model.optimiser.eta
fit!(mach, rows=train, verbosity=2)

# However, change any other parameter and training will restart from
# scratch:

model.lambda = 0.001
fit!(mach, rows=train, verbosity=2)

# Iterative models that implement warm-restart for training can be
# controlled externally (eg, using an out-of-sample stopping
# criterion). See
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/)
# for details.

# Let's train silently for a total of 50 epochs, and look at a
# prediction:

model.epochs = 50
fit!(mach, rows=train)
yhat = predict(mach, X[test,:]); # or predict(mach, rows=test)
yhat[1]

# What's going on here?

info(model).prediction_type

# **Important**: - In MLJ, a model that can predict probabilities (and
# not just point values) will do so by default.  - For most
# probabilistic predictors, the predicted object is a
# `Distributions.Distribution` object or a
# `CategoricalDistributions.UnivariateFinite` object (the case here)
# which all support the following methods: `rand`, `pdf`, `logpdf`;
# and, where appropriate: `mode`, `median` and `mean`.

# So, to obtain the probability of "Iris-virginica" in the first test
# prediction, we do

pdf(yhat[1], "Iris-virginica")

# To get the most likely observation, we do

mode(yhat[1])

# These can be broadcast over multiple predictions in the usual way:

broadcast(pdf, yhat[1:4], "Iris-versicolor")

#-

mode.(yhat[1:4])

# Or, alternatively, you can use the `predict_mode` operation instead
# of `predict`:

predict_mode(mach, X[test,:])[1:4] # or predict_mode(mach, rows=test)[1:4]

# For a more conventional matrix of probabilities you can do this:

L = levels(y)
pdf(yhat, L)[1:4, :]

# However, in a typical MLJ work-flow, this is not as useful as you
# might imagine. In particular, all probabilistic performance measures
# in MLJ expect distribution objects in their first slot:

cross_entropy(yhat, y[test]) |> mean

# To apply a deterministic measure, we first need to obtain point-estimates:

misclassification_rate(mode.(yhat), y[test])

# We note in passing that there is also a search tool for measures
# analogous to `models`:

measures()

# ### Step 4. Evaluate the model performance

# Naturally, MLJ provides boilerplate code for carrying out a model
# evaluation with a lot less fuss. Let's repeat the performance
# evaluation above and add an extra measure, `brier_score`:

evaluate!(mach, resampling=Holdout(fraction_train=0.7),
          measures=[cross_entropy, misclassification_rate, brier_score])

# Or applying cross-validation instead:

evaluate!(mach, resampling=CV(nfolds=6),
          measures=[cross_entropy, misclassification_rate, brier_score])

# Or, Monte Carlo cross-validation (cross-validation repeated
# randomized folds)

e = evaluate!(mach, resampling=CV(nfolds=6, rng=123),
              repeats=3,
              measures=[cross_entropy, misclassification_rate, brier_score])

# One can access the following properties of the output `e` of an
# evaluation: `measure`, `measurement`, `per_fold` (measurement for
# each fold) and `per_observation` (measurement per observation, if
# reported).

# We finally note that you can restrict the rows of observations from
# which train and test folds are drawn, by specifying `rows=...`. For
# example, imagining the last 30% of target observations are `missing`
# you might have a work-flow like this:

train, test = partition(eachindex(y), 0.7)
mach = machine(model, X, y)
evaluate!(mach, resampling=CV(nfolds=6),
          measures=[cross_entropy, brier_score],
          rows=train)     # cv estimate, resampling from `train`
fit!(mach, rows=train)    # re-train using all of `train` observations
predict(mach, rows=test); # and predict missing targets

# ### On learning curves

# Since our model is an iterative one, we might want to inspect the
# out-of-sample performance as a function of the iteration
# parameter. For this we can use the `learning_curve` function (which,
# incidentally can be applied to any model hyper-parameter). This
# starts by defining a one-dimensional range object for the parameter
# (more on this when we discuss tuning in Part 4):

r = range(model, :epochs, lower=1, upper=50, scale=:log10)

#-

curve = learning_curve(mach,
                       range=r,
                       resampling=Holdout(fraction_train=0.7), # (default)
                       measure=cross_entropy)

#-

using Plots
gr(size=(490,300))
plt=plot(curve.parameter_values, curve.measurements)
xlabel!(plt, "epochs")
ylabel!(plt, "cross entropy on holdout set")
savefig("learning_curve.png")
plt #!md
# ![](learning_curve.png) #md

# We will return to learning curves when we look at tuning in Part 4.

# ### Resources for Part 2

# - From the MLJ manual:
#     - [Getting Started](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/)
#     - [Model Search](https://alan-turing-institute.github.io/MLJ.jl/dev/model_search/)
#     - [Evaluating Performance](https://alan-turing-institute.github.io/MLJ.jl/dev/evaluating_model_performance/) (using `evaluate!`)
#     - [Learning Curves](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
#     - [Performance Measures](https://alan-turing-institute.github.io/MLJ.jl/dev/performance_measures/) (loss functions, scores, etc)
# - From Data Science Tutorials:
#     - [Choosing and evaluating a model](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/choosing-a-model/)
#     - [Fit, predict, transform](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/fit-and-predict/)

# ### Exercises for Part 2

# #### Exercise 4

# (a) Identify all supervised MLJ models that can be applied (without
# type coercion or one-hot encoding) to a supervised learning problem
# with input features `X4` and target `y4` defined below:

import Distributions
poisson = Distributions.Poisson

age = 18 .+ 60*rand(10);
salary = coerce(rand(["small", "big", "huge"], 10), OrderedFactor);
levels!(salary, ["small", "big", "huge"]);
small = salary[1]

#-

X4 = DataFrames.DataFrame(age=age, salary=salary)

n_devices(salary) = salary > small ? rand(poisson(1.3)) : rand(poisson(2.9))
y4 = [n_devices(row.salary) for row in eachrow(X4)]

# (b) What models can be applied if you coerce the salary to a
# `Continuous` scitype?

# #### Exercise 5 (unpack)

# After evaluating the following ...

data = (a = [1, 2, 3, 4],
        b = rand(4),
        c = rand(4),
        d = coerce(["male", "female", "female", "male"], OrderedFactor));
pretty(data)

#-

using Tables

y, X, w = unpack(data,
                 ==(:a),
                 name -> elscitype(Tables.getcolumn(data, name)) == Continuous);

# ...attempt to guess the evaluations of the following:

y

#-

pretty(X)

#-

w

# #### Exercise 6 (first steps in modeling Horse Colic)

# Here is the Horse Colic data introduced in Part 1, together with the
# type coercions we performed there:

using UrlDownload, CSV
csv_file = urldownload("https://raw.githubusercontent.com/ablaom/"*
                   "MachineLearningInJulia2020/"*
                   "for-MLJ-version-0.16/data/horse.csv");
horse = DataFrames.DataFrame(csv_file); # convert to data frame
coerce!(horse, autotype(horse));
coerce!(horse, Count => Continuous);
coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
schema(horse)

# (a) Suppose we want to use predict the `:outcome` variable, based on
# the remaining variables that are `Continuous` (one-hot encoding
# categorical variables is discussed later in Part 3) *while ignoring
# the others*.  Extract from the `horse` data set (defined in Part 1)
# appropriate input features `X` and target variable `y`. (Do not,
# however, randomize the observations.)

# (b) Create a 70:30 `train`/`test` split of the data and train a
# `LogisticClassifier` model, from the `MLJLinearModels` package, on
# the `train` rows. Use `lambda=100` and default values for the
# other hyper-parameters. (Although one would normally standardize
# (whiten) the continuous features for this model, do not do so here.)
# After training:

# - (i) Recalling that a logistic classifier (aka logistic regressor) is
#   a linear-based model learning a *vector* of coefficients for each
#   feature (one coefficient for each target class), use the
#   `fitted_params` method to find this vector of coefficients in the
#   case of the `:pulse` feature. (You can convert a vector of pairs `v =
#   [x1 => y1, x2 => y2, ...]` into a dictionary with `Dict(v)`.)

# - (ii) Evaluate the `cross_entropy` performance on the `test`
#   observations.

# - ☆(iii) In how many `test` observations does the predicted
#   probability of the observed class exceed 50%?

# - (iv) Find the `misclassification_rate` in the `test`
#   set. (*Hint.* As this measure is deterministic, you will either
#   need to broadcast `mode` or use `predict_mode` instead of
#   `predict`.)

# (c) Instead use a `RandomForestClassifier` model from the
#     `DecisionTree` package and:
#
# - (i) Generate an appropriate learning curve to convince yourself
#   that out-of-sample estimates of the `cross_entropy` loss do not
#   substantially improve for `n_trees > 50`. Use default values for
#   all other hyper-parameters, and feel free to use all available
#   data to generate the curve.

# - (ii) Fix `n_trees=90` and use `evaluate!` to obtain a 9-fold
#   cross-validation estimate of the `cross_entropy`, restricting
#   sub-sampling to the `train` observations.

# - (iii) Now use *all* available data but set
#   `resampling=Holdout(fraction_train=0.7)` to obtain a score you can
#   compare with the `KNNClassifier` in part (b)(iii). Which model is
#   better?

#