text:
using Knet, Test, Base.Iterators, Printf, LinearAlgebra, Random,  IterTools, CUDA
Knet.atype() = KnetArray{Float32}
include("/scratch/users/bberkay14/workfolder/language_to_logic/seq2tree/seq2tree/atis/tree.jl")
include("/scratch/users/bberkay14/workfolder/language_to_logic/seq2tree/seq2tree/atis/util.jl")
Knet.atype() = KnetArray{Float32}
using JLD2
using ArgParse
using Distributions

struct Memory; w; end

struct Attention; wquery; wattn; scale; end

struct Embed; w; end

function Embed(opt, vocabsize::Int, embedsize::Int)

    return Embed(param(rand(Uniform(-opt["init_weight"], opt["init_weight"]),embedsize,vocabsize)))
end

function (l::Embed)(x)

    return l.w[:,x]
end

struct Linear; w; b; end

function Linear(opt, inputsize::Int, outputsize::Int)

    return Linear(param(rand(Uniform(-opt["init_weight"], opt["init_weight"]),outputsize,inputsize)),param0(outputsize))
end

function (l::Linear)(x)

    return (mmul(l.w , x) .+ l.b)
end

struct L2L
    enc_lstm
    dec_lstm
    memory
    attention
    srcembed
    tgtembed
    dropout
    projection
end

function L2L(src_input_size, tgt_input_size, opt)
    srcembed = Embed(opt, src_input_size, opt["rnn_size"])
    tgtembed = Embed(opt, tgt_input_size, opt["rnn_size"])
    enc_lstm = RNN(opt["rnn_size"], opt["rnn_size"]; dropout=0, rnnType=:lstm)
    #input feeding deneme için aşağıdaki kodu aç
    #dec_lstm = RNN(opt["rnn_size"], 3*opt["rnn_size"]; dropout=0, rnnType=:lstm)
    dec_lstm = RNN(opt["rnn_size"], 2*opt["rnn_size"]; dropout=0, rnnType=:lstm)
    memory =1
    attention = Attention(1, param(rand(Uniform(-opt["init_weight"], opt["init_weight"]),opt["rnn_size"],2*opt["rnn_size"])), param(1) )

    dropout = opt["dropoutrec"]
    projection = Linear(opt, opt["rnn_size"],  tgt_input_size)
    return L2L(enc_lstm, dec_lstm, memory, attention, srcembed, tgtembed, dropout, projection)
end

function (m::Memory)(x)
    values = permutedims(x, [1, 3, 2])

    keys = mmul(m.w, values)
    return keys, values
end
# You can use the following helper function for scaling and linear transformations of 3-D tensors:
mmul(w,x) = (w == 1 ? x : w == 0 ? 0 : reshape(w * reshape(x,size(x,1),:), (:, size(x)[2:end]...)))

function encode(l::L2L, src, prev_c, prev_h)
    l.enc_lstm.h, l.enc_lstm.c = prev_c, prev_h
    src_embedding = l.srcembed(src)

    src_embedding = dropout(src_embedding,l.dropout; drop=true)
    src_encoded = l.enc_lstm(src_embedding)

    #s.decoder.h = permutedims(reshape(permutedims(src_encoded[:,:,end], (2,1)), (size(src_encoded[:,:,end], 2),Int64(size(src_encoded[:,:,end], 1) /2) ,2)), (2,1,3))
    l.dec_lstm.h = l.enc_lstm.h
    l.dec_lstm.c = l.enc_lstm.c
    return l.enc_lstm.c, l.enc_lstm.h
end

function (a::Attention)(cell, mem)
    query = mmul(a.wquery, cell)
    keys, vals = mem
    #attention_scores = bmm( permutedims(query, (3,1,2)) , keys)
    attention_scores = bmm( reshape(query, (size(query,3),size(query,1),size(query,2))) , keys)
    scaled_attention = a.scale .* attention_scores
    scaled_attention = softmax(scaled_attention; dims=2)
    context_vector = bmm(vals, permutedims(scaled_attention, (2,1,3)))

    concatenated_vectors = vcat(cell , permutedims(context_vector, (1,3,2)) )

    final_vector = mmul(a.wattn , concatenated_vectors)

    return final_vector
end

#function decode(l::L2L, tgt, mem, prev, parent_h)
function decode(l::L2L, tgt,  parent_h)
    tgt_embedding = l.tgtembed(tgt)
    #input feeding istiyorsan aşağıdaki kodu aç
    #input = vcat(tgt_embedding, parent_h,prev)
    input = vcat(tgt_embedding, parent_h)
    tgt_embedding = dropout(tgt_embedding, l.dropout; drop=true)
    decoder_output = l.dec_lstm(input)
    #attention_vector = l.attention(decoder_output, mem)
    #return attention_vector
    return l.dec_lstm.c, l.dec_lstm.h
end

function attention_decode(l::L2L, src_encoded, decoder_output)
    mem = l.memory(src_encoded)
    attention_vector = l.attention(decoder_output, mem)
    projected = l.projection(attention_vector)
    return projected
end

function eval_training(opt, train_loader,l::L2L, using_gpu, word_manager, form_manager)
    #encoder_optimizer.zero_grad()
    #decoder_optimizer.zero_grad()
    #attention_decoder_optimizer.zero_grad()
    enc_batch, enc_len_batch, dec_tree_batch = random_batch(train_loader)
    enc_max_len =  size(enc_batch, 2)

    #enc_outputs = zeros((enc_batch.size(0), enc_max_len, encoder.hidden_size), requires_grad=True)
    enc_outputs = param(size(enc_batch,1), enc_max_len, opt["rnn_size"]; init=zeros)
    #param((size(enc_batch,1), enc_max_len, opt["rnn_size"]); init=zeros)
    #if using_gpu:
    #    enc_outputs = enc_outputs.cuda()
    #enc_s = {}
    enc_s = Dict()
    for j in range(0,stop=opt["enc_seq_length"] )
        enc_s[j] = Dict()
    end

    dec_s = Dict()
    for i in range(0,stop=opt["dec_seq_length"] )
        dec_s[i] = Dict()
        for j in range(0,stop=opt["dec_seq_length"])
            dec_s[i][j] = Dict()
	end
    end

    for i in range(1, stop=2)
        #param a dimension belirtirken, dimensionlar knette farklı olabilir BAK
        enc_s[0][i] = param( opt["rnn_size"], opt["batch_size"]; init=zeros)

    end

    for i in range(1, stop=enc_max_len )
	println("Here")
	println(enc_s[i])
	println(length(enc_s))
	println(length(enc_s[i]))
	println(enc_batch)
	println(i)
        enc_s[i][1], enc_s[i][2] = encode(l, enc_batch[:,i], enc_s[i-1][1], enc_s[i-1][2])
        enc_outputs[:, i, :] = enc_s[i][2]
    end

    # tree decode
    queue_tree = Dict()
    for i in range(1, stop=opt["batch_size"])
        queue_tree[i] = []

        push!(queue_tree[i],Dict("tree"=>dec_tree_batch[i], "parent"=>0, "child_index"=>1))
    end
    loss = 0
    cur_index, max_index = 1,1
    dec_batch = Dict()
    #print(queue_tree[1][0]["tree"].to_string());exit()
    while (cur_index <= max_index)
        #print(cur_index)
        # build dec_batch for cur_index
        max_w_len = -1
        batch_w_list = []
        for i in range(1, stop=opt["batch_size"])
            w_list = []
            if (cur_index <= size(queue_tree[i],1))
                t = queue_tree[i][cur_index ]["tree"]
                for ic in range(1, stop=t.num_children)
                    #print("children ")
                    #print(ic +1)
                    if t.children[ic] isa Tree
                        push!(w_list, 4)
                        push!(queue_tree[i], Dict("tree"=>t.children[ic], "parent"=>cur_index, "child_index"=>ic))
                    else
                        push!(w_list, t.children[ic])
                    end
		end
                if size(queue_tree[i],1) > max_index
                    max_index = size(queue_tree[i],1)
                end
            end
            if length(w_list) > max_w_len
                max_w_len = length(w_list)
            end
            push!(batch_w_list, w_list)
        end
        #dec_batch[cur_index] = torch.zeros((opt.batch_size, max_w_len + 2), dtype=torch.long)
        #dec_batch[cur_index] = param(opt["batch_size"], max_w_len + 2; init=zeros)
	dec_batch[cur_index] = Int64.(zeros(opt["batch_size"], max_w_len + 2))
        for i in range(1, stop=opt["batch_size"])
            w_list = batch_w_list[i]
            if length(w_list) > 0
                for j in range(1, stop=length(w_list))
                    dec_batch[cur_index][i][j+1] = w_list[j]
                end
                # add <S>, <E>
                if cur_index == 1
                    dec_batch[cur_index][i][1] = 1
                else
                    dec_batch[cur_index][i][1] = get_symbol_idx(form_manager, '(')
                end
                dec_batch[cur_index][i][length(w_list) + 2] = 2
            end
        end

        for j in range(1, stop=2)
            #dec_s[cur_index][0][j] = torch.zeros((opt.batch_size, opt.rnn_size), dtype=torch.float, requires_grad=True)
            dec_s[cur_index][0][j] = param(opt["batch_size"], opt["rnn_size"] ; init=zeros)
            if using_gpu
                #dec_s[cur_index][0][j] = dec_s[cur_index][0][j].cuda()
                dec_s[cur_index][0][j] = Knet.KnetArray(dec_s[cur_index][0][j])
            end
        end

        if cur_index == 1
            for i in range(1, stop=opt["batch_size"])
                dec_s[1][0][1][i, :] = enc_s[enc_len_batch[i]][1][i, :]
                dec_s[1][0][2][i, :] = enc_s[enc_len_batch[i]][2][i, :]
            end

        else
            for i in range(1, stop=opt["batch_size"])
                if (cur_index <= size(queue_tree[i],1))
                    par_index = queue_tree[i][cur_index ]["parent"]
                    child_index = queue_tree[i][cur_index ]["child_index"]
                    dec_s[cur_index][0][1][i-1,:] = dec_s[par_index][child_index][1][i-1,:]
                    dec_s[cur_index][0][2][i-1,:] = dec_s[par_index][child_index][2][i-1,:]
                end
            end
        end
        gold_string = " "
        parent_h = dec_s[cur_index][0][2]
        for i in range(1, stop=size(dec_batch[cur_index], 2) - 1)

            dec_s[cur_index][i][1], dec_s[cur_index][i+1][2] = decode(l, dec_batch[cur_index][:,i],  parent_h)
            pred = attention_decode(l, enc_outputs, dec_s[cur_index][i][2])
            loss += nll(pred, dec_batch[cur_index][:,i+1]; dims=1, average=false)

        end
        cur_index = cur_index + 1
    end

    loss = loss / opt["batch_size"]
    return loss
end

function main(opt)
    Random.seed!(opt["seed"])
    @load  (string(opt["data_dir"]) * "/map.jld2") managers
    #------------
    word_manager, form_manager = managers
    using_gpu = false

    #L2L(src_input_size, tgt_input_size, opt, output_size)
    model = L2L(word_manager.vocab_size, form_manager.vocab_size, opt)

    ##-- load data
    train_loader = MinibatchLoader(opt, "train", using_gpu)

    if !isdir(opt["checkpoint_dir"])
        mkdir(opt["checkpoint_dir"])
    end

    ##-- start training
    step = 0
    epoch = 0
    optim_state = Dict("learningRate" => opt["learning_rate"], "alpha" =>  opt["decay_rate"])
    # default to rmsprop

    bestloss, bestmodel = eval_training(opt, train_loader,model,using_gpu, word_manager, form_manager), deepcopy(model)
    #if opt["opt_method"] == 0
    #progress!(adam(model, traindata), seconds=seconds) do y
    progress!(rmsprop(model, traindata; lr=optim_state["learningRate"], rho=optim_state["alpha"]), seconds=seconds) do y
        #devloss = loss(model, cdev)
        devloss = eval_training(opt, train_loader,model,using_gpu, word_manager, form_manager)
        if devloss < bestloss
            bestloss, bestmodel = devloss, deepcopy(model)
        end
        #println(stderr)
        (dev=devloss, tst=tstloss, mem=Float32(CUDA.usage[]))
    end

    iterations = opt["max_epochs"] * train_loader.num_batch
    #start_time = time.time()
    restarted = false
    # TODO revert back after tests
    #iterations = 2
    for i in range(1,stop=iterations)
        #epoch = i // train_loader.num_batch
        epoch = fld(i,train_loader.num_batch)
        train_loss = eval_training(opt, train_loader, encoder, decoder, attention_decoder, encoder_optimizer, decoder_optimizer, attention_decoder_optimizer, criterion, using_gpu, word_manager, form_manager)

        progress!(rmsprop(model, traindata; lr=optim_state["learningRate"], rho=optim_state["alpha"]), seconds=seconds) do y
            #devloss = loss(model, cdev)
            devloss = eval_training(opt, train_loader,model,using_gpu, word_manager, form_manager)
            if devloss < bestloss
                bestloss, bestmodel = devloss, deepcopy(model)
            end
        end
        #exponential learning rate decay
        if opt["opt_method"] == 0
            if mod(i , train_loader.num_batch) == 0 & opt["learning_rate_decay"] < 1
                if epoch >= opt["learning_rate_decay_after"]
                    decay_factor = opt["learning_rate_decay"]
                    optim_state["learningRate"] = optim_state["learningRate"] * decay_factor #decay it

                end
            end
        end
        if (epoch == opt["restart"]) &  !restarted
            restarted = True
            optim_state["learningRate"] = opt["learning_rate"]

        end

        #on last iteration
        if i == iterations
            checkpoint = Dict()

            checkpoint["lang2logic_model"] = model
            checkpoint["opt"] = opt
            checkpoint["i"] = i
            checkpoint["epoch"] = epoch

            Knet.save(opt["checkpoint_dir"] *"/model_seq2seq.jld2", "checkpoint", checkpoint)
        end
        if train_loss != train_loss
            print("loss is NaN.  This usually indicates a bug.")
            break
        end
    end
    save && Knet.save("attn-$(Int(time_ns())).jld2", "model", bestmodel)
    #bleu && Main.bleu(bestmodel,dev)
