text:
# Dynamical Ising solver
#
# Here, the dynamical variables are regarded from the perspective of the
# feasible configurations. Therefore, the Ising states are {-1, 1}^N,
# and the period of the main functions is [-2, 2).

# The main functions are F(v) and f = F'(v)
#
# F(v) is defined by
#
# C_G(v) = \sum_{(m,n) \in E} F(v_m - v_n), when v_m ∈ {-1, 1}
#
# So that F(0) = 0 (the same partition), F(±2) = -1. Hence, the period above.

# TODO:
#   0. [DONE] Fix rounding
#   1. Make energy methods for correct energy evaluations
#      Support relaxed objective functions
#   2. Enable LUT-defined methods
#   3. Add the weight support
#   4. Implement logging
#   5. Fix the chain of types
#   6. Admit differenttial equations

module Dice

using Base:Integer
using Arpack
using Graphs
using SparseArrays

export Model,
    get_connected, get_initial,
    sine, triangular,
    cut, get_best_cut, get_best_configuration, extract_configuration,
    number_to_conf,
    propagate, roundup,
    test_branch, scan_vicinity, scan_for_best_configuration,
    conf_decay,
    local_search, local_twosearch

# The main type describing the Model
# For a compact description of simulation scenarios and controlling the
# module behavior
const silence_default = 3
mutable struct Model
    graph::SimpleGraph
    method::Function # obsolete, phase out
    # coupling function
    coupling::Function
    # method_energy::Function
    # Defines the magnitude of the exchange terms
    # default = 1/Nvertices
    scale::Float64
    # The magnitude of the anisotropy term
    # default = 0
    Ks::Float64
    # anisotropy function
    # default = coupling(x, -x)
    anisotropy::Function
    extended::Bool # Whether there is an extension (remove?)
    # the inverse level of verbosity
    # default = silence_default
    silence::Integer
    Model() = new()
end

# Several explicit (legacy?) constructors
Model(graph, coupling) =
    begin
        M = Model()
        M.graph = graph
        M.coupling = coupling
        M.method = M.coupling
        # The default interpretation of anisotropy
        M.anisotropy = (x) -> M.coupling(x, -x)
        M.scale = 1/nv(graph)
        M.Ks = 0
        M.silence = silence_default
        M.extended = false
        return M;
    end

Model(graph, coupling, scale) =
    begin
        M = Model(graph, coupling)
        M.scale = scale
        M.Ks = 0
        M.silence = silence_default
        M.extended = false
        return M;
    end

Model(graph, coupling, scale, Ks) =
    begin
        M = Model(graph, coupling, scale)
        M.Ks = Ks
        return M;
    end

############################################################
#
### Internal service functions
#
############################################################

function message(model::Model, out, importance=1)
    # must be replaced by a proper logging functionality
    if importance > model.silence
        println("$out ($importance)")
    end
end

const DEBUG = 1
function debug_msg(out)
    if DEBUG > 0
        println("$out (DEBUG = $DEBUG)")
    end
end

############################################################
#
### File support
#
############################################################

"""
Read graph in the "reduced" Matrixmarket format.

Usage:
    G = loadDumpedGraph(filename::String)::SimpleGraph

INPUT:
    filename::String name of the file

OUTPUT:
    SimpleGraph object

The reduced format:
    Header: lines starting with %
    |V| |E|
    u v w_{u,v}
"""
function loadDumpedGraph(filename::String)::SimpleGraph
    # TODO simple validation

    mtxFile = open(filename, "r")
    tokens = split(lowercase(readline(mtxFile)))
    # Skip commenting header
    while tokens[1][1] == '%'
        global tokens
        tokens = split(lowercase(readline(mtxFile)))
    end

    Nodes, Edges = map((x) -> parse(Int64, x), tokens)
    G = SimpleGraph(Nodes)

    for line in readlines(mtxFile)
        global tokens
        global countEdges
        tokens = split(lowercase(line))
        # we drop the weight for now
        u, v = map((x) -> parse(Int64, x), tokens[1:2])
        add_edge!(G, u, v)
    end
    close(mtxFile)

    return G
end

"""
    dumpGraph(Graph::SimpleGraph, filename::String)

Writes the sparce adjacency matrix of `Graph` to file `filename`.

NOTE: the current version writes the (0,1) matrix only.

The adjacency matrix is written in the reduced `MatrixMarket` format
format. The first line contains two numbers: `|V| |E|`
Next lines contain three numbers: `u v A_{u, v}`
where `u` and `v` are the nodes numbers and `A_{u, v}` is the edge weight.
"""
function dumpGraph(G, filename)
    weight = 1
    open(filename, "w") do out
        println(out, nv(G), " ", ne(G))
        for edge in edges(G)
            println(out, edge.src, " ", edge.dst, " ", weight)
        end
    end
end

############################################################
#
### A library of coupling functions (TODO: detach?)
#
############################################################

function none_method(v)
    return 0 * v  # for the type stability thing
end

function none_method(v1, v2)
    return 0 * v1  # for the type stability thing
end

const Pi2 = pi / 2
const Pi4 = pi / 4

"""
    cosine(v)

The coupling energy inducing the XY-model (rank-2 relaxation).
"""
function cosine(v)
    return cos.(Pi2 .* v)
end

function sine(v)
    return sin.(Pi2 .* v)
end

function sine(v1, v2)
    return sine(v1 - v2)
end

function linear(v1, v2)
    return v1 - v2
end

function dtri(v)
    # the derivative of the triangular function
    p = Int.(floor.(v ./ 2 .+ 1 / 2))
    parity = rem.(abs.(p), 2)
    return -2 .* parity .+ 1
end

function triangular(v)
    # Odd piece-wise linear defined within the period [-2, 2] by
    # phi(-2) = phi(0) = phi(2) = 0
    # - phi(-1) = phi(1) = 1
    # This works for both scalars and vectors
    p = Int.(floor.(v ./ 2 .+ 1 / 2))
    parity = rem.(abs.(p), 2)
    return (v .- p .* 2) .* (-2 .* parity .+ 1)
end

function triangular_mod(v)
    # Odd piece-wise linear defined within the period [-2, 2] by
    # phi(-2) = phi(0) = phi(2) = 0
    # - phi(-1) = phi(1) = 1
    # this is slightly slower than triangular(v) above because of the mod
    # but looks easier to understand
    vr = mod.(v .+ 1, 4) .- 1
    return vr .+ min.(0 .* vr, -vr .+ 1).*2
end

function triangular(v1, v2)
    return triangular(v1 - v2)
end

const pwDELTA = 0.1
# const pwPERIOD = 4.0

function piecewise_fixed(v)
    Delta = pwDELTA
    vbar = mod.(v .+ 2, 4) .- 2
    s = sign.(vbar)

    ind = sign.(s .* vbar .- 1)
    out = Delta .* ind .+ 0.5

    return s .* out
end

function piecewise_generic(v, Delta = 0.1)
    # This version is slow but generic
    vbar = mod.(v .+ 2, 4) .- 2
    s = sign.(vbar)

    ind = sign.(s .* vbar .- 1)
    out = Delta .* ind .+ 0.5

    return s .* out
end

function piecewise_local(v, Delta = 0.1)
    # this version relies on the small variation of v
    # TODO: make it decent
    vabs = abs(v)
    out =
        if vabs < 1 + Delta
            sign(v)/(1+Delta)*vabs
        else
            if vabs < 3 - Delta
                -sign(v)/(1 - Delta)*(vabs - 2)
            else
                if vabs < 5 + Delta
                    sign(v)/(1+Delta)*(vabs - 4)
                else
                    # add one more step?
                    piecewise_generic(v, Delta)
                end
            end
        end
    return out
 end

function piecewise(v1, v2, Delta = 0.1)
    return piecewise_local(v1 - v2, Delta)
end

function square(v)
    return dtri(v .+ 1)
end

function square(v1, v2)
    return square(v1 - v2)
end

function bilinear(v1, v2)
    return -dtri(v1) .* triangular(v2)
end

function bisine(v1, v2)
    return -Pi2 .* cos(Pi2 .* v1) .* triangular(Pi2 .* v2)
end

function squarishk(v, k=10)
    return tanh.(k .* triangular(v))
end

function squarish(v1, v2, k=10)
    return squarishk(v1 - v2, k)
end

# this really should be a model parameter
const stW = 0.55 * 2

function skew_triangular(v)
    c1 = 1 / (stW * (stW - 2))

    vbar = mod.(v .+ 2, 4) .- 2
    s = sign.(vbar)
    svbar = s .* vbar
    ind = sign.(svbar .- stW)

    mid = (svbar .* (stW - 1) .- stW) .* c1
    Delta = (svbar .- stW) .* c1

    out = Delta .* ind .+ mid
    return s .* out
end

function skew_triangular(v1, v2)
    return skew_triangular(v1 - v2)
end

#
# Continuous representation models
#

function continuous_model_1(v1, v2)
    return -dtri(v1).*triangular(v2)
end

function dtri_cont(v, s)
    # scalar version only
    vreduced = mod(v + 2, 4) - 2
    vabs = abs(vreduced)
    m = 2/(2 - s)
    out =
        if vabs < 1 - s
            m
        elseif vabs <= 1 + s
            (1 - vabs)*m/s
        else
            -m
        end
    return out
end

function continuous_model_1_cont(v1, v2, s = 0.1)
    # change rate in Model 1 with continuous derivative
    # s is the half-width of the transitional region between -1 and 1
    # v1 is presumed to be scalar
    # noticeably slower
    return -dtri_cont(v1, s).*triangular(v2)
end

# For model 2 (differential, GW-continuation, whatnot), the triangular
# function is different. It is an even function with period 2 determined
# by phi_2(0) = phi_2(2) = 0, phi(1) = 1
function continuous_model_2(v1, v2)
    dv = v1 .- v2
    vabs = abs.(dv)
    parity = -rem.(Int.(floor.(vabs)), 2).*2 .+ 1
    return sign.(dv).*parity
end

############################################################
#
### Analysis methods
#
############################################################

"""
    roundup(V)

Return `V` folded into [-2, 2].
"""
function roundup(V)
    return mod.(V .+ 2, 4) .- 2
end

"""
    HammingD(s1, s2)

Evaluate the Hamming distance between binary strings `s1` and `s2`
"""
function HammingD(s1, s2)
    count = 0
    for i in 1:length(s1)
        if s1[i] != s2[i]
            count += 1
        end
    end
    return count
end

"""
    EuclidD(V1, V2)

Evaluate the Euclidean distance between two distributions

INPUT:
  `V1`, `V2` - two arrays on the graph vertices

OUTPUT:
      Sum_v (V_1(v) - V_2(v))^2
"""
function EuclidD(V1, V2)

    return sum((V1 .- V2).^2)
end

function av_dispersion(graph, V)
    # Evaluates average difference between the dynamical variables
    # at adjacent nodes
    #
    # INPUT:
    #   graph
    #   V - the array of dynamical variables
    #
    # OUTPUT:
    #   Σ_e |V(a) - V(b)|/|graph|

    return sum([abs(V[edge.src] - V[edge.dst]) for edge in edges(graph)]) / nv(graph)
end

function c_variance(V, intervals)
    # OBSOLETE (see NOTE)
    # Calculates first two momenta of the pieces of data that fall
    # within given M intervals provided in `intervals`
    #
    # INPUT:
    #   V - array of data
    #   intervals - M x 2 array of intervals boundaries
    #
    # OUTPUT:
    #   M x 3 array with the number of points, mean and variance of data inside the respective intervals
    #
    # NOTE: it cannot process a single interval and does not treat folding intervals

    out = zeros(size(intervals)[1], size(intervals)[2] + 1)

    for i in 1:(size(out)[1])
        ari = filter(t -> intervals[i,1] < t < intervals[i,2], V)
        ni = length(ari)
        av = sum(ari) / ni
        var = sum((ari .- av).^2) / ni
        out[i, 1] = ni
        out[i, 2] = av
        out[i, 3] = sqrt(var)
    end

    return out
end

"""
    cluster_variance(V, interval)

Calculate first two momenta of the pieces of data that fall within
the provided interval

INPUT:
      V - array of data
      interval - (vmin, vmax)

OUTPUT:
      Array[3] = [number of points, mean, variance]
"""
function cluster_variance(V, interval)

    ari = filter(t -> interval[1] < t < interval[2], V)
    ni = length(ari)
    av = sum(ari) / ni
    var = sum((ari .- av).^2) / ni

    return [ni, av, sqrt(var)]
end

# function H_Ising(graph, conf)
#     # Evaluates the Ising energy for the given graph and configuration
#     #
#     # INPUT:
#     #   graph - Graphs object
#     #   conf - configuration array with elemnts \pm 1
#     #
#     # OUTPUT:
#     #   conf * A * conf /2 energy

#     en = 0
#     for edge in edges(graph)
#         en += conf[edge.src]*conf[edge.dst]
#     end
#     return en
# end

"""
    cut(graph, conf)

Evaluate the (weighted) cut value for the given graph and binary configuration

INPUT:
    graph - Graphs object
    conf - binary configuration array with elemnts ± 1
OUTPUT:
    sum_e (1 - e1. e.2)/2
"""
function cut(graph, conf)

    if nv(graph) != length(conf)
        num = conf_to_number(conf)
        println("ERROR: The configuration $num and the vertex set have different size")
    end

    out = 0
    for edge in edges(graph)
        out += (1 - conf[edge.src] * conf[edge.dst]) / 2
    end
    return out
end

function get_rate(VFull)
    # Evaluates the magnitude of variations (discrete time derivative)
    # in a 2D array VFull[time, N].
    # Returns a 1D array of magnitudes with the last element duplicated
    # to keep the same size as the number of time points in VFull

    out = [sum((VFull[:,i + 1] - VFull[:,i]).^2) for i in 1:(size(VFull)[2] - 1)]
    return sqrt.([out; out[end]])
