text:
#using XfromProjections
#using XfromProjections.curve_utils
#using XfromProjections.snake_forward
using LinearAlgebra
using IterTools

include("../estimate_angles.jl")
include("../utils.jl")
include("../../phantoms/sperm_phantom.jl")
include("../snake_forward.jl")
include("../snake.jl")
include("../curve_utils.jl")

using Plots
using Colors
using Random

Random.seed!(0)

function reconstruct(template, r, angles, bins, projection, max_iter, w, s, degree, k, length)
    num_points = size(template, 1)
    centerline_points = recon2d_tail(deepcopy(template),r,angles,bins,projection,max_iter, s, w, degree)
    residual = parallel_forward(get_outline(centerline_points, r)[1], angles, bins) - projection
    centerline_points = keep_best_parts(residual, centerline_points, ang, bins, k, num_points, length, projection, r)
    centerline_points = recon2d_tail(deepcopy(centerline_points),r,[ang],bins,projection,1000, s, w, degree)
    return centerline_points
end

function plot_update(curve, residual, name)
    #if norm(residual) < tolerance
    @info "plotted"
    label_txt = @sprintf "%s_%f" name norm(residual)
    plot!(curve[:,1], curve[:,2], aspect_ratio=:equal, label=label_txt, linewidth=2)
    #end
end

cwd = @__DIR__
savepath = normpath(joinpath(@__DIR__, "result"))
!isdir(savepath) && mkdir(savepath)
r(s) = 1.0
detmin, detmax = -36.5, 36.5
grid = collect(detmin:0.1:detmax)
images, tracks = get_sperm_phantom(301,r,grid)

bin_width = 0.125
bins = collect(detmin:bin_width:detmax)

ang = 0.0
angles, max_iter, stepsize = [ang], 10000, 0.1
tail_length = curve_lengths(tracks[end])[end]
num_points = 30

frames2reconstruct = collect(101:10:300)
reconstructions = zeros(num_points,2,length(frames2reconstruct)+1)
#Add actual track at the end so rand sperm works and we can compare timewise
centerline_points = tracks[frames2reconstruct[end]+10]
t = curve_lengths(centerline_points)
spl = ParametricSpline(t,centerline_points',k=1, s=0.0)
tspl = range(0, t[end], length=num_points)
reconstructions[:,:,end] = spl(tspl)'
tolerance = 2.0
k=8
s=0.0
degree=1
for (iter, frame_nr) in Base.Iterators.reverse(enumerate(frames2reconstruct))
    @info iter frame_nr

    #Get projection
    @info "making forward projection for frame: " frame_nr
    outline, normals = get_outline(tracks[frame_nr], r)
    projection = parallel_forward(outline, [ang], bins)

    head = tracks[frame_nr][1,:]
    #Add noise
    @info "adding gaussian noise at level 0.01"
    rho = 0.01
    e = randn(size(projection));
    e = rho*norm(projection)*e/norm(e);
    projection = projection + e;

    #Get ground truth of same length as recon
    t = curve_lengths(tracks[frame_nr])
    spl = ParametricSpline(t,tracks[frame_nr]', k=1, s=0.0)
    tspl = range(0, t[end], length=num_points)
    gt = spl(tspl)'

    angles = [ang]

    w_u = ones(num_points*2+2).*stepsize
    w_u[num_points+4:end] .= 0.0
    w_u[[1,2]] .=0.0
    #TEST
    w_u[num_points+1] = 0.0
    w_u[num_points+3] = 0.0

    w_l = ones(num_points*2+2).*stepsize
    w_l[1:num_points] .= 0.0
    w_l[end] = 0.0
    #TEST
    w_l[num_points+1] = 0.0
    w_l[num_points+3] = 0.0

    recon1_ok = false
    recon2_ok = false

    heatmap(grid, grid, Gray.(images[:,:,frame_nr]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=true)
    cd(savepath)
    @info "setting up template"
    rebinned_bins, rebinned_projection = rebin(projection[:,1], bins, 39)
    #template = estimate_curve(rebinned_projection, rebinned_bins, angles, head, r, num_points)
    template = get_straight_template(projection[:,1], r, [0.0 0.0], ang, num_points,bins)

    best_recon = deepcopy(template)
    best_cost = cost(template, ang, bins, r, projection[:,1])
    plot_update(template,best_cost, "template")

    @info "calculating initial reconstruction"
    #Reconstruct with weights only on one side

    #recon1 = recon2d_tail(deepcopy(template),r,[ang],bins,projection,1, 0.0, w_u, 1, doplot=true)

    recon1 = reconstruct(template, r, angles, bins, projection, max_iter, w_u, s, degree, k, tail_length)
    recon2 = reconstruct(template, r, angles, bins, projection, max_iter, w_l, s, degree, k, tail_length)

    best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)
    plot_update(recon1, residual1, "recon1")
    plot_update(recon2, residual2, "recon2")
    plot_update(best_recon, best_cost, "initial")

    for i = 1:1
        # if best_cost < 1.0
        #     break
        # end
        @info "checking if any parts could need mirroring"
        initial1 = deepcopy(recon1)
        initial2 = deepcopy(recon2)

        for flip_pt=1:num_points
            recon1_flipped = flip(initial1,flip_pt,ang)
            recon2_flipped = flip(initial2,flip_pt,ang)

            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1_flipped, recon2_flipped, ang, bins, projection, best_recon, tail_length, r)

            #mirror and reconstruct with weights on both sides
            recon1 = reconstruct(recon1_flipped, r, angles, bins, projection, 100, w_u, s, degree, k, tail_length)
            recon2 = reconstruct(recon1_flipped, r, angles, bins, projection, 100, w_l, s, degree, k, tail_length)
            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)

            recon1 = reconstruct(recon2_flipped, r, angles, bins, projection, 100, w_u, s, degree, k, tail_length)
            recon2 = reconstruct(recon2_flipped, r, angles, bins, projection, 100, w_l, s, degree, k, tail_length)
            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)
        end
        plot_update(best_recon, best_cost, "flip1")
        # if best_cost < 1.0
        #     break
        # end

        for (flip_i,flip_j) in subsets(1:num_points, Val{2}())
            recon1_flipped = flip(initial1,flip_i,ang)
            recon1_flipped = flip(recon1_flipped,flip_j,ang)
            recon2_flipped = flip(initial2,flip_i,ang)
            recon2_flipped = flip(recon2_flipped,flip_j,ang)

            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1_flipped, recon2_flipped, ang, bins, projection, best_recon, tail_length, r)

            recon1 = reconstruct(recon1_flipped, r, angles, bins, projection, 100, w_u, s, degree, k, tail_length)
            recon2 = reconstruct(recon1_flipped, r, angles, bins, projection, 100, w_l, s, degree, k, tail_length)
            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)

            recon1 = reconstruct(recon2_flipped, r, angles, bins, projection, 100, w_u, s, degree, k, tail_length)
            recon2 = reconstruct(recon2_flipped, r, angles, bins, projection, 100, w_l, s, degree, k, tail_length)
            best_cost, best_recon[:,:], residual1, residual2 = try_improvement(best_cost, recon1, recon2, ang, bins, projection, best_recon, tail_length, r)
        end
        plot_update(best_recon, best_cost, "flip2")
    end
    reconstructions[:,:,iter] = best_recon
    #plot!(best_recon[:,1], best_recon[:,2], aspect_ratio=:equal, label=best_cost, linewidth=2)
    savefig(@sprintf "heuristic_result_%d" frame_nr)
    heatmap(grid, grid, Gray.(images[:,:,frame_nr]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    plot!(best_recon[:,1], best_recon[:,2], aspect_ratio=:equal, linewidth=2)
    mirror = flip(best_recon,1,ang)
    plot!(mirror[:,1], mirror[:,2], aspect_ratio=:equal, linewidth=2)
    savefig(@sprintf "result_%d" frame_nr)
end

global ps = AbstractPlot[]
for (iter, frame_nr) in enumerate(frames2reconstruct)
    best_recon = reconstructions[:,:,iter]
    mirror = flip(best_recon,1,ang)

    l = @layout [a b c]
    p = heatmap(grid, grid, Gray.(images[:,:,frame_nr]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    plot!(best_recon[:,1],best_recon[:,2], aspect_ratio=:equal, linewidth=5)
    plot!(mirror[:,1], mirror[:,2], aspect_ratio=:equal, linewidth=5)
    push!(ps,p)
    if length(ps) == 3
        plot(ps[1], ps[2], ps[3], layout = l, size=(2000,600), linewidth=5)
        savefig(@sprintf "result_all_%d" frame_nr)
        global ps = AbstractPlot[]
    end
end
cd(cwd)


text:
using Revise, CSV, DataFrames
using Tables

using DecisionTree
using Statistics

file  = "/Users/aa25desh/TimeSerise/GunPoint/GunPoint_TRAIN.ts"
df = DataFrame(CSV.File(file))
mat = Tables.matrix(df)
X = float.(mat[:, 1:149])
y = [Meta.parse(mat[i, 150]).args[3] for i in range(1, stop=49)]

function RandomForestClassifierTS(X, y; n_trees::Int=200, min_interval::Int=3)
    transform_xt = InvFeatureGen(X, n_trees=n_trees, min_interval=min_interval)
    model = DecisionTreeClassifier(pruning_purity_threshold=0.67)
    forest = Array{DecisionTreeClassifier,1}()
    for i in range(1, stop=n_trees)
        mdl = deepcopy(model)
        fit!(mdl, transform_xt[i], y)
        push!(forest, mdl)
    end
    forest
end

function InvFeatureGen(X; n_trees::Int=200, min_interval::Int=3)
    n_samps, series_length = size(X)
    transform_xt = Array{Array{Float64,2},1}()
    n_intervals = floor(Int, sqrt(series_length))
    intervals = zeros(Int, n_trees, 3*n_intervals, 2)
    for i in range(1, stop = n_trees)
       transformed_x = Array{Float64,2}(undef, 3*n_intervals, n_samps)
       for j in range(1, stop = n_intervals)
           intervals[i,j,1] = rand(1:(series_length - min_interval))
           len = rand(1:(series_length - intervals[i,j,1]))
           if len < min_interval
               len = min_interval
           end
           intervals[i,j,2] = intervals[i,j,1] + len
           x = Array(1:len+1)
           Y = X[:, intervals[i,j,1]:intervals[i,j,2]]
           means = mean(Y, dims=2)
           stds =  std(Y, dims=2)
           slope = (mean(transpose(x).*Y, dims=2) -
                    mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
           transformed_x[3*j-2,:] =  means
           transformed_x[3*j-1,:] =  stds
           transformed_x[3*j,:]   =  slope
       end
           push!(transform_xt, transpose(transformed_x))
    end
    return transform_xt
end

features = InvFeatureGen(X)
forest = RandomForestClassifierTS(X, y)

file1  = "/Users/aa25desh/TimeSerise/GunPoint/GunPoint_TEST.ts"
df1 = DataFrame(CSV.File(file1))
mat1 = Tables.matrix(df1)

X1 = float.(mat1[:, 1:149])
y1 = [Meta.parse(mat1[i, 150]).args[3] for i in range(1, stop=149)]

function predict_ff(forest, features)
    vv = zeros(Float64,2)
    for i in range(1, stop=200)
        vv += collect(proba_predict(forest,  (features[i])))
    end
    return vv
end

function proba_predict(forest, X)
    a = 0
    for tree in forest
        if predict(tree, X) == 1
            a += 1
        end
    end
    (a, length(forest)-a)
end

function InvFeatures(X, n_trees::Int=200, min_interval::Int=3)
    transform_xt = []
    series_length, = size(X)
    n_intervals = floor(Int, sqrt(series_length))
    intervals = zeros(Int, n_intervals, 2)
    for i in range(1, stop = n_trees)
        transformed_x = Array{Float64,1}(undef, 3*n_intervals)
        for j in range(1, stop = n_intervals)
            intervals[j,1] = rand(1:(series_length - min_interval))
            len = rand(1:(series_length - intervals[j,1]))
            if len < min_interval
                len = min_interval
            end
            intervals[j,2] = intervals[j,1] + len
            x = Array(1:len+1)
            Y = X[intervals[j,1]:intervals[j,2]]
            means = mean(Y)
            stds  = std(Y)
            slope = (mean(x.*Y) -
                     mean(x)*mean(Y)) / (mean(x.*x) - mean(x)^2)
            transformed_x[3*j-2] = means
            transformed_x[3*j-1] = stds
            transformed_x[3*j] = slope
        end
        push!(transform_xt, vcat(transformed_x))
    end
    return transform_xt
end

for i in range(1, stop=149)
    yy = InvFeatures(X1[i,:])
    a, b = predict_ff(forest, yy)
    c = a > b ? 1 : 2
    if y1[i] == c
        print("1")
    else
        print("0")
    end
end

#=
function predict_f(forest, X) #takes only single instance input of Type "transform_xt"
    a, b = proba_predict(forest, X)
    a > b ? 1 : 2
end
=#

ext: jl


text:
using Test, Sockets, LanguageServer, CSTParser, SymbolServer, SymbolServer.Pkg, StaticLint, JSON
using LanguageServer: Document, get_text, get_offset, get_line_offsets, get_position_at, get_open_in_editor, set_open_in_editor, is_workspace_file, applytextdocumentchanges
import JSONRPC
const LS = LanguageServer
const Range = LanguageServer.Range

# TODO Replace this with a proper mock endpoint
JSONRPC.send(::Nothing, ::Any, ::Any) = nothing
function settestdoc(text)
    empty!(server._documents)
    LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem("testdoc", "julia", 0, text)), server, nothing)

    doc = LanguageServer.getdocument(server, LanguageServer.URI2("testdoc"))
    LanguageServer.parse_all(doc, server)
    doc
end

@testset "LanguageServer" begin

    @testset "document" begin
        include("test_document.jl")
    end
    @testset "communication" begin
        include("test_communication.jl")
    end
    @testset "intellisense" begin
        include("test_intellisense.jl")

        server = LanguageServerInstance(IOBuffer(), IOBuffer(), dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH))
        server.runlinter = true
        server.jr_endpoint = nothing
        LanguageServer.initialize_request(init_request, server, nothing)

        @testset "completions" begin
            include("requests/completions.jl")
        end
        @testset "actions" begin
            include("requests/actions.jl")
        end
        @testset "features" begin
            include("requests/features.jl")
        end
        @testset "hover" begin
            include("requests/hover.jl")
        end
        @testset "textdocument" begin
            include("requests/textdocument.jl")
        end
        @testset "misc" begin
            include("requests/misc.jl")
        end
    end
    @testset "edit" begin
        include("test_edit.jl")
    end
    @testset "actions" begin
        include("test_actions.jl")
    end
    @testset "paths" begin
        include("test_paths.jl")
    end

end


text:
# Test for load-order-specific issues with WebIO and IJulia
# when running outside of Jupyter.
#
# See: https://github.com/JuliaGizmos/WebIO.jl/issues/117
using WebIO
using IJulia
@testset "IJulia initialization" begin
    @test :ijulia in WebIO.providers_initialised
end

# Test the demo jupyter notebook and make sure all of its cells execute
# without error.
using NBInclude

@testset "Demo notebook" begin
    @nbinclude(joinpath(@__DIR__, "assets", "jupyter-test.ipynb"))
end


text:
using Revise
using StaticArrays
using ZigZagBoomerang
using LinearAlgebra
using SparseArrays
using Random

Random.seed!(1)

const d = 2
const 𝕏 = SArray{Tuple{d},Float64,1,d}

n = 5

Γ = sprand(n, n, 0.1).*[0.25*SMatrix{d,d}(randn(4)) for i in 1:n, j in 1:n]
Γ = Γ + Γ' + Diagonal(fill(SMatrix{d,d}(1.0I), n))

∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x)

B
t0 = 0.0
x0 = randn(𝕏, n)
θ0 = [randn(𝕏) for i in 1:n]

μ = 0*x0
c = [50.0 for i in 1:n]
σ = [SMatrix{d,d}(1.0I) for i in 1:n]
Z = ZigZag(Γ, μ, σ; λref=0.05, ρ=0.8)
T = 2000.0

@time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)
xs = last.(collect(discretize(trace, 0.01)))