end

start = time()
s = ArgParseSettings()
@add_arg_table s begin
    "--gpuid"
        arg_type = Int
        help = "which gpu to use. -1 = use CPU"
        default = 0
    "--data_dir"
        help = "data path"
        arg_type = String
        default = "../data/"
    "--seed"
        help = "torch manual random number generator seed"
        arg_type = Int
        default = 123
    "--checkpoint_dir"
        help = "filename to autosave the checkpont to. Will be inside checkpoint_dir/"
        arg_type = String
        default = "checkpoint_dir"
    "--savefile"
        help = "max vocab size"
        arg_type = String
        default = "save"
    "--print_every"
        help = "how many steps/minibatches between printing out the loss"
        arg_type = Int
        default = 2000
    "--rnn_size"
        help = "size of LSTM internal state"
        arg_type = Int
        default = 200
    "--num_layers"
        help = "number of layers in the LSTM"
        arg_type = Int
        default = 1
    "--dropout"
        help = "dropout for regularization, used after each RNN hidden layer. 0 = no dropout"
        arg_type = Float64
        default = 0.3
    "--dropoutrec"
        help = "dropout for regularization, used after each c_i. 0 = no dropout"
        arg_type = Float64
        default = 0.3
    "--enc_seq_length"
        help = "number of timesteps to unroll for"
        arg_type = Int
        default = 60
    "--dec_seq_length"
        help = "number of timesteps to unroll for"
        arg_type = Int
        default = 220
    "--batch_size"
        help = "number of sequences to train on in parallel"
        arg_type = Int
        default = 20
    "--max_epochs"
        help = "number of full passes through the training data"
        arg_type = Int
        default = 130
    "--opt_method"
        help = "optimization method: 0-rmsprop 1-sgd"
        arg_type = Int
        default = 0
    "--learning_rate"
        help = "learning rate"
        arg_type = Float64
        default = 0.007
    "--init_weight"
        help = "initailization weight"
        arg_type = Float64
        default = 0.08
    "--learning_rate_decay"
        help = "learning rate decay"
        arg_type = Float64
        default = 0.98
    "--learning_rate_decay_after"
        help = "in number of epochs, when to start decaying the learning rate"
        arg_type = Int
        default = 5
    "--restart"
        help = "in number of epochs, when to restart the optimization"
        arg_type = Int
        default = -1
    "--decay_rate"
        help = "decay rate for rmsprop"
        arg_type = Float64
        default = 0.95
    "--grad_clip"
        help = "clip gradients at this value"
        arg_type = Int
        default = 5
end

args = parse_args(s)
main(args)
endtime = time()
time_passed = (endtime - start) / 60
print("total time: "* string(time_passed) * "minutes")