end

"""
    get_best_rounding(graph, V)

Find the best rounding (inspired by the CirCut algorithm)
Return the cut, the configuration, and the threshold (t_c)

INPUT:
  graph
  V - array with the voltage distribution assumed to be within [-2, 2]

OUTPUT:
  (bestcut, bestconf, bestbnd)
          where
              bestcut - the best cut found
              bestconf - rounded configuration
              bestleft - the left boundary of the rounding interval
"""
function get_best_rounding(graph, V)

    Nvert = nv(graph)
    # the width of the rounding interval
    d = 2

    sorting = sortperm(V)
    vvalues = V[sorting]
    push!(vvalues, 2) # a formal end point

    left = -2
    bestleft = left

    start = 1
    stop = findfirst(t -> t > 0, vvalues)
    if isnothing(stop) # there's no such element, all points are in [-2, 0]
        stop = Nvert + 1
    end
    runconf = round_configuration(V, bestleft)
    runcut = cut(graph, runconf)
    bestcut = runcut

    while left < 0
        # which side of the rounding interval meets the next point
        # if vvalues[start] - left <= vvalues[stop] - (left + d) then left
        if vvalues[start] <= vvalues[stop] - d
            flipped = sorting[start] # the true index of the flipped spin
            left = vvalues[start]
            start += 1
        else
            if stop > Nvert # trying to flip the formal point = the end
                break
            end
            flipped = sorting[stop]
            left = vvalues[stop] - d
            stop += 1
        end

        # Now, we evaluate the cut variation
        for j in neighbors(graph, flipped)
            runcut += runconf[flipped]*runconf[j]
        end
        # And update the configuration
        runconf[flipped] *= -1
        if runcut > bestcut
            bestcut = runcut
            bestleft = left
        end
    end
    return (bestcut, round_configuration(V, bestleft), bestleft)
end

function get_best_configuration(graph, V)
    # Finds the best cut following the CirCut algorithm
    # Returns the cut and the configuration
    #
    # INPUT:
    #   graph
    #   V - array with the voltage distribution assumed to be within [-2, 2]
    #
    # OUTPUT:
    #   (bestcut, bestconf)
    #           where
    #               bestcut - the best cut found
    #               bestconf - rounded configuration

    (becu, beco, beth) = get_best_rounding(graph, V)
    return (becu, beco)
end

function get_best_cut(graph, V)
    # Finds the best cut following the CirCut algorithm
    # Returns the cut
    #
    # INPUT:
    #   graph
    #   V - array with the voltage distribution assumed to be within [-2, 2]
    #
    # OUTPUT:
    #   bestcut - the best cut found

    (becu, beco, beth) = get_best_rounding(graph, V)
    return becu
end

function get_best_cut_traced(graph, V)
    # Produce the sequence of cut variations obtained while moving
    # the rounding center
    #
    # INPUT:
    #   graph
    #   V - array with the voltage distribution assumed to be within [-2, 2]
    #
    # OUTPUT:
    #   (DeltaC, rk) - tuple of two arrays
    #         DeltaC - the array of the cut variations
    #         rk      - the array of rounding centers

    Nvert = nv(graph)
    d = 2

    DeltaC = []
    rk = []

    sorting = sortperm(V)
    vvalues = V[sorting]
    push!(vvalues, 2) # a formal end point

    left = -2
    bestleft = left

    start = 1
    stop = findfirst(t -> t > 0, vvalues)
    if isnothing(stop) # there's no such element, all points are in [-2, 0]
        stop = Nvert + 1
    end
    runconf = round_configuration(V, bestleft)

    while left < 0
        if vvalues[start] <= vvalues[stop] - d
            flipped = sorting[start] # the true index of the flipped spin
            left = vvalues[start]
            start += 1
        else
            if stop > Nvert # trying to flip the formal point = the end
                break
            end
            flipped = sorting[stop]
            left = vvalues[stop] - d
            stop += 1
        end

        dC = 0
        for j in neighbors(graph, flipped)
            dC += runconf[flipped]*runconf[j]
        end
        runconf[flipped] *= -1
        push!(DeltaC, dC)
        push!(rk, left)
    end
    return (DeltaC, rk)
end

####################################################
#
###          Service methods
#
####################################################

# Patchy Bernoulli generator
function randspin(p=0.5)
    s = rand()
    return s < p ? 1 : -1
end

function randvector(len::Integer, p=0.5)
    # A Bernoulli sequence of length len
    # Can be used for generating random binary distributions
    return [randspin(p) for i in 1:len]
end

function randnode(nvert)
    # Returns a random number in [1, nvert]
    return rand(tuple(1:nvert...))
end

# Generator of transformations producing all strings with the fixed
# Hamming distance from the given one
#
# USAGE: hamseq[depth](length)
# `depth` is the Hamming distance
# `length` is the bitwise length of the string
#
# OUTPUT is [depth, C]-array of flip indices,
# where C is the number of strings at the given HD
#
# NOTE: Currently `depth` is limited by 5
# TODO: make a universal version
include("hf.jl")

function flipconf(conf, flip)
    # Changes configuration conf according to flips in the index array flip
    #
    # INPUT:
    #   conf - {+1, -1}^N array containing the original string
    #   flip - array with indices where conf should be modified
    #
    # OUTPUT:
    #   a string at the H-distance length(flip) from conf
    #
    # Q: isn't this conf[flip] .*= -1?

    for ind in flip
        conf[ind] *= -1
    end
    return conf
end

function majority_flip!(graph, conf, node)
    # Flips conf[node] to be of the opposite sign to the majority of its neighbors

    flip_flag = false
    tot = 0
    for j in neighbors(graph, node)
        tot += conf[node] * conf[j]
    end
    if tot > 0
        conf[node] *= -1
        flip_flag = true
    end
    return flip_flag
end

function majority_twoflip!(graph, conf, cut_edge)
    # Flips a cut pair if the edges adjacent to the cut edge
    # form the wrong majority
    # Preserves the node-majority
    flip_flag = false
    tot = 0
    for i in neighbors(graph, cut_edge.src)
        tot += conf[cut_edge.src] * conf[i]
    end
    for i in neighbors(graph, cut_edge.dst)
        tot += conf[cut_edge.dst] * conf[i]
    end

    if tot > -2
        conf[cut_edge.src] *= -1
        conf[cut_edge.dst] *= -1
        flip_flag = true
    end
    return flip_flag
end

function local_search(graph, conf)
    # Eliminates vertices breaking the majority rule
    # Attention, it changes conf
    # TODO It's ideologically off and probably should be phased out
    nonstop = true
    while nonstop
        nonstop = false
        for node in vertices(graph)
            nonstop |= majority_flip!(graph, conf, node)
        end
    end
    return conf
end

function local_search!(graph, conf)
    # Eliminates vertices breaking the majority rule
    # Changes the configuration in place and returns the number of passes
    nonstop = true
    count = 0
    while nonstop
        count += 1
        nonstop = false
        for node in vertices(graph)
            nonstop |= majority_flip!(graph, conf, node)
        end
    end
    return count
end

function local_twosearch(graph, conf)
    # Eliminates pairs breaking the edge-majority rule
    nonstop = true
    while nonstop
        nonstop = false
        for link in edges(graph)
            if conf[link.src] * conf[link.dst] < 1
                nonstop |= majority_twoflip!(graph, conf, link)
            end
        end
    end
    return conf
end

function local_twosearch!(graph, conf)
    # Eliminates pairs breaking the edge-majority rule
    # Changes the configuration in place and returns the number of passes
    nonstop = true
    count = 0
    while nonstop
        count += 1
        nonstop = false
        for link in edges(graph)
            if conf[link.src] * conf[link.dst] < 1
                nonstop |= majority_twoflip!(graph, conf, link)
            end
        end
    end
    return count
end

"""
    number_to_conf(number, length)

Return `number`-th configuration of a model with `length` spins

This is `number` in the binary representation
padded with leading zeros to length
"""
function number_to_conf(number, length)
    preconf = digits(number, base=2, pad=length) |> reverse
    return 2 .* preconf .- 1
end

function conf_to_number(conf)
    # Converts conf as a binary number to its decimal representation

    # TODO: This function is rarely needed, hence the dumb code.
    #        I'm not even sure it works correctly
    #        It must be checked against big-endian/little-endian
    binconf = (conf .+ 1) ./ 2

    out = 0
    for i = 1:length(binconf)
        out += binconf[i] * 2^(i - 1)
    end

    return out
end

"""
    round_configuration(V::Array, leftstop)

Maps the supplied configuration to a binary one.
The values in (leftstop + 2] are mapped to +1, the rest are mapped to -1.

INPUT:
    V - data array
    leftstop - the left boundary of the rounding interval

OUTPUT:
    conf - a binary array {-1, 1}^N, where N = length(V)

NOTE:
The function looks for points inside (L + 2], where L is chosen such
that there is no wrapping of the rounding interval. This is determined
by whether leftstop + 2 <= 2 or not. Thus, if leftstop < 0, then
L = leftstop, otherwise L = leftstop - 2 and the points inside (L + 2]
are mapped to -1, while the rest are mapped to 1.

NOTE #2:
The rounding interval (which bound is open and which is closed) is changed
compared to the legacy version.
"""
function round_configuration(V::Array, leftstop)
    # the width of the rounding interval
    width = 2
    L = leftstop

    conf = ones(length(V))
    # We don't want to deal with wrapping
    if L < 0
        conf = -conf
    else
        L -= 2
    end
    inds = L .< V .<= L + width
    conf[inds] .*= -1
    return conf
end

function extract_configuration(V::Array, center)
    # Binarizes V according to the center of the +-interval
    # In the modular form, the mapping looks like
    #
    # V ∈ [center - P/4, center + P/4) -> C_1
    #
    # where P=4 is the period
    #
    #
    # INPUT:
    #   V - data array (is presumed to be rounded and within [-2, 2])
    #   center - the rounding center
    #
    # OUTPUT:
    #   size(V) array with elements + 1 and -1 depending on the relation of
    #           the respective V elements with threshold

    width = 1 # half-width of the rounding interval

    if abs(center) <= 1
        inds = center - width .<= V .< center + width
    else
        return -extract_configuration(V, center - 2 * sign(center))
    end
    out = 2 .* inds .- 1
    return out
end

function get_connected(Nvert, prob)
    # Generate a connected Erdos-Renyi graph with `Nvert` vertices and
    # `prob` density of edges
    # More precisely. On the set of edges of a complete graph
    # K_{Nvert}, we have a (Bernoulli process) function F, which takes
    # values 1 and 0 with probabilities prob and 1 - prob, respectively.
    # The output graph is a connected subgraph of K_{Nvert} with
    # only edges where F = 1 kept in the set of edges.
    cnct = false
    G = Graph()
    while !cnct
        G = erdos_renyi(Nvert, prob)
        cnct = is_connected(G)
    end
    return G
end

function get_initial(Nvert::Integer, (vmin, vmax))
    # Generate random vector with Nvert components uniformly distributed
    # in the (vmin, vmax) interval

    bot, top = minmax(vmin, vmax)
    mag = top - bot

    return mag .* rand(Float64, Nvert) .+ bot
end

##########################################################
#
###   Dynamics methods
#
##########################################################

function step_rate(graph::SimpleGraph, method::Function, V::Array, Ks::Float64)
    # Evaluates ΔV for a single step
    # This version presumes that there is a single method for coupling and
    # anisotropy and that there is only easy-axis anisotropy
    #
    # This is a default function (notice how the anisotropy is evaluated
    # using method(V, -V)). It should be used with care as the
    # implemented idea of a "single method" is flawed. A true single
    # method would imply the presence of a field (a fixed
    # vertex). Here, however, the "single method" is used for
    # anisotropy, which presumes a tensor.
    #
    # INPUT:
    #    graph - unweighted graph carrying V's
    #    method - method to evaluate different contributions
    #   V(1:|graph|) - current distribution of dynamical variables
    #    Ks - anisotropy constant
    #
    # OUTPUT:
    #   ΔV(1:|graph|) - array of increments

    #    out = Ks .* method(V, -V) # refactoring to scalars
    out  = zeros(length(V))
    for node in vertices(graph)
        Vnode = V[node]
        out[node] =  Ks .* method(Vnode, -Vnode)

        for neib in neighbors(graph, node)
            out[node] += method(Vnode, V[neib])
        end
    end
    return out
end

function step_rate_extended(graph::SimpleGraph, coupling::Function, V::Array,
                            Ks::Float64, anisotropy::Function, R::Float64)
    # Evaluate ΔV and ΔR in the extended system
    #
    # INPUT:
    #    graph - unweighted graph carrying V's
    #    method - method to evaluate different contributions
    #   V(1:|graph|) - current distribution of dynamical variables
    #    Ks - anisotropy constant
    #
    # OUTPUT:
    #   ΔV(1:|graph|) - array of increments

    #    out = Ks .* method(V, -V) # refactoring to scalars
    out  = zeros(length(V))
    r_out = 0
    for node in vertices(graph)
        Vnode = V[node]
        an = Ks * anisotropy(Vnode - R)
        out[node] = an
        r_out += -an
        for neib in neighbors(graph, node)
            out[node] += coupling(Vnode, V[neib])
        end
    end
    return (out, r_out)
end

function trajectories(graph, duration, scale, method::Function, Ks, Vini)
    # Advances the graph duration - 1 steps forward
    # This is the verbose version, which returns the full dynamics
    #
    # scale - parameter to tweak the dynamics
    # duration - how many time points to evaluate
    # V0 - the initial conditions
    #
    # OUTPUT:
    #   VFull = [V(0) V(1) ... V(duration-1)]

    VFull = Vini
    V = Vini

    tran = 1:(duration - 1)
    for tau in tran
        ΔV = scale .* step_rate(graph, method, V, Ks)
        V += ΔV
        VFull = [VFull V]
    end

    return VFull