using Makie
using Colors
using GoldenSequences
cs = map(x->RGB(x...), (Iterators.take(GoldenSequence(3), n)))

p1 = scatter(Point2f0.(x0), markersize=0.1, color=cs)
#for x in xs
#    scatter!(p1, Point2f0.(x), markersize=0.1)
#end
for i in 1:n
    scatter!(p1, ((Point2f0∘getindex).(xs, i)), color=cs[i], markersize=0.01)

end
p1


text:
using PyPlot

#we can overload the default plot method
import PyPlot.bar
function bar{T <: Any,M <: Any}(h::ErrorHistogram{T,1,M}; kwargs...)
    es = edges(h)[1]
    cs = contents(h)
    PyPlot.bar(es[1:end-1], cs, width=diff(es); kwargs...)
end

import PyPlot.plot
function plot{T <: Any,M <: Any}(h::ErrorHistogram{T,1,M}; kwargs...)
    es = edges(h)[1]
    cs = contents(h)
    PyPlot.plot(Base.midpoints(es), cs; kwargs...)
end

import PyPlot.errorbar
function errorbar{T <: Any,M <: Any}(h::ErrorHistogram{T,1,M}; kwargs...)
    es = edges(h)[1]
    cs = contents(h)
    errs = errors(h)
    PyPlot.errorbar(Base.midpoints(es), cs, errs; kwargs...)
end


text:
function updatela!(system::CVSystem,n::Int64)
    # mitral valve state
    if system.heart.la.P[n+1] > system.heart.lv.P[n+1]
        MV = "open"
    else
        MV = "closed"
    end

    # left atrial V, Q (P updated elsewhere)
    system.heart.la.V[n+2] = (system.heart.la.V[n+1]+
        system.solverparams.h*(system.lungs.Qv[n+1,2]-
        system.heart.la.Q[n+1]));
    if MV == "open"
        system.heart.la.Q[n+2] = ((1 - system.heart.la.R/
            system.heart.la.L*system.solverparams.h)*
            system.heart.la.Q[n+1] + system.solverparams.h/
            system.heart.la.L*(system.heart.la.P[n+1]-
            system.heart.lv.P[n+1]));
    else
        system.heart.la.Q[n+2] = 0;
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
import ForwardDiff
import AbstractFFTs

# original PR by mcabbott: https://github.com/JuliaDiff/ForwardDiff.jl/pull/495

ForwardDiff.value(x::Complex{<:ForwardDiff.Dual}) = Complex(x.re.value, x.im.value)

ForwardDiff.partials(x::Complex{<:ForwardDiff.Dual}, n::Int) =
    Complex(ForwardDiff.partials(x.re, n), ForwardDiff.partials(x.im, n))

ForwardDiff.npartials(x::Complex{<:ForwardDiff.Dual{T,V,N}}) where {T,V,N} = N
ForwardDiff.npartials(::Type{<:Complex{<:ForwardDiff.Dual{T,V,N}}}) where {T,V,N} = N

ForwardDiff.tagtype(x::Complex{<:ForwardDiff.Dual{T,V,N}}) where {T,V,N} = T
ForwardDiff.tagtype(::Type{<:Complex{<:ForwardDiff.Dual{T,V,N}}}) where {T,V,N} = T

# AbstractFFTs.complexfloat(x::AbstractArray{<:ForwardDiff.Dual}) = float.(x .+ 0im)
AbstractFFTs.complexfloat(x::AbstractArray{<:ForwardDiff.Dual}) = AbstractFFTs.complexfloat.(x)
AbstractFFTs.complexfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = convert(ForwardDiff.Dual{T,float(V),N}, d) + 0im

AbstractFFTs.realfloat(x::AbstractArray{<:ForwardDiff.Dual}) = AbstractFFTs.realfloat.(x)
AbstractFFTs.realfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = convert(ForwardDiff.Dual{T,float(V),N}, d)

for plan in [:plan_fft, :plan_ifft, :plan_bfft]
    @eval begin
        AbstractFFTs.$plan(x::AbstractArray{<:ForwardDiff.Dual}, region=1:ndims(x); kwargs...) =
            AbstractFFTs.$plan(ForwardDiff.value.(x) .+ 0im, region; kwargs...)

        AbstractFFTs.$plan(x::AbstractArray{<:Complex{<:ForwardDiff.Dual}}, region=1:ndims(x); kwargs...) =
            AbstractFFTs.$plan(ForwardDiff.value.(x), region; kwargs...)
    end
end

# rfft only accepts real arrays
AbstractFFTs.plan_rfft(x::AbstractArray{<:ForwardDiff.Dual}, region=1:ndims(x); kwargs...) =
    AbstractFFTs.plan_rfft(ForwardDiff.value.(x), region; kwargs...)