text:
module Analyse
    using CSV
    using DataFrames
    using Distributions
    using JSON3
    using LightGraphs, MetaGraphs
    using Plots
    using Plots.PlotMeasures
    using PyCall
    using SparseArrays

    import Dates: now, format
    import StatsBase: countmap
    """
        DTG()

    return DTG string for logging
    """
    function DTG()
        return "$(format(now(), "HH:MM:SS"))"
    end

    """
        loaddata(userfile::Union{String,Array{String,1}}, msgfiles::Union{String,Array{String,1}})

    From a (list of) user and message CSV file(s), get the dataframes. Works on the raw files downloaded from the Twitter misinformation report.

    ### kwargs
    * :usrfields : `Array{Symbol,1}` - the columns of the user dataframe you want to use (default: all)
    * :msgfields : `Array{Symbol,1}` - the columns of the message dataframe you want to use (default: all)

    ### Examples
    ```Julia
    loaddata("/path/to/files/userfile.csv",
                ["/path/to/files/messagefile_1.csv",
                    "/path/to/files/messagefile_1.csv"])

    ```

    ### Notes
    * Makes use of [`DataFrame!`](https://juliadata.github.io/DataFrames.jl/stable/lib/functions/#DataFrames.DataFrame!)
    combined with [`CSV.File`](https://juliadata.github.io/CSV.jl/stable/#CSV.File) with a selection of columns for increased performance.
    Column types can also be passed for an additional speedup (see also [`CSV.File`](https://juliadata.github.io/CSV.jl/stable/#CSV.File) kwargs)
    * If multiple treads are available, they will be used to read the CSV files in parallel.

    See also: ['dropdataframes!'](@ref), ['setgraph!'](@ref), ['describe!'](@ref)
    """
    function loaddata(userfiles::Array{String,1}, msgfiles::Array{String,1}; kwargs...)
        @debug "$(DTG()) - loading desinfo data files"
        # default settings
        defusrfilter = [:userid, :user_screen_name]
        defmsgfilter = [:tweetid, :userid, :is_retweet, :retweet_tweetid, :retweet_userid, :in_reply_to_tweetid, :in_reply_to_userid]

        # Load up CSV data with only selected columns
        usrdf = vcat([DataFrame(CSV.File(f; select=get(kwargs,:usrfields, defusrfilter),
                                            limit=get(kwargs,:limit, nothing)))
                        for f in userfiles]...)
        msgdf = vcat([DataFrame(CSV.File(f; select=get(kwargs,:msgfields, defmsgfilter),
                                            limit=get(kwargs,:limit, nothing),
                                            truestrings=["true", "True"], falsestrings=["false", "False"], # added for stablity
                                            types=Dict(:is_retweet => Bool)))
                        for f in msgfiles ]...)
        @debug "$(DTG()) - finished desinfo loading data files"
        return usrdf, msgdf
    end

    # additional methods for flexibility
    loaddata(userfile::String, msfgile::String; kwargs...) = loaddata([userfile], [msfgile]; kwargs...)
    loaddata(userfile::String, msfgiles::Array{String,1}; kwargs...) = loaddata([userfile], msfgiles; kwargs...)
    loaddata(userfiles::Array{String,1}, msfgile::String; kwargs...) = loaddata(userfiles, [msfgile]; kwargs...)

    """
        loaddata(path::String; kwargs...)

    Load data files from a path containing user and message CSV files.

    ### kwargs:
    * `:usrdata` : `String` - pattern used to recognize user data file(s) (default: "users_csv_hashed")
    * `:msgdata` : `String` - pattern used to recognize message data file(s) (default: "tweets_csv_hashed")
    * `:usrfields` : `Array{Symbol,1}` - the fields of the user dataframe you want to use (default: all)
    * `:msgfields` : `Array{Symbol,1}` - the fields of the message dataframe you want to use (default: all)
    * `:usrfilter` : `Array{String,1}` - an additional pattern that should found in the message data files (e.g. "_2020", default:[""])
    * `:msgfilter` : `Array{String,1}` - an additional pattern that should found in the message data files (e.g. "_2020", default:[""])

    See also: ['dropdataframes!'](@ref), ['setgraph!'](@ref), ['describe!'](@ref)

    """
    function loaddata(path::String; kwargs...)
        usrfiles = joinpath.(path, filter(x -> occursin(get(kwargs, :usrdata, "users_csv_hashed"), x) &&
                                                    !occursin(".jsonl",x) &&
                                                any([occursin(t,x) for t in get(kwargs, :usrfilter, [""]) ]),
                                readdir(path)))
        msgfiles = joinpath.(path, filter(x -> occursin(get(kwargs, :usrdata, "tweets_csv_hashed"),x) &&
                                                    !occursin(".jsonl",x) &&
                                                any([occursin(t,x) for t in get(kwargs, :msgfilter, [""]) ]),
                                readdir(path)))
        @warn usrfiles, msgfiles
        return loaddata(usrfiles, msgfiles; kwargs...)
    end

    """
        get_externals(path::String)

    From a twitter misinformation report path, extract and store the data
    that still needs to be recovered (retweets/replies) using an additional tool
    in order to be able to generate an interaction network.
    """
    function get_externals(path::String)
        # export path
        fname = filter!(x->length(x)>0,split(path,"/"))[end]
        # load data
        @debug "$(DTG()) - loading data"
        udf,mdf = loaddata(path,
                            usrfields=[:userid, :user_screen_name],
                            msgfields=[:tweetid, :userid, :is_retweet, :retweet_userid,
                                       :retweet_tweetid, :in_reply_to_tweetid, :in_reply_to_userid])
        # get flagged users
        @debug "$(DTG()) - getting flagged users"
        fu = Set(udf.userid)
        # find retweets of non-flagged users
        @debug "$(DTG()) - getting retweets"
        retweets = unique(filter(row-> ismissing(row.retweet_userid), dropmissing(mdf, [:retweet_tweetid])).retweet_tweetid)
        # find replies to non-flagged users
        @debug "$(DTG()) - getting replies"
        replies = unique(filter(row-> row.in_reply_to_userid ∉ fu, dropmissing(mdf, [:in_reply_to_tweetid])).in_reply_to_tweetid)
        # write away the results
        @debug "$(DTG()) - writing results"
        open(joinpath(path, "$(fname) retweets.csv"), "w") do io
            for i in retweets
                print(io,i,"\n")
            end
        end

        open(joinpath(path, "$(fname) replies.csv"), "w") do io
            for i in replies
                print(io,i,"\n")
            end
        end
        @debug "$(DTG()) - get_externals finished"
    end

    """
        readjsonl(p::String; kind::Symbol=:retweet)

    Read external tweets from jsonl file into dataframe holding user id, user description & post id

    # Arguments:
    - `p::String`: the jsonl file contains the additional tweets that are required for the interaction graph
    - `kind:::Symbol`: kind of interacion you want (this is should occur in the filename)

    Verify this works as intended: IT DOES NOT:
    currently reads a folder and we don't want that => make it work for a single file
    """
    function readjsonl(p::String; kind::Symbol=:retweet)
        f = joinpath(p, filter(x-> occursin(isequal(kind, :retweet) ? "retweet" : "reply", x) && occursin("jsonl", x),readdir(p))[1])
        @debug "$(DTG()) - reading $(f)"
        msgids = Vector{Int64}()
        userids = Vector{String}() # user ids
        userdes = Vector{String}() # user description

        for line in eachline(f)
            res = JSON3.read(line)
            push!(msgids, parse(Int64, res["id_str"]))
            push!(userids, res["user"]["id_str"])
            push!(userdes, res["user"]["screen_name"])
            @debug "$(DTG()) - msg $(res["id_str"]) / $(res["id"]) from user ID $(res["user"]["id_str"]) ($(res["user"]["screen_name"]))" maxlog=1
        end
        return DataFrame([:msgid => msgids, :userid => userids, :user_screen_name => userdes])
    end

    """
    readjsonlnew(p::String)

    Read external tweets from jsonl file into dataframe holding: msgid, userid & user_screen_name

    This external file is the result of "hydrating tweets that were identified in the initial data set of the Twitter Information Operations dataset
    """
    function readjsonlnew(p::String)
        @debug "$(DTG()) - reading $(p)"
        msgids = Vector{Int64}()
        userids = Vector{String}() # user ids
        userdes = Vector{String}() # user description

        for line in eachline(p)
            res = JSON3.read(line)
            push!(msgids, parse(Int64, res["id_str"])) # Stored as integers
            push!(userids, res["user"]["id_str"])
            push!(userdes, res["user"]["screen_name"])
            @debug "$(DTG()) - msg $(res["id_str"]) / $(res["id"]) from user ID $(res["user"]["id_str"]) ($(res["user"]["screen_name"]))" maxlog=1
        end

        return DataFrame([:msgid => msgids, :userid => userids, :user_screen_name => userdes])
    end

    """
        directedretweetgraph(msgdf::DataFrame,msgmap::Dict, Uinfo::Dict, RUinfo::Dict)

    For a Twitter Misinformation Dataset, build directed & weighted graph from retweet interactions

    # Arguments
    - msgdf: DataFrame containing the parsed messages from Twitter Information Operations source
    - msgmap: Dict that maps messageid to a userid (linking external message to a user)
    - Uinfo: Dict that maps nodeid (integer value starting at 1) to a dictionary holding userid, user_screen_name & flagged status
    - RUinfo: Dict that maps userid (string) to nodeid (integer)

    # Method
    FOR all messages DO

        IF original message is coming from a flagged user
            add edge user > flagged user (info flow)
        ELSEIF retweeted message is coming from an external user
            IF the message has been recovered
                add edge external users > flagged user (info flow)
            ELSE
                skip the message ("weakness to method")

    # Notes
    if graph is established for multiple interaction types, the number of non-interacting nodes is likely to be higher.

    # See also
    ['directedreplygraph'](@ref), ['readjsonlnew'](@ref), ['loaddata'](@ref)
    """
    function directedretweetgraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict)
        # generate the graph, total number of users infered from Uinfo dict
        G = MetaDiGraph(length(Uinfo))
        # add nodes + node info to the graph
        for (node, atts) in Uinfo
            set_props!(G, node, atts)
        end
        # get edges from a mapping that holds the edge and its weight
        edges = Dict{Edge, Int}()
        for row in eachrow(filter(row -> row.is_retweet, msgdf)) # filter row retweets only
            src = nothing
            # determine source node
            if ismissing(row.retweet_userid) # if this is absent, we are dealing with an external tweet
                if haskey(msgmap, row.retweet_tweetid) # if it was recovered, we have a matching user
                    src = RUinfo[msgmap[row.retweet_tweetid]] # return node id from user id
                end
            else
                src = RUinfo[row.retweet_userid]
            end
            if isnothing(src)
                continue
            else
                dst = RUinfo[row.userid]
                e = Edge(src, dst)
                edges[e] = get!(edges, e, 0) + 1
            end
        end
        # add edges to the graph
        for (e,w) in edges
            add_edge!(G, e)
            set_prop!(G,e, :weight, w)
        end

        return G
    end

    """
        directedreplygraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict, flagged_users::Set)

    For a Twitter Misinformation Dataset, build directed & weighted graph from reply interactions

    # Arguments
    - msgdf: DataFrame containing the parsed messages from Twitter Information Operations source
    - msgmap: Dict that maps messageid to a userid (linking external message to a user)
    - Uinfo: Dict that maps nodeid (integer value starting at 1) to a dictionary holding userid, user_screen_name & flagged status
    - RUinfo: Dict that maps userid (string) to nodeid (integer)
    - flagged_users: Set of users that are flagged

    # Method
    FOR all messages DO

        IF original message is coming from a flagged user
            add edge user > flagged user (info flow)
        ELSEIF replies message is coming from an external user
            IF the message has been recovered
                add edge external users > flagged user (info flow)
            ELSE
                skip the message ("weakness to method")

    # Examples
    ```julia
    nothing
    ```

    # Notes
    if graph is established for multiple interaction types, the number of non-interacting nodes is likely to be higher.

    # See also
    ['directedretweetgraph'](@ref), ['readjsonlnew'](@ref), ['loaddata'](@ref)

    """
    function directedreplygraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict, flagged_users::Set)
        # generate the graph
        G = MetaDiGraph(length(Uinfo))
        ## add nodes & node info
        for (node, atts) in Uinfo
            set_props!(G, node, atts)
        end
        # get edges from a mapping that holds the edge and its weight
        edges = Dict{Edge, Int}()
        # filter replt rows only
        for row in eachrow(dropmissing(msgdf, [:in_reply_to_tweetid] ))
            src = nothing
            if row.in_reply_to_userid ∈ flagged_users # if user is flagged: get its id directly
                src = RUinfo[row.in_reply_to_userid]
            else # if it is an external user and recovered the twee: get id
                if haskey(msgmap, row.in_reply_to_tweetid)
                    src = RUinfo[msgmap[row.in_reply_to_tweetid]]
                end
            end
            if isnothing(src)
                continue
            else
                dst = RUinfo[row.userid]
                e = Edge(src, dst)
                edges[e] = get!(edges, e, 0) + 1
            end
        end
        # add edges to graph
        for (e,w) in edges
            add_edge!(G, e)
            set_prop!(G,e, :weight, w)
        end

        return G
    end

    """
        directedbipartiteretweetgraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict,)

    Generate directed bipartite retweet graph. Layer 1 = users, layer 2 = posts

    # Arguments:
    - `msgdf::DataFrame`: dataframe containing the message data in the dataset
    - `msgmap::Dict`: mapping of post id to user id
    - `Uinfo::Dict`: mapping of node if on graph to dict (twitter user id, user description, flagged status)
    - `RUinfo::Dict`: reverse mapping of Uinfo; maps twitter user id to node id

    """
    function directedbipartiteretweetgraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict)
        # generate the graph
        G = MetaDiGraph(length(Uinfo))
        ## add nodess & info
        # user nodes
        Nu = length(Uinfo)
        for (node, atts) in Uinfo
            set_props!(G, node, atts)
        end
        # post info is a mapping from post id to node id
        Pinfo = Dict{}()
        # posts from flagged users:
        p_1 = unique(filter(row -> !ismissing(row.retweet_userid) , msgdf).retweet_tweetid)
        # posts from non-flagged users that are recovered by hydration
        p_2 = unique(filter(row ->  ismissing(row.retweet_userid) &&
                                    !ismissing(row.retweet_tweetid) &&
                                    haskey(msgmap, row.retweet_tweetid),
                            msgdf).retweet_tweetid)
        for post in union(p_1, p_2)
            # add new node
            add_vertex!(G)
            @debug "number of vertices: $(nv(G))" maxlog=10
            # set settings
            set_props!(G, nv(G), Dict(:userid => post, :Label => "post", :flagged => false))
            Pinfo[post] = nv(G)
        end
        # reverse post mapping leading to:  node_id > post_id
        RPinfo = Dict(val => key for (key, val) in Pinfo)
        Np = length(Pinfo)

        ## add edges
        for row in eachrow(filter(row -> row.is_retweet, msgdf)) # filter row retweets only
            # edge from author to msg
            if ismissing(row.retweet_userid)
                if haskey(msgmap, row.retweet_tweetid) # in case of recovered tweet
                    # add edge from author to message (author_of_tweet_id > tweet_id)
                    aut = RUinfo[msgmap[row.retweet_tweetid]]
                    msg = Pinfo[row.retweet_tweetid]
                    add_edge!(G, aut, msg)
                    @debug "adding authorship from user $(msgmap[row.retweet_tweetid]) ($(src)) to msg $(row.retweet_tweetid) ($(dst))" maxlog=1
                    # add edge from message to spreader (tweet_id > userid)
                    add_edge!(G, msg, RUinfo[row.userid])
                end
            else
                # add edge from author to message (retweet_userid > tweet_id)
                aut = RUinfo[row.retweet_userid]
                msg = Pinfo[row.retweet_tweetid]
                add_edge!(G, aut, msg)
                # add edge from message to spreader
                add_edge!(G, msg, RUinfo[row.userid])
            end
        end

        set_props!(G, Dict(:Nu => Nu, :Np => Np, :kind => :retweet))

        return Pinfo, RPinfo, G
    end

    """
        directedbipartitereplygraph(msgmap::Dict, Uinfo::Dict, RUinfo::Dict, msgdf::DataFrame)

    Generate directed bipartite reply graph. Layer 1 = users, layer 2 = posts

    # Arguments:
    - `msgmap::Dict`: mapping of post id to user id
    - `Uinfo::Dict`: mapping of node if on graph to dict (twitter user id, user description, flagged status)
    - `RUinfo::Dict`: reverse mapping of Uinfo; maps twitter user id to node id
    - `msgdf::DataFrame`: dataframe containing the message data in the dataset
    - `flagged_users::Set`: set of the flagged users
    """
    function directedbipartitereplygraph(msgdf::DataFrame, msgmap::Dict, Uinfo::Dict, RUinfo::Dict, flagged_users::Set)
        # filter replies only
        mdf = dropmissing(msgdf, :in_reply_to_tweetid)
        G = MetaDiGraph(length(Uinfo))
        ## add nodes + info
        # user nodes
        Nu = length(Uinfo)
        for (node, atts) in Uinfo
            set_props!(G, node, atts)
        end
        # message mapping (post id => node id)
        Pinfo = Dict{}()
        # posts from flagged users
        p_1 = unique(filter(row -> row.in_reply_to_userid ∈ flagged_users, mdf).in_reply_to_tweetid)
        # posts from non-flagged users that are recovered by hydration
        p_2 = unique(filter(row -> row.in_reply_to_userid ∉ flagged_users && haskey(msgmap, row.in_reply_to_tweetid), mdf).in_reply_to_tweetid)
        # add post vertices
        for post in union(p_1, p_2)
                add_vertex!(G)
                @debug "number of vertices: $(nv(G))" maxlog=10
                set_props!(G, nv(G), Dict(:userid => post, :Label => "post", :flagged => false))
                Pinfo[post] = nv(G)
        end
        RPinfo = Dict(val => key for (key, val) in Pinfo)
        Np = length(Pinfo)

        ## add edges
        for row in eachrow(mdf)
            if row.in_reply_to_userid ∈ flagged_users
                # case of original message by flagged user
                aut = RUinfo[row.in_reply_to_userid]
                msg = Pinfo[row.in_reply_to_tweetid]
                add_edge!(G, aut, msg)
                @debug "adding flagged authorship from user $(aut) to msg $(msg)" maxlog=1
                # add edge from message to spreader
                add_edge!(G, msg, RUinfo[row.userid])
            else
                if haskey(msgmap, row.in_reply_to_tweetid)
                    # case of reply to recovered message from external user
                    aut = RUinfo[msgmap[row.in_reply_to_tweetid]]
                    msg = Pinfo[row.in_reply_to_tweetid]
                    add_edge!(G, aut, msg)
                    @debug "adding authorship from user $(aut) to msg $(msg)" maxlog=1
                    # add edge from message to spreader
                    add_edge!(G, msg, RUinfo[row.userid])
                end
            end
        end

        set_props!(G, Dict(:Nu => Nu, :Np => Np, :kind => :reply))

        return Pinfo, RPinfo, G
    end

    """
        grapher(path::String)

    Generate interaction graphs from a Twitter misinformation data file. Considers both retweets and replies. In the interaction graphs,
    users from both interactions will be included in order to allow multi-layer community detection based on the Leiden method.

    """
    function grapher(msgfile::Union{String,Array{String,1}}, userfile::Union{String,Array{String,1}}; α::Float64=0.05)
        #@assert occursin("tweets", msgfile)
        #@assert occursin("user", userfile)
        @debug "$(DTG()) - building graphs for $(msgfile)"
        # get raw data as dataframe
        usrdf, msgdf = loaddata(userfile, msgfile) # check this is working as intended

        # get external messages as dataframe
        @debug "$(DTG()) - loading external retweets"
        rtdf = isa(msgfile, String) ? readjsonlnew(msgfile*"_retweets.jsonl") : vcat([readjsonlnew(file*"_retweets.jsonl") for file in msgfile]...)
        @debug "$(DTG()) - loading external replies"
        rpdf = isa(msgfile, String) ? readjsonlnew(msgfile*"_replies.jsonl") : vcat([readjsonlnew(file*"_replies.jsonl") for file in msgfile]...)

        # establish nodes in graph
        Uinfo = Dict{Int64, Dict}()
        # add flagged users
        for row in eachrow(unique(usrdf, :userid))
            Uinfo[length(Uinfo) + 1] = Dict(:userid => row.userid, :Label => row.user_screen_name, :flagged => true)
        end
        # get other users from retweets
        for row in eachrow(unique(rtdf, :userid))
            Uinfo[length(Uinfo) + 1] = Dict(:userid => row.userid, :Label => row.user_screen_name, :flagged => false)
        end
        # reverse user mapping
        RUinfo = Dict(val[:userid] => key for (key, val) in Uinfo)

        # add other users from replies
        @debug "starting with replies"
        for row in eachrow(unique(rpdf, :userid))
            haskey(RUinfo, row.userid) && continue # check if user already known
            # if not, add them
            @debug "adding $(row.userid) to graph with node number $(length(Uinfo) + 1)" maxlog=10
            Uinfo[length(Uinfo) + 1] = Dict(:userid => row.userid, :Label => row.user_screen_name, :flagged => false)
            RUinfo[row.userid] = length(Uinfo)
        end

        # map messages id > user id
        ## map additional messages to user
        uniquemsg_rt = unique(rtdf, :msgid)
        msgmap_rt = Dict(zip(uniquemsg_rt.msgid, uniquemsg_rt.userid))
        uniquemsg_rp = unique(rpdf, :msgid)
        msgmap_rp = Dict(zip(uniquemsg_rp.msgid, uniquemsg_rp.userid))

        # generate graphs
        G_rt = directedretweetgraph(msgdf, msgmap_rt, Uinfo, RUinfo)
        G_rp = directedreplygraph(msgdf, msgmap_rp, Uinfo, RUinfo, Set(usrdf.userid))
        G_rt_bp = directedbipartiteretweetgraph(msgdf, msgmap_rt, Uinfo, RUinfo)[3]
        G_rp_bp = directedbipartitereplygraph(msgdf, msgmap_rp, Uinfo, RUinfo, Set(usrdf.userid))[3]
        G_rt_bp_proj = projectbipartite(G_rt_bp, α=α)[4]
        G_rp_bp_proj = projectbipartite(G_rp_bp, α=α)[4]

        @debug "retweet edge reduction $(round(ne(G_rt_bp_proj) / ne(G_rt) * 100))%"
        @debug "reply edge reduction $(round(ne(G_rp_bp_proj) / ne(G_rp) * 100))%"

        #return round(ne(G_rt_bp_proj) / ne(G_rt) * 100, digits=2), round(ne(G_rp_bp_proj) / ne(G_rp) * 100, digits=2), G_rt_bp_proj
        return G_rt, G_rp, G_rt_bp, G_rp_bp, G_rt_bp_proj, G_rp_bp_proj, Uinfo, RUinfo

    end

    """
    writeout to a edge and node file that can be read by gephi
    """
    function graphwriter(G::Union{MetaGraph, MetaDiGraph})
        io = IOBuffer()
        ## PART VERTICES
        println(io,"Id; Label; flagged")
        for v in vertices(G)
            println(io,"$(v); $(props(G, v)[:Label]); $(props(G, v)[:flagged])")
        end
        # write to file
        open("nodelist.csv", "w") do f
            print(f, String(take!(io)))
        end

        ## PART EDGES
        println(io,"Source; Target; Weight")
        for e in edges(G)
            println(io, "$(e.src); $(e.dst); $(props(G, e.src, e.dst)[:weight])")
        end
        open("edgelist.csv", "w") do f
            print(f, String(take!(io)))
        end

    end
    """
        graphbuilder(path::String; kind::Symbol=:retweet, α::Float64=0.001)

    Main function for graph generation. Reads path, loads up the data & builds MetaGraphs

    # Arguments
    - `path::String`: folder holding the data
    - `kind`::Symbol: type of interaction graph you want to build
    - `α::Float64`: limit for statistical significance of V-motif

    returns:
        - directed weighted interaction graph
        - directed bipartite interaction graph
        - projected binary interaction graph
        - observed V-motifs
        - expected V-motifs
        - p-values for projected user-interactions
    """
    function graphbuilder(path::String; kind::Symbol=:retweet, α::Float64=1e-2)
        @debug "$(DTG()) - building $(kind) graph"
        ## load up data from desinfo dataset
        usrdf, msgdf = loaddata(path,   usrfields=[:userid, :user_screen_name],
                                        msgfields=[:tweetid, :userid, :is_retweet, :retweet_userid, :retweet_tweetid, :in_reply_to_tweetid, :in_reply_to_userid])

        ## original messages per kind (as downloaded by hydrator) as DataFrame
        rtsdf = readjsonl(path, kind=kind)

        ## extract the user information
        Uinfo = Dict{Int64, Dict}()
        # get flagged users
        for row in eachrow(unique(usrdf, :userid))
            Uinfo[length(Uinfo) + 1] = Dict(:userid => row.userid, :Label => row.user_screen_name, :flagged => true)
        end
        # get other users
        for row in eachrow(unique(rtsdf, :userid))
            Uinfo[length(Uinfo) + 1] = Dict(:userid => row.userid, :Label => row.user_screen_name, :flagged => false)
        end
        # reverse user information
        RUinfo = Dict(val[:userid] => key for (key, val) in Uinfo)

        ## map additional messages to user
        uniquemsg = unique(rtsdf, :msgid)
        msgmap = Dict(zip(uniquemsg.msgid, uniquemsg.userid)) ##  maybe use integer value from the start???

        ## generate directed interaction network (retweet/reply)
        Gint = isequal(kind, :retweet) ? directedretweetgraph(msgmap, Uinfo, RUinfo, msgdf) : directedreplygraph(msgmap, Uinfo, RUinfo, msgdf, Set(usrdf.userid))

        ## generate directed bipartite usr-msg network (retweet/reply)
        Gbip = isequal(kind, :retweet) ? directedbipartiteretweetgraph(msgmap, Uinfo, RUinfo, msgdf)[3] :  directedbipartitereplygraph(msgmap, Uinfo, RUinfo, msgdf, Set(usrdf.userid))[3]

        ## generate projected user network
        V_star, V_exp, p_vals, Gproj = projectbipartite(Gbip, α=α)

        return Gint, Gbip, Gproj, V_star, V_exp, p_vals

    end

    "helper function to calculate the p-value of the poisson-binomial approximation"
    pval(v_expected, v_observed) = 1 - cdf(Poisson(v_expected), v_observed)

    "function to map bipartite user-post to user-user"
    function projectbipartite(G::MetaDiGraph; α::Float64=0.001)
        @debug "$(DTG()) - projecting graph with α = $(α)"
        # degrees of all users:
        k_λ_in =  indegree(G, 1 : get_prop(G, :Nu))
        k_λ_out = outdegree(G, 1 : get_prop(G, :Nu))
        # adjacency matrix:
        A = LightGraphs.LinAlg.adjacency_matrix(G)
        # real V motif matrix (sparse representation)
        Nu = get_prop(G, :Nu)
        V_star = A[1:Nu,:]*A[:,1:Nu]
        # expected V-motif matrix (sparse representation)
        # get locations and values
        I, J, V_star_val = findnz(V_star)
        V_exp_val = k_λ_out[I] .* k_λ_in[J] / get_prop(G, :Np)
        #expected_v_motifs(V_star, k_λ_out, k_λ_in, get_prop(G, :Np))#k_λ_out * k_λ_in' /get_prop(G, :Np)
        # p-values (sparse)
        p_vals = pval.(V_exp_val, V_star_val)
        #spzeros(Float64,get_prop(G, :Nu),get_prop(G, :Nu))#pval.(V_exp, V_star)
        # FDR procedure Benjamini–Hochberg
        #I, J, V = findnz(p_vals)
        k = sortperm(p_vals) # costly (1e6 edges: 0.09s; 1e7 edges: 1.57s; 1e8 elements: 23s)
        m = length(k)
        #@warn V_star[3121, 348], V_exp[3121, 348], p_vals[3121, 348]
        pvalfilter = p_vals .<= k/m*α # probleem!
        # generate graph from significant sparse edges
        A_filt = sparse(I[pvalfilter], J[pvalfilter], ones(Int,sum(pvalfilter)), Nu, Nu)
        G_filt = MetaDiGraph(A_filt)
        for i = 1:Nu
            set_props!(G_filt, i, props(G, i))
        end
        set_prop!(G_filt, :Nu, Nu)

        return V_star, sparse(I,J,V_exp_val), sparse(I,J,p_vals), G_filt
    end

    """
        Coherence_checks(G:MetaDiGraph)

    check if bipartite graph is built up as it should
    """
    function coherence_checks(GBip::MetaDiGraph)
        # degrees of all posts:
        k_γ_in =  indegree(GBip, get_prop(GBip, :Nu) + 1 : get_prop(GBip,:Nu) + get_prop(GBip, :Np))
        k_γ_out = outdegree(GBip, get_prop(GBip, :Nu) + 1 : get_prop(GBip,:Nu) + get_prop(GBip, :Np))
        # degrees of all users:
        k_λ_in =  indegree(GBip, 1 : get_prop(GBip, :Nu))
        k_λ_out = outdegree(GBip, 1 : get_prop(GBip, :Nu))
        @info """Running coherence checks for the graph\n
            $(get_prop(GBip, :Nu)) user nodes, $(get_prop(GBip, :Np)) msg nodes

            - indegree of posts equal to 1: $(all(isequal.(k_γ_in, 1)) ? "OK" : "PROBLEM!")
            - sum of in/out degrees of users and posts equal to one another: $(isequal(sum(map(sum, [k_γ_in, k_γ_out])), sum(map(sum, [k_λ_in, k_λ_out]))) ? "OK" : "PROBLEM")

            """
    end

    """
    degreeplotter(G::MetaDiGraph; fname::String=pwd())

    Generate illustration of in- and outdegrees and their distributions.
    illustrations are stored in the `img` subfolder of `path`
    """
    function degreeplotter(G::MetaDiGraph; path::String=pwd())
        # outpath and settings
        isdir(joinpath(path,"img")) ? mkdir(joinpath(path,"img")) : nothing
        outpath = joinpath(path,"img")
        fname = filter!(x->length(x)>0,split(dpath,"/"))[end]

        ## all nodes, showing colors for tagged data
        # indegree
        @debug "plotting indegree"
        Nu = nv(G)
        y = indegree(G,1:Nu)
        c = map(x-> get_prop(G, x ,:flagged) == true ? :red : :green, 1:Nu)
        pin =  scatter(1:Nu ,y, color=c, marker=:circle, label="", size=(2000, 1200))
        title!("projected graph indegree")
        xlabel!("node number")
        ylabel!("indegree")
        savefig(pin, joinpath(outpath,"$(fname) indegree.png"))
        # outdegree

        @debug "plotting outdegree"
        y = outdegree(G,1:Nu)
        c = map(x-> get_prop(G, x ,:flagged) == true ? :red : :green, 1:Nu)
        pout =  scatter(1:Nu, y, color=c, marker=:circle, label="", size=(2000, 1200))
        title!("projected graph outdegree")
        xlabel!("node number")
        ylabel!("outdegree")
        savefig(pout, joinpath(outpath,"$(fname) outdegree.png"))

        ## Degree distributions (global and per type)
        function minimapper(v)
            res = countmap(v)
            x = [p.first for p in res]
            f = [p.second for p in res] / length(v)
            return x,f
        end

        flaggedfilter = map(x-> get_prop(G, x ,:flagged) , 1:get_prop(G,:Nu))

        # indegree
        @debug "plotting indegree distribution"
        ppin = scatter(minimapper(indegree(G))..., label="global", yscale=:log10)
        scatter!(minimapper(indegree(G, findall(flaggedfilter)))..., label="flagged")
        scatter!(minimapper(indegree(G, findall(x-> iszero(x),  flaggedfilter)))..., label="not-flagged")
        title!("indegree distribution")
        xlabel!("indegree")
        ylabel!("frequency")
        savefig(ppin, joinpath(outpath,"$(fname) indegree distribution.png"))

        # outdegree
        @debug "plotting outdegree distribution"
        ppout = scatter(minimapper(outdegree(G)), label="global", yscale=:log10)
        scatter!(minimapper(outdegree(G, findall(flaggedfilter)))..., label="flagged")
        scatter!(minimapper(outdegree(G, findall(x-> iszero(x),  flaggedfilter)))..., label="not-flagged")
        title!("outdegree distribution")
        xlabel!("outdegree")
        ylabel!("frequency")
        savefig(ppout, joinpath(outpath,"img","$(fname) outdegree distribution.png"))

    end

    """
    function that does everything (including LaTeX table)
    """
    function global_loader(globalpath::String)
        @info "starting up"
        @assert isdir(globalpath)
        # buffer to write to
        io = IOBuffer()
        println(io, """\\begin{table*}[]\n\\begin{centering}\n\\begin{array}{lrlrlllr}\n\\hline""")
        println(io, "\\text{Dataset} & N_{\\text{retweets}} & \\%_{\\text{recovered}} & N_{\\text{replies}} & \\%_{\\text{recovered}} & \\%_{\\text{matched edges, retweet}} & \\%_{\\text{matched edges, reply}} \\\\ \\hline")
        # read folders
        for subfolder in filter(x-> isdir(joinpath(globalpath,x)), readdir(globalpath))
            # get name/ref
            @info "working on $subfolder"
            # get relevant files
            msgfiles = filter(x->x[end-3:end] == ".csv" && occursin("tweets_csv_hashed",x) , readdir(joinpath(globalpath, subfolder)))
            usrfiles = filter(x->x[end-3:end] == ".csv" && occursin("users_csv_hashed",x) , readdir(joinpath(globalpath, subfolder)))
            @info "\n\tmsgfiles: $(msgfiles)\n\tuserfiles: $(usrfiles)"

            # read log information
            tot_rt = 0; rec_rt = 0;
            tot_rp = 0; rec_rp = 0;
            for msg in msgfiles
                # read the log
                s = readchomp(joinpath(globalpath, subfolder, "$(msg)_retweets.log")) # for retweets
                c = parse.(Int,match(r"Total: (\d+), collected (\d+)", s).captures)
                tot_rt += c[1]
                rec_rt += c[2]

                s = readchomp(joinpath(globalpath, subfolder, "$(msg)_replies.log")) # for replies
                c = parse.(Int,match(r"Total: (\d+), collected (\d+)", s).captures)
                tot_rp += c[1]
                rec_rp += c[2]
                # extract value
            end

            # build the graphs
            pct_rtedges = 0.
            pct_rpedges = 0.
            gproj = nothing
            try
                pct_rtedges, pct_rpedges, gproj = grapher(joinpath.(globalpath, subfolder ,msgfiles),
                                               joinpath.(globalpath, subfolder, usrfiles), α=0.05)
            catch
                @warn "problem with $(subfolder)"
            end
            @info "$(subfolder), $(round(rec_rt / tot_rt * 100, digits=2)), $(round(rec_rp / tot_rp * 100, digits=2)) "

            println(io, "\\text{$(subfolder)} & $(tot_rt) & $(round(rec_rt / tot_rt * 100, digits=2)) & $(tot_rp) & $(round(rec_rp / tot_rp * 100, digits=2)) & $(pct_rtedges) & $(pct_rpedges) & $(!isnothing(gproj) ? ne(gproj) : 0.0) \\\\")

        end

        println(io, "\\end{array}\n\\par\\end{centering}\n\\end{table*}")
        # write result to file
        open(joinpath(globalpath,"overview.log"), "w") do f
            write(f, String(take!(io)))
        end

    end