end

function trajectories(model::Model, duration, Vini)
    # Advances the graph `(duration - 1)` steps forward
    # This is the verbose version, which returns the full dynamics
    #
    # model - the model description
    # duration - how many time points to evaluate
    # V0 - the initial conditions
    #
    # OUTPUT:    #   VFull = [V(0) V(1) ... V(duration-1)]

    return trajectories(model.graph, duration, model.scale, model.coupling,
                        model.Ks, Vini)
end

function trajectories_extended(graph, duration, scale, pair_method::Function, Vini,
                               Ks, one_method::Function, Rini)
    # Advance the graph `duration - 1` steps forward
    # This is the verbose version, which returns the full dynamics
    #
    #
    # OUTPUT:
    #   (VFull, RFull)
    #           VFull = [V(0) V(1) ... V(duration-1)]
    #           RFull - the final state of the added degree of freedom

    VFull = Vini
    V = Vini
    RFull = [Rini]
    R = Rini

    for tau in 1:(duration - 1)
        (ΔV, ΔR) = step_rate_extended(graph, pair_method, V, Ks, one_method, R)
         V += scale .* ΔV
        R += scale * ΔR
        VFull = [VFull V]
        RFull = [RFull R]
    end

    return (VFull, RFull)
end

function trajectories_extended(model::Model, duration, Vini, Rini)
    # Advances the graph `(duration - 1)` steps forward
    # This is the verbose version, which returns the full dynamics
    #
    # model - the model description
    # duration - how many time points to evaluate
    # V0 - the initial conditions
    #
    # OUTPUT:
    #   VFull = [V(0) V(1) ... V(duration-1)]

    return trajectories_extended(model.graph, duration, model.scale, model.coupling,
                        Vini, model.Ks, model.anisotropy, Rini)
end

function propagate(graph, duration, scale, method::Function, Ks, Vini)
    # Advances the graph duration - 1 steps forward
    # This is the short version, which returns only the final state vector
    #
    # scale - parameter to tweak the dynamics
    # duration - how many time points to evaluate
    # Vini - the initial conditions
    #
    # OUTPUT:
    #   [V[1] .. V[nv(graph)] at t = duration - 1
    #
    # NOTE: the order of parameters changed in 0.2.0

    V = Vini

    for tau in 1:(duration - 1)
        ΔV = scale .* step_rate(graph, method, V, Ks)
        V += ΔV
    end

    return V
end

function propagate(model::Model, duration, Vini)
    # Advances the model::Model duration - 1 steps forward
    # Returns only the final state vector
    #
    # model - the model description
    # duration - how many time points to evaluate
    # Vini - the initial conditions
    #
    # OUTPUT:
    #   [V[1] .. V[nv(graph)] at t = duration - 1

    return propagate(model.graph, duration, model.scale, model.coupling,
                     model.Ks, Vini)
end

"""
    propagate_extended(graph, duration, scale, pair_method, Vini,
                              Ks, one_method, Rini)

Propagate extended G + 1 model with the explicit specification.

INPUT:
      graph - the model description
      duration - the number of time points
      scale - the length of the single time step
      pair_method - coupling
      Vini - the initial configuration of G
      Ks - the anisotropy costant
      one_method - the anisotropy method
      Rini - the initial state of the additional degree of freedom

OUTPUT:
      (Vfinal, Rfinal) -
              Vfinal - the final configuration of G
              Rfinal - the final state of the additional degree of freedom
"""
function propagate_extended(graph, duration, scale, pair_method, Vini, Ks, one_method, Rini)

    V = Vini
    R = Rini

    for tau in 1:(duration - 1)
        (ΔV, ΔR) = step_rate_extended(graph, pair_method, V, Ks, one_method, R)
        V += scale .* ΔV
        R += scale * ΔR
    end
    return (V, R)
end

"""
    propagate_extended(model::Model, duration, Vini, Rini)

Propagate extended G + 1 model.

INPUT:
      model - the model description
      duration - the number of time points
      Vini - the initial configuration of G
      Rini - the initial state of the additional degree of freedom

OUTPUT:
      (Vfinal, Rfinal) -
              Vfinal - the final configuration of G
              Rfinal - the final state of the additional degree of freedom
"""
function propagate_extended(model::Model, duration, Vini, Rini)
    graph = model.graph
    Ks = model.Ks
    pair_method = model.coupling
    one_method = model.anisotropy
    scale = model.scale

    return propagate_extended(graph, duration, scale, pair_method, Vini,
                              Ks, one_method, Rini)
end

# function energy(graph, method::Function, V::Array)
#     # Evaluates the coupling energy corresponding to V
#     # Note: this is essential that this energy evaluates only the coupling energy without
#     #       any anisotropic terms
#     #
#     #####            NOTE: It's broken as of now! Need to pass the correct method
#     #
#     # INPUT:
#     #   method here is the energy of the elementary one-edge graph. Of course, this is
#     #          not the same method as in the propagation functions (that woould be the minus
#     #          gradient of the method passed to this function)

#     en = 0
#     for edge in edges(graph)
#         en += cosine(V[edge.src] - V[edge.dst])
#     end

#     return en
# end

############################################################
#
### Simulation
#
############################################################

function conf_decay(graph, conf::Array, listlen=3)
    # Evaluate the instability of conf
    #
    # INPUT:
    #   graph
    #   conf - a string with configuration encoded in
    #   listlen - how many eigenvalues to be returned
    #
    # OUTPUT:
    # array[1:listlen] of the largest eigenvalues

    NVert = nv(graph)
    if NVert !== length(conf)
        throw(Error("The configuration legnth does not match the size of the graph"))
    end

    L = laplacian_matrix(graph)

    D = zeros(NVert, 1) # to collect sums of rows

    for (i, j) in zip(findnz(L)...)
        L[i,j] *= conf[i] * conf[j]
        if i != j
            D[i] += L[i,j]
        end
    end

    # Correct diagonal elements of the cut Hessian
    for i in 1:NVert
        L[i,i] = -D[i]
    end

    out = eigs(L, nev=listlen, which=:LR)[1]
    return out[1:listlen]
end

function conf_decay_states(graph, conf::Array, listlen=3)
    # Evaluates the eigenvalue and eigenvector of the most unstable excitation in configuration conf
    #
    # INPUT:
    #   graph
    #   conf - a string with configuration encoded in
    #   listlen - how many eigenvalues to be returned
    #
    # OUTPUT:
    # lambda, v - eigenvalue and eigenvector

    NVert = nv(graph)
    if NVert !== length(conf)
        throw(Error("The configuration length does not match the size of the graph"))
    end

    L = laplacian_matrix(graph)

    D = zeros(NVert, 1) # to collect sums of rows
    for (i, j) in zip(findnz(L)...)
        L[i,j] *= conf[i] * conf[j]
        if i != j
            D[i] += L[i,j]
        end
    end

    # Correct diagonal elements of the cut Hessian
    for i in 1:NVert
        L[i,i] = -D[i]
    end

    eigvalue, eigvector = eigs(L, nev=1, which=:LR)
    return (eigvalue, eigvector)
end

function scan_for_best_configuration(model::Model, Vc::Array,
                                     domain::Float64, tmax::Integer, Ninitial::Integer)
    # Scans a vicinity of Vc with size given by domain in the "Monte Carlo style"
    #
    # Ninitial number of random initial conditions with individual amplitudes
    # varying between ±domain

    G = model.graph
    Nvert = nv(G)

    (bcut, bconf) = get_best_configuration(G, Vc)

    for i in 1:Ninitial
        local cucu::Integer

        Vi = domain .* get_initial(Nvert, (-1, 1))
        Vi .+= Vc;
        Vi[rand((1:Nvert))] *= -1.1  # a random node is flipped as a perturbation
        # Vi .-= sum(Vi)/Nvert

        # VF = propagateAdaptively(model, tmax, Vi);
        VF = propagate(model, tmax, Vi);

        (cucu, cuco) = get_best_configuration(G, roundup(VF))

        if cut(G, cuco) != cucu
            println("INTERNAL ERROR! Cuts are inconsistent!")
        end

        if cucu > bcut
            # println("Improvement by $(cucu - bcut)")
            bcut = cucu
            bconf = cuco
            Vc = bconf
        end
    end

    return (bcut, bconf)
end

function test_branch(model, Vstart, domain, tmax, Ni, max_depth)
    # Tests branch starting from Vstart
    # INPUT:
    #   model
    #   Vstart
    #   domain      - the size of the vicinity to take the initial states from
    #   tmax        - for how long propagate the particular initial conditions
    #   Ni          - the number of trials
    #   max_depth   - the maximal depth

    local bcut::Integer, bconf::Array{Integer}, stepcount::Integer, nextflag::Bool

    Nvert = nv(model.graph)

    nextflag = true

    (bcut, bconf) = get_best_configuration(model.graph, Vstart)

    Vcentral = Vstart
    stepcount = 0
    while nextflag
        local cucu::Integer

        stepcount += 1

        (cucut, cuconf) = scan_for_best_configuration(model, Vcentral, domain, tmax, Ni)
        message(model, "Prelocal cut = $cucut", 2)
        cuconf = local_search(model.graph, cuconf)
        message(model, "Post-local cut = $(cut(model.graph, cuconf))", 2)

        cuconf = local_twosearch(model.graph, cuconf)
        cucut = cut(model.graph, cuconf)
        message(model, "Post-local-2 cut = $cucut", 2)
        if cucut > bcut
            bcut = cucut
            message(model,
                    "Step $stepcount arrived at $bcut with the displacement $(HammingD(cuconf, bconf))", 2)
            bconf = cuconf
            Vcentral = cuconf

            # domain *= 0.9
            Ni += 10
        else
            nextflag = false
        end

        if stepcount > max_depth
            message(model, "Exceeded maximal depth", 2)
            nextflag = false
        end
    end
    return (bcut, bconf)
end

end # end of module Dice

max_stars_repo_name: merement/Dice

max_issues_repo_name: merement/Dice

max_forks_repo_name: merement/Dice


text:
"""
Structure implementing lists of intervals.
These structures are extremely useful when working with stream graphs.
"""
mutable struct Intervals
    list::Array{Tuple{Float64,Float64},1}
end

Intervals()=Intervals([])

"""
Equality between two lists of intervals
NOTE: Use approx here because of rounding errors in intervals.
"""
==(i1::Intervals,i2::Intervals)=clean(i1).list≈clean(i2).list

"""
In operator.
Returns true if the given element is inside one of the intervals of the list.
"""
function ∈(x::Float64, i::Intervals)
    for ii in i.list
        if ii[1] <= x <= ii[2]
            return true
        end
    end
    return false
end

"""Returns the number of intervals in the list."""
count(i::Intervals) = length(i.list)

"""Compute the length of a list of intervals as the sum of the interval lengths."""
length(i::Intervals)=sum([ii[2]-ii[1] for ii in i.list])

"""Merge a new interval within a list of intervals."""
merge(i::Intervals, x::Tuple{Float64,Float64})=union(i, Intervals([x]))

"""
Push a new interval at the end of the list.
WARNING: No sorting is done here...
"""
function push(i::Intervals, x::Tuple{Float64,Float64})
    push!(i.list, x)
    return i
end

"""Intersection between two list of intervals."""
function ∩(i1::Intervals, i2::Intervals)
    i1=clean(i1)
    i2=clean(i2)
    l = []
    components = vcat(i1.list,i2.list)
    if length(components)==0
        return Intervals()
    end
    sort!(components)
    push!(l, popfirst!(components))
    while length(components)>0
        modif = false
        c = popfirst!(components)
        if length(l) == 0
            push!(l,c)
            continue
        end
        if c[2] < l[end][2]
            push!(components,(c[2],l[end][2]))
            l[end] = (l[end][1],c[2])
            modif = true
        end
        if l[end][2] >= c[1] >= l[end][1]
            l[end] = (c[1], l[end][2])
            modif = true
        end
        if !modif
            pop!(l)
        else
            sort!(components)
        end
        push!(l,c)
    end
    pop!(l)
    return clean(Intervals(l))
end

"""Union between two lists of intervals."""
function ∪(i1::Intervals, i2::Intervals)
    l = []
    components = union(i1.list,i2.list)
    sort!(components)
    for c in components
        if length(l) == 0 || c[1] > l[end][2]
            push!(l,c)
        elseif c[2] > l[end][2]
            l[end] = (l[end][1], c[2])
        end
    end
    return clean(Intervals(l))
end

function clean(i::Intervals)
    l = []
    sort!(i.list)
    for c in i.list
        if length(l) == 0
            push!(l,c)
            continue
        end
        if c[1] == l[end][2]
            l[end] = (l[end][1],c[2])
            continue
        end
        push!(l,c)
    end
    return Intervals(l)
end

"""Is subset returns True if i2 is included within i1."""
⊆(i1::Intervals,i2::Intervals)=(i1 ∩ i2)==clean(i1)
⊆(i1::Tuple{Float64,Float64},i2::Intervals)=Intervals([i1]) ⊆ i2

ext: jl


text:
"""
    CompositeGP{Targs}

A GP derived from other GPs via an affine transformation. Specification given by `args`.
"""
struct CompositeGP{Targs} <: AbstractGP
    args::Targs
    n::Int
    gpc::GPC
    function CompositeGP{Targs}(args::Targs, gpc::GPC) where {Targs}
        gp = new{Targs}(args, next_index(gpc), gpc)
        gpc.n += 1
        return gp
    end
end
CompositeGP(args::Targs, gpc::GPC) where {Targs} = CompositeGP{Targs}(args, gpc)