for plan in [:plan_irfft, :plan_brfft]  # these take an extra argument, only when complex?
    @eval begin
        AbstractFFTs.$plan(x::AbstractArray{<:ForwardDiff.Dual}, region=1:ndims(x); kwargs...) =
            AbstractFFTs.$plan(ForwardDiff.value.(x) .+ 0im, region; kwargs...)

        AbstractFFTs.$plan(x::AbstractArray{<:Complex{<:ForwardDiff.Dual}}, d::Integer, region=1:ndims(x); kwargs...) =
            AbstractFFTs.$plan(ForwardDiff.value.(x), d, region; kwargs...)
    end
end

for P in [:Plan, :ScaledPlan]  # need ScaledPlan to avoid ambiguities
    @eval begin
        Base.:*(p::AbstractFFTs.$P, x::AbstractArray{<:ForwardDiff.Dual}) =
            _apply_plan(p, x)

        Base.:*(p::AbstractFFTs.$P, x::AbstractArray{<:Complex{<:ForwardDiff.Dual}}) =
            _apply_plan(p, x)

        LinearAlgebra.mul!(Y::AbstractArray, p::AbstractFFTs.$P, X::AbstractArray{<:ForwardDiff.Dual}) =
            (Y .= _apply_plan(p, X))

        LinearAlgebra.mul!(Y::AbstractArray, p::AbstractFFTs.$P, X::AbstractArray{<:Complex{<:ForwardDiff.Dual}}) =
            (Y .= _apply_plan(p, X))
    end
end

LinearAlgebra.mul!(Y::AbstractArray{<:Complex{<:ForwardDiff.Dual}}, p::AbstractFFTs.ScaledPlan{T,P,<:ForwardDiff.Dual}, X::AbstractArray{<:ComplexF64}) where {T,P} =
    (Y .= _apply_plan(p, X))

function _apply_plan(p::AbstractFFTs.Plan, x::AbstractArray{<:Complex{<:ForwardDiff.Dual{T}}}) where T
    # TODO do we want x::AbstractArray{<:ForwardDiff.Dual{T}} too?
    xtil = p * ForwardDiff.value.(x)
    dxtils = ntuple(ForwardDiff.npartials(eltype(x))) do n
        p * ForwardDiff.partials.(x, n)
    end
    map(xtil, dxtils...) do val, parts...
        Complex(
            ForwardDiff.Dual{T}(real(val), map(real, parts)),
            ForwardDiff.Dual{T}(imag(val), map(imag, parts)),
        )
    end
end

function _apply_plan(p::AbstractFFTs.ScaledPlan{T,P,<:ForwardDiff.Dual}, x::AbstractArray) where {T,P}
    _apply_plan(p.p, p.scale * x) # for when p.scale is Dual, need out-of-place
end

# this is to avoid method ambiguities between these two:
#   _apply_plan(p::AbstractFFTs.Plan, x::AbstractArray{<:Complex{<:ForwardDiff.Dual{T}}}) where T
#   _apply_plan(p::AbstractFFTs.ScaledPlan{T,P,<:ForwardDiff.Dual}, x::AbstractArray) where {T,P}
function _apply_plan(p::AbstractFFTs.ScaledPlan{T,P,<:ForwardDiff.Dual}, x::AbstractArray{<:Complex{<:ForwardDiff.Dual{Tg}}}) where {T,P,Tg}
    _apply_plan(p.p, p.scale * x)
end

# DFTK setup specific

next_working_fft_size(::Type{<:ForwardDiff.Dual}, size::Int) = size

_fftw_flags(::Type{<:ForwardDiff.Dual}) = FFTW.MEASURE | FFTW.UNALIGNED

function build_fft_plans(T::Type{<:Union{ForwardDiff.Dual,Complex{<:ForwardDiff.Dual}}}, fft_size)
    tmp = Array{complex(T)}(undef, fft_size...) # TODO think about other Array types
    opFFT  = FFTW.plan_fft(tmp, flags=_fftw_flags(T))
    opBFFT = FFTW.plan_bfft(tmp, flags=_fftw_flags(T))

    ipFFT  = DummyInplace{typeof(opFFT)}(opFFT)
    ipBFFT = DummyInplace{typeof(opBFFT)}(opBFFT)
    # backward by inverting and stripping off normalizations
    ipFFT, opFFT, ipBFFT, opBFFT
end

# PlaneWaveBasis{<:Dual} contains dual-scaled fft, which means that the result f_fourier
# must be able to hold complex dual numbers even if f_real is not dual
function r_to_G(basis::PlaneWaveBasis{T}, f_real::AbstractArray) where {T<:ForwardDiff.Dual}
    f_fourier = similar(f_real, complex(T))
    @assert length(size(f_real)) ∈ (3, 4)
    # this exploits trailing index convention
    for iσ = 1:size(f_real, 4)
        @views r_to_G!(f_fourier[:, :, :, iσ], basis, f_real[:, :, :, iσ])
    end
    f_fourier
end

# determine symmetry operations only from primal lattice values
function spglib_get_symmetry(lattice::Matrix{<:ForwardDiff.Dual}, atoms, magnetic_moments=[]; kwargs...)
    spglib_get_symmetry(ForwardDiff.value.(lattice), atoms, magnetic_moments; kwargs...)
end

function _is_well_conditioned(A::AbstractArray{<:ForwardDiff.Dual}; kwargs...)
    _is_well_conditioned(ForwardDiff.value.(A); kwargs...)
end

# other workarounds

# problem: ForwardDiff of norm of SVector gives NaN derivative at zero
# https://github.com/JuliaMolSim/DFTK.jl/issues/443#issuecomment-864930410
# solution: follow ChainRules custom frule for norm
# https://github.com/JuliaDiff/ChainRules.jl/blob/52a0eeadf8d19bff491f224517b7b064ce1ba378/src/rulesets/LinearAlgebra/norm.jl#L5
# TODO delete, once forward diff AD tools use ChainRules natively
function LinearAlgebra.norm(x::SVector{S,<:ForwardDiff.Dual{Tg,T,N}}) where {S,Tg,T,N}
    x_value = ForwardDiff.value.(x)
    y = norm(x_value)
    dy = ntuple(j->real(dot(x_value, ForwardDiff.partials.(x,j))) * pinv(y), N)
    ForwardDiff.Dual{Tg}(y, dy)
end

# problem: the derivative of 1/(1+exp(x)) = -exp(x) / (1+exp(x))^2.
# When x is too large, exp(x) = Inf and this is a NaN.
function Smearing.occupation(S::Smearing.FermiDirac, d::ForwardDiff.Dual{T}) where {T}
    x = ForwardDiff.value(d)
    if exp(x) > floatmax(typeof(x)) / 1e3
        ∂occ = -zero(x)
    else
        ∂occ = -exp(x) / (1 + exp(x))^2
    end
    ForwardDiff.Dual{T}(Smearing.occupation(S, x), ∂occ * ForwardDiff.partials(d))
end


text:
using Test
using MDCDL
using FFTW
using Random

@testset "Multiscale" begin
    include("testsetGenerator.jl")
    # include("randomInit.jl")

    ccsd1D = cnsoltValidConfigSet1D()
    ccsd2D = cnsoltValidConfigSet2D()
    ccsd3D = cnsoltValidConfigSet3D()

    ccsd = [ ccsd1D, ccsd2D, ccsd3D ]

    rcsd1D = rnsoltValidConfigSet1D()
    rcsd2D = rnsoltValidConfigSet2D()
    rcsd3D = rnsoltValidConfigSet3D()

    rcsd = [ rcsd1D, rcsd2D, rcsd3D ]

    Random.seed!(3923528829)

    @testset "ParallelFB" begin
        maxDims = 2
        for d in 1:maxDims, dt in [ Float64, Complex{Float64}]
            cfgs = vec([ (crdf.I, crord.I .- 1, nch, lv) for crdf in CartesianIndices(tuple(fill(4,d)...)), crord in CartesianIndices(tuple(fill(2+1,d)...)), nch in 2:10, lv in 1:3 ])

            subcfgs = randsubseq(cfgs, 30 / length(cfgs))

            for (df, ord, nch, lv) in subcfgs
                pfb = ParallelFB(dt, df, ord, nch)

                szFilter = df .* (ord .+ 1)
                afs = [ rand(dt, szFilter) for p in 1:sum(nch) ]
                sfs = [ rand(dt, szFilter) for p in 1:sum(nch) ]
                pfb.analysisFilters .= afs
                pfb.synthesisFilters .= sfs

                mspfb = Multiscale(pfb, lv)

                szx = rand(1:4) .* (df.^lv) .* (ord .+ 1)
                x = rand(dt, szx)

                y = analyze(mspfb, x)

                myfilter(A, h) = begin
                    ha = zero(A)
                    # ha[colon.(1,size(h))...] = h
                    ha[UnitRange.(1, size(h))...] = h
                    if dt <: Real
                        real(ifft(fft(A).*fft(ha)))
                    else
                        ifft(fft(A).*fft(ha))
                    end
                end
                offset = df .- 1

                # check analyzer
                subx = x
                for idx = 1:lv-1
                    sys = map(afs) do af
                        fx = myfilter(subx, af)
                        dwfx = downsample(fx, df, offset)
                        circshift(dwfx, (-1 .* fld.(ord,2)))
                    end
                    if idx < lv
                        @test all(y[idx] .≈ sys[2:end])
                        subx = sys[1]
                    else
                        @test all(y[idx] .≈ sys)
                    end
                end

                rx = synthesize(mspfb, y)

                # check synthesizer
                suby = y[lv]
                for idx = lv:-1:1
                    fy = map(suby, sfs) do  yp, fp
                        upyp = MDCDL.upsample(yp, df)
                        myfilter( upyp, fp )
                    end
                    subrxs = sum(fy)
                    subrx = circshift(subrxs, -1 .* df .* cld.(ord,2))
                    if idx > 1
                        suby = [subrx, y[idx-1]... ]
                    else
                        @test size(rx) == size(subrx)
                        @test subrx ≈ rx
                    end
                end
            end
        end
    end

    # @testset "CNSOLT" begin
    #     # output mode options for analyzer
    #     oms = [ :normal, :augumented ]
    #     for d in 1:length(ccsd), (df, ord, nch) in ccsd[d], lv in 1:3
    #         szx = (df.^lv) .* (ord .+ 1)
    #         nsolt = Cnsolt(df, ord, nch)
    #         rand!(nsolt)
    #         msnsolt = Multiscale(nsolt, lv)
    #
    #         x = rand(Complex{Float64}, szx...)
    #
    #         y = analyze(msnsolt, x)
    #         rx = synthesize(msnsolt, y)
    #
    #         @test rx ≈ x
    #
    #         foreach(oms) do om
    #             y = analyze(nsolt, x; shape = om)
    #             rx = synthesize(nsolt, y)
    #
    #             @test size(x) == size(rx)
    #             @test rx ≈ x
    #         end
    #     end
    # end
    # @testset "RNSOLT" begin
    #     # output mode options for analyzer
    #     oms = [ :normal, :augumented ]
    #     for d in 1:length(ccsd), (df, ord, nch) in rcsd[d], lv in 1:3
    #         szx = (df.^lv) .* (ord .+ 1)
    #         nsolt = Rnsolt(df, ord, nch)
    #         rand!(nsolt)
    #         msnsolt = Multiscale(nsolt, lv)
    #
    #         x = rand(Float64, szx...)
    #
    #         y = analyze(msnsolt, x)
    #         rx = synthesize(msnsolt, y)
    #
    #         @test rx ≈ x
    #
    #         foreach(oms) do om
    #             y = analyze(nsolt, x; shape = om)
    #             rx = synthesize(nsolt, y)
    #
    #             @test size(x) == size(rx)
    #             @test rx ≈ x
    #         end
    #     end
    # end