end

ext: jl


text:
Zygote.@nograd inner_eltype, inner_ndims, inner_axes, inner_size, inner_length

Zygote.@adjoint function slice(A::AbstractArray, alongs...)
    slice(A, alongs...), Δ -> (align(Δ, alongs...), map(_ -> nothing, alongs)...)
end

Zygote.@adjoint function align(A::NestedArray, alongs...)
    align(A, alongs...), Δ -> (slice(Δ, alongs...), map(_ -> nothing, alongs)...)
end

# flatten(A::NestedArray{V,M,N}) is equivalent to:
# align(A, [True() for _=1:M]..., [False() for _=1:N]...)
Zygote.@adjoint function flatten(A::NestedArray{V,M}) where {V,M}
    flatten(A), Δ -> (slice(Δ, Val(M)), )
end
Zygote.@adjoint flatten(A::AbstractArray) = A, Δ -> (Δ, )


text:
using Documenter, NonlinearOptics

makedocs(format = :html,
         sitename = "NonlinearOptics.jl",
         pages = Any[
             "Home" => "index.md"
         ]
)


text:
using Images
export write_grayscale_image, load_grayscale_image

# saves an image represented as a 2D or 3D array as a grayscale image
# color images are converted to grayscale
function write_grayscale_image(X, filename::String)
  img = convert(Images.Image, X)
  gray = convert(Images.Image{Images.Gray}, img)
  Images.imwrite(gray, filename)
