text:
# These functions represent combinations that are covered and uncovered.
# It is a data structure to facilitate greedy approaches to combination coverage.
using Combinatorics: combinations
import Base: eltype

###############################################################
# MatrixCoverage
###############################################################

"""
Represents covered tuples using a two-dimensional matrix.
Each column is another tuple to cover. Each row is a parameter.
A zero means this parameter is not part of the tuple.
The initial matrix represents all tuples to cover.
The `remain` integer is the number of tuples left to cover.
As tuples are covered, they are swapped to the end of the
uncovered tuples, and the `remain` value is decremented.

The arity is an integer array, the same length as the number of
parameters. Each member of the array is the number of possible
values each parameter can take.

This representation should be slow for lookup of tuples and
take lots of memory, but it is a clear representation with which
to understand the interface to the data structure.
"""
mutable struct MatrixCoverage{T <: Integer}
    allc::Array{T, 2}
    remain::Integer
    arity::Array{T, 1}
end

"""
The element type of the MatrixCoverage class is determined
by the element type of the arity because an element large
enough to hold the arity is what defines the minimum possible
element type.
"""
eltype(mc::MatrixCoverage) = eltype(mc.arity)

"""
all_combinations_matrix(arity, n_way)

Create matrix coverage that includes all n-way combinations
of parameters which can take `arity` values, where `arity`
is a vector of integers to represent the number of possible
values for each parameter, in order.
"""
function all_combinations_matrix(arity, n_way)
    allc = all_combinations(arity, n_way)
    # Every combination is nonzero.
    @assert sum(sum(allc, dims = 1) == 0) == 0
    remain = size(allc, 2)
    MatrixCoverage(allc, remain, arity)
end

"""
    one_parameter_combinations(arity, n_way)

Generates all combinations that are nonzero for the last parameter.
This is for in-parameter-order generation, where we need
only those tuples that end with this column being nonzero.

The construction method is to leave out the given parameter
and construct all `n_way` - 1 tuples. Then copy and paste that
once for each possible value of the given parameter.
"""
function one_parameter_combinations_matrix(arity, n_way)
    comb = one_parameter_combinations(arity, n_way)
    MatrixCoverage(comb, size(comb, 2), arity)
end

"""
How many tuples have not been covered yet.
"""
function remaining_uncovered(mc::MatrixCoverage)
    mc.remain
end

"""
The number of parameters for this design.
"""
function parameter_cnt(mc::MatrixCoverage)
    length(mc.arity)
end

function combination_histogram(allc, arity)
    height = maximum(arity)
    hist = zeros(Int, height, length(arity))
    for colh_idx in 1:size(allc, 2)
        for rowh_idx in 1:height
            if allc[rowh_idx, colh_idx] > 0
                hist[allc[rowh_idx, colh_idx], rowh_idx] += 1
            end
        end
    end
    hist
end

function most_to_cover(allc, col_cnt)
   argmax(vec(sum(allc[:, 1:col_cnt] .!= 0, dims = 2)))
end

"""
    coverage_by_parameter(allc::MatrixCoverage)

Given a coverage matrix, return how many times each parameter
participates in an uncovered tuple. The return value is a vector
of integers, where each value is the number of times that parameter
appears in any uncovered tuple.
"""
function coverage_by_parameter(mc::MatrixCoverage)
    vec(sum(mc.allc[:, 1:mc.remain] .!= 0, dims = 2))
end

"""
    coverage_by_value(allc, param_idx)

Given a particular parameter, look at that column of the coverage matrix
and return a histogram of how many times each value appears in
that column. We want to know which parameter has the most tuples
to cover so that we can address it first.
"""
function coverage_by_value(mc, param_idx)
    hist = zeros(Int, mc.arity[param_idx] + 1)
    for col_idx in 1:mc.remain
        hist[mc.allc[param_idx, col_idx] + 1] += 1
    end
    hist[2:end]
end

function most_common_value(allc, col_cnt, arity, param_idx)
    hist = zeros(Int, arity[param_idx] + 1)
    for col_idx in 1:col_cnt
        hist[allc[param_idx, col_idx] + 1] += 1
    end
    argmax(hist[2:end])
end

"""
    most_matches_existing(mc::MatrixCoverage, existing, param_idx)

The `existing` vector is a partial test case. Each nonzero entry in this
test case is considered as decided. This function then determines, for
the next parameter, chosen by param_idx, which value of param_idx, between
1 and its arity, would cover the most tuples. It returns a histogram of
how many uncovered tuples could be covered, given the existing choices
and a particular value of this parameter.

This returns all tuples that the existing case could _exactly_ cover
if it didn't have missing values. That means the existing vector
either has the same number of entries as the tuple or it has
fewer entries than the tuple.

A value of the `param_idx` parameter matches with a tuple if its
existing choices match at least some part of the tuple and don't
disagree with any part of the tuple.
"""
function most_matches_existing(mc::MatrixCoverage, existing, param_idx)
    @assert existing[param_idx] == 0
    param_cnt = length(mc.arity)
    # The match count can only be as large as n_way - 1 because the
    # column with the parameter is non-zero and will be zero in existing.
    # If n-way is 3 and we only know one parameter so far, that's another limit
    # to the possible match size.
    max_known = sum(existing .!= 0)
    # params_known = min(sum(existing .!= 0), n_way - 1)
    hist = zeros(Int, mc.arity[param_idx])
    for col_idx in 1:mc.remain
        # The given parameter column is part of this match.
        if mc.allc[param_idx, col_idx] != 0
            match_cnt = 0
            n_way = 0
            for match_idx in 1:param_cnt
                if mc.allc[match_idx, col_idx] != 0
                    n_way += 1
                end
                if existing[match_idx] != 0 && mc.allc[match_idx, col_idx] == existing[match_idx]
                    match_cnt += 1
                end
            end
            if match_cnt == min(max_known, n_way)
                hist[mc.allc[param_idx, col_idx]] += 1
            end
        end
    end
    hist
end

"""
The logic of tuple comparison is oddly complicated.
We're going to write this out the first round in order to know
what our tools are.

If you take a case and a tuple to cover, then there are
five states for any pair of values:
               case tuple
    ignores    0    0
    skips      a    0
    misses     0    b
    matches    a == b
    mismatch   a != b

It's always one of those five. In this language,

covers = all(matches | irrelevant)
"""

ignores(a, b) = a == 0 && b == 0
skips(a, b) = a != 0 && b == 0
misses(a, b) = a == 0 && b != 0
matches(a, b) = a != 0 && b != 0 && a == b
mismatch(a, b) = a != 0 && b != 0 && a != b

"""
Example matches.
      crossed      incomplete cover
case  [0 1 0 2]    [1 0 0 3]  [1 1 0 2]
tuple [1 0 3 0]    [1 1 0 3]  [1 0 0 2]
"""
function case_compatible_with_tuple(case, tuple)
    !any(mismatch(a, b) for (a, b) in zip(case, tuple))
end

"""
Example matches.
      incomplete cover
case  [1 0 0 3]  [1 1 0 2]
tuple [1 1 0 3]  [1 0 0 2]
"""
function case_partial_cover(case, tuple)
    (sum(matches(a, b) for (a, b) in zip(case, tuple)) > 0 &&
     !any(mismatch(a, b) for (a, b) in zip(case, tuple)))
end

"""
Example matches.
case  [1 1 0 2]
tuple [1 0 0 2]
"""
function case_covers_tuple(case, tuple)
    all(matches(a, b) || skips(a, b) || ignores(a, b) for (a, b) in zip(case, tuple))
end

"""
Given an entry in the test set that has missing values, which are
zeros, find any matches that could be created by setting those
missing values. Return a new version of the entry.
"""
function matches_from_missing(mc::MatrixCoverage, entry, missing_param, matcher = case_compatible_with_tuple)
    param_cnt = parameter_cnt(mc)
    hist = zeros(eltype(mc), mc.arity[missing_param])
    for tuple_idx in 1:mc.remain
        if mc.allc[missing_param, tuple_idx] != 0
            # case_covers_tuple - variation.
            if matcher(entry, mc.allc[:, tuple_idx])
                hist[mc.allc[missing_param, tuple_idx]] += 1
            end
        end  # No matches unless the particular column is nonzero.
    end
    hist
end

"""
    first_match_for_parameter(mc::MatrixCoverage, param_idx)

Find the first uncovered tuple for a particular parameter value.
"""
function first_match_for_parameter(mc::MatrixCoverage, param_idx)
    for col_idx in 1:mc.remain
        if mc.allc[param_idx, col_idx] != 0
            return mc.allc[:, col_idx]
        end  # else keep looking
    end
    return zeros(eltype(mc), length(mc.arity))
end

"""
    fill_consistent_matches(mc::MatrixCoverage, entry)

Given an entry that's partially decided, fill in every covering tuple
that matches the existing values, in any order.
"""
function fill_consistent_matches(mc::MatrixCoverage, entry, matcher = case_compatible_with_tuple)
    param_cnt = parameter_cnt(mc)
    for col_idx in 1:mc.remain
        if matcher(entry, mc.allc[:, col_idx])
            for copy_idx in 1:param_cnt
                if mc.allc[copy_idx, col_idx] != 0
                    entry[copy_idx] = mc.allc[copy_idx, col_idx]
                end
            end
        end
    end
    entry
end

"""
    add_coverage!(allc, entry)

The coverage matrix, `allc`, has `row_cnt` entries that are considered
uncovered, and the rest have been covered already. This adds a new entry,
which is a complete set of parameter choices. For every parameter that's
covered, this moves those to the end. It returns a new `row_cnt` which is
the number of initial uncovered rows.
"""
function add_coverage!(mc::MatrixCoverage, entry)
    param_cnt = length(entry)

    covers = zeros(Int, mc.remain)
    cover_cnt = 0
    # Find matches before reordering them to the end.
    for col_idx in 1:mc.remain
        if case_covers_tuple(entry, mc.allc[:, col_idx])
            cover_cnt += 1
            covers[cover_cnt] = col_idx
        end
    end
    # Given the matches, we can swap them to the end of the matrix.
    # Work from the end in case a match is near row_cnt.
    for cover_idx in cover_cnt:-1:1
        if mc.remain > 1
            save = mc.allc[:, mc.remain]
            mc.allc[:, mc.remain] = mc.allc[:, covers[cover_idx]]
            mc.allc[:, covers[cover_idx]] = save
            mc.remain -= 1
        else
            mc.remain -= 1
        end
    end
    mc.remain
end

"""
    match_score(allc, entry)

If we were to add this entry to the trials, how many uncovered
n-tuples would it now cover? `allc` is the matrix of n-tuples.
`row_cnt` is the first set of rows, representing uncovered tuples.
`n_way` is the length of each tuple. Entry is a set of putative
parameter values. Returns an integer number of newly-covered tuples.

This is a place to look into alternative algorithms. We could weight
the match score by increasing the score for greater wayness.
"""
function match_score(mc::MatrixCoverage, entry)
    param_cnt = length(entry)
    cover_cnt = 0
    for col_idx in 1:mc.remain
        param_match_cnt = 0
        n_way = 0
        for match_idx in 1:param_cnt
            if mc.allc[match_idx, col_idx] == entry[match_idx]
                param_match_cnt += 1
            end
            if mc.allc[match_idx, col_idx] != 0
                n_way += 1
            end
        end
        if param_match_cnt == n_way
            cover_cnt += 1
        end
    end
    cover_cnt
end