end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
include("types.jl")

"Ordered timestamps from blotter associative collection."
vtblotter(blotter::Blotter) = sort!(collect(keys(blotter)))

"""
Amount `Vector{Int64)` and price `Vector{Float64)` from blotter
in chronological order (returns vector tuple).
"""
function vapblotter(blotter::Blotter)
  vt = vtblotter(blotter)
  nt = length(vt)
  vqty = Array(Int64, nt)
  vprc = zeros(nt)
  for i = 1:nt
    ti = vt[i]
    bi = blotter[ti]
    @inbounds vqty[i] = bi[1]
    @inbounds vprc[i] = bi[2]
  end
  return vqty, vprc
end

"Print a text line from string vector."
function printvecstring(io, vstring::Vector,
                        separator::Char, quotemark::Char)
  p = length(vstring)
  for j in 1:p
    print(io, quotemark)
    print(io, vstring[j])
    print(io, quotemark)
    if j < p
      print(io, separator)
    else
      print(io, '\n')
    end
  end
  return
end

"Print blotter transactions. Resembles DataFrames.printtable."
function printblotter(io::IO, blotter::Blotter;
                      dtformat::AbstractString = "yyyy-mm-ddTHH:MM:SS",
                      separator::Char = ',', quotemark::Char = '"')
  # ordered timestamps
  vt = vtblotter(blotter)
  nt = length(vt)

  # column names: print header
  cnames = ["Timestamp", "Amount", "Price"]
  printvecstring(io, cnames, separator, quotemark)

  # print contents
  for ti in vt
    # timestamp
    print(io, quotemark)
    print(io, Dates.format(ti, dtformat))
    print(io, quotemark)
    print(io, separator)
    # amount and price
    bi = blotter[ti]
    print(io, bi[1])
    print(io, separator)
    print(io, bi[2])
    print(io, '\n')
  end

  return
end

"Write blotter transactions to file."
function writeblotter(filename::AbstractString, blotter::Blotter;
                      dtformat::AbstractString = "yyy-mm-ddTHH:MM:SS",
                      separator::Char = ',', quotemark::Char = '"')
  open(filename, "w") do io
    printblotter(io, blotter,
                 dtformat=dtformat, separator=separator, quotemark=quotemark)
  end
  return
end

"""
Trade analysis for `blotter` provided as
`DateTime => (Qty::Int64, FillPrice::Float64)` assoc. collection.
Input `metrics` specifies what to calculate (PnL included already - others).
Returns: tuple ( DateTime (ordered) array , assoc. collection of perf metrics ).
Basic transaction info is also included (quantity, fill price).

**CAUTION**: PnL and drawdown are calculated here based on the transaction blotter
only, not the price history. Hence, price swing effects while holding
an open position are not showing up in the results. Use `orderhandling!`
output if performance metrics over the whole price history are needed
(as typically done when analyzing PnL and drawdown).
"""
function tradeperf(blotter::Blotter,
                   metrics::Vector{Symbol} = [:DDown])
  perfm = Dict{Symbol,Vector{Float64}}()
  ### TODO (later): accociative collections syntax changes in Julia 0.4

  # timestamps in order
  vt = vtblotter(blotter)
  nt = length(vt)

  # basic info for each transaction
  vqty, vprc = vapblotter(blotter)
  perfm[:Qty] = vqty
  perfm[:FillPrice] = vprc

  # PnL and return-based drawdown (always calculating those metrics)
  vpl = zeros(nt); vdd = zeros(nt)
  acsum = 0; pnlcsum = 0.0; pnlcmax = 0.0
  for i = 2:nt
    bprev = blotter[vt[i-1]]

    # lagged cumulative position held up to the current timestamp
    @inbounds acsum += bprev[1]

    # cumulative profit/loss at the current timestamp
    pnlcsum += acsum * (blotter[vt[i]][2] - bprev[2])
    @inbounds vpl[i] = pnlcsum

    # return-based drawdown
    if pnlcsum > pnlcmax # positive or negative
      pnlcmax = pnlcsum
    end
    @inbounds vdd[i] = pnlcsum - pnlcmax
  end
  perfm[:PnL] = vpl

  # return-based drawdown
  if findfirst(metrics, :DDown) > 0
    perfm[:DDown] = vdd
  end

  return vt, perfm
end

"Cumulative position, profit/loss, last fill price for blotter."
function apnlcum(blotter::Blotter)
  # timestamps in order
  vt = vtblotter(blotter)
  nt = length(vt)
  if nt == 0
    # no transactions yet
    return 0, 0.0, 0.0
  end

  acsum = 0; pnlcsum = 0.0
  for i = 2:nt
    bprev = blotter[vt[i-1]]

    # lagged cumulative position held up to the current timestamp
    @inbounds acsum += bprev[1]

    # cumulative profit/loss at the current timestamp
    pnlcsum += acsum * (blotter[vt[i]][2] - bprev[2])
  end

  blast = blotter[vt[nt]]
  return acsum + blast[1], pnlcsum, blast[2]
end

"""
Final profit/loss for `blotter` provided as
`DateTime => (Qty::Int64, FillPrice::Float64)` assoc. collection.
Faster verision (minimizing memory allocation) to be used
in e.g. parameter optimization workflow.

Returns: final profit/loss `Float64` scalar.
"""
:tradepnlfinal

"Based on blotter only, ending at the last transaction timestamp."
tradepnlfinal(blotter::Blotter) = apnlcum(blotter)[2]

"Adding current price as the last timestamp."
function tradepnlfinal(blotter::Blotter, pnow::Float64)
  # up to pcur
  acsumb, pnlb, pblast = apnlcum(blotter)

  # with pcur
  return pnlb + acsumb * (pnow - pblast)
end

"Performance metrics helper function for use in foldp."
function tradeperffold(perfprev::@compat(Tuple{Float64,Float64}), statusnow::@compat(Tuple{Bool, Float64}))
  # current max. PnL
  pnlprev = perfprev[1]
  pnlnow = statusnow[2]
  pnlcummax = pnlnow > pnlprev ? pnlnow : pnlprev

  # current return-based drawdown
  ddnow = abs(pnlnow - pnlcummax)

  # maximum drawdown
  ddprev = perfprev[2]
  ddmax = ddnow > ddprev ? ddnow : ddprev

  return pnlcummax, ddmax
end

"""
Selected performance metrics from `runtrading!` signal output.

Output tuple-signal components:

* `Float64` cumulative maximum PnL;
* `Float64` maximum drawdown over the entire trading session hisotry.

NOTE: Use this function only if needed, otherwise save resources; it is
not required for running the trading session.
"""
tradeperfcurr(s_status::Signal{@compat(Tuple{Bool, Float64})}) = foldp(tradeperffold, (0.0, 0.0), s_status)

"""
Selected metrics for completed trades out of transactions blotter.

Return tuple contains:

- `Vector{Float64}` profit/loss for each completed trade;
- `Int64` number of winning trades;
- `Float64` average winning trade profit;
- `Int64` number of loosing trades;
- `Float64` average loosing trade loss.
"""
function vtradespnl(blotter::Blotter)
  # timestamps in order
  vt = vtblotter(blotter)
  nt = length(vt)
  if nt == 0
    # no transactions
    return Array(Float64, 0), 0, 0.0, 0, 0.0
  end

  vtrpnl = Array(Float64, 0)
  acsum = 0; pnlcsum = 0.0
  pnlcsum_trprevend = 0.0
  ntrpos = 0; ntrneg = 0
  profcum_trpos = 0.0; losscum_trneg = 0.0
  for i = 2:nt
    bprev = blotter[vt[i-1]]; bnow = blotter[vt[i]]
    # lagged cumulative position held up to the current timestamp
    @inbounds acsum += bprev[1]
    # cumulative profit/loss at the current timestamp
    pnlcsum += acsum * (bnow[2] - bprev[2])

    # new position based on the current transaction amount
    posi = acsum + bnow[1]
    if posi == 0 || sign(posi) == -sign(acsum)
      # trade completed at the current timestep
      pnli = pnlcsum - pnlcsum_trprevend
      pnlcsum_trprevend = pnlcsum
      push!(vtrpnl, pnli)
      if pnli > 0.0
        # winning trade
        ntrpos += 1
        profcum_trpos += pnli
      else
        # loosing trade
        ntrneg += 1
        losscum_trneg += abs(pnli)
      end
    end
  end

  if ntrpos < 1
    avwin = 0.0
  else
    avwin = profcum_trpos / ntrpos
  end
  if ntrneg < 1
    avloss = 0.0
  else
    avloss = losscum_trneg / ntrneg
  end

  return vtrpnl, ntrpos, avwin, ntrneg, avloss
end

"""
Pessimistic return on margin `marg`.

For `pror = true` returns pessimistic rate of return (can be thought of as
a more realistic profit factor).

Optional `nbest_remove` argument specifies how many best wins
are dropped (to increase PROM strictness).

**CAUTION** `-Inf` returned if `nbest_remove` exceeds the
total number of trades.
"""
:perf_prom

"Pessimistic rate of return with extreme case handling."
function perf_pror_auxil(ppos::Float64, pneg::Float64)
  if abs(pneg) < eps() && abs(ppos) < eps()
    # 0/0 situation: no profit anyway
    return 0.0
  end
  # with finite ppos +/-Inf is acceptable if pneg=0
  return ppos/pneg