end

# loads an image and converts it to grayscale
# returns an 2d array with each value from 0 (black) to 1 (white)
function load_grayscale_image(filename::String)
  img = Images.imread(filename)
  gray = convert(Images.Image{Images.Gray{Float64}}, img)
  gray_mat = convert(Array, gray)
  gray_arr = convert(Array{Float64}, gray_mat)
  return gray_arr
end


text:
using Distributed
addprocs(6)
@everywhere include("reproducibility.jl")
include("reproducibility.jl")
##

@everywhere using QSWalk
@everywhere using LightGraphs
@everywhere using NPZ
@everywhere using LinearAlgebra, SparseArrays
data_folder_name = "convergence_qsw"

@everywhere function connected_ER(n::Int, p::Float64; is_directed::Bool=false)
    g = erdos_renyi(n, p, is_directed=is_directed)
    while !is_weakly_connected(g)
        g = erdos_renyi(n, p, is_directed=is_directed)
    end
    g
end

function distances(g, targets::Vector{Int})
    [minimum(length.([a_star(g, v, t) for t=targets])) for v=1:nv(g)]
end

## #########################################################################
# ##################### structure observance ###############################
# ##########################################################################

ωs = 0.5:0.05:1.
t_inf = 10_000
repeating = 50
n = 15

println("Path LQSW")
# path_graph LQSW
filename = "$data_folder_name/local_structure_observance_path_$n.npz"
if isfile(filename)
    println("$filename exists -- skipped")
else
    data = zeros(length(ωs), n)
    a = diagm(1 => ones(n-1))
    l = local_lind(a)
    h = diagm(1 => ones(n-1), -1 => ones(n-1))
    init_state = sparse(diagm(fill(1/n, n)))
    for (ωind, ω) = enumerate(ωs)
        s = evolve_generator(h, l, ω)
        data[ωind, :] = abs.(diag(evolve(s, init_state, t_inf)))
    end
    npzwrite(filename, data)
end

# BA LQSW
for m0=[1, 2]
    println("BA $m0 LQSW")
    local filename = "$data_folder_name/local_structure_observance_$n-$m0.npz"
    if isfile(filename)
        println("$filename exists -- skipped")
        continue
    end

    local data = zeros(repeating, length(ωs)+1, n) # +1 for storing distances
    for m=1:repeating
        @show m
        g = barabasi_albert(n, m0, complete=true, is_directed=true)

        l = local_lind(adjacency_matrix(g, dir=:in))
        h = adjacency_matrix(Graph(g))
        init_state = sparse(diagm(fill(1/n, n)))
        for (ωind, ω) = enumerate(ωs)
            s = evolve_generator(h, l, ω)
            data[m, ωind, :] = abs.(diag(evolve(s, init_state, t_inf)))
        end
        data[m, end, :] = distances(g, collect(1:m0))
    end
    npzwrite(filename, data)
end

# path_graph GQSW
println("Path GQSW")
filename = "$data_folder_name/global_structure_observance_path_$n.npz"
if isfile(filename)
    println("$filename exists -- skipped")
else
    data = zeros(length(ωs), n)
    l = diagm(1 => ones(n-1))
    h = diagm(1 => ones(n-1), -1 => ones(n-1))
    init_state = sparse(diagm(fill(1/n, n)))
    for (ωind, ω) = enumerate(ωs)
        s = evolve_generator(h, [l], ω)
        data[ωind, :] = abs.(diag(evolve(s, init_state, t_inf)))
    end
    npzwrite(filename, data)
end

# BA GQSW
for m0=[1, 2]
    println("BA $m0 GQSW")

    local filename = "$data_folder_name/global_structure_observance_$n-$m0.npz"
    if isfile(filename)
        println("$filename exists -- skipped")
        continue
    end

    local data = zeros(repeating, length(ωs)+1, n) # +1 for storing distances
    for m=1:repeating
        @show m
        g = barabasi_albert(n, m0, complete=true, is_directed=true)

        l = [adjacency_matrix(g, dir=:in)]
        h = adjacency_matrix(Graph(g))
        init_state = sparse(diagm(fill(1/n, n)))
        for (ωind, ω) = enumerate(ωs)
            s = evolve_generator(h, l, ω)
            data[m, ωind, :] = abs.(diag(evolve(s, init_state, t_inf)))
        end
        data[m, end, :] = distances(g, collect(1:m0))
    end
    npzwrite(filename, data)
end

# path_graph NGQSW
println("Path NGQSW")
n = 15
filename = "$data_folder_name/nonmoral_structure_observance_path_$n.npz"
if isfile(filename)
    println("$filename exists -- skipped")
else
    local data = zeros(length(ωs), n)
    a = diagm(1 => ones(n-1))
    l, vset = nm_lind(a)
    h = nm_glob_ham(a)
    hrot = nm_loc_ham(vset)
    init_state = nm_init(vlist(vset), vset)
    for (ωind, ω) = enumerate(ωs)
        s = evolve_generator(h, [l], hrot, ω)
        data[ωind, :] = abs.(nm_measurement(evolve(s, init_state, t_inf), vset))
    end
    npzwrite(filename, data)
end

# BA NGQSW
function random_rot(kmax::Int)
    result = Dict{Int,AbstractMatrix{ComplexF64}}()
    for k=1:kmax
        h = rand(k,k)
        result[k] = h + transpose(h)
        h = rand(k,k)
        result[k] += 1im*(h - transpose(h))
    end
    result
end

for m0=[1, 2]
    println("BA $m0 NGQSW")
    local filename = "$data_folder_name/nonmoral_structure_observance_$n-$m0.npz"
    if isfile(filename)
        println("$filename exists -- skipped")
        continue
    end

    local data = zeros(repeating, length(ωs)+1, n) # +1 for storing distances
    for m=1:repeating
        @show m
        g = barabasi_albert(n, m0, complete=true, is_directed=true)
        a = adjacency_matrix(g, dir=:in)
        l, vset = nm_lind(a)
        h = nm_glob_ham(a)
        hrot = nm_loc_ham(vset, random_rot(n))
        init_state = nm_init(vlist(vset), vset)
        data[m, end, :] = distances(g, collect(1:m0))
        for (ωind, ω) = enumerate(ωs)
            s = evolve_generator(h, [l], hrot, ω)
            data[m, ωind, :] =  abs.(nm_measurement(evolve(s, init_state, t_inf), vset))
        end
    end
    npzwrite(filename, data)
end

## #########################################################################
# ############################# NMGQSW probs ###############################
# ##########################################################################
println("Convergence in probs")

n = 12
models = []
push!(models, (n -> barabasi_albert(n, 1, is_directed=true), "ba_1_directed"))
push!(models, (n -> barabasi_albert(n, 1), "ba_1"))
push!(models, (n -> barabasi_albert(n, 3, complete=true), "ba_3"))
push!(models, (n -> barabasi_albert(n, 3, complete=true, is_directed=true), "ba_3_directed"))
push!(models, (n -> connected_ER(n, 0.4, is_directed=true), "er_0.4_directed"))
push!(models, (n -> connected_ER(n, 0.4), "er_0.4"))

t_step = 100
t_count = 100
repeating = 500

@everywhere function random_s(g)
    a = adjacency_matrix(g, dir=:in)
    l, vset = nm_lind(a)
    h = nm_glob_ham(a)
    hrot = nm_loc_ham(vset)
    s = Matrix(evolve_generator(h, [l], hrot))
    s, vset
end

@everywhere function data_gen(m, generator, n, t_count, t_step)
    @show m

    s, vset = random_s(generator(n))
    @time v = eigvals(s)
    list_imaginary = sort(imag.(collect(filter(x->abs(real(x)) < 1e-10 && abs(imag(x)) > 1e-10 , v))), by = x->abs(x))
    while length(list_imaginary) == 0
        s, vset = random_s(generator(n))
        @time v = eigvals(s)
        list_imaginary = sort(imag.(collect(filter(x->abs(real(x)) < 1e-10 && abs(imag(x)) > 1e-10 , v))), by = x->abs(x))
    end
    @show size(s, 1)
    flush(stdout)

    @time begin
        u = evolve_operator(s, t_step)
        rho = nm_init(vlist(vset), vset)
        probs = [nm_measurement(rho, vset)]
        @time for _ = 1:t_count
            rho = evolve(u, rho)
            push!(probs, nm_measurement(rho, vset))
        end
    end
    flush(stdout)
    hcat(probs...)
end

for (generator, gen_label) = reverse(models)
    println("########################### $gen_label ###########################")
    data = pmap(m -> data_gen(m, generator, n, t_count, t_step), 1:repeating)

    npzwrite("$data_folder_name/nonmoral_probs_$gen_label-$n.npz", cat(data..., dims=3))
end


text:
# port of https://github.com/Blosc/c-blosc/blob/3a668dcc9f61ad22b5c0a0ab45fe8dad387277fd/hdf5/blosc_filter.c (copyright 2010 Francesc Alted, license: MIT/expat)

import Blosc

const FILTER_BLOSC = API.H5Z_filter_t(32001) # Filter ID registered with the HDF Group for Blosc
const FILTER_BLOSC_VERSION = 2
const blosc_name = "blosc"

function blosc_set_local(dcpl::API.hid_t, htype::API.hid_t, space::API.hid_t)
    blosc_flags = Ref{Cuint}()
    blosc_values = Vector{Cuint}(undef,8)
    blosc_nelements = Ref{Csize_t}(length(blosc_values))
    blosc_chunkdims = Vector{API.hsize_t}(undef,32)

    API.h5p_get_filter_by_id(dcpl, FILTER_BLOSC, blosc_flags, blosc_nelements, blosc_values, 0, C_NULL, C_NULL)
    flags = blosc_flags[]

    nelements = max(blosc_nelements[], 4) # First 4 slots reserved

    # Set Blosc info in first two slots
    blosc_values[1] = FILTER_BLOSC_VERSION
    blosc_values[2] = Blosc.VERSION_FORMAT

    ndims = API.h5p_get_chunk(dcpl, 32, blosc_chunkdims)
    chunksize = prod(resize!(blosc_chunkdims, ndims))
    if ndims < 0 || ndims > 32 || chunksize > Blosc.MAX_BUFFERSIZE
        return API.herr_t(-1)
    end

    htypesize = API.h5t_get_size(htype)
    if API.h5t_get_class(htype) == API.H5T_ARRAY
        hsuper = API.h5t_get_super(htype)
        basetypesize = API.h5t_get_size(hsuper)
        API.h5t_close(hsuper)
    else
        basetypesize = htypesize
    end

    # Limit large typesizes (they are pretty inefficient to shuffle
    # and, in addition, Blosc does not handle typesizes larger than
    # blocksizes).
    if basetypesize > Blosc.MAX_TYPESIZE
        basetypesize = 1
    end
    blosc_values[3] = basetypesize
    blosc_values[4] = chunksize * htypesize # size of the chunk

    API.h5p_modify_filter(dcpl, FILTER_BLOSC, flags, nelements, blosc_values)

    return API.herr_t(1)
end