"""
    remove_combinations!(mc::MatrixCoverage, disallow)

Remove all tuples that are disallowed by the `disallow` function.
Input to the function is a vector of parameter indexes.
The output to the function is whether they are disallowed,
as a boolean. This reduces the size of total tuples in the
coverage matrix. It doesn't move them to covered tuples. It
deletes them.
"""
function remove_combinations!(mc::MatrixCoverage, disallow)
    allow_cnt = 0
    allowed = zeros(Int, size(mc.allc, 2))
    for i in 1:size(mc.allc, 2)
        if !disallow(mc.allc[:, i])
            allow_cnt += 1
            allowed[allow_cnt] = i
        # else disallowed
        end
    end
    mc.allc = mc.allc[:, allowed[1:allow_cnt]]
    mc.remain = allow_cnt
end

"""
Builds a coverage set where subsets of variables have different wayness.

The all-pairs algorithm is 2-way. All-triples is 3-way. This function
takes a `base_wayness`, which would be 2 or 3 for those examples.
Then it takes a dictionary that specifies higher wayness for subsets
of parameters.

The wayness is a dictionary from an integer, the wayness, to a set
of lists of indices that should have that wayness together. The `base_wayness`
is the `n_way` for the rest of the variables. It should be less than
the other sets of waynesses.
"""
function multi_way_coverage(arity, wayness, base_wayness)
    orders = sort(collect(keys(wayness)), rev = true)
    @assert minimum(orders) > base_wayness
    @assert maximum(orders) <= length(arity)

    param_cnt = length(arity)
    order_combos = Any[]
    for order in orders
        # The parameter set is a list of parameter indices.
        for param_set in wayness[order]
            high_combos = all_combinations(arity[param_set], order)
            widened = zeros(eltype(arity), param_cnt, size(high_combos, 2))
            widened[param_set, :] = high_combos
            push!(order_combos, widened)
        end
        # We could remove duplicates of higher orders.
    end
    push!(order_combos, all_combinations(arity, base_wayness))
    hcat(order_combos...)
end


text:
#last Commits on Sep 14, 2018
struct GError
    domain::UInt32
    code::Cint
    message::Ptr{UInt8}
end #ok
@make_gvalue(GError, Ptr{GError}, :boxed, (:g_error, :libgobject)) #GV
convert(::Type{GError}, err::Ptr{GError}) = GError(err)

GError(err::Ptr{GError}) = unsafe_load(err)
function GError(f::Function)
    err = mutable(Ptr{GError})
    err.x = C_NULL
    if !f(err) || err[] != C_NULL
        gerror = GError(err[])
        emsg = bytestring(gerror.message)
        ccall((:g_clear_error, libglib), Nothing, (Ptr{Ptr{GError}},), err)
        error(emsg)
    end
end

ext: jl

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
function codepart(cnt,prevc)
    codepart=cnt>1 ? string(cnt) : ""
    return codepart * prevc
end

function encode(s)
    # prevc=s
    prevc = length(s)>0 ? s[1] : s
    coded = ""
    cnt=0

    for c in s
        if c == prevc
            cnt = cnt + 1
        else
            # codepart=cnt>1 ? string(cnt) : ""
            # print(prevc)

            # print(codepart(cnt,prevc))

            coded = coded * codepart(cnt,prevc)
            # prevc = ""
            cnt = 1
        end
        prevc = c
    end

    return coded*codepart(cnt,prevc)
end

function decode(s)
    concat = "1"
    decoded = ""

    for c in s
        isnumb = isnumeric(c)

        if isnumb
            # collect new chars as part of number
            concat!="1" ? concat = concat * c : concat=c
        else
            cnt = parse(Int64, concat)
            decoded = decoded * repeat(c, cnt)
            concat = "1"
        end
    end
    return decoded
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
#= This is one of the famous Interview Problems. You are given
an array of numbers, treat this array as stone heights. When the rainfalls
the water get's filled inside the stone gaps. You need to calculate the amount
of water which will get trapped inside the stones.
The code calculates the answer in the best time complexity.
=#

## Function

function rainWater(n, array)
    left = zeros(n)
    left[1] = array[1]
    for i = 2:n
        left[i] = max(left[i-1], array[i])
    end
    right = zeros(n)
    right[n] = array[n]
    for i = n-1:-1:1
        right[i] = max(right[i+1], array[i])
    end
    total = 0
    for i = 1:n
        total += min(left[i], right[i]) - array[i]
    end
    return total
end

## Input

println("Enter the length of array!")
n = readline()
n = parse(Int64, n)
array = Int64[]
println("Enter the array !")
for i = 1:n
    temp = readline()
    temp = parse(Int64, temp)
    push!(array, temp)
end

## Calling the function

rainWater(n, array)

#=
Sample Test Case:
Input:
    n = 3
    array = 2 0 2
Output:
    2

Time complexity: O ( N )
=#

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
using CachedCalls
using Test

descr = "subtr"
f = -
args = (5,2)

cc = CachedCall(descr, f, args)
cc2 = CachedCall(f, args)

@test call!(cc) == call!(cc2) == 3

@test is_executed(cc)
@test func_args(cc) == args
@test func(cc) == f
@test result(cc) == 3
@test isequal(description(cc), descr)
#test show (uncalled)
@test isa(repr(cc), AbstractString)
uncall!(cc)
@test !is_executed(cc)

#test show (called)
@test isa(repr(cc), AbstractString)

##### test recursive args #####
ccr = CachedCall(descr, f, (4, cc))
@test call!(ccr) == 1 # 4-3
@test is_executed(ccr)
@test func_args(ccr) == (4, cc)
@test is_executed(cc) #the arg should be executed
@test func(ccr) == f
@test result(ccr) == 1
@test isa(repr(ccr), AbstractString)
uncall!(ccr)
@test !is_executed(ccr)
@test !is_executed(cc) #the arg should be unexecuted
@test isa(repr(ccr), AbstractString)


text:
module CommunalHelperCustomSummitGem

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/CustomSummitGem" CustomSummitGem(x::Integer, y::Integer,
	index::Integer=0)

const placements = Ahorn.PlacementDict(
	"Summit Gem (Communal Helper)" => Ahorn.EntityPlacement(
		CustomSummitGem
	)
)

const sprites = ["collectables/summitgems/$i/gem00" for i = 0:7]

function getSprite(index)
	if index > length(sprites)
		return sprites[end]
	end
	return sprites[index]
end

# positive numbers only
function getClampedIndex(entity::CustomSummitGem)
	index = Int(get(entity.data, "index", 0))
	if index < 0
		index = 0
	end
	entity.data["index"] = index
	return index
end

function Ahorn.selection(entity::CustomSummitGem)
	x, y = Ahorn.position(entity)

	return Ahorn.getSpriteRectangle(getSprite(getClampedIndex(entity) + 1) , x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSummitGem) = Ahorn.drawSprite(ctx, getSprite(getClampedIndex(entity) + 1), 0, 0)

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
function gen(data,mva)
    data["pmax"]=mva
    data["pmin"]=mva
    return data
end

function converterdc_ne(conv,nde,mbase)
    #NOTE change
    conv["cost"]=nde["xfm_off"].costs.cpx_p+nde["xfm_off"].costs.cpx_i+nde["conv_off"].costs.cpx+nde["plat"].costs.cpx+nde["cableWF"].costs.cpx_p+nde["cableWF"].costs.cpx_i;
    conv["Pacmax"]=nde["conv_off"].mva/mbase
    conv["Pacmin"]=-1*nde["conv_off"].mva/mbase
    conv["Qacmax"]=nde["conv_off"].mva/mbase
    conv["Qacmin"]=-1*nde["conv_off"].mva/mbase
    return conv
end

function convdc_ne(conv,nde,mbase)
    #NOTE change
    conv["cost"]=nde["dc"]["xfm"].costs.cpx_p+nde["dc"]["xfm"].costs.cpx_i+nde["dc"]["conv"].costs.cpx
    #conv["cost"]=0
    if (haskey(nde["dc"],"plat")); conv["cost"]=conv["cost"]+nde["dc"]["plat"].costs.cpx;end
    #if (haskey(nde["dc"],"plat")); conv["cost"]=0;end
    #########
    conv["Pacmax"]=nde["mva"]/mbase
    conv["Pacmin"]=-1*nde["mva"]/mbase
    conv["Qacmax"]=nde["mva"]/mbase
    conv["Qacmin"]=-1*nde["mva"]/mbase
    return conv
end

function branchdc_on(brc,dcss,mbase)
    brc["cost"]=brc["cost"]+dcss["xfm_on"].costs.cpx_p+dcss["xfm_on"].costs.cpx_i+dcss["conv_on"].costs.cpx
    return brc
end

function converterdc(conv,nde,mbase)
    conv["Pacmax"]=nde.mva/mbase
    conv["Pacmin"]=-1*nde.mva/mbase
    conv["Qacmax"]=nde.mva/mbase
    conv["Qacmin"]=-1*nde.mva/mbase
    return conv
end

function branchdc_ne(b_ne,dc_can,mbase)
    b_ne["rateA"]=(dc_can.num*dc_can.elec.mva)/mbase
    b_ne["rateB"]=(dc_can.num*dc_can.elec.mva)/mbase
    b_ne["rateC"]=(dc_can.num*dc_can.elec.mva)/mbase
    b_ne["r"]=dc_can.elec.ohm

    #NOTE change
    b_ne["cost"]=dc_can.costs.cpx_p+dc_can.costs.cpx_i
    #b_ne["cost"]=0
    #######
    return b_ne
end

function ne_branch(ne_b,ac_can,mbase)
    ne_b["mva"]=ac_can.num*ac_can.elec.mva
    ne_b["rate_a"]=ac_can.num*ac_can.elec.mva/mbase
    ne_b["rate_b"]=ac_can.num*ac_can.elec.mva/mbase
    ne_b["rate_c"]=ac_can.num*ac_can.elec.mva/mbase
    ne_b["br_r"]=ac_can.elec.ohm
    ne_b["br_x"]=ac_can.elec.xl
    #data["ne_branch"][datai]["construction_cost"]=hoa["candidate_lines"]["ac"][hoai].costs.cpx_p+hoa["candidate_lines"]["ac"][hoai].costs.cpx_i+hoa["candidate_lines"]["ac"][hoai].costs.qc+hoa["candidate_lines"]["ac"][hoai].costs.sg
    #ne_b["construction_cost"]=(ac_can.costs.cpx_p+ac_can.costs.cpx_i)/mbase

    #NOTE change
    ne_b["construction_cost"]=ac_can.costs.cpx_p+ac_can.costs.cpx_i
    #ne_b["construction_cost"]=0
    ##########
    return ne_b
end

function branch_xfo(ne_b,ne_xfo,mbase)
    #ne_b["construction_cost"]=ne_b["construction_cost"]+(ne_xfo.costs.cpx_p+ne_xfo.costs.cpx_i)/mbase

    #NOTE change
    ne_b["construction_cost"]=ne_b["construction_cost"]+ne_xfo.costs.cpx_p+ne_xfo.costs.cpx_i
    #ne_b["construction_cost"]=ne_b["construction_cost"]+0
    ##########
    return ne_b
end

function branch_plat(ne_b,ne_plat,mbase)
    #ne_b["construction_cost"]=ne_b["construction_cost"]+(ne_plat.costs.cpx)/mbase

    #NOTE change
    ne_b["construction_cost"]=ne_b["construction_cost"]+ne_plat.costs.cpx
    #ne_b["construction_cost"]=ne_b["construction_cost"]+0
    ##########
    return ne_b
end