end

"From blotter using completed trades."
function perf_prom(blotter::Blotter;
                   pror::Bool = false,
                   marg::Float64 = 1.0,
                   nbest_remove::Int64 = 0)
  marg = abs(marg) # positive denominator
  vtrpnl, ntrpos, avwin, ntrneg, avloss = vtradespnl(blotter)

  if nbest_remove == 0
    ppos = avwin*(ntrpos - sqrt(float(ntrpos)))
    pneg = avloss*(ntrneg + sqrt(float(ntrneg)))
    if pror
      return perf_pror_auxil(ppos, pneg)
    else
      return (ppos - pneg) / marg
    end
  end
  return perf_prom(vtrpnl, marg = marg, pror = pror, nbest_remove = nbest_remove)
end

"From profit/loss vector of completed trades."
function perf_prom(vtrpnl::Vector{Float64};
                   marg::Float64 = 1.0,
                   pror::Bool = false,
                   nbest_remove::Int64 = 0)
  ntr = length(vtrpnl)
  marg = abs(marg) # positive denominator
  if ntr < 1
    # no trades
    return 0.0
  end

  if nbest_remove > 0
    if ntr <= nbest_remove
      return -Inf
    end
    pnllim = sort(vtrpnl, rev=true)[nbest_remove]
  else
    # no trades are dropped: not limiting PnL
    pnllim = Inf
  end

  ntrpos = 0; ntrneg = 0
  profcum_trpos = 0.0; losscum_trneg = 0.0
  for pnli in vtrpnl
    # take into account PnL limit after best trades removal
    if 0.0 < pnli < pnllim
      ntrpos += 1
      profcum_trpos += pnli
    end
    if pnli < 0.0
      ntrneg += 1
      losscum_trneg += abs(pnli)
    end
  end

  if ntrpos < 1
    avwin = 0.0
  else
    avwin = profcum_trpos / ntrpos
  end
  if ntrneg < 1
    avloss = 0.0
  else
    avloss = losscum_trneg / ntrneg
  end

  ppos = avwin*(ntrpos - sqrt(ntrpos))
  pneg = avloss*(ntrneg + sqrt(ntrneg))
  if pror
    return perf_pror_auxil(ppos, pneg)
  end
  return (ppos - pneg) / marg
end


text:
using Dictionaries
using Tables
using Test
using TypedTables
using SplitApplyCombine
using Adapt
using Tables

# Here because Julia v1.0 doesn't allow type definitions inside of testsets:
struct TestArrayConverter end
Adapt.adapt_storage(::TestArrayConverter, xs::AbstractArray) = convert(Array, xs)

include("properties.jl")
include("Table.jl")
include("FlexTable.jl")
include("DictTable.jl")


text:
using MAT
using Random
using StatsBase
using SparseArrays
include("../include/HyperLocal/HyperLocal.jl")
include("../include/HyperLocal/SparseCardHyperLocal.jl")
include("../src/hypergraph-clustering-utils.jl")

## Load in the data
include("load_stackoverflow.jl")

## Generate seed sets
# NOTE: original seed sets are saved in Seeds folder

# include("generate_stackoverflow_seeds.jl")

## Run Hyperlocal+SparseCard with alternative hyperedge cut penalties
types = ["clique", "x9", "sqrt"]
lnum = length(labels)

for trial = 1:10
    SS = matread("Seeds/Seed_sets_stack45_$trial.mat")
    seed_sets = SS["seed_sets"]

    # Run for the three different alternate penalties
    for tt = 1:3
        epsilon = 1.0
        sparsityeps1 = 1.0
        sparsityeps2 = 0.1
        sparsityeps3 = 0.01
        type = types[tt]
        output = "Output/SCHL_$(type)_45_$trial.mat"

        if tt == 1
            EdgesW = clique_expansion_EdgesW(order)
        elseif tt == 2
            fun = x->x.^(.9)
            EdgesW = generic_EdgesW(order,fun)
        else
            funsqrt = x->sqrt.(x)
            EdgesW = generic_EdgesW(order,funsqrt)
        end

        sparse1_stats = zeros(lnum,6)
        sparse2_stats = zeros(lnum,6)
        sparse3_stats = zeros(lnum,6)

        for lab = 1:length(LabelNames)
            T = findnz(LabelMatrix[:,lab])[1]
            nT = length(T)
            println("$lab \t $nT \t $(LabelNames[label])")

            # Generate a new seed set
            seednum = round(Int64,nT/20)
            grownum = round(Int64,min(nT*2,n/2-seednum))
            Rstart = findall(x->x>0,seed_sets[:,lab])
            OneHop = get_immediate_neighbors(H,Ht,Rstart)
            Rmore = BestNeighbors(H,d,Rstart,OneHop,grownum)
            R = union(Rmore,Rstart)
            Rs = findall(x->in(x,Rstart),R)     # Force seed nodes to be in output set

            # Run HyperLocal + SparseCard with 3 different sparsity epsilons
            sparsities = [1.0, 0.1, 0.01]
            for sss = 1:3
                sparsityeps = sparsities[sss]
                s = time()
                S, lcond = SparseCardHyperLocal(H,Ht,EdgesW,Edges,sparsityeps,order,d,R,epsilon,Rs)
                timer = time()-s
                pr, re, f1 = PRF(T,S)
                nS =  length(setdiff(S,R))
                if sss == 1
                    sparse1_stats[lab,:] = [pr, re, f1, length(S), nS, timer]
                elseif sss == 2
                    sparse2_stats[lab,:] = [pr, re, f1, length(S), nS, timer]
                else
                    sparse3_stats[lab,:] = [pr, re, f1, length(S), nS, timer]
                end
            end
            # matwrite(output,Dict("sparse1_stats"=>sparse1_stats,"sparse2_stats"=>sparse2_stats,"sparse3_stats"=>sparse3_stats))
        end
    end

    # Run Hyperlocal with standard delta linear threshold penalty
    delta = 5000.0
    epsilon = 1.0
    delta_stats = zeros(lnum,6)
    output = "Output/deltalinear_45_$trial.mat"
    for lab = 1:length(labels)
        label = labels[lab]
        T = findnz(LabelMatrix[:,label])[1]
        nT = length(T)
        println("$lab \t $nT \t $label")

        # Generate a new seed set
        seednum = round(Int64,nT/20)
        grownum = round(Int64,min(nT*2,n/2-seednum))
        Rstart = findall(x->x>0,seed_sets[:,lab])
        OneHop = get_immediate_neighbors(H,Ht,Rstart)
        Rmore = BestNeighbors(H,d,Rstart,OneHop,grownum)
        R = union(Rmore,Rstart)
        Rs = findall(x->in(x,Rstart),R)     # Force seed nodes to be in output set

        # Run HyperLocal with delta = 5000.0
        s = time()
        S, lcond = HyperLocal(H,Ht,order,d,R,epsilon,delta,Rs,true)
        timer = time()-s
        pr, re, f1 = PRF(T,S)
        nS =  length(setdiff(S,R))
        delta_stats[lab,:] = [pr, re, f1, length(S), nS, timer]
        # matwrite(output,Dict("delta_stats"=>delta_stats))
    end

end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
struct Ybus{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF64}
    data::SparseArrays.SparseMatrixCSC{ComplexF64, Int64}
    axes::Ax
    lookup::L
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int64},
    b::Line,
    num_bus::Dict{Int64, Int64},
)

    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_l = (1 / (get_r(b) + get_x(b) * 1im))
    Y11 = Y_l + (1im * get_b(b).from)

    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = -Y_l
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_l + (1im * get_b(b).to)
    ybus[bus_to_no, bus_to_no] += Y22

    return

end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int64},
    b::Transformer2W,
    num_bus::Dict{Int64, Int64},
)

    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    Y11 = Y_t
    b = get_primaryshunt(b)

    ybus[bus_from_no, bus_from_no] += Y11
    ybus[bus_from_no, bus_to_no] += -Y_t
    ybus[bus_to_no, bus_from_no] += -Y_t
    ybus[bus_to_no, bus_to_no] += Y_t + (1im * b)

    return

end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int64},
    b::TapTransformer,
    num_bus::Dict{Int64, Int64},
)

    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    c = 1 / get_tap(b)
    b = get_primaryshunt(b)

    Y11 = (Y_t * c^2)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t * c)
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)

    return

end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int64},
    b::PhaseShiftingTransformer,
    num_bus::Dict{Int64, Int64},
)

    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    tap = (get_tap(b) * exp(get_α(b) * 1im))
    c_tap = (get_tap(b) * exp(-1 * get_α(b) * 1im))
    b = get_primaryshunt(b)

    Y11 = (Y_t / abs(tap)^2)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t / c_tap)
    ybus[bus_from_no, bus_to_no] += Y12
    Y21 = (-Y_t / tap)
    ybus[bus_to_no, bus_from_no] += Y21
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)

    return

end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int64},
    fa::FixedAdmittance,
    num_bus::Dict{Int64, Int64},
)

    bus = get_bus(fa)
    bus_no = num_bus[get_number(bus)]

    ybus[bus_no, bus_no] += fa.Y

    return

end

function _buildybus(branches, nodes, fixed_admittances)

    buscount = length(nodes)
    num_bus = Dict{Int64, Int64}()

    for (ix, b) in enumerate(nodes)
        num_bus[get_number(b)] = ix
    end

    ybus = SparseArrays.spzeros(ComplexF64, buscount, buscount)

    for (ix, b) in enumerate(branches)
        if get_name(b) == "init"
            throw(DataFormatError("The data in Branch is invalid"))
        end

        _ybus!(ybus, b, num_bus)

    end

    for fa in fixed_admittances
        _ybus!(ybus, fa, num_bus)
    end

    return ybus

end

function Ybus(branches, nodes)

    #Get axis names
    bus_ax = [get_number(bus) for bus in nodes]
    axes = (bus_ax, bus_ax)
    look_up = (_make_ax_ref(bus_ax), _make_ax_ref(bus_ax))

    ybus = _buildybus(branches, nodes, Vector{FixedAdmittance}())

    return Ybus(ybus, axes, look_up)

end

function Ybus(sys::System)
    branches = get_components(ACBranch, sys)
    nodes = sort(collect(get_components(Bus, sys)), by = x -> x.number)
    fixed_admittances = get_components(FixedAdmittance, sys)

    # Get axis names
    bus_ax = [get_number(bus) for bus in nodes]
    axes = (bus_ax, bus_ax)
    look_up = (_make_ax_ref(bus_ax), _make_ax_ref(bus_ax))

    ybus = _buildybus(branches, nodes, fixed_admittances)

    return Ybus(ybus, axes, look_up)
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using PyPlot
using StableRNGs