function blosc_filter(flags::Cuint, cd_nelmts::Csize_t,
                      cd_values::Ptr{Cuint}, nbytes::Csize_t,
                      buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Cvoid}})
    typesize = unsafe_load(cd_values, 3) # The datatype size
    outbuf_size = unsafe_load(cd_values, 4)
    # Compression level:
    clevel = cd_nelmts >= 5 ? unsafe_load(cd_values, 5) : Cuint(5)
    # Do shuffle:
    doshuffle = cd_nelmts >= 6 ? unsafe_load(cd_values, 6) != 0 : true

    if (flags & API.H5Z_FLAG_REVERSE) == 0 # compressing
        # Allocate an output buffer exactly as long as the input data; if
        # the result is larger, we simply return 0. The filter is flagged
        # as optional, so HDF5 marks the chunk as uncompressed and proceeds.
        outbuf_size = unsafe_load(buf_size)
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return Csize_t(0)

        compname = if cd_nelmts >= 7
            compcode = unsafe_load(cd_values, 7)
            Blosc.compname(compcode)
        else
            "blosclz"
        end
        Blosc.set_compressor(compname)
        status = Blosc.blosc_compress(clevel, doshuffle, typesize, nbytes,
                                      unsafe_load(buf), outbuf, nbytes)
        status < 0 && (Libc.free(outbuf); return Csize_t(0))
    else # decompressing
        # Extract the exact outbuf_size from the buffer header.
        #
        # NOTE: the guess value got from "cd_values" corresponds to the
        # uncompressed chunk size but it should not be used in a general
        # cases since other filters in the pipeline can modify the buffer
        # size.
        outbuf_size, cbytes, blocksize = Blosc.cbuffer_sizes(unsafe_load(buf))
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return Csize_t(0)
        status = Blosc.blosc_decompress(unsafe_load(buf), outbuf, outbuf_size)
        status <= 0 && (Libc.free(outbuf); return Csize_t(0))
    end

    if status != 0
        Libc.free(unsafe_load(buf))
        unsafe_store!(buf, outbuf)
        unsafe_store!(buf_size, outbuf_size)
        return Csize_t(status) # size of compressed/decompressed data
    end
    Libc.free(outbuf)
    return Csize_t(0)
end