mean_vector(f::CompositeGP, x::AV) = mean_vector(f.args, x)

cov(f::CompositeGP, x::AV) = cov(f.args, x)
cov_diag(f::CompositeGP, x::AV) = cov_diag(f.args, x)

cov(f::CompositeGP, x::AV, x′::AV) = cov(f.args, x, x′)
cov_diag(f::CompositeGP, x::AV, x′::AV) = cov_diag(f.args, x, x′)

function cov(f::AbstractGP, f′::AbstractGP, x::AV, x′::AV)
    @assert f.gpc === f′.gpc
    if f.n === f′.n
        return cov(f.args, x, x′)
    elseif f isa GP && f.n > f′.n || f′ isa GP && f′.n > f.n
        return zeros(length(x), length(x′))
    elseif f.n >= f′.n
        return cov(f.args, f′, x, x′)
    else
        return cov(f, f′.args, x, x′)
    end
end
function cov_diag(f::AbstractGP, f′::AbstractGP, x::AV, x′::AV)
    @assert f.gpc === f′.gpc
    if f.n === f′.n
        return cov_diag(f.args, x, x′)
    elseif f isa GP && f.n > f′.n || f′ isa GP && f′.n > f.n
        return zeros(length(x))
    elseif f.n >= f′.n
        return cov_diag(f.args, f′, x, x′)
    else
        return cov_diag(f, f′.args, x, x′)
    end
end


text:
const CSTACK1 = Array{Bool,1}()
const CSTACK2 = Array{Bool,1}()
const CSTACK3 = Array{Bool,1}()

macro if2(condition1, condition2, alltrue)
    if !(length(CSTACK1) == length(CSTACK2) == length(CSTACK3))
        throw("prior if2 statement error: must have if2, elseif1, elseif2, and elseifneither")
    end
    ifcond1 = eval(condition1)
    ifcond2 = eval(condition2)
    if ifcond1 && ifcond2
        eval(alltrue)
        push!(CSTACK1, false)
        push!(CSTACK2, false)
        push!(CSTACK3, false)
    elseif ifcond1
        push!(CSTACK1, true)
        push!(CSTACK2, false)
        push!(CSTACK3, false)
    elseif ifcond2
        push!(CSTACK1, false)
        push!(CSTACK2, true)
        push!(CSTACK3, false)
    else
        push!(CSTACK1, false)
        push!(CSTACK2, false)
        push!(CSTACK3, true)
    end
end

macro elseif1(block)
    quote
        if pop!(CSTACK1)
            $block
        end
    end
end

macro elseif2(block)
    quote
        if pop!(CSTACK2)
            $block
        end
    end
end

macro elseifneither(block)
    quote
        if pop!(CSTACK3)
            $block
        end
    end
end

# Testing of code starts here

@if2(2 != 4, 3 != 7, begin x = "all"; println(x) end)

@elseif1(begin println("one") end)

@elseif2(begin println("two") end)

@elseifneither(begin println("neither") end)

@if2(2 != 4, 3 == 7, println("all"))

@elseif1(begin println("one") end)

@elseif2(begin println("two") end)

@elseifneither(begin println("neither") end)

@if2(2 == 4, 3 != 7, begin x = "all"; println(x) end)

@elseif1(begin println("one") end)

@elseif2(begin println("two") end)

@elseifneither(begin println("neither") end)

@if2(2 == 4, 3 == 7, println("last all") end)

@elseif1(begin println("one") end)

@elseif2(begin println("two") end)

@elseifneither(begin println("neither") end)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Fairly generic data structures.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-06-06

module Structs

export Identifier, Identities, id

immutable Identifier{T,U<:Integer}
    d::Dict{T,U}
end

Identifier(T::Type=Any,U::Type=UInt) = Identifier(Dict{T,U}())

import Base: length
length(c::Identifier) = length(c.d)

id(c::Identifier, o) = get!(c.d, o, length(c) + 1)

immutable Identities{T,U<:Integer}
    d::Dict{T,U}
    default::U
end

function Identities{T,U<:Integer}(c::Identifier{T,U})
    Identities(c.d, zero(U))
end

function Identities{T,U<:Integer}(c::Identifier{T,U}, default::U)
    Identities(c.d, default)
end

import Base: length
length(c::Identities) = length(c.d)

id(c::Identities, o) = get(c.d, o, c.default)

end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
## BortnikRayTrace.jl
# Todd Anderson
# 5 October 2020
#
# Ray-tracing of electromagnetic waves from upper ionosphere to plasmasphere using Haselgrove equations, following method in Jacob Bortnik's thesis.  Includes main ODE solver function calls.
#
# Ray-tracing functions:
# a.  dr/dt   = (1/μ²)*(ρᵣ - μ*(δμ/δρᵣ))
# b.  dθ/dt   = (1/rμ²)*(ρ_θ - μ*(δμ/δρ_θ))
# c.  dϕ/dt   = (1/rμ²sinθ)*(ρ_ϕ - μ*(δμ/δρ_ϕ))
# d.  dρᵣ/dt  = (1/μ)*(δμ/δr) + ρ_θ*dθ/dt + ρ_ϕ*dϕ/dt*sinθ
# e.  dρ_θ/dt = (1/r)*((1/μ)*(δμ/δθ) - ρ_θ*dr/dt + r*ρ_ϕ*dϕ/dt*cosθ)
# f.  dρ_ϕ/dt = (1/r*sinθ)*((1/μ)*(δμ/δϕ) - ρ_ϕ*dr/dt*sinθ - r*ρ_ϕ*dθ/dt*cosθ)
#
# Assume all rays fixed to meridional plane -> dϕ/dt, ρ_ϕ, dρ_ϕ/dt = 0
#   -> ignore eqns (c) and (f), new forms for (d) and (e):
#   d2. dρᵣ/dt  = (1/μ)*(δμ/δr) + ρ_θ*dθ/dt
#   e2. dρ_θ/dt = (1/r)*((1/μ)*(δμ/δθ) - ρ_θ*dr/dt)
# Functions:
# 1. refractive_index: calculate
#
# TODO April 16 2021
# 1. Calculate ψ, μ from components ρ_r, ρ_λ, ρ_ϕ
# 	-> most likely do this in refractive_index, or in separate magnetic field/plasma environment functions
# 	-> calculation of ψ requires magnetic field direction, from dip calculation or lookup
# 	--> need three-component B field, not just dip angle w.r.t ground
# 2. Check for latitude/colatitude consistency: should be using colatitude (θ) for ray-tracing

using DifferentialEquations
using LinearAlgebra
using BenchmarkTools
using LSODA
using Sundials
using Plots

# physical constants
c = 2.99792458e8	# speed of light in m/s
re = 6.3712e6      	# radius of Earth in meters
B0 = 3.12e-5   		# magnitude of B field at Earth's equator surface in Tesla
e = 1.602e-19 		# electric charge in Coulombs
me = 9.1093e-31		# electron rest mass
mp = 1.6726219e-27	# proton rest mass
eps = 8.854e-12 	# permittivity of free space in mks

# plasmasphere parameters
Kp_max = 3                             	# maximum Kp index
Lppi = 5.6 - 0.46*Kp_max                # plasmapause inner limit in L
d = 0                              		# day number
R̄ = 90                                  # 13-month average sunspot number
mlt = 2                             	# magnetic local time
Lppo = initialize_plasmasphere(Kp_max, Lppi, d, R̄, mlt)

b = function magnetic_field(r,θ,ϕ)
	# get magnetic field components at any location in magnetic local time
	# centered dipole field model (i.e. no correction for solar forcing or higher-order geomagnetic model terms): valid for L ≈ 2 - 5
	λ = π/2 - θ								# convert colatitude to latitude
	Br = -2*B0*(re/r)^3*sin(λ)              # radial compoenent of B-field
    Bλ = B0*(re/r)^3*cos(λ)                 # latitude component of B-field
	Bθ = -1.0*Bλ							# colatitude component of B-field
	Bϕ = 0.0
    Bmag = B0*(re/r)^3*sqrt(1 + 3*sin(λ)^2)  # magnitude of B-field == sqrt(Br^2 + Bλ^2)

	b = [Br,Bθ,Bϕ]

end

Lppo = function initialize_plasmasphere(Kp_max, Lppi, d, R̄, mlt)
	# pre-solve plasmasphere from Carpenter and Anderson 1992 model
	r_range = [re:1000:(10*re);]
	λ_range = 0
	L = r_range./(re*cos(λ_range)^2) # range of L-values where Lppo will be searched for

	Kp_max = 3              		# for example
	Lppi = 5.6 - 0.46*Kp_max        # plasmapause inner limit in L
	#Lppo = 5.75                    # placeholder
	d = 0                           # day number
	R̄ = 90							 # 13-month average sunspot number
	mlt = 2                         # magnetic local time

	ne_Lppi = 10^((-0.3145*Lppi + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-Lppi)/1.5))

	# pre-solve individual components in order to find Lppo
	ne_plasma_1 = @. 10 .^((-0.3145.*L + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2 .-L)./1.5))
	ne_plasma_2 = ne_Lppi*10 .^(-1.0.*(L .- Lppi)./0.1)
	ne_plasma_3 = (5800 + 300*mlt).*L.^(-4.5) + (1 .- exp.((2 .-L)./10))

	min_index = findmin(abs.(ne_plasma_2 - ne_plasma_3))[2]
	Lppo = L[min_index]

end

u = function refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
    # 1. set up medium
    # 1a. magnetic field

	Bvec = magnetic_field(r,θ,ϕ)
	Br, Bθ, Bϕ = Bvec
	Bmag = sqrt(Br^2 + Bθ^2 + Bϕ^2)

    # 1b. plasma
    # assumptions: cold plasma, 2.25 < L < 8, midnight sector, based on Carpenter and Anderson 1992
	#Lppo = initialize_plasmasphere(Kp_max, Lppi, d, R̄, mlt) # include this line if any Lppo variables change during integration (they shouldn't)

    L = r/(re*cos(π/2 - θ)^2)
	ne_iono = (1.8e5.*exp.(-4.183119.*((r./re).-1.0471))) # cm^-3

    if L <= Lppi
        log_ne = (-0.3145*L + 3.9043) + (0.15*(cos((2*π*(d+9)))/365 - 0.5*cos((4*π*(d+9)))/365) + 0.00127*R̄ - 0.0635)*exp((2-L)/1.5)
		ne_plasma = 10^(log_ne)
    elseif Lppi < L <= Lppo
        ne_plasma = ne_Lppi*10^((Lppi-L)/0.1)
    elseif Lppo < L
        ne_plasma = (5800 + 300*mlt)*L^(-4.5) + (1 - exp((2-L)/10))
    else
        ne_plasma = 0.0
    end

	n_e = (ne_plasma + ne_iono)*1e6 	# add plasmasphere + ionosphere, convert from cm^-3 to m^-3
    n_p = n_e

    # electron plasma angular frequency squared
    ω_e2 = (n_e*(e^2))/(eps*me)
	# proton plasma angular frequency squared
    ω_p2 = (n_p*(e^2))/(eps*mp)

	# electron angular gyrofrequency
    Ω_e = (e*Bmag)/(me)
	# proton angular gyrofrequency
    Ω_p = (e*Bmag)/(mp)

	# convert wave frequency to gyrofrequency
	ω = 2*π*freq

	# 2. define dispersion relation
	# 2a. calculate ψ from ρ_r, ρ_θ, ρ_ϕ

	μvec = [ρ_r, ρ_θ, ρ_ϕ]
	μmag = sqrt(ρ_r^2 + ρ_θ^2 + ρ_ϕ^2)
	cos_ψ = (Bvec ⋅ μvec)/(Bmag*μmag)
	ψ = acos(cos_ψ)

	# define R, L, P, D, S,
	# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
	R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
	#println("R = ",R)

	# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
	L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
	#println("L = ",L)

	# P ≡ 1 - Σ(ωₖ²/ω²)
	P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
	#println("P = ",P)

	# D ≡ ½(R-L); S ≡ ½(R+L)
	D = (R - L)/2.0
	S = (R + L)/2.0

	# define parts of dispersion relation
	# Aμ⁴ - Bμ² + C = 0
	# A = Ssin²ψ + Pcos²ψ
	A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0

	# B = RLsin²ψ + PS(1 + cos²ψ)
	B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0)

	# C  = PRL
	C = P*R*L

	# solve the dispersion relation for μ:
	# μ² = (B +- F)/2A
	# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
	F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
	F = sqrt(F2)

	# Rice 1997: typical solution to dispersion relation with quadratic formula
	μ2_minus = (B - F)/(2.0*A)
	μ2_plus = (B + F)/(2.0*A)

	μ_plus = try
		sqrt(abs(μ2_plus)) # abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			# println("DomainError in μ_plus!")
			# println("B=", B)
			# println("F=", F)
			# println("A=", A)
		else
			# println(err)
		end
	end

	# Electron whistler case: ψ = 0, μ² = R; this is the μ_plus case
	#μ = μ_minus		# EMIC case
	μ = μ_plus		# electron whistler case

	# Find dA/dψ, dB/dψ, dC/dψ, dμ/dψ
	dAdψ = 2.0*(S-P)*sin(ψ)*cos(ψ)
	dBdψ = 2.0*(R*L-P*S)*sin(ψ)*cos(ψ)
    dCdψ = 0.0
    #dμdψ = ((μ^4.0)*dAdψ-(μ^2.0)*dBdψ+dCdψ)/(4.0*A*(μ^3.0)-2.0*B*μ)

	dFdψ = 1/(2.0*F)*((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ))
	#dFdψ = sqrt(abs((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)))
	dμdψ = 1/(2.0*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B + F)/(2*A^2))
	# NOTE 11/24/2020: from Rice 1997; choose '+' for (B +- F) and (dBdψ +- dFdψ)

	#DEBUG check values
	#println("μ = ",μ)
	#println("dμdψ = ", dμdψ)

	# non-mutating output
	u = [μ,dμdψ,ψ,Bvec]