MLJ.color_off() # hide
LGBMRegressor = @load LGBMRegressor

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module MLJSerialization

# export IterationControl controls:
export Save

include("machines.jl")
include("controls.jl")

end # module


text:
include("utils.jl")

@testset "Curve warping action" begin
    N = 20
    r2 = Euclidean(2)
    knots = range(0.0, 1.0, length=N)
    dcsr2 = DiscretizedCurves(r2, knots)
    M = CurveWarpingSpace(knots)
    G = CurveWarpingGroup(M)
    A_left = CurveWarpingAction(dcsr2, G)

    #@test repr(A_left) == "CurveWarpingAction($(repr(dcsr2)), $(repr(G)))"

    @test g_manifold(A_left) == dcsr2
    @test base_group(A_left) == G
    @test isa(A_left, AbstractGroupAction{LeftAction})
    @test base_manifold(G) == M

    a_funs = [x -> x^2, x -> x^3, x -> sqrt(x)]
    a_pts = [FunManifolds.make_warping(M, map(f, knots)) for f in a_funs]
    m_pts = [
        discretize(dcsr2, t -> project(r2, @SVector [t^2, sin(t)])),
        discretize(dcsr2, t -> project(r2, @SVector [sin(t^2), 1])),
        discretize(dcsr2, t -> project(r2, @SVector [t, cos(t)])),
    ]

    test_action(
        A_left,
        a_pts,
        m_pts;
        atol=2 / N,
        test_optimal_alignment=false,
        test_diff=false,
        test_mutating_group=false,
        test_mutating_action=false,
        test_switch_direction=false,
    )
end

@testset "Curve SRSF warping action" begin
    N = 100
    r2 = Euclidean(2)
    knots = range(0.0, 1.0, length=N)
    dcsr2 = DiscretizedCurves(r2, knots)
    M = CurveWarpingSRSFSpace(knots)
    M2 = CurveWarpingSpace(knots)
    G = CurveWarpingSRSFGroup(M)
    p = [0.0, 0.0]
    dcsr2t = DiscretizedCurves(TangentSpaceAtPoint(r2, p), knots)
    A_left = CurveWarpingSRSFAction(dcsr2t, p, G)

    #@test repr(A_left) == "CurveWarpingSRSFAction($(repr(dcsr2)), $(repr(G)))"

    @test g_manifold(A_left) == dcsr2t
    @test base_group(A_left) == G
    @test isa(A_left, AbstractGroupAction{LeftAction})
    @test base_manifold(G) == M

    a_funs = [x -> x^2, x -> x^3, x -> sqrt(x)]
    a_pts = [srsf(M2, FunManifolds.make_warping(M2, map(f, knots))) for f in a_funs]
    m_pts = [
        tsrvf(
            dcsr2,
            discretize(dcsr2, t -> project(r2, @SVector [t^2, sin(t)])),
            p,
            FunManifolds.ProjectedDifferenceBackend(nothing),
        ),
        tsrvf(
            dcsr2,
            discretize(dcsr2, t -> project(r2, @SVector [sin(t^2), 1])),
            p,
            FunManifolds.ProjectedDifferenceBackend(nothing),
        ),
        tsrvf(
            dcsr2,
            discretize(dcsr2, t -> project(r2, @SVector [t, cos(t)])),
            p,
            FunManifolds.ProjectedDifferenceBackend(nothing),
        ),
    ]

    test_action(
        A_left,
        a_pts,
        m_pts;
        atol=3 / sqrt(N),
        atol_ident_compose=2 * sqrt(eps(Float64)),
        test_optimal_alignment=true,
        test_diff=false,
        test_mutating_group=false,
        test_mutating_action=false,
        test_switch_direction=false,
    )
end

@testset "Mean amplitude" begin
    grid_small = 0.0:0.125:1.0
    r1 = Euclidean(1)
    r2 = Euclidean(2)
    f1 = t -> begin
        [if t < 0.5
            0.0
        elseif t < 0.75
            t - 0.5
        else
            1.0 - t
        end]
    end
    f2 = t -> begin
        [if t < 0.25
            0.0
        elseif t < 0.5
            t - 0.25
        elseif t < 0.75
            0.75 - t
        else
            0.0
        end]
    end
    f3 = t -> begin
        [if t < 0.25
            0.0
        elseif t < 0.5
            -20 * (t - 0.25)
        elseif t < 0.75
            -20 * (0.75 - t)
        else
            0.0
        end]
    end
    dcsr1 = DiscretizedCurves(r1, grid_small)
    M = CurveWarpingSRSFSpace(grid_small)
    M2 = CurveWarpingSpace(grid_small)
    G = CurveWarpingSRSFGroup(M)
    p = [0.0]
    dcsr1t = DiscretizedCurves(TangentSpaceAtPoint(r1, p), grid_small)
    A_left = CurveWarpingSRSFAction(dcsr1t, p, G)

    f1qc = tsrvf(
        dcsr1,
        discretize(dcsr1, f1),
        p,
        FunManifolds.ProjectedDifferenceBackend(nothing),
    )
    f2qc = tsrvf(
        dcsr1,
        discretize(dcsr1, f2),
        p,
        FunManifolds.ProjectedDifferenceBackend(nothing),
    )
    f3qc = tsrvf(
        dcsr1,
        discretize(dcsr1, f3),
        p,
        FunManifolds.ProjectedDifferenceBackend(nothing),
    )
    μq = FunManifolds.karcher_mean_amplitude(A_left, [f1qc, f2qc])
    for pt in μq
        @test abs(pt) < 2.0
    end
    (p̃, γs, aligned_ps) = FunManifolds.phase_amplitude_separation(A_left, [f1qc, f2qc])
    @test is_point(dcsr1t, p̃)
    for γ in γs
        @test is_point(M, γ)
    end
    for p in aligned_ps
        @test is_point(dcsr1t, p)
    end
end


text:
#
println("Ed) Photoabsorption of Ne^+: AverageSCA model.")
using JLD
#
if  false
    # Calculation of the Ne^+  2s^2 2p^5; 2s --> 2p, 3p, 4p, 5p excitation cross sections
    setDefaults("print summary: open", "Ne-plus-photoexcitation.sum")
    setDefaults("method: continuum, asymptotic Coulomb")                                                      ## "method: continuum, Galerkin"
    grid = Radial.Grid(Radial.Grid(false); rnt = 3.0e-6, h = 2.0e-2, hp = 3.0e-2, NoPoints=1110)
    name = "Photoabsorption calculations for Ne^+ for energies = [1, 4] a.u."
    wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(10.), grid=grid, approach=Cascade.AverageSCA(),
                               scheme=Cascade.PhotonExcitationScheme([PhotoExc()], [E1], 0.5, 4.0, 1, [Shell("2s"), Shell("2p")],
                                                                     [Shell("2s"), Shell("2p"), Shell("3p"), Shell("4p"), Shell("5p")]),
                               initialConfigs=[Configuration("1s^2 2s^2 2p^5")] )
    println(wa)
    wb = perform(wa; output=true)
    setDefaults("print summary: close", "")
    #
elseif false
    # Calculation of the Ne^+  2s^2 2p^5; photo-ionization cross sections
    setDefaults("print summary: open", "Ne-plus-photoionization.sum")

    name = "Photoionization of Si- "
    grid = Radial.Grid(Radial.Grid(false); rnt = 3.0e-6, h = 2.0e-2, hp = 3.0e-2, NoPoints=1110)
    wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(10.), grid=grid, approach=Cascade.AverageSCA(),
                               scheme=Cascade.PhotonIonizationScheme([Photo()], 1, [1.0, 3.0, 5.0]),
                               initialConfigs=[Configuration("1s^2 2s^2 2p^5")] )
    println(wa)
    wb = perform(wa; output=true)
    setDefaults("print summary: close", "")

elseif false
    # Simulation of the Ne^+  2s^2 2p^5; photo-ionization cross sections
    setDefaults("print summary: open", "Ne-plus-photoabsorption-simulation.sum")

    data = [JLD.load("zzz-cascade-ionizing-computations-2020-07-10T19.jld"), JLD.load("zzz-cascade-excitation-computations-2020-07-10T19.jld")]
    ##x JLD2.@load "zzz-cascade-ionizing-computations-2020-06-23T09.jld"
    ##x resIon  = results
    ##x JLD2.@load "zzz-cascade-excitation-computations-2020-06-23T09.jld"
    ##x resExc = results
    ##x data = [resExc]

    name = "Simulation of photoabsorption for Ne^+"

    wc   = Cascade.Simulation(Cascade.Simulation(), name=name, property=Cascade.PhotoAbsorption(),
                              settings=Cascade.SimulationSettings(0., 0., 1., 4., 3., [(1, 1.0)]), computationData=data )
    println(wc)
    wd = perform(wc; output=true)
    setDefaults("print summary: close", "")
#
elseif false
    # Example for Atomic cascade computation (symmetry, 2021)
    # Calculation of the Mg  1s 2s^2 2p^6 3s^2 decay cascade after 1s photoionization
    setDefaults("print summary: open", "Mg-1s-decay.sum")

    decayScheme = Cascade.StepwiseDecayScheme([Auger(), Radiative()], 3, Dict{Int64,Float64}(), 0., Shell[], Shell[])

    grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)
    name = "Computation of Mg  1s 2s^2 2p^6 3s^2 decay cascade after 1s photoionization"
    wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(12.), grid=grid,
                               approach=Cascade.AverageSCA(), scheme=decayScheme,
                               ## initialConfigs=[Configuration("1s 2s^2 2p^6 3s^2")] )
                               initialConfigs=[Configuration("1s 2s^2 2p^6 3s^2 3p")] )
    @show wa
    wb = perform(wa; output=true)
    setDefaults("print summary: close", "")
#
elseif true
    # Simulation of the Mg  1s 2s^2 2p^6 3s^2 decay cascade after 1s photoionization
    setDefaults("print summary: open", "zzz-Cascade-simulation.sum")

    ## simulationSettings = Cascade.SimulationSettings(0., 0., 0., 0., 0., [(1, 1.0)])
    ## simulationSettings = Cascade.SimulationSettings(0., 0., 0., 0., 0., [(4, 1.0)])
    simulationSettings = Cascade.SimulationSettings(0., 0., 0., 0., 0., [(1, 0.25), (2, 0.25), (3, 0.25), (4, 0.25)])

    data = [JLD.load("zzz-cascade-decay-computations-2021-03-14T19.jld")]
    name = "Simulation of Mg  1s 2s^2 2p^6 3s^2 decay cascade after 1s photoionization"

    wc   = Cascade.Simulation(Cascade.Simulation(), name=name, property=Cascade.IonDistribution(),
                              settings=simulationSettings, computationData=data )
    @show wc
    wd = perform(wc; output=true)
    setDefaults("print summary: close", "")
    #