# register the Blosc filter function with HDF5
function register_blosc()
    c_blosc_set_local = @cfunction(blosc_set_local, API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    c_blosc_filter = @cfunction(blosc_filter, Csize_t,
                                (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                 Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))
    API.h5z_register(API.H5Z_class_t(API.H5Z_CLASS_T_VERS, FILTER_BLOSC, 1, 1, pointer(blosc_name), C_NULL, c_blosc_set_local, c_blosc_filter))

    return nothing
end

"""
    BloscFilter(;level=5, shuffle=true, compressor="blosclz")

The Blosc compression filter, using [Blosc.jl](https://github.com/JuliaIO/Blosc.jl). Options:

 - `level`: compression level
 - `shuffle`: whether to shuffle data before compressing (this option should be used instead of the [`Shuffle`](@ref) filter)
 - `compressor`: the compression algorithm. Call `Blosc.compressors()` for the available compressors.

# External links
- [What Is Blosc?](https://www.blosc.org/pages/blosc-in-depth/)
"""
struct BloscFilter <: Filter
    blosc_version::Cuint
    version_format::Cuint
    typesize::Cuint
    bufsize::Cuint
    level::Cuint
    shuffle::Cuint
    compcode::Cuint
end

function BloscFilter(;level=5, shuffle=true, compressor="blosclz")
    compcode = Blosc.compcode(compressor)
    BloscFilter(0,0,0,0,level,shuffle,compcode)
end

function Base.show(io::IO, blosc::BloscFilter)
    print(io, BloscFilter,
          "(level=", Int(blosc.level),
          ",shuffle=", blosc.shuffle!=0,
          ",compressor=", Blosc.compname(blosc.compcode),
          ")")
end

filterid(::Type{BloscFilter}) = FILTER_BLOSC
FILTERS[FILTER_BLOSC] = BloscFilter

function Base.push!(f::FilterPipeline, blosc::BloscFilter)
    0 <= blosc.level <= 9 || throw(ArgumentError("blosc compression $(blosc.level) not in [0,9]"))
    ref = Ref(blosc)
    GC.@preserve ref begin
        API.h5p_set_filter(f.plist, filterid(BloscFilter), API.H5Z_FLAG_OPTIONAL, div(sizeof(BloscFilter), sizeof(Cuint)), pointer_from_objref(ref))
    end
    return f
end


text:
module CUBLAS

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

import CUDAnative

using ..CuArrays
using ..CuArrays: active_context, unsafe_free!
using LinearAlgebra

using CEnum

const libcublas = if Sys.iswindows()
    # no ccall by soname, we need the filename
    # NOTE: we discover the full path here, while only the wordsize and toolkit versions
    #       would have been enough to construct "cublas64_10.dll"
    toolkit = find_toolkit()
    path = find_cuda_library("cublas", toolkit)
    if path === nothing
        error("Could not find libcublas")
    end
    const libcublas = basename(path)
else
    # ccall by soname; CuArrays.__init__ will have populated Libdl.DL_LOAD_PATH
    const libcublas = "libcublas"
end

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

const _handles = Dict{CuContext,cublasHandle_t}()
const _xt_handles = Dict{CuContext,cublasXtHandle_t}()
const _handle = Ref{cublasHandle_t}(C_NULL)
const _xt_handle = Ref{cublasXtHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        CUDAnative.maybe_initialize("CUBLAS")
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cublasCreate_v2()

            # enable tensor math mode if our device supports it, and fast math is enabled
            dev = CUDAdrv.device(context)
            if Base.JLOptions().fast_math == 1 && CUDAdrv.capability(dev) >= v"7.0" && version() >= v"9"
                cublasSetMathMode(CUBLAS_TENSOR_OP_MATH, handle)
            end

            atexit(()->CUDAdrv.isvalid(context) && cublasDestroy_v2(handle))
            handle
        end
    end

    return _handle[]
end

function xt_handle()
    if _xt_handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _xt_handle[] = get!(_xt_handles, active_context[]) do
            context = active_context[]
            handle = cublasXtCreate()
            devs = convert.(Cint, CUDAdrv.devices())
            cublasXtDeviceSelect(handle, length(devs), devs)
            atexit(()->CUDAdrv.isvalid(context) && cublasXtDestroy(handle))
            handle
        end
    end
    return _xt_handle[]
end

end


text:
function add_ref_el(el::XMLElement, param::MyXMLElement)
    ref_el = new_element(name(param))
    set_attribute(ref_el, bn.IDREF, get_id(param))
    add_child(el, ref_el)
end

################################################################################
## Parameter
################################################################################
import Base: size

mutable struct Parameter <: MyXMLElement
    el::XMLOrNothing
    val::T where T <: AbstractVector{Float64}
    dim::Int
    lower::Float64
    upper::Float64
    id::String
end

function Parameter(val::Vector{Float64}, id::String; lower::Float64 = NaN)
    return Parameter(nothing, val, length(val), lower, NaN, id)
end

function Parameter(val::Vector{<:Real})
    return Parameter(val, "")
end

function make_parameter(;id::String = "",
        value::AbstractArray{T} = Vector{Int8}(undef, 0),
        lower::Float64 = NaN, upper::Float64 = NaN,
        dim::Int = length(value)) where T <: Real

    el = new_element(bn.PARAMETER)
    if id != ""
        set_attribute(el, bn.ID, id)
    end
    if length(value) != 0
        set_attribute(el, bn.VALUE, join(value, ' '))
    end
    if !isnan(lower)
        set_attribute(el, bn.LOWER, lower)
    end
    if !isnan(upper)
        set_attribute(el, bn.UPPER, upper)
    end
    if dim != length(value)
        set_attribute(el, bn.DIMENSION, dim)
    end

    return el
end

function name(::Parameter)
    return bn.PARAMETER
end

function add_parameter(pel::XMLElement; id::String = "",
                       value::AbstractArray{<:Real} = Float64[],
                       lower::Float64 = NaN, upper::Float64 = NaN,
                       dim::Int = length(value)) where T <: Real

    el = make_parameter(id = id, value = value, lower = lower, upper = upper,
                        dim = dim)
    add_child(pel, el)
    return el
end

function make_xml(param::Parameter)
    el = make_parameter(id = param.id, value = param.val, lower = param.lower,
                        upper = param.upper, dim = param.dim)
    param.el = el
    return el
end

function set_value(p::Parameter, val::AbstractVector{Float64})
    if length(val) != length(p.val)
        error("Incompatible dimensions.")
    end
    p.val .= val
end

import Base: length
function length(p::Parameter)
    return p.dim
end

function get_loggables(p::Parameter)
    make_xml(p)
    return [p.el]
end

function get_hmc_parameter(param::Parameter)
    make_xml(param)
    return param.el
end

################################################################################
## MatrixParameter
################################################################################

mutable struct MatrixParameter <: MyXMLElement
    el::XMLOrNothing
    mat::AbstractMatrix{Float64}
    id::String
    ids::Vector{String}
end

function MatrixParameter(mat::AbstractMatrix{Float64}, id::String, ids::Vector{String})
    m, n = size(mat)
    if length(ids) != m
        error("The matrix has $m rows, but $(length(ids)) row ids were provided.")
    end

    return MatrixParameter(nothing, mat, id, ids)
end

function MatrixParameter(mat::AbstractMatrix{Float64}, id::String)
    m = size(mat, 1)
    return MatrixParameter(mat, id, ["" for i = 1:m])
end

function name(::MatrixParameter)
    return bn.MATRIX_PARAMETER
end

function make_xml(mp::MatrixParameter)
    el = new_element(bn.MATRIX_PARAMETER)
    set_attribute(el, bn.ID, mp.id)
    n = size(mp.mat, 1)
    for i = 1:n
        add_parameter(el, id = mp.ids[i], value = @view mp.mat[i, :])
    end
    mp.el = el
    return el
end

function set_mat!(mp::MatrixParameter, mat::S) where S <: AbstractMatrix{Float64}
    if size(mat) != size(mp.mat)
        error("Supplied matrix has dimensions $(size(mat)), while matrix " *
              "parameter has dimensions $(size(mp.mat))")
    end

    mp.mat = mat
end

function size(p::MatrixParameter)
    return size(p.mat)
end

function size(p::MatrixParameter, dim::Int)
    return size(p.mat, dim)
end

function get_loggables(mp::MatrixParameter)
    make_xml(mp)
    return [mp.el]
end

################################################################################
## DiagonalParameter
################################################################################

mutable struct DiagonalMatrixParameter <: MyXMLElement
    el::XMLOrNothing
    vals::AbstractArray{Float64}
    lower::Float64
    upper::Float64
    id::String
end

function DiagonalMatrixParameter(vals::AbstractVector{Float64}, id::String)
    return DiagonalMatrixParameter(nothing, vals, NaN, NaN, id)
end

function DiagonalMatrixParameter(vals::AbstractVector{Float64})
    return DiagonalMatrixParameter(nothing, vals, NaN, NaN, "")
end

function make_xml(dmp::DiagonalMatrixParameter)
    el = new_element(bn.DIAGONAL_MATRIX)
    if dmp.id != ""
        set_id!(el, dmp.id)
    end
    p_el = make_parameter(value = dmp.vals, lower=dmp.lower, upper=dmp.upper)
    add_child(el, p_el)
    dmp.el = el
    return el
end

function size(d::DiagonalMatrixParameter)
    p = length(d.vals)
    return (p, p)
end

function size(d::DiagonalMatrixParameter, dim::Int)
    p = length(d.vals)
    return p
end

################################################################################
## Multiplicative parameter
################################################################################

mutable struct MultiplicativeParameter <: MyXMLElement
    el::XMLOrNothing
    param::Parameter
    id::String

    function MultiplicativeParameter(param::Parameter, id::String)
        return new(nothing, param, id)
    end
end

function make_xml(mp::MultiplicativeParameter)
    el = new_element(bn.MULTIPLICATIVE_PARAMETER)
    set_id!(el, mp.id)

    make_xml(mp.param)
    add_ref_el(el, mp.param.el)

    mp.el = el
    return el
end

function name(::MultiplicativeParameter)
    return bn.MULTIPLICATIVE_PARAMETER
end

################################################################################
## Masked parameter
################################################################################

mutable struct MaskedParameter <: MyXMLElement
    el::XMLOrNothing
    param::MyXMLElement
    mask::Vector{Float64}

    function MaskedParameter(param::MyXMLElement, mask::Vector{Float64})
        return new(nothing, param, mask)
    end
end

function make_xml(mp::MaskedParameter)
    el = new_element(bn.MASKED_PARAMETER)
    add_ref_el(el, mp.param)
    m_el = new_child(el, bn.MASK)
    add_child(m_el, make_xml(Parameter(mp.mask)))
    mp.el = el
    return el
end

################################################################################
## compound parameter
################################################################################

mutable struct CompoundParameter <: MyXMLElement
    el::XMLOrNothing
    parameters::Vector{<:MyXMLElement}
    id::String

    function CompoundParameter(parameters::Vector{<:MyXMLElement}, id::String)
        return new(nothing, parameters, id)
    end
end

function name(::CompoundParameter)
    return bn.COMPOUND_PARAMETER
end

function make_xml(cp::CompoundParameter)
    el = new_element(name(cp))
    for param in cp.parameters
        p_el = make_xml(param)
        add_child(el, p_el)
    end

    set_attribute(el, bn.ID, cp.id)
    cp.el = el
    return el
end

function get_loggables(cp::CompoundParameter)
    return make_xml(cp)
end


text:
abstract type AbstractArrayStream{S,T,N,I} <: AbstractArray{T,N} end

# IO interactions
stream(s::AbstractArrayStream) = error("stream is not defined for $(typeof(s))")
Base.seek(s::AbstractArrayeStream, n::Integer) = seek(stream(s), n)
Base.position(s::AbstractArrayStream) = position(stream(s))
Base.skip(s::AbstractArrayStream, n::Integer) = skip(stream(s), n)
Base.close(s::AbstractArrayStream) = close(stream(s))
Base.eof(s::AbstractArrayStream) = eof(stream(s))
Base.isopen(s::AbstractArrayStream) = isopen(stream(s))

# Field accessors
ownstream(s::AbstractArrayStream) = getfield(s, :ownstream)
streamindices(s::AbstractArrayStream) = getfield(s, :streamindices)
needswap(s::AbstractArrayStream) = getfield(s :needswap)

Base.size(s::AbstractArrayStream{S,T,N}) where {S,T,N} = (S.parameters...,)
Base.size(s::AbstractArrayStream{S,T,N}, i::Int) where {S,T,N} = S.parameters[i]
Base.ndims(s::AbstractArrayStream{S,T,N}) where {S,T,N} = N
Base.length(s::AbstractArrayStream{S,T,N}) where {S,T,N} = prod(S.parameters)
Base.eltype(s::AbstractArrayStream{S,T,N}) where {S,T,N} = T

Base.IndexStyle(::Type{AbstractArrayStream}) = IndexLinear

Base.firstindex(s::AbstractArrayStream) = streamindices(s)[1]
Base.lastindex(s::AbstractArrayStream) = streamindices(s)[end]

"""
    getindex(AbstractArrayStream, I)

While an `ArrayStream` may be indexed, this will only result in performance gains
if the indices represent a significantly smaller portion of the original `ArrayStream`.
This is because an indexed ArrayStream has to perform a `seek(ArrayStream, index)`
step at every index before reading. This is in contrast to reading in a single
chunk of memory, which is considerably faster.
"""
Base.getindex(s::AbstractArrayStream{S,T,N,IndexLinear}, i::Int) where {S,T,N} =
    ArrayStreams(stream(s), streamindices(s)[i], ownstream(s), needswap(s))

function Base.getindex(s::AbstractArrayStream{S,T,N,L}, I::Vararg{Int, N}) where {S,T,N,L}
    si = s.streamindices[LinearIndices(size(s))[I...]]
    ArrayStreams{Tuple{1},T,1,1}(stream(s), s.streamindices[LinearIndices(size(s))[I...]], ownstream(s), needswap(s))
end

# leverage SubArrays which determine whether fast indexing is possible
function Base.getindex(s::AbstractArrayStream{S,T,N,IndexLinear}, I...) where {S,T,N}
    li = view(LinearIndices(s), I...)
    # fast linear indexing, so just need first and last
    if typeof(li).parameters[5]
        ArrayStream{Tuple{size(li)...},T,ndims(li)}(stream(s), streamindices(s)[first(li):last(li)], ownstream(s), needswap(s))
    else
        ArrayStream{Tuple{size(li)...},T,ndims(li)}(stream(s), streamindices(s)[li...], ownstream(s), needswap(s))
    end
end

Base.setindex!(s::AbstractArrayStream, I) = error("setindex is not defined for $(typeof(s)). Use write to change stream content.")

function Base.show(io::IO, s::AbstractArrayStream)
    println(summary(s))
    println("ownstream: ", ownstream(s))
    println("needswap: ", needswap(s))
    println("streamindices: \n", streamindices(s))
end


text:
module JLScriptManager

using CoreExtended
using FileManager.FileSource

export JLComponent

""" TODO """
abstract type JLComponent <: AbstractObjectReference end

""" TODO """
type JLInvalidComponent <: JLComponent end

""" TODO """
const JL_INVALID_COMPONENT = JLInvalidComponent()

""" TODO """
type JLStateListComponent <: JLComponent
	isInitalized	::Bool
	isRunning			::Bool
	isTerminated	::Bool

	JLStateListComponent() = new(false,false,false)
end

export JLScriptFunction

""" TODO """
type JLScriptFunction
	func::Function
	JLScriptFunction(f::Function) = new(stabilize(f))
end

export JLScript

""" TODO """
type JLScript
  id        ::Symbol
	mod				::Module
  source    ::FileSource

	args			::Tuple # arguments passed on main()
  cache			::Dict{Symbol, Any}
  listener	::Dict{Symbol, Function}
	extern		::Dict{Symbol, Function}
	state			::JLComponent
	objref		::JLComponent
	funcs			::JLComponent
	events		::JLComponent
end

""" TODO """
JLScript(id::Symbol) = JLScript(id,FileSource())

""" TODO """
function JLScript(id::Symbol, source::FileSource)
  this=JLScript(id,Module(id),source,(),Dict(),Dict(),Dict(),JLStateListComponent(),JL_INVALID_COMPONENT,JL_INVALID_COMPONENT,JL_INVALID_COMPONENT)
  JLSCRIPTS[id]=this
  this
end

JLSCRIPTS = Dict{Symbol,JLScript}()

""" TODO """
loop(f::Function) = for (k,s) in JLSCRIPTS f(s) end

""" TODO """
listen(this::JLScript, k::Symbol, f::Function) = (this.listener[k]=f)

""" TODO """
function run(this::JLScript, args...)
	debug("run $(this.id)")
	#Module(:__anon__)
	result = execute(this)
	if !result[1] result = nothing
	else result = result[2]
	end
  result
end

""" TODO """
(this::JLScript)(s::Symbol, args...) = CoreExtended.execute(this.mod,s,args...) #@eval $f($args...)

""" TODO """
exists(this::JLScript, s::Symbol) = CoreExtended.exists(this.mod,s)

""" TODO """
function execute(this::JLScript, compile_args=[], args...)
  debug("execute $(this.id)")
	result = compile(this, compile_args...)
	if result[1] result = (true, execute(result[2], args...)) end
	result
end

""" TODO """
execute(f::JLScriptFunction, args...) = (debug("execute function("*string(args...)*")"); f.func(args...))

""" TODO """
function compile(this::JLScript, args...)
	debug("compile $(this.id)")
	result = (false, nothing)

	catchException(function()
		#if length(args)>0
			result = eval(this, args...)
		#else
		#	result = evalfile(this.source.path)
		#end

		if isa(result, Function) result = JLScriptFunction(result)
		else warn("Main Function not found.")
		end

		result = (true, result)
	end)

	result
end

""" TODO """
function cleanCode(code::String)

	# problem using match(): only "function (name)" will be detected!
	# for x in eachmatch(r"function \s(\w)", code) println(x.captures[1]) end

	# TODO merge code lines
	code=Base.replace(code, "\r", "")
	#code=Base.replace(code, r"\#\=(?(?=\=\#)then|else)*\=\#", "")
	code=Base.replace(code, r"(?s)(?<=\#\=)(.*?)(?=\=\#)", "")
	code=Base.replace(code, r"(\#[^\n]+)", "")
	code=Base.replace(code, r"\n\s+\n", "\n")
	code=Base.replace(code, r"\n+", "\n")
	code=Base.replace(code, r"\s*(,)\s*", s"\1")
	code=Base.replace(code, r"(\()\s*", s"\1")
	code=Base.replace(code, r"\s*(\))", s"\1")
	code=Base.replace(code, r"\n", ";")
end

""" TODO """
function eval(this::JLScript, args...)

	this.mod = Module(this.id) #reset modul

	debug("eval $(this.id)")

	# parse
	#################################################

	code = ""
	open(this.source.path) do f	code = readstring(f) end
	code=cleanCode(code)
	code=parse(code)

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

	funcs = ""
	funcsext = "" #"_(x) = (args...)->execute(x,args...);"

	xTypBody = ""
	xTypCall = ""
	fTypBody = ""
	fTypCall = ""
	eTypBody = ""
	eTypCall = ""

	def = (name) -> "$name::Function;"
	dec = (name) -> "(args...)->execute($name,args...)"

	for (name,f) in this.extern
		funcsext *= "global $name = this.extern[:$name];"
		#xTypBody *= def(name)
		#xTypCall *= (xTypCall != "" ? "," : "") * dec(name)
	end

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

	for x in code.args
		if x.head == Symbol("function")
			name=Symbol(Base.replace(string(x.args[1]), r"\(.*", ""))

			if ismatch(r"^On\w+", string(name))
				# event functions
				eTypBody *= def(name)
				eTypCall *= (eTypCall != "" ? "," : "") * dec(name)
			else
				# other functions
				fTypBody *= def(name)
				fTypCall *= (fTypCall != "" ? "," : "") * dec(name)
			end

		end
	end

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

	imports = ""
	imports *= "using CoreExtended;"
	imports *= "using JLScriptManager;"

	typs = ""
	typ = string(this.mod) #*"_JLFunctionListComponent"
	#if createExtFuncList typs *= "type "*typ*"_EXTFUNCTIONS <: JLComponent;"*xTypBody*";end;" end
	typs *= "type "*typ*"_FUNCTIONS <: JLComponent;"*fTypBody*";end;"
	typs *= "type "*typ*"_EVENTS <: JLComponent;"*eTypBody*";end;"

	mainFunc = "function main(args...);"
	mainFunc *= "this.args=args;"
	mainFunc *= funcsext
	#if createExtFuncList mainFunc *= "this.extfuncs = "*typ*"_EXTFUNCTIONS("*xTypCall*");" end
	mainFunc *= "this.funcs = "*typ*"_FUNCTIONS("*fTypCall*");"
	mainFunc *= "this.events = "*typ*"_EVENTS("*eTypCall*");"
	mainFunc *= "end;"

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

	ex=Expr(:toplevel,
		:(const ARGS = $(this.args)),
		:(eval(x) = Main.Core.eval($(this.mod),x)),
		#:(eval(m,x) = Main.Core.eval(m,x)),
		:(this=$this),
		args...,
		:(Main.Base.include($(this.source.path))), #which is faster?
		#code, #which is faster?
		parse(imports*typs*funcs*mainFunc),
	)

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

	#println("EXTFUNCS:\n", Base.replace(extfuncs, ",(", "\n("),"\n")
	#println("FUNCS: ", fTypBody)
	#println("EVENTS: ", eTypBody)
	#dump(ex)

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

	eval(this.mod, ex)
end

end # ScriptManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
module TestAveraging

using DataTools
using Test
using Transducers: Filter, Map

reduce_bs1(args...; kw...) = reduce(args...; basesize = 1, kw...)

@testset for fold in [foldl, reduce_bs1, reduce]
    @test fold(averaging, Filter(isodd), 1:10) == 5
    @test fold(
        oncol(a = averaging, b = averaging),
        Map(identity),
        [(a = 1, b = 2), (a = 2, b = 3)],
    ) == (a = 1.5, b = 2.5)
end

end  # module


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
# Random.seed!(233)

function spatial_fem_stiffness(hmat,m,n,h)
    spatial_fem_stiffness_ = load_op_and_grad("./build/libSpatialFemStiffness","spatial_fem_stiffness", multiple=true)
    hmat,m,n,h = convert_to_tensor([hmat,m,n,h], [Float64,Int32,Int32,Float64])
    spatial_fem_stiffness_(hmat,m,n,h)
end

# TODO: specify your input parameters
m = n = 20
h = 0.1
D = diagm(0=>[1,1,0.5])
Ds = zeros(4*m*n, 3, 3)
for i = 1:m
    for j = 1:n
        Ds[4*((j-1)*m+i-1)+1,:,:] = D
        Ds[4*((j-1)*m+i-1)+2,:,:] = D
        Ds[4*((j-1)*m+i-1)+3,:,:] = D
        Ds[4*((j-1)*m+i-1)+4,:,:] = D
    end
end
u = spatial_fem_stiffness(Ds,m,n,h)
sess = Session(); init(sess)
run(sess, u);

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m_)
    return sum(spatial_fem_stiffness(m_,m,n,h)[3]^4)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(size(Ds)...))
v_ = rand(4*m*n,3,3)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1e-4 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
# struct ParentsBeforeChildrenIterator
#     ops::Vector{Operation}
#     visited::Vector{Bool}
# end
# function iterate(pbci::ParentsBeforeChildrenIterator)
#     for (i,op) ∈ enumerate(pbci.ops)
#         if iszero(length(parents(op)))
#             pbci.visited[i] = true
#             return op, pbci
#         end
#     end
#     nothing
# end
# function iterate()

# end

function dependent_outer_reducts(ls::LoopSet, op)
    for i ∈ ls.outer_reductions
        search_tree(parents(operations(ls)[i]), name(op)) && return true
    end
    false
end

function isnopidentity(ls::LoopSet, op::Operation, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, u₂max::Int)
    parents_op = parents(op)
    if iscompute(op) && instruction(op).instr === :identity && isone(length(parents_op)) && name(first(parents_op)) === name(op)
        loopistiled = u₂max ≠ -1
        # mvar, u₁unrolledsym, u₂unrolledsym = variable_name_and_unrolled(op, u₁loop, u₂loop, vectorized, Core.ifelse(isu₂unrolled(op), u₂max, -1))
        # parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loop, u₂loop, u₂max)
        # @show  (u₁unrolledsym, first(parents_u₁syms)), (isu₂unrolled(op), parents_u₂syms[1])
        # @show op parents(op) isu₁unrolled(op), isu₁unrolled(only(parents(op)))
        # if (u₁unrolledsym == first(parents_u₁syms)) && (isu₂unrolled(op) == parents_u₂syms[1])
        opp = only(parents_op)
        if (isu₁unrolled(op) == isu₁unrolled(opp)) & (isu₂unrolled(op) == isu₂unrolled(opp))
            #TODO: identifer(first(parents_op)) ∉ ls.outer_reductions is going to miss a lot of cases
            #Should probably replace that with `DVec` (demoting Vec) types, that demote to scalar.
            #TODO: document (after finding out...) why only checking `isvectorized(first(parents_op))` -- why not `any(isvectorized, parents_op)`???
            if (isvectorized(opp) && !isvectorized(op)) && !dependent_outer_reducts(ls, op)
                op.instruction = reduction_to_scalar(instruction(opp))
                op.mangledvariable = gensym(op.mangledvariable)
                false
            else
                true
            end
        else
            false
        end
    else
        false
    end
end

function set_upstream_family!(adal::Vector{T}, op::Operation, val::T, ld::Vector{Symbol}, id::Int) where {T}
    adal[identifier(op)] == val && return # must already have been set
    if ld != loopdependencies(op) || id == identifier(op)
        (adal[identifier(op)] = val)
    end
    for opp ∈ parents(op)
        identifier(opp) == identifier(op) && continue
        set_upstream_family!(adal, opp, val, ld, id)
    end
end

function addoptoorder!(
    ls::LoopSet, included_vars::Vector{Bool}, place_after_loop::Vector{Bool}, op::Operation,
    loopsym::Symbol, _n::Int, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, u₂max::Int
)
    lo = ls.loop_order
    id = identifier(op)
    included_vars[id] && return nothing
    loopsym ∈ loopdependencies(op) || return nothing
    for opp ∈ parents(op) # ensure parents are added first
        addoptoorder!(ls, included_vars, place_after_loop, opp, loopsym, _n, u₁loop, u₂loop, vectorized, u₂max)
    end
    included_vars[id] && return nothing
    included_vars[id] = true
    isunrolled = (isu₁unrolled(op)) + 1
    istiled = isu₂unrolled(op) + 1
    # optype = Int(op.node_type) + 1
    after_loop = place_after_loop[id] + 1
    if !isloopvalue(op)
        isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max) || push!(lo[isunrolled,istiled,after_loop,_n], op)
        # if istiled
        #     isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max, u₂max) || push!(lo[isunrolled,2,after_loop,_n], op)
        # else
        #     isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max, nothing) || push!(lo[isunrolled,1,after_loop,_n], op)
        # end
    end
    # @show op, after_loop
    # isloopvalue(op) || push!(lo[isunrolled,istiled,after_loop,_n], op)
    # all(opp -> iszero(length(reduceddependencies(opp))), parents(op)) &&
    set_upstream_family!(place_after_loop, op, false, loopdependencies(op), identifier(op)) # parents that have already been included are not moved, so no need to check included_vars to filter
    nothing
end

function fillorder!(ls::LoopSet, order::Vector{Symbol}, u₁loop::Symbol, u₂loop::Symbol, u₂max::Int, vectorized::Symbol)
    lo = ls.loop_order
    resize!(lo, length(ls.loopsymbols))
    ro = lo.loopnames # reverse order; will have same order as lo
    nloops = length(order)
    ops = operations(ls)
    nops = length(ops)
    included_vars = fill!(resize!(ls.included_vars, nops), false)
    place_after_loop = fill!(resize!(ls.place_after_loop, nops), true)
    # to go inside out, we just have to include all those not-yet included depending on the current sym
    empty!(lo)
    for _n ∈ 1:nloops
        n = 1 + nloops - _n
        ro[_n] = loopsym = order[n]
        #loopsym = order[n]
        for op ∈ ops
            addoptoorder!( ls, included_vars, place_after_loop, op, loopsym, _n, u₁loop, u₂loop, vectorized, u₂max )
        end
    end
end


text:
"""
`KernelSum(kernels::Array{Kernel};weights::Array{Real}=ones(length(kernels)))`
Create a positive weighted sum of kernels.
One can also use the operator `+`
```
k1 = SqExponentialKernel()
k2 = LinearKernel()
k = KernelSum([k1,k2])
kernelmatrix(k,X) == kernelmatrix(k1,X).+kernelmatrix(k2,X)
kernelmatrix(k,X) == kernelmatrix(k1+k2,X)
kweighted = 0.5*k1 + 2.0*k2
```
"""
struct KernelSum{T,Tr} <: Kernel{T,Tr}
    kernels::Vector{Kernel}
    weights::Vector{Real}
    function KernelSum{T,Tr}(kernels::AbstractVector{<:Kernel},weights::AbstractVector{<:Real}) where {T,Tr}
        new{T,Tr}(kernels,weights)
    end
end

function KernelSum(kernels::AbstractVector{<:Kernel}; weights::AbstractVector{<:Real}=ones(Float64,length(kernels)))
    @assert length(kernels)==length(weights) "Weights and kernel vector should be of the same length"
    @assert all(weights.>=0) "All weights should be positive"
    KernelSum{eltype(kernels),Transform}(kernels,weights)
end

params(k::KernelSum) = (k.weights,params.(k.kernels))
opt_params(k::KernelSum) = (k.weights,opt_params.(k.kernels))
duplicate(k::KernelSum,θ) = KernelSum(duplicate.(k.kernels,θ[end]),weights=first(θ))

Base.:+(k1::Kernel,k2::Kernel) = KernelSum([k1,k2],weights=[1.0,1.0])
Base.:+(k1::KernelSum,k2::KernelSum) = KernelSum(vcat(k1.kernels,k2.kernels),weights=vcat(k1.weights,k2.weights))
Base.:+(k::Kernel,ks::KernelSum) = KernelSum(vcat(k,ks.kernels),weights=vcat(1.0,ks.weights))
Base.:+(ks::KernelSum,k::Kernel) = KernelSum(vcat(ks.kernels,k),weights=vcat(ks.weights,1.0))
Base.:*(w::Real,k::Kernel) = KernelSum([k],weights=[w]) #TODO add tests
Base.:*(w::Real,k::KernelSum) = KernelSum(k.kernels,weights=w*k.weights) #TODO add tests

Base.length(k::KernelSum) = length(k.kernels)
metric(k::KernelSum) = metric.(k.kernels)
transform(k::KernelSum) = transform.(k.kernels)
transform(k::KernelSum,x::AbstractVecOrMat) = transform.(k.kernels,[x])
transform(k::KernelSum,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim)

function kernelmatrix(
    κ::KernelSum,
    X::AbstractMatrix;
    obsdim::Int=defaultobs)
    sum(κ.weights[i]*kernelmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
end

function kernelmatrix(
    κ::KernelSum,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int=defaultobs)
    sum(κ.weights[i]*_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
end

function kerneldiagmatrix(
    κ::KernelSum,
    X::AbstractMatrix;
    obsdim::Int=defaultobs)
    sum(κ.weights[i]*kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
end


text:
"""
    ContinuousSystemModel{T, N}

An abstract type that represents a set of model parameters of type T for a continous dynamical system,
consisting of `N` variables.

As a minimum, such types always include a field `dt` (giving the sampling time step, which is the
`dt` argument to `DynamicalSystems.trajectory`) and a field `ui` (giving the initial condition). Hence,
define a continous model system by (1) its parameters, (2) its initial condition, (3) its sampling rate.

Concrete subtypes must implement the following methods:

- **[`ContinuousDynamicalSystem(x::ContinuousSystemModel)`](@ref)**:
    Create a `ContinuousDynamicalSystem` from the model.
- **[`get_dt`](@ref)**: Returns the time step.
- **[`get_ui`](@ref)**: Returns the initial condition.
- **[`get_nvars`](@ref)**: Returns the number of variables.
- **[`interaction_matrix`](@ref)**: Returns the interaction matrix for the model (respecting the
    coupling parameter for that particular model).

Concrete subtypes *may* implement the following methods:

- **[`rand`](@ref)**. If implemented for a continuous system model of type `SM`,
    returns an instance of the model with randomised parameters.
"""
abstract type ContinuousSystemModel{T, N} <: AbstractSystemModel{T, N} end

"""
    ContinuousDynamicalSystem(x::ContinuousSystemModel) -> ContinuousDynamicalSystem

Convert a system model to a `ContinuousDynamicalSystem` by connecting
its equations of motion to it parameters.
"""
ContinuousDynamicalSystem(x::ContinuousSystemModel)


text:
"""
    VarDims

A type for calculating conversions between dimensions of the dynamical
system and the augmented system (continuation problem).

| Bifurcation | Variables  | Equations  | `eigvec_dim` | `eigvec_eltype` |
| ----------- | ---------- | ---------- | ------------ | --------------- |
| Saddle-node | ``2n + 2`` | ``2n + 1`` | ``n``        | `<: Real`       |
| Hopf        | ``3n + 3`` | ``3n + 2`` | ``2n``       | `<: Complex`    |

where ``n`` is `ds_dim`, the number of variables is `dim(dom(H))`, and
the number of equations is `dim(codom(H))`.

# Fields
- `ds_dim::Int`: Dimension of the dynamical system.
- `eigvec_dim::Int`: Number of reals required to represent a eigenvector.
- `eigvec_eltype::Type`
"""
struct VarDims
    ds_dim::Int
    eigvec_dim::Int
    eigvec_eltype::Type
end

"""
    augsys_dim(d::VarDims) :: Int

Number of variables `dim(dom(H))` of the augmented system.
"""
augsys_dim(d::VarDims) =
    d.ds_dim + d.eigvec_dim + 2 + (d.eigvec_eltype <: Complex)
eigvec_range(d::VarDims) = d.ds_dim .+ (1:d.eigvec_dim)
eigval_index(d::VarDims) = d.ds_dim + d.eigvec_dim + 1

function dims_from_augsys(as_dim::Int, T::Type)
    if T <: Real
        ds_dim, r = divrem(as_dim - 2, 2)
        eigvec_dim = ds_dim
    else
        @assert T <: Complex
        ds_dim, r = divrem(as_dim - 3, 3)
        eigvec_dim = 2 * ds_dim
    end
    return _VarDims(ds_dim, eigvec_dim, T, as_dim, r)
end

function _VarDims(ds_dim, eigvec_dim, T, as_dim, r)
    @assert r == 0
    d = VarDims(ds_dim, eigvec_dim, T)
    @assert augsys_dim(d) == as_dim
    return d
end

function dims_from_augsys(as_dom_dim::Int, ::SaddleNodeCont)
    ds_dim, r = divrem(as_dom_dim - 2, 2)
    eigvec_dim = ds_dim
    eigvec_eltype = Real  # TODO: don't
    return _VarDims(ds_dim, eigvec_dim, eigvec_eltype, as_dom_dim, r)
end

function dims_from_augsys(as_dom_dim::Int, ::HopfCont)
    ds_dim, r = divrem(as_dom_dim - 3, 3)
    eigvec_dim = 2 * ds_dim
    eigvec_eltype = Complex  # TODO: don't
    return _VarDims(ds_dim, eigvec_dim, eigvec_eltype, as_dom_dim, r)
end

dims_from_augsys(J::AbstractMatrix, ckind::ContinuationKind) =
    dims_from_augsys(size(J, 2), ckind)

# ds_dim(x::AbstractVector) = ds_dim(length(x))
# ds_dim(var_dim) = var_dim ÷ 2 - 1
eq_dim(N) = 2N + 1
var_dim(N) = 2N + 2

cat_outputs(H1, H2, H3) = vcat(H1, H2, H3)
cat_outputs(H1::SVector, H2::SVector, H3::Number) :: SVector =
    vcat(H1, H2, SVector(H3))

function _as_reals(S, T)
    values = [:($part(v[$i])) for i in 1:S for part in [:real, :imag]]
    quote
        SVector{$(2S), $T}($(values...))
    end
end

as_reals(v::AbstractVector{<: Real}) = v
as_reals(v::Vector{<: Complex{T}}) where T = reinterpret(T, v)

# manually doing this since reinterpret does not work with SVector:
@generated as_reals(v::SVector{S, Complex{T}}) where {S, T} =
    _as_reals(S, T)

as_reals(s::Real) = SVector(s)
as_reals(s::Complex) = SVector(real(s), imag(s))