end

# Derivatives of phase refractive index w.r.t. r, λ, ϕ, ρ_r, ρ_λ, ρ_ϕ, freq
# calculate derivative of μ w.r.t. r
dμdr = function ddr(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ,dμdψ = u			# unpack
	dr = 1.0e-11		# r step size

	μ_l = refractive_index(r-dr/2.0,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_l = μ
	μ_r = refractive_index(r+dr/2.0,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_r = μ

	dμdr = (μ_r[1] - μ_l[1])/dr
end

dμdθ = function ddθ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ,dμdψ = u			# unpack
	dθ = 1.0e-11		# θ step size

	μ_l = refractive_index(r,θ-dθ/2.0,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_l = μ
	#println("μ_l = ", μ)
	μ_r = refractive_index(r,θ+dθ/2.0,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_r = μ
	#println("μ_r = ", μ)

	dμdθ = (μ_r[1] - μ_l[1])/dθ
end

dμdϕ = function ddϕ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ,dμdψ = u			# unpack
	dϕ = 1.0e-11		# χ step size

	μ_l = refractive_index(r,θ,ϕ-dϕ/2.0,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_l = μ
	μ_r = refractive_index(r,θ,ϕ+dϕ/2.0,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ_r = μ

	dμdϕ = (μ_r[1] - μ_l[1])/dϕ
end

dμdρr = function ddρr(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dρr = 1.0e-11		# ρ_r step size

	# Kimura 1966:
	# dμ/dρ_k = (dμ/dψ)(dψ/dρ_k) = (dμ/dψ)((ρ_k*cosψ - μcos(α_Bk)/(μ²sinψ))
	# where k = r, λ, ϕ; and α_Bk = angle between the magnetic field vector and k

	μ, dμdψ, ψ, Bvec = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq) #DEBUG: BoundsError at this location 04/21/2021

	ρrvec = [ρ_r, 0.0, 0.0]
	Bmag = sqrt(Bvec[1]^2 + Bvec[2]^2 + Bvec[3]^2)
	ρrmag = abs(ρ_r)

	cos_αr = (Bvec ⋅ ρrvec)/(Bmag*ρrmag)

	dμdρr = dμdψ*((ρ_r*cos(ψ) - μ*cos_αr)/(μ^2 * sin(ψ)))
end

dμdρθ = function ddρλ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dρθ = 1.0e-11		# ρ_λ step size

	# Kimura 1966:
	# dμ/dρ_k = (dμ/dψ)(dψ/dρ_k) = (dμ/dψ)((ρ_k*cosψ - μcos(α_Bk)/(μ²sinψ))
	# where k = r, λ, ϕ; and α_Bk = angle between the magnetic field vector and k
	μ, dμdψ, ψ, Bvec = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq) # unpack

	ρθvec = [0.0, ρ_θ, 0.0]
	Bmag = sqrt(Bvec[1]^2 + Bvec[2]^2 + Bvec[3]^2)
	ρθmag = abs(ρ_θ)

	cos_αθ = (Bvec ⋅ ρθvec)/(Bmag*ρθmag)

	dμdρθ = dμdψ*((ρ_θ*cos(ψ) - μ*cos_αθ)/(μ^2 * sin(ψ)))
end

dμdρϕ = function ddρϕ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq) #TODO λ → θ
	dρϕ = 1.0e-11		# ρ_ϕ step size

	# Kimura 1966:
	# dμ/dρ_k = (dμ/dψ)(dψ/dρ_k) = (dμ/dψ)((ρ_k*cosψ - μcos(α_Bk)/(μ²sinψ))
	# where k = r, λ, ϕ; and α_Bk = angle between the magnetic field vector and k
	μ, dμdψ, ψ, Bvec = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq) # unpack

	ρϕvec = [0.0, 0.0, ρ_ϕ]
	Bmag = sqrt(Bvec[1]^2 + Bvec[2]^2 + Bvec[3]^2)
	ρϕmag = abs(ρ_ϕ)

	cos_αϕ = (Bvec ⋅ ρϕvec)/(Bmag*ρϕmag)

	dμdρϕ = dμdψ*((ρ_ϕ*cos(ψ) - μ*cos_αϕ)/(μ^2 * sin(ψ)))
end

dμdf = function ddf(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	#μ,dμdψ = u			# unpack
	df = 1.0e-5			# f step size

	μ_l = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq-df/2.0)
	#μ_l = μ
	μ_r = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq+df/2.0)
	#μ_r = μ

	dμdf = (μ_r[1] - μ_l[1])/df
end

# calculate derivatives w.r.t. time
function haselgrove!(du,u,p,t)
# 	du = drdt, dθdt, dϕdt, dρrdt, dρθdt, dρϕdt
# 	u = r, θ, ϕ, ρ_r, ρ_θ, ρ_ϕ, T
#  	p = freq
# 	t = time

	r, θ, ϕ, ρ_r, ρ_θ, ρ_ϕ, T = u 	# unpack u
	freq = p[1] 						# unpack p

	dμdr = ddr(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdθ = ddθ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdϕ = ddϕ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdρr = ddρr(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdρθ = ddρλ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdρϕ = ddρϕ(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	dμdf = ddf(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)

	v = refractive_index(r,θ,ϕ,ρ_r,ρ_θ,ρ_ϕ,freq)
	μ = v[1]
	dμdψ = v[2]

	# TODO: ensure ρ_k have appropriate initial conditions!

    du[1] = 1/(μ^2)*(ρ_r-μ*dμdρr) 							# r: radial distance (m)
    du[2] = 1/(r*μ^2)*(ρ_θ-μ*dμdρθ)							# θ: colatitude (radians)
    du[3] = 1/(r*μ^2*sin(θ))*(ρ_ϕ-μ*dμdρϕ)					# ϕ: longitude (radians)
	du[4] = (1/μ)*dμdr + ρ_θ*du[2] + ρ_ϕ*du[3]*sin(θ)		# ρ_r: r-component of μ
	du[5] = (1/r)*((1/μ)*dμdθ - ρ_θ*du[1] + r*ρ_ϕ*du[3]*cos(θ)) # ρ_θ: θ-component of μ
	du[6] = (1/(r*sin(θ)))*((1/μ)*dμdϕ - ρ_ϕ*du[1]*sin(θ) - r*ρ_ϕ*du[2]*cos(θ)) # ρ_ϕ: ϕ-component of μ
	du[7] = (1/c)*(1 + (freq/μ)*dμdf) # T: group delay time (s)

end

## callbacks!
# re_cb: terminate when ray intersects Earth surface
function re_term_condition(u,t,integrator)
	u[1] - re
end

function terminate_affect!(integrator)
	terminate!(integrator)
end

re_cb = ContinuousCallback(re_term_condition,terminate_affect!)

# TODO: adapt saving callback to new ray-tracing equations
# #saveμ_cb: save μ, dμ values at each timestep
# function save_func(u,t,integrator)
# 	r, θ, ϕ, freq = u
# 	v = refractive_index(r, θ, χ, f)
# 	dip = atan(2.0*tan(λ))
# 	ψ = π/2.0 + dip + χ
# 	vcat(v,dip,ψ)#,r,λ)
# end

# saved_μ = SavedValues(Float64, Array{Float64,1})
# saveμ_cb = SavingCallback(save_func, saved_μ)

cb = CallbackSet(re_cb) #, saveμ_cb)

## solve ODEs and plot
u0 = [re+1.0e+6, 1.0*pi/4, 0.0, 1.0, 1.0, 0.0, 0.0]	# r0, θ0, ϕ0, ρ_r0, ρ_θ0, ρ_ϕ0, T0
p = [1000.0]										# f0
tspan = (0.0,5.0e+9)

hasel_prob = ODEProblem(haselgrove!,u0,tspan,p)
hasel_soln = solve(hasel_prob,  callback=cb)#, reltol=1e-7, dtmax = 1e6, dtmin = 1e-8)

## scratch
# plasma model:
# superposition of ionosphere and plasmasphere. ionosphere is modeled as isotropic exponential decay; plasmasphere model from Carpenter and Anderson 1992, extened to L = 1
r = [re:1000:(10*re);]
λ = 0
L = r./(re*cos(λ)^2)

Kp_max = 3                             # for example
Lppi = 5.6 - 0.46*Kp_max                # plasmapause inner limit in L
Lppo = 5.75                             # placeholder
d = 0                                   # day number
R̄ = 90                                  # 13-month average sunspot number
t = 2                                   # magnetic local time

ne_iono = (1.8e5.*exp.(-4.183119.*((r./re).-1.0471))) # cm^-3

ne_Lppi = 10^((-0.3145*Lppi + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-Lppi)/1.5))

# individual components
ne_plasma_1 = @. 10 .^((-0.3145.*L + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2 .-L)./1.5))
ne_plasma_2 = ne_Lppi*10 .^(-1.0.*(L .- Lppi)./0.1)
ne_plasma_3 = (5800 + 300*t).*L.^(-4.5) + (1 .- exp.((2 .-L)./10))

plot(L, ne_plasma_1, yaxis=:log, label="plasmasphere")
plot!(L, ne_plasma_2, yaxis=:log, label="plasmapause inner limit")
plot!(L, ne_plasma_3, yaxis=:log, label="plasmasphere outer limit")
plot!(L,ne_iono, yaxis=:log, label="ionosphere")

min_index = findmin(abs.(ne_plasma_2 - ne_plasma_3))[2]
Lppo = L[min_index]

ne_plasma = ones(length(L))

for i = 1:length(L)
    if L[i] <= Lppi
        log_ne = (-0.3145*L[i] + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-L[i])/1.5)
        ne_plasma[i] = 10^(log_ne)
    elseif Lppi < L[i] <= Lppo
        ne_plasma[i] = ne_Lppi*10^((Lppi-L[i])/0.1)
    elseif Lppo < L[i]
        ne_plasma[i] = (5800 + 300*t)*L[i]^(-4.5) + (1 - exp((2-L[i])/10))
    else
        ne_plasma[i] = 0.0
    end
end

ne_total = ne_iono + ne_plasma

using Plots
plot(L, ne_iono, yaxis=:log, ylims=[1,10^6], label="ionosphere")
plot!(L, ne_plasma, yaxis=:log, label="plasmasphere")
plot!(L, ne_total, yaxis=:log, label="aggregate n_e")

##
r = re*2
θ = π/4
ϕ = 0.0
Bvec = magnetic_field(r,θ,ϕ)
Br, Bθ, Bϕ = Bvec
Bmag = sqrt(Br^2 + Bθ^2 + Bϕ^2)

# 1b. plasma
# assumptions: cold plasma, 2.25 < L < 8, midnight sector, based on Carpenter and Anderson 1992
#Lppo = initialize_plasmasphere(Kp_max, Lppi, d, R̄, mlt) # include this line if any Lppo variables change during integration (they shouldn't)

L = r/(re*cos(π/2 - θ)^2)
ne_iono = (1.8e5.*exp.(-4.183119.*((r./re).-1.0471))) # cm^-3

if L <= Lppi
	log_ne = (-0.3145*L + 3.9043) + (0.15*(cos((2*π*(d+9)))/365 - 0.5*cos((4*π*(d+9)))/365) + 0.00127*R̄ - 0.0635)*exp((2-L)/1.5)
	ne_plasma = 10^(log_ne)
elseif Lppi < L <= Lppo
	ne_plasma = ne_Lppi*10^((Lppi-L)/0.1)
elseif Lppo < L
	ne_plasma = (5800 + 300*mlt)*L^(-4.5) + (1 - exp((2-L)/10))
else
	ne_plasma = 0.0
end

n_e = (ne_plasma + ne_iono)*1e6 	# add plasmasphere + ionosphere, convert from cm^-3 to m^-3
n_p = n_e

# electron plasma angular frequency squared
ω_e2 = (n_e*(e^2))/(eps*me)
# proton plasma angular frequency squared
ω_p2 = (n_p*(e^2))/(eps*mp)

# electron angular gyrofrequency
Ω_e = (e*Bmag)/(me)
# proton angular gyrofrequency
Ω_p = (e*Bmag)/(mp)

# convert wave frequency to gyrofrequency
ω = 2*π*freq

# 2. define dispersion relation
# 2a. calculate ψ from ρ_r, ρ_θ, ρ_ϕ
ρ_r = 1.0
ρ_θ = 1.0
ρ_ϕ = 0.0
μvec = [ρ_r, ρ_θ, ρ_ϕ]
μmag = sqrt(ρ_r^2 + ρ_θ^2 + ρ_ϕ^2)
cos_ψ = (B ⋅ μvec)/(Bmag*μmag)
ψ = acos(cos_ψ)

# define R, L, P, D, S,
# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
#println("R = ",R)

# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
#println("L = ",L)

# P ≡ 1 - Σ(ωₖ²/ω²)
P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
#println("P = ",P)

# D ≡ ½(R-L); S ≡ ½(R+L)
D = (R - L)/2.0
S = (R + L)/2.0

# define parts of dispersion relation
# Aμ⁴ - Bμ² + C = 0
# A = Ssin²ψ + Pcos²ψ
A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0

# B = RLsin²ψ + PS(1 + cos²ψ)
B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0)

# C  = PRL
C = P*R*L

# solve the dispersion relation for μ:
# μ² = (B +- F)/2A
# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
F = sqrt(F2)

# Rice 1997: typical solution to dispersion relation with quadratic formula
μ2_minus = (B - F)/(2.0*A)
μ2_plus = (B + F)/(2.0*A)

μ_plus = try
	sqrt(abs(μ2_plus)) # abs() is not physical! for test only

catch err
	if isa(err,DomainError)
		# println("DomainError in μ_plus!")
		# println("B=", B)
		# println("F=", F)
		# println("A=", A)
	else
		# println(err)
	end
end

# Electron whistler case: ψ = 0, μ² = R; this is the μ_plus case
#μ = μ_minus		# EMIC case
μ = μ_plus		# electron whistler case

# Find dA/dψ, dB/dψ, dC/dψ, dμ/dψ
dAdψ = 2.0*(S-P)*sin(ψ)*cos(ψ)
dBdψ = 2.0*(R*L-P*S)*sin(ψ)*cos(ψ)
dCdψ = 0.0
#dμdψ = ((μ^4.0)*dAdψ-(μ^2.0)*dBdψ+dCdψ)/(4.0*A*(μ^3.0)-2.0*B*μ)

dFdψ = 1/(2.0*F)*((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ))
#dFdψ = sqrt(abs((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)))
dμdψ = 1/(2.0*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B + F)/(2*A^2))
# NOTE 11/24/2020: from Rice 1997; choose '+' for (B +- F) and (dBdψ +- dFdψ)

#DEBUG check values
#println("μ = ",μ)
#println("dμdψ = ", dμdψ)

# non-mutating output
u = [μ,dμdψ,ψ,Bvec]

# try unpacking u

μ1, dμdψ1, ψ1, B1 = u

B1[1]

max_stars_repo_name: andersontodds/ray-trace

max_issues_repo_name: andersontodds/ray-trace

max_forks_repo_name: andersontodds/ray-trace


text:
using DependencyParser
using Test
using DataStructures
using Flux

# FACTORIES

include("factories/dependency_tree.jl")
include("factories/sentence.jl")
include("factories/configuration.jl")
include("factories/gold_state.jl")
include("factories/model.jl")
include("factories/dep_parser.jl")

@testset "DependencyParser.jl" begin
    include("dependecy_parsing/arc_eager/moves.jl")
    include("dependecy_parsing/arc_eager/gold_state.jl")
    include("dependecy_parsing/connlu.jl")
    include("dependecy_parsing/model.jl")
    include("dependecy_parsing/dep_parser.jl")
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
#!/usr/bin/julia
# julia examples from LJR

f(x) = x^2 * cos(x)*tan(x)*besselj(x, 1.5)

@time println(f(2.5))
@time println(f(2.5))

A = rand(3, 3)
A * A
A .* A

A = rand(100, 10)
b = rand(100)
@time A\b
@time A\b

x = 10
println(2x)

x = y = z = 10
println( :(x=y=z=10) )
println( typeof( :(x=y=z=10) ))

println( :(x^2x) )
println( :(x^2*x) )
println( :(2x) )

println( macroexpand( :(@evalpoly x 2 3) ) )

println( @evalpoly x 1 2 3 )
println( 1 + 2x + 3x^2 )

function f(n)
   sum(0:3:n) + sum(0:5:n) - sum(0:15:n)
end

@time println( f(Int128(2)^64))
@time println( f(Int128(2)^64))

function nextfib(n)
   a, b = zero(n), one(n)
   while b < n
      a, b = b, a+b
   end
   return b
end

@time println( nextfib( Int128(2)^120 ))
@time println( nextfib( Int128(2)^120 ))

@code_llvm nextfib( Int128(2)^120 )

nheads = @parallel (+) for i=1:200000000
   Int(rand(Bool))
end
println(nheads)

x = 0
@simd for i=1:20000000
   x += Int(rand(Bool))
end
println(x)

a = SharedArray(Float64, 10)
@parallel for i=1:10
   a[i] = i
end
println(a)

#using Threads
#
#println(Threads.nthreads())
#Threads.@threads for i = 1:10
#   a[i] = Threads.threadid()
#end
#println(a)

#using ArrayFire
#a = rand(AFArray{Float32}, 100, 100)
#@time a^log(a)
#@time a^log(a)

x1 = rand(1000000);
x2 = rand(1000000);
@time dot(x1, x2);

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using Test, CacheServers

@testset "core APIs" begin
    include("Core.jl")
end

@testset "default" begin
    include("Default.jl")
end

@testset "get global server" begin
    s = get_server(DefaultServer, Int, Dict{Any, Any})
    @test s isa DefaultServer{Int, Dict{Any, Any}}
    cs = get_server(DefaultServer, Int, Dict{Any, Any})

    @test cs === s
end


text:
module ConfigLH

using Pkg, Printf

export add_to_path!, remove_from_path!, shared_dir, test_dir

include("computer.jl")
include("package.jl")

"""
add_to_path!

Add a directory to the LOAD_PATH without duplication
"""
function add_to_path!(newDir :: String)
	if !in(newDir, LOAD_PATH)
		push!(LOAD_PATH, newDir);
	end
	return nothing
end

"""
remove_from_path!

Remove a directory to the LOAD_PATH, if it exists
"""
function remove_from_path!(newDir :: String)
	idx = findfirst(x -> x == newDir, LOAD_PATH);
	if !isnothing(idx)
		deleteat!(LOAD_PATH, idx);
	end
	return nothing
end

include("projectLH.jl")

end # module


text:
function DEP(stochasticprogram::StochasticProgram{2}, structure::DeterministicEquivalent; optimizer = nothing)
    # Ensure stochastic program has been generated at this point
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    if optimizer == nothing
        return structure.model
    end
    structure_model = copy(structure.model)
    set_optimizer(structure_model, optimizer)
    return structure_model
end


text:
""" [Brown2002](@cite) """
function ARM_SGP_z(::Type{FT}) where {FT}
    z_in = FT[0.0, 50.0, 350.0, 650.0, 700.0, 1300.0, 2500.0, 5500.0 ] #LES z is in meters
    return z_in
end
""" [Brown2002](@cite) """
function ARM_SGP_θ_liq_ice(::Type{FT}) where {FT}
    z_in = ARM_SGP_z(FT)
    θ_liq_ice_in = FT[299.0, 301.5, 302.5, 303.53, 303.7, 307.13, 314.0, 343.2] # K
    profile = Dierckx.Spline1D(z_in, θ_liq_ice_in; k = 1)
    return profile
end
""" [Brown2002](@cite) """
ARM_SGP_u(::Type{FT}) where {FT} = z -> FT(10)

""" [Brown2002](@cite) """
function ARM_SGP_q_tot(::Type{FT}) where {FT}
    z_in = ARM_SGP_z(FT)
    r_in = FT[15.2,15.17,14.98,14.8,14.7,13.5,3.0,3.0] ./ 1000 # qt should be in kg/kg
    q_tot_in = r_in ./ (1 .+ r_in)
    profile = Dierckx.Spline1D(z_in, q_tot_in; k = 1)
    return profile
end
""" [Brown2002](@cite) """
function ARM_SGP_tke(::Type{FT}) where {FT}
    return z -> if z <= 2500.0
        FT(1) - z / 3000
    else
        FT(0)
    end
end

""" [Brown2002](@cite) """
function ARM_SGP_time(::Type{FT}) where {FT}
    t_in = FT[0.0, 3.0, 6.0, 9.0, 12.0, 14.5] .* 3600 #LES time is in sec
    return t_in
end

""" [Brown2002](@cite) """
function ARM_SGP_dTdt(::Type{FT}) where {FT}
    t_in = ARM_SGP_time(FT)
    # Advective forcing for theta [K/h] converted to [K/sec]
    AT_in = FT[0.0, 0.0, 0.0, -0.08, -0.016, -0.016] ./ 3600
    # Radiative forcing for theta [K/h] converted to [K/sec]
    RT_in = FT[-0.125, 0.0, 0.0, 0.0, 0.0, -0.1] ./ 3600
    dTdt_A = Dierckx.Spline1D(t_in, AT_in; k = 1)
    dTdt_R = Dierckx.Spline1D(t_in, RT_in; k = 1)
    return (t, z) -> if z <= 1000.0
        dTdt_A(t)+dTdt_R(t)
    elseif z > 1000.0 && z <= 2000.0
        (dTdt_A(t)+dTdt_R(t)) * (1 - (z - 1000) / 1000)
    else
        FT(0)
    end
end

""" [Brown2002](@cite) """
function ARM_SGP_dqtdt(::Type{FT}) where {FT}
    t_in = ARM_SGP_time(FT)
    # Radiative forcing for qt converted to [kg/kg/sec]
    Rqt_in = FT[0.08, 0.02, 0.04, -0.1, -0.16, -0.3] ./ 1000 ./ 3600
    dqtdt = Dierckx.Spline1D(t_in, Rqt_in; k = 1)
    return (Π, t, z) -> if z <= 1000.0
        dqtdt(t) * Π
    elseif z > 1000.0 && z <= 2000.0
        dqtdt(t) * Π * (1 - (z - 1000) / 1000)
    else
        FT(0)
    end
end

""" [Brown2002](@cite) """
function ARM_SGP_shf(::Type{FT}) where {FT}
    t_Sur_in = FT[0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5] .* 3600 #LES time is in sec
    shf = FT[-30.0, 90.0, 140.0, 140.0, 100.0, -10, -10] # W/m^2
    profile = Dierckx.Spline1D(t_Sur_in, shf; k = 1)
    return profile
end

""" [Brown2002](@cite) """
function ARM_SGP_lhf(::Type{FT}) where {FT}
    t_Sur_in = FT[0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5] .* 3600 #LES time is in sec
    lhf = FT[5.0, 250.0, 450.0, 500.0, 420.0, 180.0, 0.0] # W/m^2
    profile = Dierckx.Spline1D(t_Sur_in, lhf; k = 1)
    return profile
end


text:
using Test, PolynomialBases

if !haskey(ENV, "JULIA_PKGEVAL") # sympy is not install on https://github.com/JuliaComputing/NewPkgEval.jl
  import SymPy

  tol = 5.e-15

  for p in 0:6
      basis_sympy = LobattoLegendre(p, SymPy.Sym)
      basis_float = LobattoLegendre(p, Float64)
      @test maximum(abs.( float.(basis_sympy.nodes) - basis_float.nodes )) < tol
      @test maximum(abs.( float.(basis_sympy.weights) - basis_float.weights )) < tol
      @test maximum(abs.( float.(basis_sympy.baryweights) - basis_float.baryweights )) < tol
      @test maximum(abs.( float.(basis_sympy.D) - basis_float.D, )) < 2tol
  end
  @test_throws ArgumentError LobattoLegendre(7, SymPy.Sym)

  for p in 0:4
      basis_sympy = GaussLegendre(p, SymPy.Sym)
      basis_float = GaussLegendre(p, Float64)
      @test maximum(abs.( float.(basis_sympy.nodes) - basis_float.nodes )) < tol
      @test maximum(abs.( float.(basis_sympy.weights) - basis_float.weights )) < tol
      @test maximum(abs.( float.(basis_sympy.baryweights) - basis_float.baryweights )) < tol
      @test maximum(abs.( float.(basis_sympy.D) - basis_float.D, )) < 2tol
  end
  @test_throws ArgumentError GaussLegendre(5, SymPy.Sym)

  interpolation_matrix([-1, 1], LobattoLegendre(4, SymPy.Sym))
  interpolation_matrix([-1, 1], GaussLegendre(4, SymPy.Sym))
end


text:
using Compat
import Compat: unsafe_string, String, is_windows, is_unix

depsfile = joinpath(dirname(@__FILE__),"deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    f = open(depsfile,"w")
    if is_windows()
        path = replace(path, "\\", "\\\\")
    end
    println(f,"const xprs = \"$(path)\"")
    close(f)
end

aliases = Compat.ASCIIString["xprs"]

paths_to_try = copy(aliases)

for a in aliases
    if is_windows()
        if haskey(ENV, "XPRESSDIR")
            push!(paths_to_try, joinpath(ENV["XPRESSDIR"],"bin",string(a,".",Libdl.dlext)))
        end
    elseif is_unix()
        if haskey(ENV, "XPRESS")
            push!(paths_to_try, joinpath(ENV["XPRESS"],"lib",string("lib",a,".",Libdl.dlext)))
            push!(paths_to_try, joinpath(ENV["XPRESS"],"..","lib",string("lib",a,".",Libdl.dlext)))
        end
    end
end

println(paths_to_try)
found = false
for l in paths_to_try[2:end]
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        found = true
        write_depsfile(l)
        break
    end
end

if !found
    error("Unable to locate Xpress installation, please check your enviromental variable (XPRESSDIR in windows and XPRESS in unix). Note that this must be downloaded separately from fico.com")
end


text:
## example for README

using AutoAligns
aa = AutoAlign(align = Dict(1 => left, :default => right));
for (i, r) in zip([1, 100, 10000], ["a", "bb", "ccc"])
    print(aa, r)
    for j in 1:5
        print(aa, "  ", i + j) # padding
    end
    println(aa)
end
print(STDOUT, aa)


text:
# @recipe function f(Vol::SSD.ConeMantle{T}, n_aux_lines =0) where T
#
#     rStart = Vol.cone.r_interval.left
#     rStop = Vol.cone.r_interval.right
#     φStart = Vol.cone.φ_interval.left
#     φStop = Vol.cone.φ_interval.right
#     zStart = Vol.cone.z_interval.left
#     zStop = Vol.cone.z_interval.right
#
#     orientation = Vol.cone.orientation
#
#     if orientation == :bottom_left || orientation == :top_right
#         @series begin
#             line_3d(rStop,rStart,φStart,φStart,zStart,zStop)
#         end
#         @series begin
#             line_3d(rStop,rStart,φStop,φStop,zStart,zStop)
#         end
#         for ia in 0:n_aux_lines
#             @series begin
#                 line_3d(rStop,rStart,φStart+ia*(φStop-φStart)/(n_aux_lines+1),φStart+ia*(φStop-φStart)/(n_aux_lines+1),zStart,zStop)
#             end
#         end
#     end
# end


text:
__precompile__()

module TimedSound

using Lazy: @>>, @>, @_
using Unitful

export @>>, @>, @_

depsjl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsjl)
  include(depsjl)
else
  error("TimedSound not properly installed. "*
        "Please run\nPkg.build(\"TimedSound\")")
end

include(joinpath(@__DIR__,"units.jl"))
include(joinpath(@__DIR__,"timing.jl"))
include(joinpath(@__DIR__,"sound.jl"))
include(joinpath(@__DIR__,"stream.jl"))
include(joinpath(@__DIR__,"playback.jl"))

include(joinpath(@__DIR__,"audio.jl"))

const localunits = Unitful.basefactors
const localpromotion = Unitful.promotion
function __init__()
  merge!(Unitful.basefactors,localunits)
  merge!(Unitful.promotion, localpromotion)
end

end


text:
using HorizonSideRobots

function ugl(r)
    Flag = true
    y = true
    while !isborder(r,Sud) || Flag
        if !isborder(r,Sud)
            move!(r,Sud)
        end
        if isborder(r,Sud)
            Flag =  Obhod(r,2)
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    if !isborder(r,Nord)
        move!(r,Nord)
    else
        Obhod(r,0)
    end
end

function Obhod(r,HS)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 3) % 4))
        move!(r,HorizonSide((HS + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(HS))
        move!(r,HorizonSide(HS))
    end
    while !isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))
        move!(r,HorizonSide(HS))
    end
    for i = 1:x
        move!(r,HorizonSide((HS+1) % 4))
    end
    return Flag