end


text:
using Pkg
using GitHub
using DataFrames
using Downloads

function pkgversion()
    projecttoml_filename = joinpath(dirname(dirname(@__FILE__)),"Project.toml")
    projecttoml_parsed = Pkg.TOML.parse(read(projecttoml_filename, String))
    return VersionNumber(projecttoml_parsed["version"])
end

function readprinteddataframe(str::String)
    lines = split(str, "\n")
    if occursin("DataFrame",str)
        start = findfirst(occursin.("DataFrame", lines))
    else
        start = 1 #try first line
    end
    df = DataFrame(CSV.File(IOBuffer(str), delim="│", header=start, datarow=start+4))
    select!(df, Not([1,2,length(names(df))]))
    names!(df,Symbol.(strip.(names(df))))
    df[!,:cat] = string.(strip.(df[:cat]))
    df[!,:testname] = string.(strip.(df[:testname]))
    df[!,:res] = string.(strip.(df[:res]))
    return df
end

function getsubmittedbenchmarks(;repo::String="IanButterworth/SystemBenchmark.jl", issue::Int=8,
                    refname::String="ref.txt", transpose::Bool=true, authkey=get(ENV,"PERSONAL_ACCESS_TOKEN",nothing))

    if !isnothing(authkey)
        auth = GitHub.authenticate(authkey)
        comments = GitHub.comments(repo, issue, :issue, auth=auth)[1]
    else
        comments = GitHub.comments(repo, issue, :issue)[1]
    end
    i = 1
    ref = readbenchmark(joinpath(dirname(@__DIR__), "ref", refname))
    master_res = DataFrame(cat=["info","info"],testname=["user","datetime"],units=Union{String,Missing}[missing,missing],res=["ref","2020-05-10"])
    append!(master_res, ref)
    rename!(master_res, [:cat,:testname,:units,:ref])
    files = map(comment->map(x->string(x.match), collect(eachmatch(r"https:\/\/github\.com\/.*\.txt",comment.body, overlap=false))),comments)
    nresults = sum(length.(files))
    filter!(x->occursin(r"https:\/\/github\.com\/.*\.txt", x.body), comments)
    prog = Progress(nresults, desc = "Downloading $(nresults) results... ")
    filedict = Dict{String,String}()
    @sync for comment in comments
        resulturls = map(x->string(x.match), collect(eachmatch(r"https:\/\/github\.com\/.*\.txt",comment.body, overlap=false)))
        for resulturl in resulturls
            @async begin
                filedict[resulturl] = Downloads.download(resulturl)
                next!(prog)
            end
        end
    end
    for comment in comments
        resulturls = map(x->string(x.match), collect(eachmatch(r"https:\/\/github\.com\/.*\.txt",comment.body, overlap=false)))
        for resulturl in resulturls
            username = "@$(comment.user.login)"
            datetime = comment.updated_at
            file = filedict[resulturl]
            res = readbenchmark(file)
            if ("units" in names(res))
                "test_res" in names(res) && (res = DataFrame(cat=res.cat, testname=res.testname, units=res.units, res=res.test_res))
                res_formatted = DataFrame(cat=["info","info"],testname=["user","datetime"],units=Union{String,Missing}[missing,missing],res=[username,datetime])
                append!(res_formatted, res)
                rename!(res_formatted, [:cat,:testname,:units,Symbol("res_$i")])
                master_res = DataFrames.outerjoin(master_res, res_formatted, on = [:cat,:units,:testname], matchmissing=:equal)
            else
                "test_res" in names(res) && (res = DataFrame(cat=res.cat, testname=res.testname, res=res.test_res))
                res_formatted = DataFrame(cat=["info","info"],testname=["user","datetime"],res=[username,datetime])
                append!(res_formatted, res)
                rename!(res_formatted, [:cat,:testname,Symbol("res_$i")])
                master_res = DataFrames.outerjoin(master_res, res_formatted, on = [:cat,:testname], matchmissing=:equal)
            end
            i += 1
            next!(prog)
        end
    end
    finish!(prog)
    if transpose
        return restranspose(master_res)
    else
        return master_res
    end
end
function restranspose(res::DataFrame)
    resfilt = res[:,4:end]
    rows =  collect.(eachrow(resfilt))
    i = 1
    for c in res.cat
        if c != "info"
            rows[i] = map(x->ismissing(x) ? missing : parse(Float64,x),rows[i])
        end
        i += 1
    end
    df = DataFrame([[names(resfilt)]; rows], [:resultid; Symbol.(res.testname)])
    return df
end


text:
using MacroTools: @q
import PyCall
using MLStyle

import SymPy
using SymPy: Sym, sympy, symbols, free_symbols
import SymPy.sympy

const symfuncs = Dict()

_pow(a,b) = Base.:^(float(a),b)

function __init__()
    stats = PyCall.pyimport_conda("sympy.stats", "sympy")
    SymPy.import_from(stats)
    global stats = stats

    # for dist in [:Normal, :Cauchy, :Laplace, :Beta, :Uniform]
    #     @eval begin
    #         function Distributions.$dist(μ::Sym, σ::Sym)
    #             println("Evaluating ",$dist, "(",μ," :: Sym, ", σ, " :: Sym)")
    #             stats.$dist(:dist, μ,σ) |> SymPy.density
    #         end

    #         Distributions.$dist(μ,σ) = $dist(promote(μ,σ)...)
    #     end
    # end

    # https://discourse.julialang.org/t/pyobjects-as-keys/26521/2
    merge!(symfuncs, Dict(
        sympy.log => Base.log
      , sympy.Pow => _pow
      , sympy.Abs => Base.abs
      , sympy.Indexed => getindex
    ))

    @eval begin
    @gg function codegen(_m::Model, _args, _data)
        f = _codegen(type2model(_m))
        :($f(_args, _data))
    end
end

end

struct SymDist
    logpdf :: Function
end

logpdf(d::SymDist, x) = d.logpdf(sym(x))

Distributions.Bernoulli(p::Sym) = SymDist(y -> y * log(p) + (1-y) * log(1-p))

for dist in [:Bernoulli]
    @eval begin
        logpdf(d::$dist, x::Sym) = logpdf($dist(sym.(Distributions.params(d))...), x)

    end
end

"Half" distributions
for dist in [:Normal, :Cauchy]
    let half = Symbol(:Half, dist)
        @eval begin
            logpdf(d::$half, x::Sym) = 2 * logpdf($dist(0, sym.(d.σ)), x)
        end
    end
end

for dist in [:Normal, :Cauchy, :Laplace, :Beta, :Uniform]
    @eval begin
        function Distributions.$dist(μ::Sym, σ::Sym)
            stats.$dist(:dist, μ,σ) |> SymPy.density
        end

        Distributions.$dist(μ,σ) = $dist(promote(μ,σ)...)
    end
end

export sym
sym(s::Symbol) = sympy.IndexedBase(s, real=true)
sym(s) = Base.convert(Sym, s)
function sym(expr::Expr)
    @match expr begin
        Expr(:call, f, args...) => :($f($(map(sym,args)...)))
        :($x[$j]) => begin
            j = symbols(:j, cls=sympy.Idx)
            x = sympy.IndexedBase(x, real=true)
            return x[j]
        end
        _ => begin
                 @show expr
                 error("sym: Argument type not implemented")
             end
    end
end

# export symlogpdf
# # function symlogpdf(m::Model)
# #     m = canonical(m)

# #     result = @q begin
# #         ctx = Dict()

# #         ℓ = zero(Soss.Sym)
# #     end

# #     exprs = []

# #     for st in map(v -> findStatement(m,v), toposortvars(m))
# #         push!(exprs, symlogpdf(st))
# #     end

# #     append!(result.args, exprs)

# #     # result

# #     push!(result.args, :(ctx,ℓ))

# #     result

# #     # expandSums(ℓ)
# # end

# # function symlogpdf(st::Soss.Sample)
# #     d = st.rhs
# #     x = st.x
# #     :(ℓ += $(symlogpdf(d,x)))
# # end

# # function symlogpdf(st::Soss.Assign)
# #     val = st.rhs
# #     x = st.x
# #     :(ctx[$(QuoteNode(x))] = Expr(:$, $val))
# #     # :($x = $val)
# # end

# # function symlogpdf(d::Expr, x::Symbol)
# #     @match d begin
# #         :(iid($n,$dist)) => begin
# #             j = symbols(:j, cls=sympy.Idx)
# #             dist = sym(dist)
# #             x = sympy.IndexedBase(x)
# #             n = sym(n)
# #             :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
# #         end

# #         :(For($f, (1:$n,))) => begin
# #             n = sym(n)
# #             @match f begin
# #                 :(($j,) -> begin $lineno; $dist end) => begin
# #                             j = symbols(j) # , cls=sympy.Idx)
# #                             # @show j
# #                             dist = sym(dist)
# #                             # @show dist
# #                             x = sympy.IndexedBase(x)
# #                             return :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
# #                 end

# #                 f => begin
# #                     @show f
# #                     error("symlogpdf: bad argument")
# #                 end
# #             end

# #         end

# #         _ => :(logpdf($(sym(d)), $(sym(x))))
# #     end
# # end

export expandSums
function expandSums(s::Sym)
    s.args==() && return s
    if s.func==sympy.Sum
        s2 = expandSum(s.args...)
        # Maybe it didn't simplify and we're done
        s2 == s && return s2
        # Or not, and we recurse
        s = s2
    end
    newargs = [expandSums(t) for t in s.args]
    result = s.func(newargs...)
    # @show result - s
    # @assert sympy.simplify(result - s) == sym(0)
    return result
end

function expandSum(s::Sym, limits::Sym...)

    hasIdx(s) || return maybeSum(s, limits...)

    func = s.func
    args = s.args
    args == () && return maybeSum(s,limits...)

    if func == sympy.Add
        return func([expandSum(t, limits...) for t in args]...)
    elseif func == sympy.Mul
        return expandMulSum(args, limits...)
    else
        return sympy.Sum(s, limits...)
    end

end