function pm_dict_candidateICs(data,hoa,candidates)
    #WF AC candidate lines
    for (i,k) in enumerate(vcat(candidates[2],candidates[2]))
        if (i<=length(candidates[2]))
            data["ne_branch"][string(i)]=ne_branch(data["ne_branch"][string(i)], hoa["candidates"]["wf"]["ac"][string(k)]["cableBE"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_xfo(data["ne_branch"][string(i)],hoa["candidates"]["wf"]["ac"][string(k)]["xfm_on"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_xfo(data["ne_branch"][string(i)],hoa["candidates"]["wf"]["ac"][string(k)]["xfm_off"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_plat(data["ne_branch"][string(i)], hoa["candidates"]["wf"]["ac"][string(k)]["plat"], data["gen"]["1"]["mbase"])
        else
            data["ne_branch"][string(i)]=ne_branch(data["ne_branch"][string(i)], hoa["candidates"]["wf"]["ac"][string(k)]["cableUK"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_xfo(data["ne_branch"][string(i)],hoa["candidates"]["wf"]["ac"][string(k)]["xfm_on"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_xfo(data["ne_branch"][string(i)],hoa["candidates"]["wf"]["ac"][string(k)]["xfm_off"], data["gen"]["1"]["mbase"])
            data["ne_branch"][string(i)]=branch_plat(data["ne_branch"][string(i)], hoa["candidates"]["wf"]["ac"][string(k)]["plat"], data["gen"]["1"]["mbase"])
        end
    end

    for (i,k) in enumerate(vcat(candidates[1],candidates[1]))
        if (i<=length(candidates[2]))
            data["branchdc_ne"][string(i)]=branchdc_ne(data["branchdc_ne"][string(i)],hoa["candidates"]["ic"][string(k)]["cableBE"],data["gen"]["1"]["mbase"])
            data["branchdc_ne"][string(i)]=branchdc_on(data["branchdc_ne"][string(i)],hoa["candidates"]["ic"][string(k)],data["gen"]["1"]["mbase"])
        else
            data["branchdc_ne"][string(i)]=branchdc_ne(data["branchdc_ne"][string(i)],hoa["candidates"]["ic"][string(k)]["cableUK"],data["gen"]["1"]["mbase"])
            data["branchdc_ne"][string(i)]=branchdc_on(data["branchdc_ne"][string(i)],hoa["candidates"]["ic"][string(k)],data["gen"]["1"]["mbase"])
        end
    end

    for (i,k) in enumerate(candidates[2])
        data["convdc_ne"][string(i)]=converterdc_ne(data["convdc_ne"][string(i)],hoa["candidates"]["wf"]["dc"][string(k)],data["gen"]["1"]["mbase"])
    end
    size_set=[];for (k,v) in hoa["candidates"]["ic"]; push!(size_set,parse(Float64,k));end
    mx_key=string(maximum(size_set))
    data["convdc"]["1"]=converterdc(data["convdc"]["1"],hoa["candidates"]["ic"][mx_key]["conv_on"],data["gen"]["1"]["mbase"])
    data["convdc"]["2"]=converterdc(data["convdc"]["2"],hoa["candidates"]["ic"][mx_key]["conv_on"],data["gen"]["1"]["mbase"])

    return data
end

function additional_candidates(data,candidates)
    #AC
    c2be=deepcopy(data["ne_branch"]["1"]);c2uk=deepcopy(data["ne_branch"]["2"])
    data["ne_branch"]=Dict{String,Any}()
    #WF
    for i=1:1:length(candidates[2]); c2be["source_id"][2]=i;c2be["index"]=i; push!(data["ne_branch"],string(i)=>deepcopy(c2be)); end
    for i=length(candidates[2])+1:1:length(candidates[2])*2; c2uk["source_id"][2]=i;c2uk["index"]=i; push!(data["ne_branch"],string(i)=>deepcopy(c2uk)); end
    #DC
    #WF
    convwf=deepcopy(data["convdc_ne"]["1"])
    data["convdc_ne"]=Dict{String,Any}()
    for i=1:1:length(candidates[2]); convwf["source_id"][2]=i;convwf["index"]=i; push!(data["convdc_ne"],string(i)=>deepcopy(convwf)); end
    #IC
    c2be=deepcopy(data["branchdc_ne"]["1"]);c2uk=deepcopy(data["branchdc_ne"]["2"])
    data["branchdc_ne"]=Dict{String,Any}()
    for i=1:1:length(candidates[1]); c2be["source_id"][2]=i;c2be["index"]=i; push!(data["branchdc_ne"],string(i)=>deepcopy(c2be)); end
    for i=length(candidates[1])+1:1:length(candidates[1])*2; c2uk["source_id"][2]=i;c2uk["index"]=i; push!(data["branchdc_ne"],string(i)=>deepcopy(c2uk)); end
    return data
end

#calculates candidate equipment and sorts in struct
function hoa_datastruct_candidateICs(ic_mva,owpp_mva,ic_length,owpp_km,candidates)
    #location
    uk_pcc=xy();uk_pcc.x=ic_length;uk_pcc.y=0;
    be_pcc=xy();be_pcc.x=0;be_pcc.y=0;
    owpp_xy=xy();owpp_xy.x=owpp_km;owpp_xy.y=0
    #datastructure
    hoa=Dict{String,Any}();
    push!(hoa,"nodes"=>Dict{String,Any}());
    push!(hoa["nodes"],"ac"=>Dict{String,Any}());
    push!(hoa["nodes"],"dc"=>Dict{String,Any}());

    [push!(hoa["nodes"]["ac"],n=>Dict{String,Any}()) for n in ["1","2","3"]]
    [push!(hoa["nodes"]["dc"],n=>Dict{String,Any}()) for n in ["1","2","3","4"]]
    #position of nodes
    push!(hoa["nodes"]["ac"]["1"],"xy"=>be_pcc);
    push!(hoa["nodes"]["ac"]["2"],"xy"=>uk_pcc);
    push!(hoa["nodes"]["ac"]["3"],"xy"=>owpp_xy);
    #position of nodes
    push!(hoa["nodes"]["dc"]["1"],"xy"=>be_pcc);
    push!(hoa["nodes"]["dc"]["2"],"xy"=>uk_pcc);
    push!(hoa["nodes"]["dc"]["3"],"xy"=>owpp_xy);
    push!(hoa["nodes"]["dc"]["4"],"xy"=>owpp_xy);
    #rating of nodes
    push!(hoa["nodes"]["ac"]["1"],"mva"=>ic_mva);
    push!(hoa["nodes"]["ac"]["2"],"mva"=>ic_mva);
    push!(hoa["nodes"]["ac"]["3"],"mva"=>owpp_mva);
    #rating of nodes
    push!(hoa["nodes"]["dc"]["1"],"mva"=>ic_mva);
    push!(hoa["nodes"]["dc"]["2"],"mva"=>ic_mva);
    push!(hoa["nodes"]["dc"]["3"],"mva"=>owpp_mva);
    push!(hoa["nodes"]["dc"]["4"],"mva"=>owpp_mva);
    #lines
    push!(hoa,"candidates"=>Dict{String,Any}());
    push!(hoa["candidates"],"wf"=>Dict{String,Any}());
    push!(hoa["candidates"],"ic"=>Dict{String,Any}());
    push!(hoa["candidates"]["wf"],"ac"=>Dict{String,Any}());
    push!(hoa["candidates"]["wf"],"dc"=>Dict{String,Any}());

    [push!(hoa["candidates"]["wf"]["ac"],string(p)=>Dict{String,Any}()) for p in candidates[2]]
    [push!(hoa["candidates"]["wf"]["dc"],string(p)=>Dict{String,Any}()) for p in candidates[2]]
    [push!(hoa["candidates"]["ic"],string(p)=>Dict{String,Any}()) for p in candidates[1]]

    #wf-ac
    [push!(hoa["candidates"]["wf"]["ac"][k],"xfm_on"=>on_xfm(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    [push!(hoa["candidates"]["wf"]["ac"][k],"xfm_off"=>off_xfm(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    [push!(hoa["candidates"]["wf"]["ac"][k],"plat"=>AC_plat(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    [push!(hoa["candidates"]["wf"]["ac"][k],"xfm"=>on_xfm(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    [push!(hoa["candidates"]["wf"]["ac"][k],"cableBE"=>AC_cbl(owpp_mva*parse(Float64,k), owpp_km)) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    [push!(hoa["candidates"]["wf"]["ac"][k],"cableUK"=>AC_cbl(owpp_mva*parse(Float64,k), ic_length-owpp_km)) for (k,p) in hoa["candidates"]["wf"]["ac"]]
    #wf-dc
    [push!(hoa["candidates"]["wf"]["dc"][k],"xfm_off"=>off_xfm(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["dc"]]
    [push!(hoa["candidates"]["wf"]["dc"][k],"conv_off"=>off_conv(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["dc"]]
    [push!(hoa["candidates"]["wf"]["dc"][k],"plat"=>DC_plat(owpp_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["wf"]["dc"]]
    [push!(hoa["candidates"]["wf"]["dc"][k],"cableWF"=>DC_cbl(owpp_mva/2*parse(Float64,k), 1)) for (k,p) in hoa["candidates"]["wf"]["dc"]]
    #ic dc
    [push!(hoa["candidates"]["ic"][k],"xfm_on"=>on_xfm(ic_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["ic"]]
    [push!(hoa["candidates"]["ic"][k],"conv_on"=>on_conv(ic_mva*parse(Float64,k))) for (k,p) in hoa["candidates"]["ic"]]
    [push!(hoa["candidates"]["ic"][k],"cableBE"=>DC_cbl(ic_mva/2*parse(Float64,k), owpp_km)) for (k,p) in hoa["candidates"]["ic"]]
    [push!(hoa["candidates"]["ic"][k],"cableUK"=>DC_cbl(ic_mva/2*parse(Float64,k), ic_length-owpp_km)) for (k,p) in hoa["candidates"]["ic"]]

    return hoa
end


text:
# # [Introduction to periodic problems and plane-wave discretisations](@id periodic-problems)
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/guide/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/guide/@__NAME__.ipynb)

# In this example we want to show how DFTK can be used to solve simple one-dimensional
# periodic problems. Along the lines this notebook serves as a concise introduction into
# the underlying theory and jargon for solving periodic problems using plane-wave
# discretizations.

# ## Periodicity and lattices
# A periodic problem is characterized by being invariant to certain translations.
# For example the ``\sin`` function is periodic with periodicity ``2π``, i.e.
# ```math
#    \sin(x) = \sin(x + 2πm) \quad ∀ m ∈ \mathbb{Z},
# ```
# This is nothing else than saying that any translation by an integer multiple of ``2π``
# keeps the ``\sin`` function invariant. In a more formal one can use the
# translation operator ``T_{-2πm}`` to write this as:
# ```math
#    T_{-2πm} \, \sin(x) = \sin(x + 2πm) = \sin(x).
# ```
#
# Whenever such periodicity exists one can exploit it to save computational work.
# Consider a problem in which we want to find a function ``f : \mathbb{R} → \mathbb{R}``,
# but *a priori* the solution is known to be periodic with periodicity ``a``. As a consequence
# of said periodicity it is sufficient to determine the values of ``f`` for all ``x`` from the
# interval ``[-a/2, a/2)`` to uniquely define ``f`` on the full real axis. Naturally exploiting
# periodicity in a computational procedure thus greatly reduces the required amount of work.
#
# Let us introduce some jargon: The periodicity of our problem implies that we may define
# a **lattice**
# ```
#         -3a/2      -a/2      +a/2     +3a/2
#        ... |---------|---------|---------| ...
#                 a         a         a
# ```
# with lattice constant ``a``. Each cell of the lattice is an identical periodic image of
# any of its neighbors. For finding ``f`` it is thus sufficient to consider only the
# problem inside a **unit cell** ``[-a/2, a/2)``. In passing we note that the definition
# of the unit cell is itself only unique up to translations. A choice ``[0, a)``,
# for example, would have done just as well.
#
# ## Periodic operators and the Bloch transform
# Not only functions, but also operators can feature periodicity.
# Consider for example the **free-electron Hamiltonian**
# ```math
#     H = -\frac12 Δ.
# ```
# In free-electron model (which gives rise to this Hamiltonian) electron motion is only
# by their own kinetic energy. As this model features no potential which could make one point
# in space more preferred than another, we would expect this model to be periodic.
# If an operator is periodic with respect to a lattice such as the one defined above,
# than it commutes with all lattice translations. For the free-electron case ``H``
# one can easily show exactly that, i.e.
# ```math
#    T_{ma} H = H T_{ma} \quad  ∀ m ∈ \mathbb{Z}.
# ```
# We note in passing that the free-electron model is actually very special in the sense that
# the choice of ``a`` is completely arbitrary here. In other words ``H`` is periodic
# with respect to any translation. In general, however, periodicity is only
# attained with respect to a finite number of translations ``a`` and we will take this
# viewpoint here.
#
# **Bloch's theorem** now tells us that for periodic operators,
# the solutions to the eigenproblem
# ```math
#     H ψ_{kn} = ε_{kn} ψ_{kn}
# ```
# satisfy a factorization
# ```math
#     ψ_{kn}(x) = e^{i k⋅x} u_{kn}(x)
# ```
# into a plane wave ``e^{i k⋅x}`` and a lattice-periodic function
# ```math
#    T_{ma} u_{kn}(x) = u_{kn}(x - ma) = u_{kn}(x) \quad ∀ m ∈ \mathbb{Z}.
# ```
# In this ``n`` is a labeling integer index and ``k`` is a real number,
# whose details will be clarified in the next section.
# The index ``n`` is sometimes also called the **band index** and
# functions ``ψ_{kn}`` satisfying this factorization are also known as
# **Bloch functions** or **Bloch states**.
#
# Consider the application of ``2H = -Δ = - \frac{d^2}{d x^2}``
# to such a Bloch wave. First we notice for any function ``f``
# ```math
#    -i∇ \left( e^{i k⋅x} f \right)
#    = -i\frac{d}{dx} \left( e^{i k⋅x} f \right)
#    = k e^{i k⋅x} f -i∇ e^{i k⋅x} f = (-i∇ + k) e^{i k⋅x} f.
# ```
# Using this result twice one shows that applying ``-Δ`` yields
# ```math
# \begin{aligned}
#    -\Delta \left(e^{i k⋅x} u_{kn}(x)\right)
#    &= -i∇ ⋅ \left[-i∇ \left(u_{kn}(x) e^{i k⋅x} \right) \right] \\
#    &= -i∇ ⋅ \left[(-i∇ + k) u_{kn}(x) e^{i k⋅x} \right] \\
#    &= (-i∇ + k)^2 u_{kn}(x) e^{i k⋅x} \\
#    &= e^{i k⋅x} 2H_k u_{kn}(x),
# \end{aligned}
# ```
# where we defined
# ```math
#     H_k = \frac12 (-i∇ + k)^2.
# ```
# The action of this operator on a function ``u_{kn}`` is given by
# ```math
#     H_k u_{kn} = e^{-i k⋅x} H e^{i k⋅x} u_{kn},
# ```
# which in particular implies that
# ```math
#    H_k u_{kn} = ε_{kn} u_{kn} \quad ⇔ \quad H (e^{i k⋅x} u_{kn}) = ε_{kn} (e^{i k⋅x} u_{kn}).
# ```
# To seek the eigenpairs of ``H`` we may thus equivalently
# find the eigenpairs of *all* ``H_k``.
# The point of this is that the eigenfunctions ``u_{kn}`` of ``H_k``
# are periodic (unlike the eigenfunctions ``ψ_{kn}`` of ``H``).
# In contrast to ``ψ_{kn}`` the functions ``u_{kn}`` can thus be fully
# represented considering the eigenproblem only on the unit cell.
#
# A detailed mathematical analysis shows that the transformation from ``H``
# to the set of all ``H_k`` for a suitable set of values for ``k`` (details below)
# is actually a unitary transformation, the so-called **Bloch transform**.
# This transform brings the Hamiltonian into the symmetry-adapted basis for
# translational symmetry, which are exactly the Bloch functions.
# Similar to the case of choosing a symmetry-adapted basis for other kinds of symmetries
# (like the point group symmetry in molecules), the Bloch transform also makes
# the Hamiltonian ``H`` block-diagonal[^1]:
# ```math
#     T_B H T_B^{-1} ⟶ \left( \begin{array}{cccc} H_1&&&0 \\ &H_2\\&&H_3\\0&&&\ddots \end{array} \right)
# ```
# with each block ``H_k`` taking care of one value ``k``.
# This block-diagonal structure under the basis of Bloch functions lets us
# completely describe the spectrum of ``H`` by looking only at the spectrum
# of all ``H_k`` blocks.
#
# [^1]: Notice that block-diagonal is a bit an abuse of terms here, since the Hamiltonian
#       is not a matrix but an operator and the number of blocks is essentially infinite.
#       The mathematically precise term is that the Bloch transform reveals the fibers
#       of the Hamiltonian.
#
# ## The Brillouin zone
#
# We now consider the parameter ``k`` of the Hamiltonian blocks in detail.
#
# - As discussed ``k`` is a real number. It turns out, however, that some of
#   these ``k∈\mathbb{R}`` give rise to operators related by unitary transformations
#   (again due to translational symmetry).
# - Since such operators have the same eigenspectrum, only one version needs to be considered.
# - The smallest subset from which ``k`` is chosen is the **Brillouin zone** (BZ).
#
# - The BZ is the unit cell of the **reciprocal lattice**, which may be constructed from
#   the **real-space lattice** by a Fourier transform.
# - In our simple 1D case the reciprocal lattice is just
#   ```
#     ... |--------|--------|--------| ...
#            2π/a     2π/a     2π/a
#   ```
#   i.e. like the real-space lattice, but just with a different lattice constant
#   ``b = 2π / a``.
# - The BZ in our example is thus ``B = [-π/a, π/a)``. The members of ``B``
#   are typically called ``k``-points.
#
# ## Discretization and plane-wave basis sets
#
# With what we discussed so far the strategy to find all eigenpairs of a periodic
# Hamiltonian ``H`` thus reduces to finding the eigenpairs of all ``H_k`` with ``k ∈ B``.
# This requires *two* discretisations:
#
#   - ``B`` is infinite (and not countable). To discretize we first only pick a finite number
#     of ``k``-points. Usually this **``k``-point sampling** is done by picking ``k``-points
#     along a regular grid inside the BZ, the **``k``-grid**.
#   - Each ``H_k`` is still an infinite-dimensional operator.
#     Following a standard Ritz-Galerkin ansatz we project the operator into a finite basis
#     and diagonalize the resulting matrix.
#
# For the second step multiple types of bases are used in practice (finite differences,
# finite elements, Gaussians, ...). In DFTK we currently support only plane-wave
# discretizations.
#
# For our 1D example normalized plane waves are defined as the functions
# ```math
# e_{G}(x) = \frac{e^{i G x}}{\sqrt{a}}  \qquad G \in b\mathbb{Z}
# ```
# and typically one forms basis sets from these by specifying a
# **kinetic energy cutoff** ``E_\text{cut}``:
# ```math
# \left\{ e_{G} \, \big| \, (G + k)^2 \leq 2E_\text{cut} \right\}
# ```
#
# ## Correspondence of theory to DFTK code
#
# Before solving a few example problems numerically in DFTK, a short overview
# of the correspondence of the introduced quantities to data structures inside DFTK.
#
# - ``H`` is represented by a `Hamiltonian` object and variables for hamiltonians are usually called `ham`.
# - ``H_k`` by a `HamiltonianBlock` and variables are `hamk`.
# - ``ψ_{kn}`` is usually just called `ψ`.
#   ``u_{kn}`` is not stored (in favor of ``ψ_{kn}``).
# - ``ε_{kn}`` is called `eigenvalues`.
# - ``k``-points are represented by `Kpoint` and respective variables called `kpt`.
# - The basis of plane waves is managed by `PlaneWaveBasis` and variables usually just called `basis`.
#
# ## Solving the free-electron Hamiltonian
#
# One typical approach to get physical insight into a Hamiltonian ``H`` is to plot
# a so-called **band structure**, that is the eigenvalues of ``H_k`` versus ``k``.
# In DFTK we achieve this using the following steps:
#
# Step 1: Build the 1D lattice. DFTK is mostly tailored for 3D problems.
# Therefore quantities related to the problem space are have a fixed
# dimension 3. The convention is that for 1D / 2D problems the
# trailing entries are always zero and ignored in the computation.
# For the lattice we therefore construct a 3x3 matrix with only one entry.
using DFTK

lattice = zeros(3, 3)
lattice[1, 1] = 20.;

# Step 2: Select a model. In this case we choose a free-electron model,
# which is the same as saying that there is only a Kinetic term
# (and no potential) in the model. The `n_electrons` is dummy here.
model = Model(lattice; n_electrons=0, terms=[Kinetic()])

# Step 3: Define a plane-wave basis using this model and a cutoff ``E_\text{cut}``
# of 300 Hartree. The ``k``-point grid is given as a regular grid in the BZ
# (a so-called **Monkhorst-Pack** grid). Here we select only one ``k``-point (1x1x1),
# see the note below for some details on this choice.
basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1))

# Step 4: Plot the bands! Select a density of ``k``-points for the ``k``-grid to use
# for the bandstructure calculation, discretize the problem and diagonalize it.
# Afterwards plot the bands.

using Unitful
using UnitfulAtomic
using Plots

plot_bandstructure(basis; n_bands=6, kline_density=30)

# !!! note "Selection of k-point grids in `PlaneWaveBasis` construction"
#     You might wonder why we only selected a single ``k``-point (clearly a very crude
#     and inaccurate approximation). In this example the `kgrid` parameter specified
#     in the construction of the `PlaneWaveBasis`
#     is not actually used for plotting the bands. It is only used when solving more
#     involved models like density-functional theory (DFT) where the Hamiltonian is
#     non-linear. In these cases before plotting the bands the self-consistent field
#     equations (SCF) need to be solved first. This is typically done on
#     a different ``k``-point grid than the grid used for the bands later on.
#     In our case we don't need this extra step and therefore the `kgrid` value passed
#     to `PlaneWaveBasis` is actually arbitrary.

# ## Adding potentials
# So far so good. But free electrons are actually a little boring,
# so let's add a potential interacting with the electrons.
#
# - The modified problem we will look at consists of diagonalizing
#   ```math
#   H_k = \frac12 (-i \nabla + k)^2 + V
#   ```
#   for all ``k \in B`` with a periodic potential ``V`` interacting with the electrons.
#
# - A number of "standard" potentials are readily implemented in DFTK and
#   can be assembled using the `terms` kwarg of the model.
#   This allows to seamlessly construct
#
#   * density-functial theory (DFT) models for treating electronic structures
#     (see the [Tutorial](@ref)).
#   * Gross-Pitaevskii models for bosonic systems
#     (see [Gross-Pitaevskii equation in one dimension](@ref))
#   * even some more unusual cases like anyonic models.
#
# In this tutorial we will go a little more low-level and directly provide
# an analytic potential describing the interaction with the electrons to DFTK.
#
# First we define a custom Gaussian potential as a new "element" inside DFTK:

struct ElementGaussian <: DFTK.Element
    α  # Prefactor
    L  # Extend
end

## Some default values
ElementGaussian() = ElementGaussian(0.3, 10.0)

## Real-space representation of a Gaussian
function DFTK.local_potential_real(el::ElementGaussian, r::Real)
    -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)
end

## Fourier-space representation of the Gaussian
function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)
    ## = ∫ -α exp(-(r/L)^2 exp(-ir⋅q) dr
    -el.α * exp(- (q * el.L)^2 / 2)
end

# A single potential looks like:

using Plots
using LinearAlgebra
nucleus = ElementGaussian()
plot(r -> DFTK.local_potential_real(nucleus, norm(r)), xlims=(-50, 50))

# With this element at hand we can easily construct a setting
# where two potentials of this form are located at positions
# ``20`` and ``80`` inside the lattice ``[0, 100]``:

using LinearAlgebra

## Define the 1D lattice [0, 100]
lattice = diagm([100., 0, 0])

## Place them at 20 and 80 in *fractional coordinates*,
## that is 0.2 and 0.8, since the lattice is 100 wide.
nucleus = ElementGaussian()
atoms = [nucleus => [[0.2, 0, 0], [0.8, 0, 0]]]

## Assemble the model, discretize and build the Hamiltonian
model = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])
basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));
ham   = Hamiltonian(basis)

## Extract the total potential term of the Hamiltonian and plot it
potential = DFTK.total_local_potential(ham)[:, 1, 1]
rvecs = collect(r_vectors_cart(basis))[:, 1, 1]  # slice along the x axis
x = [r[1] for r in rvecs]                        # only keep the x coordinate
plot(x, potential, label="", xlabel="x", ylabel="V(x)")

# This potential is the sum of two "atomic" potentials (the two "Gaussian" elements).
# Due to the periodic setting we are considering interactions naturally also occur
# across the unit cell boundary (i.e. wrapping from `100` over to `0`).
# The required periodization of the atomic potential is automatically taken care,
# such that the potential is smooth across the cell boundary at `100`/`0`.
#
# With this setup, let's look at the bands:

using Unitful
using UnitfulAtomic

plot_bandstructure(basis; n_bands=6, kline_density=15)

# The bands are noticeably different.
#  - The bands no longer overlap, meaning that the spectrum of $H$ is no longer continuous
#    but has gaps.
#
#  - The two lowest bands are almost flat. This is because they represent
#    two tightly bound and localized electrons inside the two Gaussians.
#
#  - The higher the bands are in energy, the more free-electron-like they are.
#    In other words the higher the kinetic energy of the electrons, the less they feel
#    the effect of the two Gaussian potentials. As it turns out the curvature of the bands,
#    (the degree to which they are free-electron-like) is highly related to the delocalization
#    of electrons in these bands: The more curved the more delocalized. In some sense
#    "free electrons" correspond to perfect delocalization.


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bilevel_pilotables_scenarios" begin

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")

    #=
    TS: [11h]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |wind_2_1
    Pmin=0, Pmax=200    |                      |Pmin=0, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------| S1: 30
                        |         35           | S2: 70
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 0             |                      | S1: 100
      S2: 0             |                      | S2: 100

      In scenario 1,
          we can use all the limitable power (30MW),
          this leaves a residual demand of 70MW
          prod_1_1 can only supply 35MW due to RSO constraints
          => a Lol of 35MW
      In scenario 2,
          we can use all the limitable power (70MW),
          this leaves a residual demand of 30MW
          prod_1_1 can supply the required 30MW

      The TSO will need to take an action to limit production of prod_1_1 especially in scenario S1.
    =#
    function create_instance(limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "wind_2_1", PSCOPF.Networks.LIMITABLE,
                                                0., 200.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_2_1", TS[1], "S1", 30.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_2_1", TS[1], "S2", 70.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", 0.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S2", 0.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", 100.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S2", 100.)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    SortedDict{String,PSCOPF.GeneratorState}(), #gen_initial_state
                                    uncertainties, nothing, logs)

        return context
    end

    @testset "different_impositions_by_scenario" begin

        context = create_instance(35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(LINK_SCENARIOS_PILOTABLE_LEVEL=false,
                                LINK_SCENARIOS_PILOTABLE_ON=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution has slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test_broken 35. ≈ value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) #LoL is for RSO reason but is only noticed in market

        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"])

        #wind_2_1
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S1"]) > 30. - 1e-09
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S2"]) > 70. - 1e-09

        #Market chooses the levels sets the bounds for pilotable production
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 30. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"])

        #tso distributes limitable power
        @test 30. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S1"])
        @test 70. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S2"])

    end

    @testset "link_scenarios_in_tso" begin

        context = create_instance(35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(LINK_SCENARIOS_PILOTABLE_LEVEL=true,
                                LINK_SCENARIOS_PILOTABLE_ON=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution has slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test_broken 35. ≈ value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) #LoL is for RSO reason but is only noticed in market

        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1
        @test ( value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
                ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) )
        @test ( value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
                ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"]) )

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) < 1e-09
        @test 35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"])

        #wind_2_1
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S1"]) > 30. - 1e-09
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S2"]) > 70. - 1e-09

        #Market chooses the levels sets the bounds for pilotable production
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 30. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"])

        #tso distributes limitable power
        @test 30. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S1"])
        @test 70. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S2"])
    end

    #=
    This is not what is illustrated here but keep in mind !

    I'm afraid the TSO can cheat the market when
    LINK_SCENARIOS_PILOTABLE_LEVEL = false
    and
    LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET = true :
    The TSO can penalize one scenario (leaving the others free which costs less)
    but all the market scenarios will obey to the most strict scenario
    since scenarios are linked in the market.
    =#
    @testset "link_scenarios_in_market" begin

        context = create_instance(35.,)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution has slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test_broken 35. ≈ value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) #LoL is for RSO reason but is only noticed in market

        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        @test 5. ≈ value(result.lower.limitable_model.p_global_capping[TS[1],"S2"])
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]) > 35. - 1e-09

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) < 1e-09
        @test value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"]) > 35. - 1e-09

        @test ( (35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]))
                || (35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"])) )
        @test ( (200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]))
                || (200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"])) )

        #wind_2_1
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S1"]) > 30. - 1e-09
        @test 65. ≈ value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S2"])

        #Market chooses the levels sets the bounds for pilotable production
        @test ( value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
                ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"]) )
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"])

        #tso distributes limitable power
        @test 30. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S1"])
        @test 65. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S2"])

    end

    @testset "link_scenarios" begin

        context = create_instance(35.,)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET=true,
                                                        LINK_SCENARIOS_PILOTABLE_LEVEL=true,
                                                        LINK_SCENARIOS_PILOTABLE_ON=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution has slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test_broken 35. ≈ value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) #LoL is for RSO reason but is only noticed in market

        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S2"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        @test 5. ≈ value(result.lower.limitable_model.p_global_capping[TS[1],"S2"])
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S2"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1
        @test ( value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
                ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) )
        @test ( value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
                ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"]) )

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])

        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S2"]) < 1e-09
        @test 35. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S2"])

        #wind_2_1
        @test value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S1"]) > 30. - 1e-09
        @test 65. ≈ value(result.upper.limitable_model.p_limit["wind_2_1",TS[1],"S2"])

        #Market chooses the levels sets the bounds for pilotable production
        @test ( value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
                ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"]) )
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 35. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S2"])

        #tso distributes limitable power
        @test 30. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S1"])
        @test 65. ≈ value(result.upper.limitable_model.p_injected["wind_2_1",TS[1],"S2"])

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# ---
# title: 947. Most Stones Removed with Same Row or Column
# id: problem947
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Union Find
# link: <https://leetcode.com/problems/most-stones-removed-with-same-row-or-column/description/>
# hidden: true
# ---
#
# On a 2D plane, we place stones at some integer coordinate points.  Each
# coordinate point may have at most one stone.
#
# Now, a _move_ consists of removing a stone that shares a column or row with
# another stone on the grid.
#
# What is the largest possible number of moves we can make?
#
#
#
# **Example 1:**
#
#
#
#     Input: stones = [[0,0],[0,1],[1,0],[1,2],[2,1],[2,2]]
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: stones = [[0,0],[0,2],[1,1],[2,0],[2,2]]
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: stones = [[0,0]]
#     Output: 0
#
#
#
#
# **Note:**
#
#   1. `1 <= stones.length <= 1000`
#   2. `0 <= stones[i][j] < 10000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export readjats

import EzXML

function readjats(filepath::String)
    article = root(readxml(filepath))
    nodename(article) == "article" || throw("No article.")

    journal = findfirst(article, "front/journal-meta/journal-title-group/journal-title") |> nodecontent
    title = findfirst(article, "front/article-meta/title-group/article-title") |> nodecontent
    search!(title, pdf)
end

function aaa()
    nodename(xml_article) == "article" || throw("No article.")
    countelements(xml_article) < 3 && throw("#xml elements < 3")

    article = Tree("article")
    xml_front = findfirst(xml_article, "front")
    push!(article, parse_front(xml_front))

    body = findfirst(xml_article, "body")
    push!(article, parse_body(body))

    back = find(xml_article, "back")
    if !isempty(back)
        push!(article, parse_back(back[1]))
    end

    push!(article, Tree("floats-group"))
    append!(article[end], findfloats(article))
    floats = find(xml_article, "floats-group")
    if !isempty(floats)
        append!(article[end], parse_body(floats[1]).children)
    end
    isempty(article[end]) && deleteat!(article,length(article)) # no floats

    maths = findall(article, "math")
    for i = 1:length(maths)
        math = maths[i]
        mathml = root(parsexml(toxml(math)))
        normalize_mathml!(mathml)
        replace!(math, convert(Tree,mathml))
    end

    # tokenize_word!(article)
    postprocess!(article)
    nonrecursive!(article)
    article
end

function align(pdftxt::String, xmlfile::String)
    pdchars = readpdftxt(pdftxt)
    #pdchars = filter(istext, pdchars)
    #iddict = Dict{String,Int}()
    #pdids = map(x -> get!(iddict,x.c,length(iddict)+1), pdchars)

    xml = readjats(xmlfile)
    xmltree = xmltree[findfirst(c -> c.name == "body", xmltree.children)]
    tokenize!(xmltree)
    xmlchars = findall(isempty, xmltree)
    xmlids = map(c -> get!(iddict,string(c.name),length(iddict)+1), xmlchars)
    pairs = lcsmatch(pdids, xmlids)
end


text:
using ImplicitEquations
using Test

f(x,y) = y-x
g(x,y) = y+x

ImplicitEquations.GRAPH(f ⩵ 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≦ 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≪  g, -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(f ≶  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≷  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(Neq(f, g), -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(Eq(f, g), -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ⩵ g, -5, 5, -5, 5, 2^4, 2^4)


text:
function __init__reddit()
    DEPNAME = "Reddit"
    LINK = "http://snap.stanford.edu/graphsage/reddit.zip"
    DOCS = "http://snap.stanford.edu/graphsage/"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME Dataset
        Website: $DOCS
        Download Link: $LINK
        """,
        LINK,
        "25337a21540cd373e4cee3751e6600324ab6a7377ef3966bb49f57412a17ed02",
        post_fetch_method=unpack
    ))
end

"""
    Reddit(; full=true, dir=nothing)

The Reddit dataset was introduced in Ref [1].
It is a graph dataset of Reddit posts made in the month of September, 2014.
The dataset contains a single post-to-post graph, connecting posts if the same user comments on both.
The node label in this case is one of the 41 communities, or “subreddit”s, that a post belongs to.
This dataset contains 232,965 posts.
The first 20 days are used for training and the remaining days for testing (with 30% used for validation).
Each node is represented by a 602 word vector.

Use `full=false` to lod only a subsample of the dataset.

# References
[1]: [Inductive Representation Learning on Large Graphs](https://arxiv.org/abs/1706.02216)
[2]: [Benchmarks on the Reddit Dataset](https://paperswithcode.com/dataset/reddit)
"""
struct Reddit <: AbstractDataset
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
end

function Reddit(; full=true, dir=nothing)
    DATAFILES =  [
        "reddit-G.json", "reddit-G_full.json", "reddit-adjlist.txt",
        "reddit-class_map.json", "reddit-feats.npy", "reddit-id_map.json"
        ]
    DATA = joinpath.("reddit", DATAFILES)
    DEPNAME = "Reddit"

    if full
        graph_json = datafile(DEPNAME, DATA[2], dir)
    else
        graph_json = datafile(DEPNAME, DATA[1], dir)
    end

    class_map_json = datafile(DEPNAME, DATA[4], dir)
    id_map_json = datafile(DEPNAME, DATA[6], dir)
    feat_path = datafile(DEPNAME, DATA[5], dir)

    # Read the json files
    graph = open(JSON3.read, graph_json)
    class_map = open(JSON3.read, class_map_json)
    id_map = open(JSON3.read, id_map_json)

    # Metadata
    directed = graph["directed"]
    multigraph = graph["multigraph"]
    links = graph["links"]
    nodes = graph["nodes"]
    num_edges = directed ? length(links) : length(links) * 2
    num_nodes = length(nodes)
    num_graphs = length(graph["graph"]) # should be zero

    # edges
    s = get.(links, "source", nothing) .+ 1
    t = get.(links, "target", nothing) .+ 1
    if !directed
        s, t = [s; t], [t; s]
    end

    # labels
    node_keys = get.(nodes, "id", nothing)
    node_idx = [id_map[key] for key in node_keys] .+ 1

    sort_order = sortperm(node_idx)
    node_idx = node_idx[sort_order]
    labels = [class_map[key] for key in  node_keys][sort_order]
    @assert length(node_idx) == length(labels)

    # features
    features = npzread(feat_path)[node_idx, :]
    features = permutedims(features, (2, 1))

    # split
    test_mask =  get.(nodes, "test", nothing)[sort_order]
    val_mask = get.(nodes, "val", nothing)[sort_order]
    # A node should not be both test and validation
    @assert sum(val_mask .& test_mask) == 0
    train_mask = nor.(test_mask, val_mask)

    train_idx = node_idx[train_mask]
    test_idx = node_idx[test_mask]
    val_idx = node_idx[val_mask]

    split = Dict(
        "train" => train_idx,
        "test" => test_idx,
        "val"  => val_idx
    )

    metadata = Dict{String, Any}("directed" => directed, "multigraph" => multigraph,
                "num_graphs" => num_graphs, "num_edges" => num_edges, "num_nodes" => num_nodes,
                "split" => split)
    g = Graph(
        num_edges=num_edges, num_nodes=num_nodes,
        edge_index=(s, t),
        node_data= (; labels, features)
    )
    return Reddit(metadata, [g])
end

Base.length(d::Reddit) = length(d.graphs)
Base.getindex(d::Reddit) = d.graphs
Base.getindex(d::Reddit, i) = getindex(d.graphs, i)


text:
export cleanup!

@generated function cleanup!(f, args)
    if args<:Tuple
        quote
            try
                f(args...)
            finally
                for arg in args
                    release!(args)
                end
            end
        end
    else
        quote
            try
                f(args)
            finally
                release!(args)
            end
        end
    end
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# Generic functions (on top of the type-specific implementations)

## Basic functions

Base.eltype{T<:Real}(a::AbstractPDMat{T}) = T
Base.ndims(a::AbstractPDMat) = 2
Base.size(a::AbstractPDMat) = (dim(a), dim(a))
Base.size(a::AbstractPDMat, i::Integer) = 1 <= i <= 2 ? dim(a) : 1
Base.length(a::AbstractPDMat) = abs2(dim(a))

## arithmetics

pdadd!{T<:Real}(r::Matrix, a::Matrix, b::AbstractPDMat{T}) = pdadd!(r, a, b, one(T))

pdadd!(a::Matrix, b::AbstractPDMat, c) = pdadd!(a, a, b, c)
pdadd!{T<:Real}(a::Matrix, b::AbstractPDMat{T}) = pdadd!(a, a, b, one(T))

pdadd{T<:Real, S<:Real, R<:Real}(a::Matrix{T}, b::AbstractPDMat{S}, c::R) = pdadd!(similar(a, promote_type(T, S, R)), a, b, c)
pdadd{T<:Real, S<:Real}(a::Matrix{T}, b::AbstractPDMat{S}) = pdadd!(similar(a, promote_type(T, S)), a, b, one(T))

+(a::Matrix, b::AbstractPDMat) = pdadd(a, b)
+(a::AbstractPDMat, b::Matrix) = pdadd(b, a)

*{T<:Real}(a::AbstractPDMat, c::T) = a * c
*{T<:Real}(c::T, a::AbstractPDMat) = a * c
/{T<:Real}(a::AbstractPDMat, c::T) = a * inv(c)

## whiten and unwhiten

whiten!(a::AbstractPDMat, x::StridedVecOrMat) = whiten!(x, a, x)
unwhiten!(a::AbstractPDMat, x::StridedVecOrMat) = unwhiten!(x, a, x)

whiten(a::AbstractPDMat, x::StridedVecOrMat) = whiten!(similar(x), a, x)
unwhiten(a::AbstractPDMat, x::StridedVecOrMat) = unwhiten!(similar(x), a, x)

## quad

function quad{T<:Real, S<:Real}(a::AbstractPDMat{T}, x::StridedMatrix{S})
    @check_argdims dim(a) == size(x, 1)
    quad!(Array(promote_type(T, S), size(x,2)), a, x)
end

function invquad{T<:Real, S<:Real}(a::AbstractPDMat{T}, x::StridedMatrix{S})
    @check_argdims dim(a) == size(x, 1)
    invquad!(Array(promote_type(T, S), size(x,2)), a, x)
end


text:
using Memcache
using Test

function test_begin()
    mc = MemcacheClient()
    mc.debug = true

    st = stats(mc)
    @test !isempty(st["version"])

    st = version(mc)
    @test !isempty(st)

    flush_all(mc)

    st = stats(mc, "settings")
    @test st["tcpport"] == "11211"

    mc
end

function test_end(mc)
    flush_all(mc)

    st = stats(mc)
    @test parse(Int, st["get_hits"]) > 0

    @test length(stats(mc, "items")) >= 0
    @test length(stats(mc, "sizes")) >= 0
    @test parse(Int, stats(mc, "slabs")["active_slabs"]) >= 0

    # commented out to avoid julia issue #5793
    #close(mc)
    nothing
end

function test_set_get(mc)
    set(mc, "key1", "val1")
    @test get(mc, "key1") == "val1"

    set(mc, "key2", "val2")
    res = get(mc, "key1", "key2")
    @test res["key1"] == "val1"
    @test res["key2"] == "val2"

    delete(mc, "key1")
    res = get(mc, "key1", "key2")
    @test collect(keys(res)) == ["key2"]
end

function test_touch(mc)
    set(mc, "touch_key1", "touch_val1", exptime=2)
    set(mc, "touch_key2", "touch_val2", exptime=1000)
    sleep(3)
    res = get(mc, "touch_key1", "touch_key2")
    @test collect(keys(res)) == ["touch_key2"]

    touch(mc, "touch_key2", 2)
    set(mc, "touch_key1", "touch_val1", exptime=1000)
    sleep(3)
    res = get(mc, "touch_key1", "touch_key2")
    @test collect(keys(res)) == ["touch_key1"]
end

function test_incr_decr(mc)
    set(mc, "num1", 1)
    incr(mc, "num1", 5)
    @test get(mc, "num1") == 6

    decr(mc, "num1", 2)
    @test get(mc, "num1") == 4
end

function test_cas(mc)
    set(mc, "cas_key", 1)
    res = get(mc, "cas_key", cas=true)
    val,casval = res["cas_key"]
    @test val == 1

    cas(mc, "cas_key", 2, casval)
    res = get(mc, "cas_key", cas=true)
    val,casval2 = res["cas_key"]
    @test val == 2
    @test casval != casval2

    @test false == try cas(mc, "cas_key", 3, casval); true; catch; false; end
end

mutable struct mytype
    i::Int
    s::AbstractString
    k::Dict
end

function test_julia_type(mc)
    t = mytype(10, "hello", Dict(1=>"A", 2=>"B"))
    set(mc, "jul", t)
    t1 = get(mc, "jul")
    @test t.i == t1.i
    @test t.s == t1.s
    @test t.k == t1.k
end

mc = test_begin()
println("began test")

test_set_get(mc)
println("tested set get")

test_touch(mc)
println("tested touch")

test_incr_decr(mc)
println("tested incr decr")

test_cas(mc)
println("tested cas")

test_julia_type(mc)
println("tested julia type")

test_end(mc)
println("end of test")

exit()


text:
using TimeAxisArrays
using Base.Test
using Base.Dates

Base.isequal(A1::TimeAxisArray, A2::TimeAxisArray) =
    (collect(timestamps(A1)) == collect(timestamps(A2))) && isequal(A1.data, A2.data)

function isapproxequal(A1::TimeAxisArray, A2::TimeAxisArray)
    notnan = find(!isnan(A1))
    return (collect(timestamps(A1)) == collect(timestamps(A2))) &&
           (notnan == find(!isnan(A2))) &&
           isapprox(A1.data[notnan], A2.data[notnan])
end #isapproxequal

@test isa(TimeAxisArray(randn(5,2,2), 1:5, [:A, :B], [:High, :Low]), RegularTimeAxisArray)
@test isa(TimeAxisArray(randn(5,2,2), collect(1:5), [:A, :B], [:High, :Low]), TimeAxisArray)

n = now()

A1d1 = TimeAxisArray(randn(11), n-Hour(2):Hour(1):n+Hour(8))
A1d2 = TimeAxisArray(randn(5), n-Hour(2):Hour(1):n+Hour(2))
@test isa(A1d1, RegularTimeAxisArray)

data = randn(5,2)
A2d1 = TimeAxisArray(data, n-Hour(2):Hour(1):n+Hour(2),
                         Axis{:Columns}([:A, :B]))
A2d2 = TimeAxisArray(data, n-Hour(2):Hour(1):n+Hour(2), [:A, :B])
@test isa(A2d1, RegularTimeAxisArray)
@test A2d1 == A2d2

A3d = TimeAxisArray(randn(7,2,3), n-Day(1):Hour(8):n+Day(1),
              Axis{:Columns}([:A, :B]),
              Axis{:Pages}([:I, :II, :III]))
@test isa(A3d, RegularTimeAxisArray)

A4d = TimeAxisArray(randn(49,2,5,3), n-Day(1):Hour(1):n+Day(1),
              Axis{:Columns}([:A, :B]),
              Axis{:Pages}([:I, :II, :III, :IV, :V]),
              Axis{:Superpages}([:one, :two, :three]))
@test isa(A4d, RegularTimeAxisArray)

include("operations.jl")
include("io.jl")

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
############################################################################################
@testset "Flattening - Probabilistic Models - flatten type consistency" begin
    @test _checkparams(val_dist)
    _params = merge(val_dist, val_dist_nested)
    ## Iterate trough all Params in TestHelper.jl file
    for sym in eachindex(_params)
        # println(sym)
        param = _params[sym]
        θ = _get_val(param)
        constraint = _get_constraint(param)
        ## Check all flatten possibilities
        for unflat in unflattenmethods
            for flat in flattentypes
                for floattypes in outputtypes
                  # println(unflat, " ", flat, " ", floattypes)
                    flatdefault = FlattenDefault(;
                        output = floattypes,
                        flattentype = flat,
                    )
                    _flatten, _unflatten = construct_flatten(flatdefault, unflat, constraint, θ)
                    θ_flat = _flatten(θ)
                    θ_unflat = _unflatten(θ_flat)

                    #!NOTE Do not test if FlattenContinuous and empty Integer Param struct is evaluated
                    if θ_flat isa Vector{T} where {T<:AbstractFloat}
                        if unflat isa UnflattenStrict
                            @test eltype(θ_flat) == floattypes
                            @test typeof(θ_unflat) == typeof(θ)
                        end
                        if unflat isa UnflattenFlexible && θ_unflat isa Array{T} where {T<:AbstractFloat}
                            @test eltype(θ_unflat) == eltype(θ_flat)
                        end
                    end
                end
            end
        end
        @test _checkparams(param)
        @test _checkfinite(θ)
        @test _checkprior(constraint)
        _anyparam(_checksampleable(constraint))
        @test _allparam(_checksampleable(constraint))
        _transform, _inversetransform = construct_transform(constraint, θ)
        ## Type Check 1 - Flatten/Unflatten
        _flatten, _unflatten = construct_flatten(FlattenDefault(), UnflattenStrict(), constraint, θ)
        _flattenAD, _unflattenAD = construct_flatten(FlattenDefault(), UnflattenFlexible(), constraint, θ)
        θ_flat = _flatten(θ)
        θ_unflat = _unflatten(θ_flat)
        @test typeof(_unflatten(θ_flat)) == typeof(θ)
        ## Type Check 2 - Constrain/Unconstrain
        θ_unconstrained = unconstrain(_transform, θ)
        θ_constrained = constrain(_inversetransform, θ_unconstrained)
        @test typeof(θ) == typeof(θ_constrained)
        ## Type Check 3 - size of flatten(constrained) == flatten(unconstrained) for current Bijectors
        _θ_vec1 = _flatten(θ)
        _θ_vec2 = _flatten(θ_unconstrained)
        @test length(_θ_vec1) == length(_θ_vec2)
        @test typeof(_unflatten(_θ_vec1)) == typeof(_unflatten(_θ_vec2))
        ## If applicable, check if gradients for supported AD frameworks can be computed
        if length(θ_flat) > 0
            function check_AD_closure(constraint, val)
                reconstruct = ReConstructor(constraint, val)
                bij, bij⁻¹ = construct_transform(constraint, val)
                function check_AD(θₜ::AbstractVector{T}) where {T<:Real}
                    θ = unflattenAD(reconstruct, θₜ)
                    θ = constrain(bij⁻¹, θ)
                    return log_prior(constraint, θ) + log_abs_det_jac(bij, θ)
                end
            end
            check_AD = check_AD_closure(constraint, θ)
            check_AD(θ_flat)
            grad_mod_fd = ForwardDiff.gradient(check_AD, θ_flat)
            grad_mod_rd = ReverseDiff.gradient(check_AD, θ_flat)
            grad_mod_zy = Zygote.gradient(check_AD, θ_flat)[1]
            @test sum(abs.(grad_mod_fd - grad_mod_rd)) ≈ 0 atol = _TOL
            @test sum(abs.(grad_mod_fd - grad_mod_zy)) ≈ 0 atol = _TOL
        end
    end
end

############################################################################################
@testset "Flattening - Non-Probabilistic Models - flatten type consistency" begin
    @test _checkparams(val_constrained)
    _params = merge(val_constrained)
    ## Iterate trough all Params in TestHelper.jl file
    for sym in eachindex(_params)
        #println(sym)
        param = _params[sym]
        θ = _get_val(param)
        constraint = _get_constraint(param)
        for unflat in unflattenmethods
            for flat in flattentypes
                for floattypes in outputtypes
                  # println(unflat, " ", flat, " ", floattypes)
                    flatdefault = FlattenDefault(;
                        output = floattypes,
                        flattentype = flat,
                    )
                    _flatten, _unflatten = construct_flatten(flatdefault, unflat, constraint, θ)
                    θ_flat = _flatten(θ)
                    θ_unflat = _unflatten(θ_flat)

                    #!NOTE Do not test if FlattenContinuous and empty Integer Param struct is evaluated
                    if !(constraint isa Fixed) && θ_flat isa Vector{T} where {T<:AbstractFloat}
                        if unflat isa UnflattenStrict
                            @test eltype(θ_flat) == floattypes
                            @test typeof(θ_unflat) == typeof(θ)
                        end
                        if unflat isa UnflattenFlexible && θ_unflat isa Array{T} where {T<:AbstractFloat}
                            @test eltype(θ_unflat) == eltype(θ_flat)
                        end
                    end
                end
            end
        end
        @test _checkparams(param)
        @test _checkfinite(θ)
        @test !_checkprior(constraint)
        _anyparam(_checksampleable(constraint))
        @test !_allparam(_checksampleable(constraint))
        _transform, _inversetransform = construct_transform(constraint, θ)
        ## Type Check 1 - Flatten/Unflatten
        _flatten, _unflatten = construct_flatten(FlattenDefault(), UnflattenStrict(), constraint, θ)
        _flattenAD, _unflattenAD = construct_flatten(FlattenDefault(), UnflattenFlexible(), constraint, θ)
        θ_flat = _flatten(θ)
        θ_unflat = _unflatten(θ_flat)
        @test typeof(_unflatten(θ_flat)) == typeof(θ)
        ## Type Check 2 - Constrain/Unconstrain
        θ_unconstrained = unconstrain(_transform, θ)
        θ_constrained = constrain(_inversetransform, θ_unconstrained)
        @test typeof(θ) == typeof(θ_constrained)
        ## Type Check 3 - size of flatten(constrained) == flatten(unconstrained) for current Bijectors
        _θ_vec1 = _flatten(θ)
        _θ_vec2 = _flatten(θ_unconstrained)
        @test length(_θ_vec1) == length(_θ_vec2)
        @test typeof(_unflatten(_θ_vec1)) == typeof(_unflatten(_θ_vec2))
        ## If applicable, check if gradients for supported AD frameworks can be computed
        if length(θ_flat) > 0
            function check_AD_closure(constraint, val)
                reconstruct = ReConstructor(constraint, val)
                bij, bij⁻¹ = construct_transform(constraint, val)
                function check_AD(θₜ::AbstractVector{T}) where {T<:Real}
                    θ = unflattenAD(reconstruct, θₜ)
                    θ = constrain(bij⁻¹, θ)
                    return log_prior(constraint, θ) + log_abs_det_jac(bij, θ)
                end
            end
            check_AD = check_AD_closure(constraint, θ)
            check_AD(θ_flat)
            grad_mod_fd = ForwardDiff.gradient(check_AD, θ_flat)
            grad_mod_rd = ReverseDiff.gradient(check_AD, θ_flat)
            #!NOTE: Zygote would just record "Nothing" as gradient for Fixed/Unconstrained without a functor
            #grad_mod_zy = Zygote.gradient(check_AD, θ_flat)[1]
            @test sum(abs.(grad_mod_fd - grad_mod_rd)) ≈ 0 atol = _TOL
            #@test sum(abs.(grad_mod_fd - grad_mod_zy)) ≈ 0 atol = _TOL
        end
    end
end

############################################################################################
@testset "Flattening - flatten type consistency" begin
    _temp = Param(Float32(2.0), Distributions.Normal{Float32}(Float32(1.0), Float32(2.0))) #In flatten case, everything would be returned as Float64 when unflattening
    _flatten, _unflatten = construct_flatten(FlattenDefault(), UnflattenStrict(), _temp.constraint, _temp.val)
    θ_flat = _flatten(_temp.val)
    θ_unflat = _unflatten(θ_flat)
    @test typeof(_temp.val) == typeof(θ_unflat)
end

############################################################################################
@testset "Flattening - discrete parameter are not flattened" begin
    _temp = (
        a=Param(1.0, Distributions.Normal()),
        b=Param(2, Distributions.Categorical(40)),
        c=Param([3.0, 4.0], [Distributions.Normal(), Distributions.Normal(2.0, 3.0)]),
        d=Param([5, 6], [Distributions.Categorical(10), Distributions.Categorical(10)]),
    )
    _temp_val = _get_val(_temp)
    _temp_constraint = _get_constraint(_temp)
    _flatten, _unflatten = construct_flatten(FlattenDefault(), UnflattenStrict(), _temp_constraint, _temp_val)
    θ_flat = _flatten(_temp_val)
    θ_unflat = _unflatten(θ_flat)
    @test _allparam(_checksampleable(_temp_constraint))
    @test _anyparam(_checksampleable(_temp_constraint))
    @test _checkprior(_temp_constraint)

    @test length(θ_flat) == 3
    @test typeof(_unflatten(θ_flat)) == typeof(_temp_val)
end

############################################################################################
@testset "Flattening - discrete parameter flattened" begin
    df_all = FlattenDefault(Float64, FlattenAll())
    _temp = (
        a=Param(1.0, Distributions.Normal()),
        b=Param(2, Distributions.Categorical(40)),
        c=Param([3.0, 4.0], [Distributions.Normal(), Distributions.Normal(2.0, 3.0)]),
        d=Param([5, 6], [Distributions.Categorical(10), Distributions.Categorical(10)]),
    )
    _temp_val = _get_val(_temp)
    _temp_constraint = _get_constraint(_temp)
    _flatten, _unflatten = construct_flatten(df_all, _temp_constraint, _temp_val)
    _val = _flatten(_temp_val)
    θ_unflat = _unflatten(_val)
    @test length(_val) == 6
    @test typeof(θ_unflat) == typeof(_temp_val)
end


text:
struct RNNWrapper{T}
  rnn::T
end
Flux.@functor RNNWrapper

# Need to specialize for RNNWrapper.
fw(r::RNNWrapper, X::Vector{<:AbstractArray}) = begin
  Flux.reset!(r.rnn)
  [r.rnn(x) for x in X]
end

fw(r::RNNWrapper, X) = begin
  Flux.reset!(r.rnn)
  r.rnn(X)
end

fwbw(r::RNNWrapper, ps, X::Vector{<:AbstractArray}) = gradient(ps) do
  y = fw(r, X)
  sum(sum(y))
end

pb(r::RNNWrapper, ps, X::Vector{<:AbstractArray}) = pullback(ps) do
  y = fw(r, X)
  sum(sum(y))
end

function rnn_benchmark_sweep(data_creator::Function, rnn_type)
  for n in [2, 20, 200, 1000], ts in [1, 4, 16, 64]
    x, x_n = data_creator(n, ts)
    model = RNNWrapper(rnn_type(n, n))

    println("$rnn_type $x_n CPU n=$n, ts=$ts")
    run_benchmark(model, x, cuda=false)

    println("$rnn_type $x_n CUDA n=$n, ts=$ts")
    try
      run_benchmark(model, x, cuda=true)
    catch ex
      @show typeof(ex)
      if ex isa OutOfGPUMemoryError
        @warn "Not enough GPU memory to run test"
      else
        rethrow(ex)
      end
    end
  end
end

for rnn_type in [Flux.RNN, Flux.GRU, Flux.LSTM]
  rnn_benchmark_sweep(rnn_type) do n, ts
    [randn(Float32, n, n) for _ in 1:ts], "Vec"
  end
end

for rnn_type in [Flux.RNN, Flux.GRU, Flux.LSTM]
  rnn_benchmark_sweep(rnn_type) do n, ts
    randn(Float32, n, n, ts), "Block"
  end
end


text:
using ODEFilters
using ODEFilters: stack
using DifferentialEquations
using Statistics
using Plots

# ODE definition as in DifferentialEquations.jl
function fitz(u, p, t)
    a, b, c = p
    return [c*(u[1] - u[1]^3/3 + u[2])
            -(1/c)*(u[1] -  a - b*u[2])]
end
u0 = [-1.0; 1.0]
tspan = (0., 20.)
p = (0.2,0.2,3.0)
prob = ODEProblem(fitz, u0, tspan, p)

appxsol = solve(prob, abstol=1e-12, reltol=1e-12)

############################################################################################
# Solve with integrator interface and make animation
integ = init(prob, EK0(order=1), adaptive=false, dt=7e-2)

anim = @animate for i in 1:(prob.tspan[2] / integ.dt)
    step!(integ)

    plot(integ.sol; color=[1 2], ribbon_width=3,
         xlims=prob.tspan, ylims=(-2.5, 2.5),
         ylabel="u(t)", xlabel="t",
         label=["u₁(t)" "u₂(t)"],
         )

    times = integ.t:0.1:prob.tspan[2]
    post = integ.sol(times)
    plot!(times, stack(mean(post)),
          ribbon=3stack(std(post)),
          linestyle=:dot,
          color=[1 2],
          fillalpha=0.2,
          label="",
         )
    scatter!(integ.sol.t, stack(integ.sol.u), color=[1 2], label="")
    plot!(appxsol, linestyle=:dash, color=:black, label="")
    vline!([integ.t], label="", color=:black)
end

gif(anim, "./fitzhughnagumo_solve.gif")


text:
using Plots

p1=plot(sin)
p2=plot!(p1, cos)
p1 |> display
p2 |> display
plot([sin,cos]) |> display
plot(rand(10),rand(10)) |> display
scatter(rand(10),rand(10)) |> display


text:
struct ClusterWolff end
"""
Wolff single cluster update

#Arguments:
* `spins` : array of spin values
* `nearest_neighbors` : function that returns a list of nearest neighbors to index i
* `beta` : inverse temperature
* `rng`  : random number generator

Ref: U. Wolff, Phys. Rev. Lett. 62, 361 (1989)

consider for Potts:
https://stanford.edu/~kurinsky/ClassWork/Physics271_Final_Paper.pdf

consider BlumeCapel implementation

todo:
return dE?
how to specify function type for performance?
"""
function update(alg::ClusterWolff, rng::AbstractRNG, spins::Vector{Int}, nearest_neighbors::Function, betaJ::Float64)
    N = length(spins)

    # get random spin and flip it
    index_i = rand(rng, 1:N)
    s_i = spins[index_i]
    spins[index_i] *= -1

    # probability (per bond!)
    p = 1.0 - exp(-2.0 * betaJ)

    # check all bonds and already flip cluster
    proposed = copy(nearest_neighbors(index_i))
    while !isempty(proposed)
        # check proposed sites for alignment with base spin s_i and add to cluster with probability p
        index_j = pop!(proposed)
        # no need to check if site is already in cluster, because spin flip disqualifies it already
        if spins[index_j] == s_i && rand(rng) < p
            spins[index_j] *= -1
            append!(proposed, nearest_neighbors(index_j))
        end
    end
end


text:
using Revise
using ADSeismic
using ADCME
using Optim
using PyPlot
using MAT
using DelimitedFiles
using JLD2
using PyCall
using ADOPT
using Dates
matplotlib.use("Agg")
close("all")

ADCME.options.customop.verbose = false

data_dir = "data/elastic"
method = "FWI/elastic"

figure_dir = string("figure/",method)
result_dir = string("result/",method)
loss_file = joinpath(result_dir, "loss_$(Dates.now()).txt")

check_path(dir) = !ispath(dir) ? mkpath(dir) : nothing
check_path(figure_dir)
check_path(result_dir)

################### Generate synthetic data #####################
reset_default_graph()
mpi_init()

# model_name = "models/marmousi2-model-true-large.mat"
# model_name = "models/marmousi2-model-true.mat"
# model_name = "models/BP-model-true.mat"
model_name = "models/marmousi2-model-smooth.mat"

param = load_params(model_name, "MPIElastic")
compute_PML_Params!(param)

vp0 = matread(model_name)["vp"]
vs0 = matread(model_name)["vs"]
# rho0 = matread(model_name)["rho"]
##
# vp0 = matread("models/marmousi2-model-true.mat")["vp"]
vs0 = matread("models/marmousi2-model-true.mat")["vs"]
rho0 = matread("models/marmousi2-model-true.mat")["rho"]

@info "loading data..."
@load joinpath(data_dir, "data$(mpi_rank())-$(mpi_size()).jld2") Rs

vp_ = Variable(vp0)
var_size = size(vp_)
ps_ratio_ = Variable(1.75)
vp = mpi_bcast(vp_)
ps_ratio = mpi_bcast(ps_ratio_, deps=sum(vp))

if exists(matopen(model_name), "mask")
    @info "using mask..."
    mask = matread(model_name)["mask"]
    vp = mask .* vp + (1.0.-mask) .* vp0
end

ps_ratio = tanh(ps_ratio - 1.75)*0.5 + 1.75
# vs = vp / ps_ratio
vs = constant(vs0)
rho = constant(rho0)

vp = extract_local_patch(param, vp, tag = 0)
vs = extract_local_patch(param, vs, deps = sum(vp), tag = 810000)
rho = extract_local_patch(param, rho, deps = sum(vs), tag = 820000)

λ, μ, ρ = compute_lame_parameters(vp, vs, rho)

src = load_elastic_source(model_name, use_mpi=true, param=param)
receiver = load_elastic_receiver(model_name, use_mpi=true, param=param)

Rs_ = []
Vs_ = []

local_loss = constant(0.0)
# for i = 1:length(src)
for i in [3,6]

  rcvv = Rs[i]
  src_ = src[i]
  receiver_ = receiver[i]

  propagator = MPIElasticPropagatorSolver(param, src_, ρ, λ, μ; tag_offset = i*100000, dep = length(Vs_)==0 ? nothing : sum(Vs_[end]))
  MPISimulatedObservation!(propagator, receiver_)

  push!(Rs_, receiver_.rcvv)
  push!(Vs_, propagator.vy)

  if !isnothing(rcvv)
      global local_loss += sum((receiver_.rcvv - rcvv)^2)
  end

end

local_loss += sum(Vs_[end]) * 1e-20

loss = mpi_sum(local_loss) * 1e7
# g = gradients(loss, [vp_, ps_ratio_])
g = gradients(loss, vp_)

py"""
from tensorflow.core.protobuf import rewriter_config_pb2
import tensorflow as tf
config = tf.ConfigProto()
off = rewriter_config_pb2.RewriterConfig.OFF
config.graph_options.rewrite_options.memory_optimization  = off
"""
session_conf = py"config"
sess = Session(config = session_conf); init(sess)
#sess = Session(); init(sess)

# @info "running loss"
# @time x = run(sess, loss)
# @info "running gradients"
# @time x = run(sess, g)

#using ADOPT
@info "starting..."
@info "initial local loss = $(run(sess, local_loss)) at rank = $(mpi_rank())"

function calculate_loss(x)
    p = reshape(x[1:prod(var_size)], var_size)
    L = run(sess, loss, vp_=>p)

    # ps =  x[prod(var_size)+1]
    # L = run(sess, loss, vp_=>p, ps_ratio_=>ps)
    return L
end

iter = 0
iter_result = Array{Float64, 2}[]
function calculate_gradients(G, x)
    p = reshape(x[1:prod(var_size)], var_size)

    G[:] = run(sess, g, vp_=>p)

    # ps =  x[prod(var_size)+1]
    # ps = tanh(ps - 1.75)*0.5 + 1.75
    # s = p ./ ps
    # if mpi_rank() == 0
    #   @info ps, maximum(p), maximum(s)
    # end
    # G_tmp = run(sess, g, vp_=>p, ps_ratio_=>ps)
    # G[:] = append!(reshape(G_tmp[1], prod(var_size)), G_tmp[2])

    if mpi_rank() == 0
      @info "iter = ", iter
      if iter % 1 == 0
          clf()
          subplot(211)
          pcolormesh(p', shading="auto", cmap="jet")
          axis("scaled")
          colorbar()
          xlabel("x (km)")
          ylabel("z (km)")
          gca().invert_yaxis()
          title("Iteration = $iter")
        #   savefig(joinpath(figure_dir, "inv_p_$(lpad(iter,5,"0")).png"), bbox_inches="tight")
        #   clf()

          # subplot(212)
          # pcolormesh(s', shading="auto", cmap="jet")
          # axis("scaled")
          # colorbar()
          # xlabel("x (km)")
          # ylabel("z (km)")
          # gca().invert_yaxis()

        #   subplot(312)
        #   pcolormesh(r', cmap="jet")
        #   gca().invert_yaxis()
        #   colorbar()
          savefig(joinpath(figure_dir, "inv_$(lpad(iter,5,"0")).png"), bbox_inches="tight")
          writedlm(joinpath(result_dir, "inv_vp_$(lpad(iter,5,"0")).txt"), p)
          # writedlm(joinpath(result_dir, "inv_vs_$(lpad(iter,5,"0")).txt"), s)
      end
      global iter += 1
    end
end

losses = Float64[]
function step_callback(x)
    @info "loss = $x"
    push!(losses, x)
end

# initial_x = run(sess, [vp_, ps_ratio_])
# # initial_x = cat(initial_x..., dims=2)[:]
# initial_x = append!(reshape(initial_x[1], prod(var_size)), initial_x[2])
initial_x = run(sess, vp_)[:]
options = Options()
# result = ADOPT.mpi_optimize(calculate_loss, calculate_gradients, initial_x, LBFGS(), options; step_callback = step_callback)
result = ADCME.mpi_optimize(calculate_loss, calculate_gradients, initial_x)#, LBFGS(), options; step_callback = step_callback)

if mpi_size()>1
    mpi_finalize()
end


text:
using FilesLH, Test

function make_dir_test()
    @testset "Make dir" begin
        d = joinpath(test_file_dir(), "make_dir_test");
        isdir(d)  &&  rm(d);
        make_dir(d);
        @test isdir(d);
    end
end

function right_dirs_test()
    @testset "Parent dirs" begin
        @test right_dirs("abc", 1) == "abc"
        @test right_dirs("/abc", 1) == "abc"
        @test right_dirs("/abc", 2) == "/abc"
        @test right_dirs("abc/def", 1) == "def"
        @test right_dirs("abc/d/ef/", 2) == "d/ef"
        @test right_dirs("abc/d/ef/", 4) == "abc/d/ef"
    end
end

function files_in_dir_test()
    @testset "files in dir" begin
        d = @__DIR__;
        fList = files_in_dir(d);
        @test "directories_test.jl" ∈ fList
        @test "coverage/Project.toml" ∈ fList

        @test isempty(files_not_in_dir2(d, d))

        d2 = normpath(d, "..");
        @test isequal(files_not_in_dir2(d, d2), fList)
    end
end

# This test just the filtering of excluded strings.
# Should also test `dir_diff_report` directly +++
function filter_excludes_test()
    @testset "filter excludes" begin
        v = ["a", "abc", "abcdef", "cdefg"];
        v2 = copy(v);
        FilesLH.filter_excludes!(v2, "cb");
        @test v == v2

        FilesLH.filter_excludes!(v2, "ab");
        @test v2 == v[[1,4]];

        v2 = copy(v);
        FilesLH.filter_excludes!(v2, ["ab", "ba"]);
        @test v2 == v[[1,4]];
	end
end

function find_common_base_dir_test()
    @testset "Common base dir" begin
        d1 = "/a/bc/def/ghk";
        d2 = "/a/bc/defg/ghk";
        cbd = find_common_base_dir(d1, d2);
        @test cbd == "/a/bc"

        @test find_common_base_dir([d1, d2]) == cbd
    end
end

@testset "Directories" begin
    tgDir = joinpath(test_file_dir(), "clear_dir_test");
    isdir(tgDir) || mkpath(tgDir);
    tgFile = joinpath(tgDir, "test.txt");
    open(tgFile, "w") do io
        println(io, "test");
    end
    @test !is_dir_empty(tgDir)

    clear_directory(tgDir; askFirst = false);
    @test isdir(tgDir)
    @test is_dir_empty(tgDir)

    right_dirs_test();
    files_in_dir_test()
    find_common_base_dir_test()
    filter_excludes_test();
end

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