end

function count(r)
    count1 = 0
    count2 = 0
    StopCount = true
    while  !isborder(r,Ost) || StopCount
        if isborder(r,Sud)
            count1 += 1
        else
            if count1 != 0
                count2 += 1
            end
            count1 = 0
        end
        if !isborder(r,Ost)
            move!(r,Ost)
        else
            StopCount = false
            if count1 != 0
                count2 += 1
            end
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    return(count2)

end

function all(r)
    ugl(r)
    Flag2 = true
    Summ = 0
    while !isborder(r,Nord) || Flag2
        Summ += count(r)
        if !isborder(r,Nord)
            move!(r,Nord)
        else
            Flag2 = Obhod(r,0)
        end
    end
    return Summ
end
all(r)


text:
function FindData(fpath::String;query="", ignore="", hideList::Bool=false)

    # ESTABLISH "QUERY" AND "IGNORE" LISTS
    if isa(query, String)
        query = split(query,"&")
    elseif !isa(query, Array)
        println("'query' should be: (1) a string with '&' separating multiple substrings where necessary, (2) an array of strings with containing multiple substrings")
    end
    if isa(ignore, String)
        ignore = split(ignore,"&")
    elseif !isa(ignore, Array)
        println("'ignore' should be: (1) a string with '&' separating multiple substrings where necessary, (2) an array of strings with containing multiple substrings")
    end
    # FIND DIRECTORY CONTENTS
    dir = readdir(fpath);

    # INITIALISE
    dirQry = []

    # IDENTIFY AND LOCATE "QUERY" ARGUEMNTS
    for j in 1:length(query)
        if query[j] != ""
            rslt = [occursin(query[j],dir[i]) for i in 1:length(dir)]
            dirQry = vcat(dirQry,dir[rslt])
        end
    end

    # REMOVE DUPLICATES
    dirQry = unique(dirQry)

    # IDENTIFY AND REMOVE "IGNORE" ARGUEMNTS
    for j in 1:length(ignore)
        if ignore[j] != ""
            rslt = [occursin(ignore[j],dirQry[i]) for i in 1:length(dirQry)]
            dirQry = dirQry[rslt.<true]
        end
    end

    if hideList == false
        println(" ")
        [println("$(i) \t →     $(dirQry[i])") for i in 1:length(dirQry)];
        println(" ")
    end

    return dirQry
end

function LoadData(fname;pth::String="",resolution::Integer=60,from::String="",days::Integer=0,UKDS_1amCorrection::Bool=false)

    ## ESTABLISH PATH TO DATA FROM FIRST TWO ARGUMENTS
    fullPath = joinpath(pth,fname)
    ## ADD ".csv" IF MISSING
    if fullPath[end-3:end] != ".csv"
        fullPath = fullPath*".csv"
    end
    ## INITIAL LOAD OF "raw" ASSUMES HEADERS EXIST (EXCEPTION FOLLOWS)
    raw = CSV.read(fullPath)

    ## CHECK FOR EXCEPTION ⇒ IF NO HEADERS EXIST IN ORIGINAL, RELOAD "raw" WITH "header=0" (TRY/CATCH REQUIRED TO AVOID ERROR)
    header = 1
    try
        firstColName = String(names(raw)[1])
        firstColValue = parse(Float64,firstColName)
        if firstColValue .> 38000 && firstColValue < 44000
            raw = CSV.read(fullPath,header=0)
            header=0
        else
            error("1. Time vector needs to be identified / convert input file to have headers (including time or date in time vector heading)")
        end
    catch
    end

    ## IDENTIFY DATE COLUMN NAME
    if header == 0
        dateTimeCheck = [col for col in names(raw) if occursin("Column1",String(col))]
    else
        ## CHECK IF ORIGINAL FILE INCLUDES TIME VECTOR
        dateTimeCheck = [col for col in names(raw) if occursin("date",lowercase(String(col))) || occursin("time",lowercase(String(col)))]
    end

    ## WHERE PREVIOUSLY CLIPPED DATA INCLUDES EXTRA ROW (MIDNIGHT THE FOLLOWING DAY), REMOVE THIS TO AVOID BLANK DAY AT END OF MATRIX
    # println(raw[dateTimeCheck][1])
    modulusOfT1 = raw[!,dateTimeCheck][!,1][1]%1 ## DETERMINE WHETHER t[1] IS MIDNIGHT OR NOT (WHOLE NUMBER)
    spareRows = (size(raw,1))%(24*60/resolution) ## DETERMINE WHETHER THERE IS AN EXTRA ROW (MIDNIGHT THE FOLLOWING DAY)
    if modulusOfT1==0 && spareRows==1 ## IF BOTH TRUE, REMOVE LAST (SPARE) ROW
        raw = raw[1:end-1,:]
    end

    ## EXTRACT "t" AND "dat" AS ARRAYS
    if !isempty(dateTimeCheck)
        t = raw[!,dateTimeCheck][!,1]
        dataTag = [col for col in names(raw) if col!=dateTimeCheck[1]]
        dat = raw[!,dataTag][!,1]
    else
        error("2. Time vector needs to be identified / convert input file to have headers (including time or date in time vector heading)")
    end

    ## FILL GAPS WITH NANS
    tResSec = minimum(round.(unique(diff(t))*24*60*60))
    ptsPerDay = 24*60*60/tResSec
    # println("Input dat resolution: $(tResSec)sec \t($(ptsPerDay) points per day)")
    padStart = Int(round(t[1]%1*ptsPerDay))
    ind = 1 .+ padStart .+ Int.(round.(1*(t.-t[1])*ptsPerDay)/1)
    padLen = Int(ceil(ind[end]/ptsPerDay)*ptsPerDay)
    dataFillNans = fill!(Array{Float64,1}(undef,padLen),NaN)
    dataFillNans[ind] = dat

    tFill = range(floor(t[1]),stop=(ceil(t[end])),length=padLen+1)
    tFill = tFill[1:end-1]

    ## CONVERT FROM XLS DATE FORMAT TO JULIA (DateTime) FORMAT
    tFill = util.xlsDate.(tFill)

    if from != ""
        ind = 1:length(tFill)
        if isa(from,Date)
            from = DateTime(from)
        elseif isa(from,String)
            from = DateTime(parse(Int64,from[1:4]),parse(Int64,from[6:7]),parse(Int64,from[9:10]))
        end
        pos = ind[from.==tFill]
        tFill = tFill[(pos[1]):end]
        dataFillNans = dataFillNans[(pos[1]):end]

        if days != 0
            ind = 1:length(tFill)
            tTarg = DateTime(Dates.UTM(Dates.value(tFill[1])+(days)*(1000*60*60*24)))
            pos = ind[tFill.==tTarg]
            tFill = tFill[1:(pos[1]-1)]
            dataFillNans = dataFillNans[1:(pos[1]-1)]
        end
    end

    ## RESHAPE DATA TO DAILY ARRAY
    nStep = Int((24*60*60/tResSec))
    nDays = Int(length(dataFillNans)/nStep)
    dataFillNans = reshape(dataFillNans, nStep, nDays)

    ## LINEARLY INTERPOLATE GAP AT 01:00
    if UKDS_1amCorrection == true
        dataFillNans[3,:] = [mean([dataFillNans[2,c],dataFillNans[4,c]]) for c in  1:size(dataFillNans,2)]
    end

    ## EXTRACT SMARTMETER NAME
    n = split(split(fullPath,"/")[end],".")[1]

    ## CREATE STUCT
    smartMeter = SmartMeter(n, dataFillNans, tFill, Second(tResSec), Day(nDays))
    return smartMeter
end

function InitProject(;N,pth::String,query::String)
    ## INIT PROJECT STRUCT
    global proj = Project_SMSPS()
    proj.path = Dict("data"=>pth)
    proj.fullDataList = FindData(proj.path["data"]; query=query,hideList=true)
    proj.dataIndices = vcat([N]...)
    proj.verbose = Array{DataFrame}(undef,length(proj.dataIndices))

    nDwlg = length(N)
    #=
    n = fill!(Array{Int32}(undef,nDwlg),0)
    region = fill!(Array{String}(undef,nDwlg),"tbc")
    numDays = fill!(Array{Int16}(undef,nDwlg),0)
    completeness = fill!(Array{Float64}(undef,nDwlg),0)
    minVal = fill!(Array{Float64}(undef,nDwlg),0)
    maxVal = fill!(Array{Float64}(undef,nDwlg),0)
    spare1 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare2 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare3 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare4 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare5 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare6 = Array{Union{Missing,Int16}}(missing,nDwlg)
    spare7 = fill!(Array{String}(undef,nDwlg),"-")
    spare8 = fill!(Array{String}(undef,nDwlg),"-")
    =#

    proj.smartMeters = [LoadData(proj.fullDataList[i],pth=proj.path["data"],UKDS_1amCorrection=true) for i in proj.dataIndices]
    #=
    for i in 1:length(proj.smartMeters)
        # trgFile = proj.fullDataList[i]
        # proj.smartMeters = LoadData((pathData,trgFile))
        n[i] = parse(Int32,split(split(proj.smartMeters[i].id,"_")[end],".")[1])
        region[i] = split(proj.smartMeters[i].id,"_")[3]
        numDays[i] = Int(size(proj.smartMeters[i].dat)[2])
        completeness[i] = 1-sum([1 for j in proj.smartMeters[i].dat if isnan(j)])/(size(proj.smartMeters[i].dat)[1]*size(proj.smartMeters[i].dat)[2])
        minVal[i] = minimum(replace(proj.smartMeters[i].dat, NaN=>Inf))
        maxVal[i] = maximum(replace(proj.smartMeters[i].dat, NaN=>0))
    end

    proj.summaryData = DataFrames.DataFrame(DwellingID=n[:], Region=region,NumDays=numDays[:],
                    Completeness=completeness[:],MinValue=minVal[:],MaxValue=maxVal[:],
                    HtgCont=spare1,HtgProg=spare2,HtgRand=spare3,HtgIdle=spare4,HtgVoid=spare5,
                    NumProgs=spare6, HtgPrograms=spare7, RandomReoccur=spare8)

    proj.htgSchedules = DataFrame(Shape=String[], Dwelling=Int32[], ID=Int16[], HtgSPHigh=Bool[]) ## BLANK DATAFRAME FOR BELOW

    =#
    return proj
end


text:
using SymbolicUtils: Rewriters

const KEEP = typemin(Int)

function alias_eliminate_graph!(state::TransformationState)
    mm = linear_subsys_adjmat(state)
    size(mm, 1) == 0 && return nothing, mm # No linear subsystems

    @unpack graph, var_to_diff = state.structure

    ag, mm = alias_eliminate_graph!(graph, complete(var_to_diff), mm)
    return ag, mm
end

# For debug purposes
function aag_bareiss(sys::AbstractSystem)
    state = TearingState(sys)
    mm = linear_subsys_adjmat(state)
    return aag_bareiss!(state.structure.graph, complete(state.structure.var_to_diff), mm)
end

function alias_elimination(sys)
    state = TearingState(sys; quick_cancel=true)
    ag, mm = alias_eliminate_graph!(state)
    ag === nothing && return sys

    fullvars = state.fullvars
    graph = state.structure.graph

    subs = OrderedDict()
    for (v, (coeff, alias)) in pairs(ag)
        subs[fullvars[v]] = iszero(coeff) ? 0 : coeff * fullvars[alias]
    end

    dels = Set{Int}()
    eqs = collect(equations(state))
    for (ei, e) in enumerate(mm.nzrows)
        vs = 𝑠neighbors(graph, e)
        if isempty(vs)
            # remove empty equations
            push!(dels, e)
        else
            rhs = mapfoldl(+, pairs(nonzerosmap(@view mm[ei, :]))) do (var, coeff)
                iszero(coeff) && return 0
                return coeff * fullvars[var]
            end
            eqs[e] = 0 ~ rhs
        end
    end
    dels = sort(collect(dels))
    deleteat!(eqs, dels)

    dict = Dict(subs)
    for (ieq, eq) in enumerate(eqs)
        eqs[ieq] = fixpoint_sub(eq.lhs, dict) ~ fixpoint_sub(eq.rhs, dict)
    end

    newstates = []
    for j in eachindex(fullvars)
        if !(j in keys(ag))
            isdervar(state.structure, j) || push!(newstates, fullvars[j])
        end
    end

    sys = state.sys
    @set! sys.eqs = eqs
    @set! sys.states = newstates
    @set! sys.observed = [observed(sys); [lhs ~ rhs for (lhs, rhs) in pairs(subs)]]
    return invalidate_cache!(sys)
end

"""
$(SIGNATURES)

Find the first linear variable such that `𝑠neighbors(adj, i)[j]` is true given
the `constraint`.
"""
@inline function find_first_linear_variable(
        M::SparseMatrixCLIL,
        range,
        mask,
        constraint,
    )
    eadj = M.row_cols
    for i in range
        vertices = eadj[i]
        if constraint(length(vertices))
            for (j, v) in enumerate(vertices)
                (mask === nothing || mask[v]) && return (CartesianIndex(i, v), M.row_vals[i][j])
            end
        end
    end
    return nothing
end

@inline function find_first_linear_variable(
        M::AbstractMatrix,
        range,
        mask,
        constraint,
    )
    for i in range
        row = @view M[i, :]
        if constraint(count(!iszero, row))
            for (v, val) in enumerate(row)
                iszero(val) && continue
                if mask === nothing || mask[v]
                    return CartesianIndex(i, v), val
                end
            end
        end
    end
    return nothing
end

function find_masked_pivot(variables, M, k)
    r = find_first_linear_variable(M, k:size(M,1), variables, isequal(1))
    r !== nothing && return r
    r = find_first_linear_variable(M, k:size(M,1), variables, isequal(2))
    r !== nothing && return r
    r = find_first_linear_variable(M, k:size(M,1), variables, _->true)
    return r
end

"""
    AliasGraph

When eliminating variables, keeps track of which variables where eliminated in
favor of which others.

Currently only supports elimination as direct aliases (+- 1).

We represent this as a dict from eliminated variables to a (coeff, var) pair
representing the variable that it was aliased to.
"""
struct AliasGraph <: AbstractDict{Int, Pair{Int, Int}}
    aliasto::Vector{Union{Int, Nothing}}
    eliminated::Vector{Int}
    function AliasGraph(nvars::Int)
        new(fill(nothing, nvars), Int[])
    end
end

Base.length(ag::AliasGraph) = length(ag.eliminated)

function Base.getindex(ag::AliasGraph, i::Integer)
    r = ag.aliasto[i]
    r === nothing && throw(KeyError(i))
    coeff, var = (sign(r), abs(r))
    if var in keys(ag)
        # Amortized lookup. Check if since we last looked this up, our alias was
        # itself aliased. If so, just adjust the alias table.
        ac, av = ag[var]
        nc = ac * coeff
        ag.aliasto[var] = nc > 0 ? av : -av
    end
    return (coeff, var)
end

function Base.iterate(ag::AliasGraph, state...)
    r = Base.iterate(ag.eliminated, state...)
    r === nothing && return nothing
    c = ag.aliasto[r[1]]
    return (r[1] => (c == 0 ? 0 :
                     c >= 0 ? 1 :
                             -1, abs(c))), r[2]
end

function Base.setindex!(ag::AliasGraph, v::Integer, i::Integer)
    @assert v == 0
    if ag.aliasto[i] === nothing
        push!(ag.eliminated, i)
    end
    ag.aliasto[i] = 0
    return 0=>0
end

function Base.setindex!(ag::AliasGraph, p::Pair{Int, Int}, i::Integer)
    (c, v) = p
    @assert v != 0 && c in (-1, 1)
    if ag.aliasto[i] === nothing
        push!(ag.eliminated, i)
    end
    ag.aliasto[i] = c > 0 ? v : -v
    return p
end

function Base.get(ag::AliasGraph, i::Integer, default)
    i in keys(ag) || return default
    return ag[i]
end

struct AliasGraphKeySet <: AbstractSet{Int}
    ag::AliasGraph
end
Base.keys(ag::AliasGraph) = AliasGraphKeySet(ag)
Base.iterate(agk::AliasGraphKeySet, state...) = Base.iterate(agk.ag.eliminated, state...)
function Base.in(i::Int, agk::AliasGraphKeySet)
    aliasto = agk.ag.aliasto
    1 <= i <= length(aliasto) && aliasto[i] !== nothing
end

count_nonzeros(a::AbstractArray) = count(!iszero, a)

# N.B.: Ordinarily sparse vectors allow zero stored elements.
# Here we have a guarantee that they won't, so we can make this identification
count_nonzeros(a::SparseVector) = nnz(a)

function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
    mm = copy(mm_orig)
    is_linear_equations = falses(size(AsSubMatrix(mm_orig), 1))
    for e in mm_orig.nzrows
        is_linear_equations[e] = true
    end

    # Variables that are highest order differentiated cannot be states of an ODE
    is_not_potential_state = isnothing.(var_to_diff)
    is_linear_variables = copy(is_not_potential_state)
    for i in 𝑠vertices(graph); is_linear_equations[i] && continue
        for j in 𝑠neighbors(graph, i)
            is_linear_variables[j] = false
        end
    end
    solvable_variables = findall(is_linear_variables)

    function do_bareiss!(M, Mold=nothing)
        rank1 = rank2 = nothing
        pivots = Int[]
        function find_pivot(M, k)
            if rank1 === nothing
                r = find_masked_pivot(is_linear_variables, M, k)
                r !== nothing && return r
                rank1 = k - 1
            end
            if rank2 === nothing
                r = find_masked_pivot(is_not_potential_state, M, k)
                r !== nothing && return r
                rank2 = k - 1
            end
            return find_masked_pivot(nothing, M, k)
        end
        function find_and_record_pivot(M, k)
            r = find_pivot(M, k)
            r === nothing && return nothing
            push!(pivots, r[1][2])
            return r
        end
        function myswaprows!(M, i, j)
            Mold !== nothing && swaprows!(Mold, i, j)
            swaprows!(M, i, j)
        end
        bareiss_ops = ((M,i,j)->nothing, myswaprows!, bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
        rank3, = bareiss!(M, bareiss_ops; find_pivot=find_and_record_pivot)
        rank1 = something(rank1, rank3)
        rank2 = something(rank2, rank3)
        (rank1, rank2, rank3, pivots)
    end

    return mm, solvable_variables, do_bareiss!(mm, mm_orig)
end

function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
    # Step 1: Perform bareiss factorization on the adjacency matrix of the linear
    #         subsystem of the system we're interested in.
    #
    # Let `m = the number of linear equations` and `n = the number of
    # variables`.
    #
    # `do_bareiss` conceptually gives us this system:
    # rank1 | [ M₁₁  M₁₂ | M₁₃  M₁₄ ]   [v₁] = [0]
    # rank2 | [ 0    M₂₂ | M₂₃  M₂₄ ] P [v₂] = [0]
    # -------------------|------------------------
    # rank3 | [ 0    0   | M₃₃  M₃₄ ]   [v₃] = [0]
    #         [ 0    0   | 0    0   ]   [v₄] = [0]
    mm, solvable_variables, (rank1, rank2, rank3, pivots) =
        aag_bareiss!(graph, var_to_diff, mm_orig)

    # Step 2: Simplify the system using the Bareiss factorization
    ag = AliasGraph(size(mm, 2))
    for v in setdiff(solvable_variables, @view pivots[1:rank1])
        ag[v] = 0
    end

    # Kind of like the backward substitution, but we don't actually rely on it
    # being lower triangular. We eliminate a variable if there are at most 2
    # variables left after the substitution.
    diff_to_var = invview(var_to_diff)
    function lss!(ei::Integer)
        vi = pivots[ei]
        # the differentiated variable cannot be eliminated
        islowest = isnothing(diff_to_var[vi]) && isnothing(var_to_diff[vi])
        locally_structure_simplify!((@view mm[ei, :]), vi, ag, islowest)
    end

    # Step 2.1: Go backwards, collecting eliminated variables and substituting
    #         alias as we go.
    foreach(lss!, reverse(1:rank2))

    # Step 2.2: Sometimes Bareiss can make the equations more complicated.
    #         Go back and check the original matrix. If this happened,
    #         Replace the equation by the one from the original system,
    #         but be sure to also run `lss!` again, since we only ran that
    #         on the Bareiss'd matrix, not the original one.
    reduced = mapreduce(|, 1:rank2; init=false) do ei
        if count_nonzeros(@view mm_orig[ei, :]) < count_nonzeros(@view mm[ei, :])
            mm[ei, :] = @view mm_orig[ei, :]
            return lss!(ei)
        end
        return false
    end

    # Step 2.3: Iterate to convergence.
    #         N.B.: `lss!` modifies the array.
    # TODO: We know exactly what variable we eliminated. Starting over at the
    #       start is wasteful. We can lookup which equations have this variable
    #       using the graph.
    reduced && while any(lss!, 1:rank2); end

    # Step 3: Reflect our update decisions back into the graph
    for (ei, e) in enumerate(mm.nzrows)
        set_neighbors!(graph, e, mm.row_cols[ei])
    end

    return ag, mm
end

function exactdiv(a::Integer, b)
    d, r = divrem(a, b)
    @assert r == 0
    return d
end

function locally_structure_simplify!(adj_row, pivot_col, ag, may_eliminate)
    pivot_val = adj_row[pivot_col]
    iszero(pivot_val) && return false

    nirreducible = 0
    alias_candidate = 0

    # N.B.: Assumes that the non-zeros iterator is robust to modification
    # of the underlying array datastructure.
    for (var, val) in pairs(nonzerosmap(adj_row))
        # Go through every variable/coefficient in this row and apply all aliases
        # that we have so far accumulated in `ag`, updating the adj_row as
        # we go along.
        var == pivot_col && continue
        iszero(val) && continue
        alias = get(ag, var, nothing)
        if alias === nothing
            nirreducible += 1
            alias_candidate = val => var
            continue
        end
        (coeff, alias_var) = alias
        # `var = coeff * alias_var`, so we eliminate this var.
        adj_row[var] = 0
        if alias_var != 0
            # val * var = val * (coeff * alias_var) = (val * coeff) * alias_var
            val *= coeff
            # val * var + c * alias_var + ... = (val * coeff + c) * alias_var + ...
            new_coeff = (adj_row[alias_var] += val)
            if alias_var < var
                # If this adds to a coeff that was not previously accounted for,
                # and we've already passed it, make sure to count it here. We
                # need to know if there are at most 2 terms left after this
                # loop.
                #
                # We're relying on `var` being produced in sorted order here.
                nirreducible += 1
                alias_candidate = new_coeff => alias_var
            end
        end
    end

    if may_eliminate && nirreducible <= 1
        # There were only one or two terms left in the equation (including the
        # pivot variable). We can eliminate the pivot variable.
        #
        # Note that when `nirreducible <= 1`, `alias_candidate` is uniquely
        # determined.
        if alias_candidate !== 0
            d, r = divrem(alias_candidate[1], pivot_val)
            if r == 0 && (d == 1 || d == -1)
                alias_candidate = -d => alias_candidate[2]
            else
                return false
            end
        end
        ag[pivot_col] = alias_candidate
        zero!(adj_row)
        return true
    end
    return false
end

swap!(v, i, j) = v[i], v[j] = v[j], v[i]

function getcoeff(vars, coeffs, var)
    for (vj, v) in enumerate(vars)
        v == var && return coeffs[vj]
    end
    return 0
end

"""
$(SIGNATURES)

Use Kahn's algorithm to topologically sort observed equations.

Example:
```julia
julia> @variables t x(t) y(t) z(t) k(t)
(t, x(t), y(t), z(t), k(t))

julia> eqs = [
           x ~ y + z
           z ~ 2
           y ~ 2z + k
       ];

julia> ModelingToolkit.topsort_equations(eqs, [x, y, z, k])
3-element Vector{Equation}:
 Equation(z(t), 2)
 Equation(y(t), k(t) + 2z(t))
 Equation(x(t), y(t) + z(t))
```
"""
function topsort_equations(eqs, states; check=true)
    graph, assigns = observed2graph(eqs, states)
    neqs = length(eqs)
    degrees = zeros(Int, neqs)

    for 𝑠eq in 1:length(eqs); var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            # 𝑠eq => 𝑑eq
            degrees[𝑑eq] += 1
        end
    end

    q = Queue{Int}(neqs)
    for (i, d) in enumerate(degrees)
        d == 0 && enqueue!(q, i)
    end

    idx = 0
    ordered_eqs = similar(eqs, 0); sizehint!(ordered_eqs, neqs)
    while !isempty(q)
        𝑠eq = dequeue!(q)
        idx+=1
        push!(ordered_eqs, eqs[𝑠eq])
        var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            degree = degrees[𝑑eq] = degrees[𝑑eq] - 1
            degree == 0 && enqueue!(q, 𝑑eq)
        end
    end

    (check && idx != neqs) && throw(ArgumentError("The equations have at least one cycle."))

    return ordered_eqs
end

function observed2graph(eqs, states)
    graph = BipartiteGraph(length(eqs), length(states))
    v2j = Dict(states .=> 1:length(states))

    # `assigns: eq -> var`, `eq` defines `var`
    assigns = similar(eqs, Int)

    for (i, eq) in enumerate(eqs)
        lhs_j = get(v2j, eq.lhs, nothing)
        lhs_j === nothing && throw(ArgumentError("The lhs $(eq.lhs) of $eq, doesn't appear in states."))
        assigns[i] = lhs_j
        vs = vars(eq.rhs)
        for v in vs
            j = get(v2j, v, nothing)
            j !== nothing && add_edge!(graph, i, j)
        end
    end

    return graph, assigns
end

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(eqs, dict)
    sub = Base.Fix2(fixpoint_sub, dict)
    map(eq->eq.lhs ~ sub(eq.rhs), eqs)
end