# Expand a sum of a `Mul`
# Currently a greedy algorithm, may need to optimize later
function expandMulSum(factors::NTuple{N,Sym}, limits::Sym...) where {N}
    limits == () && return prod(factors)

    for fac in factors
        for lim in limits
            (ix, ixlo, ixhi) = lim.args
            if ix ∉ fac
                inSummand = prod(allbut(factors, fac))
                inSum = expandSum(inSummand, lim)
                outLims = allbut(limits, lim)
                return expandSum(fac*inSum, outLims...)
            end
        end
    end
    return maybeSum(prod(factors), limits...)
end

function atoms(s::Sym)
    result = free_symbols(s)
    union(result, map(x -> x.args, result)...)
end

import Base.in
function Base.in(j::Sym, s::Sym)
    j ∈ atoms(s)
    # for t in s.args
    #     if j==t || in(j,t)
    #         return true
    #     end
    # end
    # return false
end

hasIdx(s::Sym) = any(startswith.(getproperty.(Soss.atoms(s), :name), "_j"))

function allbut(tup, x)
    result = filter(collect(tup)) do v
        v ≠ x
    end
    tuple(result...)
end

# Force computation of sums that don't involve the index
function maybeSum(t::Sym, limits::Sym...)
    length(limits) > 0 || return t

    for lim in limits
        (ix, ixlo, ixhi) = lim.args
        ix ∈ t || return maybeSum(t * (ixhi - ixlo + 1), allbut(limits, lim)...)
    end

    return sympy.Sum(t, limits...)
end

# # integrate(exp(ℓ), (sym(:μ), -oo, oo), (sym(:logσ),-oo,oo))
export marginal
function marginal(ℓ,v)
    f = ℓ.func
    f == sympy.Add || return ℓ
    newargs = filter(t -> sym(v) in t, collect(ℓ.args))
    foldl(+,newargs)
end

# marginal(m::Model, v) = marginal(m |> symlogpdf, v)

# # We should be able to reason about a marginal from its derivative
# export dmarginal
# function dmarginal(ℓ, v)
#     @as x ℓ begin
#         marginal(x,sym(v))
#         diff(x, sym(v))
#         expand(x)
#         sympy.collect(x, sym(v))
#     end
# end

# dmarginal(m::Model, v) = dmarginal(m |> symlogpdf, v)

# logpdf(Normal(sym(:μ),sym(:σ)), :x) |> SymPy.cse

# # macro symdist(n, dist)
# #     p = Expr(:tuple,gensym.(Symbol.(:p,1:n)))
# #     @q begin
# #         $dist($p)
# #     end
# # end
# # @macroexpand @symdist(4,Normal)

function symvar(st::Sample)
    st.rhs.args[1] ∈ [:For, :iid] && return IndexedBase(st.x)
    return sym(st.x)
end

export sourceSymlogpdf
function sourceSymlogpdf()
    function(_m::Model)
        function proc(_m, st :: Assign)
            # :($(st.x) = $(st.rhs))
            x = st.x
            xname = QuoteNode(x)
            :($x = sympy.IndexedBase($xname))
        end

        function proc(_m, st :: Sample)
            @q begin
                _ℓ += symlogpdf($(st.rhs), $(st.x))
            end
        end
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        function wrap(kernel)
            q = @q begin
                _ℓ = 0.0
            end

            for x in variables(_m)
                xname = QuoteNode(x)
                push!(q.args, :($x = sympy.IndexedBase($xname)))
            end

            for st in map(v -> findStatement(_m,v), toposortvars(_m))

                typeof(st) == Sample || continue
                x = st.x
                xname = QuoteNode(x)
                rhs = st.rhs
                xsym = ifelse(rhs.args[1] ∈ [:For, :iid]
                    , :(sympy.IndexedBase($xname))
                    , :(sym($xname))
                )
                # push!(q.args, :($x = $xsym))
            end

            @q begin
                $q
                $kernel
                return _ℓ
            end
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end

For(f, θ::Sym) = For(f, (θ,))

function For(f::F, θ::NTuple{N,Sym}) where {F, N}
    For{F,NTuple{N,Sym},Sym,Sym}(f,θ)
end

# logpdf(d::For{F,N,T}, x::Array{Symbol, N}) where {F,N,T} = logpdf(d,sym.(x))

export symlogpdf
function symlogpdf(d::For{F,T,D,X}, x::Sym) where {F, N, J <: Union{Sym,Integer}, T <: NTuple{N,J}, D,  X}
    js = symbols.(Symbol.(:_j,1:N), cls=sympy.Idx)
    x = sympy.IndexedBase(x)
    result = symlogpdf(d.f(js...), x[js...])

    for k in N:-1:1
        result = sympy.Sum(result, (js[k], 1, d.θ[k]))
    end
    result
end

symlogpdf(d::iid, x::Sym) = symlogpdf(For(j -> d.dist, d.size), x)

symlogpdf(d::For{F,T,D,X}, x::Symbol) where {F,T,D,X} = symlogpdf(d,sym(x))

symlogpdf(d::Normal, x::Sym) = symlogpdf(Normal(sym(d.μ),sym(d.σ)), x)

symlogpdf(d::Cauchy, x::Sym) = symlogpdf(Cauchy(sym(d.μ),sym(d.σ)), x)

symlogpdf(d::Beta, x::Sym) = symlogpdf(Beta(sym(d.α),sym(d.β)), x)

# @generated function symlogpdf(d,x::Sym)
#     quote
#         args = propertynames(d)

#     end
# end

logpdf(d::Sym, x::Sym) = symlogpdf(d,x)

function symlogpdf(d::Sym, x::Sym)
    d.func
    result = d.pdf(x) |> log
    sympy.expand_log(result,force=true)
end

symlogpdf(d,x::Sym) = logpdf(d,x)

function symlogpdf(m::JointDistribution)
    return _symlogpdf(m.model)
end

function symlogpdf(m::Model)
    return _symlogpdf(m)
end

@gg function _symlogpdf(_m::Model)
    type2model(_m) |> canonical |> sourceSymlogpdf()
end

# # s = symlogpdf(normalModel).args[7].args[3]

# # export fexpr
# # fexpr = quote
# #     f = function(μ,σ,x)
# #         a = $(codegen(s))
# #         return a
# #     end
# # end

# # julia> i,j = sympy.symbols("i j", integer=True)
# # (i, j)

# # julia> x = sympy.IndexedBase("x")
# # x

# # julia> a = sympy.Sum(x[i], (i, 1, j))
# #   j
# #  ___
# #  ╲
# #   ╲   x[i]
# #   ╱
# #  ╱
# #  ‾‾‾
# # i = 1

# # julia> SymPy.walk_expression(a)
# # :(Sum(Indexed(IndexedBase(x), i), (:i, 1, :j)))

export tolatex
function tolatex(ℓ::SymPy.Sym)
    r = r"_j(?<num>\d+)"
    s = s"j_{\g<num>}"
    Base.replace(sympy.latex(ℓ), r => s)
end

export foldConstants
function foldConstants(s::Sym)
    s.func ==sympy.Integer && return s
    isempty(free_symbols(s)) && return Float64(SymPy.N(s))
    s.func == sympy.Sum && return sympy.Sum(foldConstants(s.args[1]), s.args[2:end]...)
    s.func == sympy.Indexed && return s
    s.func == sympy.IndexedBase && return s
    s.func == sympy.Symbol && return s
    newargs = foldConstants.(s.args)
    return s.func(newargs...)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Model And Observed Maps
#
# Here we (1) read `CBIOMES model output` and `ocean color data` from file, and (2) illustrate these arrays as global maps in `julia`.
#
# <img src="../figs/cbiomes-01.png" alt="Drawing" style="height: 100px;"/>

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Read variables from file
#
# Two-dimensional arrays for longitude (`lon`), latitude (`lat`), irradiance reflectance from a model (`drwn3_Rirr` at `wv_drwn3`), and remotely sensed reflectance from satellite data (`cci_Rrs_490` at 490nm) are read from files in the `samples/` folder.

# + {"slideshow": {"slide_type": "-"}}
dirIn="../samples/"

fld = Array{Float32,2}(undef,(720,360))
fid = open(dirIn*"lon.bin"); read!(fid,fld); lon = hton.(fld)

fld = Array{Float32,2}(undef,(720,360))
fid = open(dirIn*"lat.bin"); read!(fid,fld); lat = hton.(fld)

fld = Array{Float32,2}(undef,(720,360))
fid = open(dirIn*"cci_Rrs_490.bin"); read!(fid,fld); cci_Rrs_490 = hton.(fld)
cci_Rrs_490[findall(cci_Rrs_490.==0)].=NaN

fld = Array{Float32,3}(undef,(720,360,13))
fid = open(dirIn*"drwn3_Rirr.bin"); read!(fid,fld); drwn3_Rirr = hton.(fld)
drwn3_Rirr[findall(drwn3_Rirr.==0)].=NaN;

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Model and data wavebands
#
# Currently, the `OC-CCI` [satellite data set](https://esa-oceancolour-cci.org) provides remotely sensed reflectance at 6 wavelengths (`wv_cci` in `nm`) while the `CBIOMES-global` [ocean model](https://cbiomes.readthedocs.io/) outputs irradiance reflectance at 13 wavelengths (`wv_drwn3` in `nm`).
# -

wv_cci=[412, 443, 490, 510, 555, 670]
wv_drwn3=[400,425,450,475,500,525,550,575,600,625,650,675,700];

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Display model and data maps
#
# The `Plots.jl` package provides convenient plotting recipes like `heatmap`. the following command assumes that `Plots.jl` has aleary been installed using `julia`'s package manager (documentation available [here](https://docs.julialang.org/en/)).
# -

using Plots

# + {"slideshow": {"slide_type": "subslide"}}
ii=3
heatmap(vec(lon[:,1]),vec(lat[1,:]), transpose(drwn3_Rirr[:,:,ii]), clims=(0,0.08))
title!("irradiance reflectance from a model (at $(wv_drwn3[ii])nm)")

# + {"slideshow": {"slide_type": "subslide"}}
heatmap(vec(lon[:,1]),vec(lat[1,:]), transpose(cci_Rrs_490), clims=(0,0.01))
title!("remotely sensed reflectance from satellite (at 490nm)")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Your turn!
#
# Here is one idea in case you want to take this a bit further:
#
# - Turn recipes from `OceanColourAlgorithms.ipynb` into functions.
# - Apply these to all points in `drwn3_Rirr` using for loops or broadcast.
# - Plot the resulting map of `Rrs` to compare with the data map (`cci_Rrs_490`).

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks