text:
using Base.Math: @horner

# Compute the sine integral: ∫_0^x sin(t)/t dt,
# and the cosine integral: γ + log x + ∫_0^x (cos(t)-1)/t dt,
# using the rational approximants tabulated in:
#     A.J. MacLeod, "Rational approximations, software and test methods for
#     sine and cosine integrals," Numer. Algor. 12, pp. 259--272 (1996).
#         http://dx.doi.org/10.1007/BF02142806
#         https://link.springer.com/article/10.1007/BF02142806
#
# Note: the second zero of Ci(x) has a typo that is fixed:
#
# r1 = 3.38418 04228 51186 42639 78511 46402 in the article, but is in fact:
#
# r1 = 3.38418 04225 51186 42639 78511 46402.
#

function sinint(x::Float64)
    t = x*x
    if t ≤ 36.0
        return x * @horner(t, 1.00000000000000000000E0,
                             -0.44663998931312457298E-1,
                              0.11209146443112369449E-2,
                             -0.13276124407928422367E-4,
                              0.85118014179823463879E-7,
                             -0.29989314303147656479E-9,
                              0.55401971660186204711E-12,
                             -0.42406353433133212926E-15) /
                   @horner(t, 1.00000000000000000000E0,
                              0.10891556624243098264E-1,
                              0.59334456769186835896E-4,
                              0.21231112954641805908E-6,
                              0.54747121846510390750E-9,
                              0.10378561511331814674E-11,
                              0.13754880327250272679E-14,
                              0.10223981202236205703E-17)
    elseif t ≤ 144.0
        invt = inv(t)
        return copysign(π/2, x) - cos(x) *
                                         @horner(invt, 0.99999999962173909991E0,
                                                       0.36451060338631902917E3,
                                                       0.44218548041288440874E5,
                                                       0.22467569405961151887E7,
                                                       0.49315316723035561922E8,
                                                       0.43186795279670283193E9,
                                                       0.11847992519956804350E10,
                                                       0.45573267593795103181E9) /
                                    (x * @horner(invt, 1.00000000000000000000E0,
                                                       0.36651060273229347594E3,
                                                       0.44927569814970692777E5,
                                                       0.23285354882204041700E7,
                                                       0.53117852017228262911E8,
                                                       0.50335310667241870372E9,
                                                       0.16575285015623175410E10,
                                                       0.11746532837038341076E10)) -
                           sin(x)*invt * @horner(invt, 0.99999999920484901956E0,
                                                       0.51385504875307321394E3,
                                                       0.92293483452013810811E5,
                                                       0.74071341863359841727E7,
                                                       0.28142356162841356551E9,
                                                       0.49280890357734623984E10,
                                                       0.35524762685554302472E11,
                                                       0.79194271662085049376E11,
                                                       0.17942522624413898907E11) /
                                         @horner(invt, 1.00000000000000000000E0,
                                                       0.51985504708814870209E3,
                                                       0.95292615508125947321E5,
                                                       0.79215459679762667578E7,
                                                       0.31977567790733781460E9,
                                                       0.62273134702439012114E10,
                                                       0.54570971054996441467E11,
                                                       0.18241750166645704670E12,
                                                       0.15407148148861454434E12)
    elseif t < Inf
        invt = inv(t)
        return copysign(π/2, x) - cos(x) / x * (1.0 -
                                                @horner(invt, 0.19999999999999978257E1,
                                                              0.22206119380434958727E4,
                                                              0.84749007623988236808E6,
                                                              0.13959267954823943232E9,
                                                              0.10197205463267975592E11,
                                                              0.30229865264524075951E12,
                                                              0.27504053804288471142E13,
                                                              0.21818989704686874983E13) /
                                                @horner(invt, 1.00000000000000000000E0,
                                                              0.11223059690217167788E4,
                                                              0.43685270974851313242E6,
                                                              0.74654702140658116258E8,
                                                              0.58580034751805687471E10,
                                                              0.20157980379272098841E12,
                                                              0.26229141857684496445E13,
                                                              0.87852907334918467516E13)*invt) -
                         sin(x) * invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
                                                              0.96527746044997139158E4,
                                                              0.56077626996568834185E7,
                                                              0.15022667718927317198E10,
                                                              0.19644271064733088465E12,
                                                              0.12191368281163225043E14,
                                                              0.31924389898645609533E15,
                                                              0.25876053010027485934E16,
                                                              0.12754978896268878403E16) /
                                                @horner(invt, 1.00000000000000000000E0,
                                                              0.16287957674166143196E4,
                                                              0.96636303195787870963E6,
                                                              0.26839734750950667021E9,
                                                              0.37388510548029219241E11,
                                                              0.26028585666152144496E13,
                                                              0.85134283716950697226E14,
                                                              0.11304079361627952930E16,
                                                              0.42519841479489798424E16)*invt)
    elseif isnan(x)
        return NaN
    else
        return copysign(π/2, x)
    end
end

function cosint(x::Float64)
    t, r0, r1 = x*x, 0.616505485620716233797110404100, 3.384180422551186426397851146402
    r01, r02 = 0.6162109375, 0.29454812071623379711E-3
    r11, r12 = 3.3837890625, 0.39136005118642639785E-3
    if x < 0.0
        throw(DomainError(x, "`x` must be nonnegative."))
    elseif x ≤ 3.0
        return log(x/r0) + ((x - r01) - r02) * (x + r0) *
                                       @horner(t, -0.24607411378767540707E0,
                                                   0.72113492241301534559E-2,
                                                  -0.11867127836204767056E-3,
                                                   0.90542655466969866243E-6,
                                                  -0.34322242412444409037E-8,
                                                   0.51950683460656886834E-11) /
                                        @horner(t, 1.00000000000000000000E0,
                                                   0.12670095552700637845E-1,
                                                   0.78168450570724148921E-4,
                                                   0.29959200177005821677E-6,
                                                   0.73191677761328838216E-9,
                                                   0.94351174530907529061E-12)
    elseif x ≤ 6.0
        return log(x/r1) + ((x - r11) - r12) * (x + r1) *
                                       @horner(t, -0.15684781827145408780E0,
                                                   0.66253165609605468916E-2,
                                                  -0.12822297297864512864E-3,
                                                   0.12360964097729408891E-5,
                                                  -0.66450975112876224532E-8,
                                                   0.20326936466803159446E-10,
                                                  -0.33590883135343844613E-13,
                                                   0.23686934961435015119E-16) /
                                        @horner(t, 1.00000000000000000000E0,
                                                   0.96166044388828741188E-2,
                                                   0.45257514591257035006E-4,
                                                   0.13544922659627723233E-6,
                                                   0.27715365686570002081E-9,
                                                   0.37718676301688932926E-12,
                                                   0.27706844497155995398E-15)
    elseif x ≤ 12.0
        invt = inv(t)
        return sin(x) * @horner(invt, 0.99999999962173909991E0,
                                      0.36451060338631902917E3,
                                      0.44218548041288440874E5,
                                      0.22467569405961151887E7,
                                      0.49315316723035561922E8,
                                      0.43186795279670283193E9,
                                      0.11847992519956804350E10,
                                      0.45573267593795103181E9) /
                   (x * @horner(invt, 1.00000000000000000000E0,
                                      0.36651060273229347594E3,
                                      0.44927569814970692777E5,
                                      0.23285354882204041700E7,
                                      0.53117852017228262911E8,
                                      0.50335310667241870372E9,
                                      0.16575285015623175410E10,
                                      0.11746532837038341076E10)) -
        cos(x) * invt * @horner(invt, 0.99999999920484901956E0,
                                      0.51385504875307321394E3,
                                      0.92293483452013810811E5,
                                      0.74071341863359841727E7,
                                      0.28142356162841356551E9,
                                      0.49280890357734623984E10,
                                      0.35524762685554302472E11,
                                      0.79194271662085049376E11,
                                      0.17942522624413898907E11) /
                        @horner(invt, 1.00000000000000000000E0,
                                      0.51985504708814870209E3,
                                      0.95292615508125947321E5,
                                      0.79215459679762667578E7,
                                      0.31977567790733781460E9,
                                      0.62273134702439012114E10,
                                      0.54570971054996441467E11,
                                      0.18241750166645704670E12,
                                      0.15407148148861454434E12)
    elseif x < Inf
        invt = inv(t)
        return sin(x)/x * (1.0 - @horner(invt, 0.19999999999999978257E1,
                                               0.22206119380434958727E4,
                                               0.84749007623988236808E6,
                                               0.13959267954823943232E9,
                                               0.10197205463267975592E11,
                                               0.30229865264524075951E12,
                                               0.27504053804288471142E13,
                                               0.21818989704686874983E13) /
                                 @horner(invt, 1.00000000000000000000E0,
                                               0.11223059690217167788E4,
                                               0.43685270974851313242E6,
                                               0.74654702140658116258E8,
                                               0.58580034751805687471E10,
                                               0.20157980379272098841E12,
                                               0.26229141857684496445E13,
                                               0.87852907334918467516E13)*invt) -
            cos(x)*invt * (1.0 - @horner(invt, 0.59999999999999993089E1,
                                               0.96527746044997139158E4,
                                               0.56077626996568834185E7,
                                               0.15022667718927317198E10,
                                               0.19644271064733088465E12,
                                               0.12191368281163225043E14,
                                               0.31924389898645609533E15,
                                               0.25876053010027485934E16,
                                               0.12754978896268878403E16) /
                                 @horner(invt, 1.00000000000000000000E0,
                                               0.16287957674166143196E4,
                                               0.96636303195787870963E6,
                                               0.26839734750950667021E9,
                                               0.37388510548029219241E11,
                                               0.26028585666152144496E13,
                                               0.85134283716950697226E14,
                                               0.11304079361627952930E16,
                                               0.42519841479489798424E16)*invt)
    elseif isnan(x)
        return NaN
    else
        return 0.0
    end
end

for f in (:sinint, :cosint)
    @eval begin
        ($f)(x::Float32) = Float32(($f)(Float64(x)))
        ($f)(x::Float16) = Float16(($f)(Float64(x)))
        ($f)(x::Real) = ($f)(float(x))
        ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))
    end
end

"""
    sinint(x)

Compute the sine integral function of `x`, defined by ``\\operatorname{Si}(x) := \\int_0^x\\frac{\\sin t}{t} dt``
for real `x`.
"""
sinint

"""
    cosint(x)

Compute the cosine integral function of `x`, defined by ``\\operatorname{Ci}(x) := \\gamma + \\log x + \\int_0^x \\frac{\\cos t - 1}{t} dt``
for real `x > 0`, where ``\\gamma`` is the Euler-Mascheroni constant.
"""
cosint


text:
"""
Describes a hydropower unit.

Hydropower is both a purpose and a possible dam output. The optimisation process could ignore the purpose side of
hydropower, and only consider it as a possible (constant) output; but optimising this purpose needs the dedicated
penstock. As a consequence, the
"""
immutable HydropowerUnit
  id::Symbol
  max_discharge::Float64
  max_power::Float64
  efficiency::Float64

  # Private constructor, without units nor defaut values.
  function HydropowerUnit(id::Symbol, max_discharge::Float64, max_power::Float64, efficiency::Float64)
    if efficiency < 0. || efficiency > 1.
      error("Wrong efficiency: " * string(efficiency) * " not between 0.0 and 1.0")
    end

    new(id, max_discharge, max_power, efficiency)
  end
end

# Public constructor, with units and keyword arguments.
HydropowerUnit(;id::Symbol=:hydropower,
                max_discharge::Union{Discharge, Float64}=error("Missing maximum discharge"),
                max_power::Union{Power, Float64}=error("Missing maximum power not specified"),
                efficiency::Float64=error("Missing efficiency")
          ) = HydropowerUnit(id, _from_unitful(max_discharge), _from_unitful(max_power), efficiency)
HydropowerUnit(max_discharge::Union{Discharge, Float64}, max_power::Union{Power, Float64}, efficiency::Float64) =
  HydropowerUnit(max_discharge=max_discharge, max_power=max_power, efficiency=efficiency)

copy(hp::HydropowerUnit; id::Symbol=hp.id,
     max_discharge::Union{Discharge, Float64}=hp.max_discharge, max_power::Union{Power, Float64}=hp.max_power,
     efficiency::Float64=hp.efficiency) =
  # Always call the internal constructor!
  Reservoir(id, _from_unitful(max_discharge), _from_unitful(max_power), efficiency)

getId(hp::HydropowerUnit) = hp.id
getPurposeId(hp::HydropowerUnit) = symbol(getId(hp), "Purpose")
getDamOutputId(hp::HydropowerUnit) = symbol(getId(hp), "DamOutput")
getMaximumDischarge(hp::HydropowerUnit) = hp.max_discharge
getMaximumPower(hp::HydropowerUnit) = hp.max_power
getEfficiency(hp::HydropowerUnit) = hp.efficiency

id(hp::HydropowerUnit) = getId(hp)
purposeId(hp::HydropowerUnit) = getPurposeId(hp)
damOutputId(hp::HydropowerUnit) = getDamOutputId(hp)
maxDischarge(hp::HydropowerUnit) = getMaximumDischarge(hp)
maxPower(hp::HydropowerUnit) = getMaximumPower(hp)
efficiency(hp::HydropowerUnit) = getEfficiency(hp)

ext: jl


text:
using DataFrames, CSV, Gadfly
df = CSV.read("/Users/Matthieu/Dropbox/Github/FixedEffectModels.jl/benchmark/benchmark.csv")

df.R = df.R ./ df.Julia
df.Stata = df.Stata ./ df.Julia
df.Julia = df.Julia ./ df.Julia

mdf = melt(df[!, [:Command, :Julia, :R, :Stata]], :Command)
mdf = rename(mdf, :variable => :Language)
p = plot(mdf, x = "Language", y = "value", color = "Command", Guide.ylabel("Time (Ratio to Julia)"), Guide.xlabel("Model"), Guide.yticks(ticks= [1, 5, 10, 15]))

draw(PNG("/Users/Matthieu/Dropbox/Github/FixedEffectModels.jl/benchmark/fixedeffectmodels_benchmark.png", 8inch, 5inch, dpi=300), p)


text:
println("# VCM TESTS...")

# run demo code at first
println(" * try demo...")

D = 40
K = 2
N = 100

C = zeros(D, K)
C[1:20, 1] = 1
C[21:40, 2] = 1

cc = zeros(1,D)
cc[1:20] = 1
cc[21:40] = 2

mu_x = 0.0;
sigma_x = 1.0;
mu_g = 0.0;
sigma_g = 1.0;
mu_noise = 0.0;
sigma_noise = 0.1;

# set G and X (factor loadings, and latent factors) of test data
G = mu_g + sigma_g * randn(D, K)
X = mu_x + sigma_x * randn(K, N)

# set noise
ϵ = mu_noise + sigma_noise * randn(D,N);

# test data
Y = (G .* C) * X + ϵ;

train(VCM(), Gibbs(), IncrementalInitialisation(), Y);
@time train(VCM(), Gibbs(), IncrementalInitialisation(), Y);

println("# FINISHED")


text:
importall CVModule

type CVTimer
    tfd::Float64
    tc::Float64
    ts::Float64
    tr::Float64
    tt::Float64
    function CVTimer()
        this = new();
        tfd = 0;
        tc = 0;
        ts = 0;
        tr = 0;
        tt = 0;
        return this
    end
end

function main()

# options
# filename = "arterytree.csv"; # default artery data file for new simulation
filename = "run1.mat"; # filename for restarting simulation
rstflag = "yes" # restarting from scratch or previous simulation
hemoflag = "no" # 10% hemorrhage from left femoral artery
saveflag = "yes" # save solution to .mat file
coupleflag = "no" # coupling to 3D liver tissue model
timeflag = "yes" # solver timing

# build solution struct
system = CVModule.buildall(filename;numbeatstotal=1,restart=rstflag,injury=hemoflag);

# timers
times = CVTimer();

# collect all IDs of terminal/non-terminal branches
terms = Int64[];
splits = Int64[];
for i = 1:length(system.branches.ID)
    if isempty(system.branches.children[i])
        push!(terms,i);
    else
        push!(splits,i);
    end
end

n = system.solverparams.nstart;

# coupling to 3D liver model
if coupleflag == "yes"
    ctx=ZMQ.Context();
    sender=ZMQ.Socket(ctx,ZMQ.REQ);
    ZMQ.connect(sender,"tcp://127.0.0.1:5555");
end

tic();
# solver loop
ret1 = Array{Float64}(system.solverparams.JL-2); # placeholder arrays for interior 1d points
ret2 = Array{Float64}(system.solverparams.JL-2);
while system.solverparams.numbeats < system.solverparams.numbeatstotal
    tic();
    for i = 1:length(system.branches.ID)
        CVModule.predictorfluxes!(system.branches.Fp[i],system.branches.Q[i][:,n+1],
            system.branches.A[i][:,n+1],system.branches.beta[i][end],system.solverparams.rho,
            system.solverparams.acols,system.solverparams.qcols);
        CVModule.predictorstep!(system.branches.Aforward[i],system.branches.Qforward[i],
            system.branches.Abackward[i],system.branches.Qbackward[i],
            system.branches.A[i][:,n+1],system.branches.Q[i][:,n+1],
            system.solverparams.h,system.branches.k[i],system.branches.Fp[i],
            system.solverparams.colsint,system.solverparams.acolspre,system.solverparams.qcolspre,
            system.solverparams.diffusioncoeff,system.solverparams.mu,system.solverparams.rho);
        CVModule.correctorfluxes!(system.branches.Fbarforward[i],system.branches.Fbarbackward[i],
            system.branches.Qforward[i],system.branches.Aforward[i],system.branches.Qbackward[i],
            system.branches.Abackward[i],system.branches.beta[i][end],system.solverparams.rho,
            system.solverparams.acolscor,system.solverparams.qcolscor);
        CVModule.correctorstep!(ret1,ret2,
            system.branches.A[i][:,n+1],system.branches.Q[i][:,n+1],system.branches.Fbarforward[i],
            system.branches.Fbarbackward[i],system.branches.Qforward[i],system.branches.Aforward[i],
            system.branches.Qbackward[i],system.branches.Abackward[i],system.solverparams.colsint,
            system.solverparams.acolscor,system.solverparams.qcolscor,system.solverparams.h,system.branches.k[i],
            system.solverparams.diffusioncoeff,system.solverparams.mu,system.solverparams.rho);
        system.branches.A[i][system.solverparams.colsint,n+2] .= ret1;
        system.branches.Q[i][system.solverparams.colsint,n+2] .= ret2;
    end
    times.tfd += toq();
    tic();
    if hemoflag == "no"
        CVModule.applyendbcs!(system,n,terms);
    elseif hemoflag == "yes"
        CVModule.applyendbcs!(system,n,terms,hemoflag);
    end
    times.tc += toq();
    tic();
    for i = 1:length(splits)
        children = system.branches.children[splits[i]];
        if length(children) == 1
            f = CVModule.fsingle;
            J = CVModule.Jsingle;
        elseif length(children) == 2
            f = CVModule.fdouble;
            J = CVModule.Jdouble;
        elseif length(children) == 3
            f = CVModule.ftriple;
            J = CVModule.Jtriple;
        elseif length(children) == 4
            f = CVModule.fquad;
            J = CVModule.Jquad;
        end
        W = zeros(length(children)+1);
        beta = zeros(length(children)+1);
        A0 = zeros(length(children)+1);
        c0 = zeros(length(children)+1);
        k = zeros(length(children)+1);
        Q = zeros(length(children)+1);
        Qo = zeros(length(children)+1);
        A = zeros(length(children)+1);
        Ao = zeros(length(children)+1);
        beta[1] = system.branches.beta[splits[i]][end];
        A0[1] = system.branches.A0[splits[i]][end];
        c0[1] = system.branches.c0[splits[i]][end];
        k[1] = system.branches.k[splits[i]][end];
        Q[1] = system.branches.Q[splits[i]][system.solverparams.JL,n+1];
        Qo[1] = system.branches.Q[splits[i]][system.solverparams.JL-1,n+1];
        A[1] = system.branches.A[splits[i]][system.solverparams.JL,n+1];
        Ao[1] = system.branches.A[splits[i]][system.solverparams.JL-1,n+1];
        for j = 1:length(children)
            beta[j+1] = system.branches.beta[children[j]][end];
            A0[j+1] = system.branches.A0[children[j]][end];
            c0[j+1] = system.branches.c0[children[j]][end];
            k[j+1] = system.branches.k[children[j]][end];
            Q[j+1] = system.branches.Q[children[j]][1,n+1];
            Qo[j+1] = system.branches.Q[children[j]][2,n+1];
            A[j+1] = system.branches.A[children[j]][1,n+1];
            Ao[j+1] = system.branches.A[children[j]][2,n+1];
        end
        CVModule.splitinvariants!(W,Q,A,Qo,Ao,beta,c0,system.solverparams.rho,
            k,system.solverparams.h,system.solverparams.diffusioncoeff,
            system.solverparams.mu,children);
        system.branches.W1[splits[i]] = W[1];
        system.branches.W2[children] .= W[2:end];
        Qnew = zeros(length(children)+1);
        Anew = zeros(length(children)+1);
        iters = Int64[0];
        if hemoflag == "no"
            CVModule.solvesplits!(iters,Qnew,Anew,children,Q,A,W,beta,A0,c0,
                system.solverparams.rho,f,J,system.solverparams.maxiter,
                system.solverparams.maxval,system.solverparams.epsJ,
                system.solverparams.epsN)
        elseif hemoflag == "yes"
            curloss = [0.];
            CVModule.solvesplits!(iters,Qnew,Anew,children,Q,A,W,beta,A0,c0,
                system.solverparams.rho,f,J,system.solverparams.maxiter,
                system.solverparams.maxval,system.solverparams.epsJ,
                system.solverparams.epsN,hemoflag,system.hemo.injured[splits[i]],
                system.solverparams.h,system.hemo.Ph,curloss);
            system.hemo.Vloss += curloss[1];
            if system.hemo.injured[splits[i]] == true
                system.branches.W2[splits[i]] = W[1];
                system.branches.W1[children] .= W[2:end];
            end
            CVModule.applytourniquet!(system,n); # turn off to allow continual bleeding
        end
        system.solverparams.totaliter += iters[1];
        system.branches.Q[splits[i]][system.solverparams.JL,n+2] = Qnew[1];
        system.branches.A[splits[i]][system.solverparams.JL,n+2] = Anew[1];
        for j = 1:length(children)
            system.branches.Q[children[j]][1,n+2] = Qnew[j+1];
            system.branches.A[children[j]][1,n+2] = Anew[j+1];
        end
    end
    times.ts += toq();
    tic();
    CVModule.arterialpressure!(system,n);
    CVModule.regulateall!(system,n,terms);
    if coupleflag == "yes"
        CVModule.senddata(system,n,sender);
    end
    n+=1;
    times.tr += toq();
end
times.tt += toc();

CVModule.updatevolumes!(system,n);

if coupleflag == "yes"
    ZMQ.close(sender)
    ZMQ.close(ctx)
end

if saveflag == "yes"
    file = MAT.matopen("run2.mat", "w")
    write(file, "system", system)
    close(file)
end

if timeflag == "yes"
    return system, n, times
elseif timeflag == "no"
    return system, n
end

end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module ClimateEasy

# Main file for the ClimateEasy module that provides a range of useful
# commands written in Julia that allows one to easily perform commonly-needed
# tasks such as:
#   - Converting dates and times to string format (e.g. YYYYMMDD)
#   - Generating vectors of dates
#   - Decomposing datestrings into their relevant components
#   - Plotting on maps (requires CartoPy in conda package manager)

## Modules Used
using Dates, DelimitedFiles, Printf, Statistics
using PyCall, PyPlot

## Adding relevant Python Functions
const ccrs  = PyNULL()
const pyhdf = PyNULL()

function __init__()
    copy!(ccrs,pyimport_conda("cartopy","cartopy"))
    copy!(pyhdf,pyimport_conda("pyhdf","pyhdf","conda-forge"))
end

## Exporting functions
export
       yr2str, mo2str, dy2str, hr2str, mi2str, yrmo2str, ymd2str, ymdhm2str,
       yrmo2dir, ymd2dir, yrdy2dir,
       regionload, regioninfodisplay, regionshortname, regionfullname, regionparent,
       regionbounds, regionpoint, regiongrid, regionisglobe,
       ispointinregion, isgridinregion,
       regionextract, regionextractpoint, regionextractgrid,
       from0360to180, from180to0360,
       convert2hourly, convert23hourly, convert2daily,
       convert2regionmean, convert2zonalmean, convert2meridionalmean,
       mollweide, robinson, coastlines,
       h4read

include("dnt2str.jl")
include("regions.jl")
include("processing.jl")
include("hdf4.jl")
include("mapping.jl")
include("misc.jl")

end # module


text:
using Compat.Printf
using Compat.LinearAlgebra

include("MBHFMM2D.jl")

using PyPlot

## parameters

lambda = 1.0e-2

ifcharge = false
ifdipole = false
ifquad = false
ifoct = true
ifpot = true
ifgrad = true
ifhess = true

## define source and target points

ns = 10000
nt = 2000

# sources on an ellipse, targets inside

src = zeros(Float64,2,ns)
rnorm = zeros(Float64,2,ns)
dsdt = zeros(Float64,ns)
targ = zeros(Float64,2,nt)

charge = zeros(Float64,1)
dipstr = zeros(Float64,1)
dipvec = zeros(Float64,2,1)
quadstr = zeros(Float64,1)
quadvec = zeros(Float64,3,1)
octstr = zeros(Float64,ns) # octopole strength
octvec = zeros(Float64,4,ns) # octopole "direction"

# storage for solution

pottarg1 = zeros(Float64,nt)
gradtarg1 = zeros(Float64,2,nt)
hesstarg1 = zeros(Float64,3,nt)
pottarg2 = zeros(Float64,nt)
gradtarg2 = zeros(Float64,2,nt)
hesstarg2 = zeros(Float64,3,nt)

center = [6.1;-1.3]

h = 2.0*pi/ns
a = 1.123
b = 0.7

for i = 1:ns
    t = (i-1)*h
    src[:,i] = center + [a*cos(t);b*sin(t)]
    dx = -a*sin(t)
    dy = b*cos(t)
    dsdt[i] = sqrt(dx^2+dy^2)
    rnorm[1,i] = dy/dsdt[i]
    rnorm[2,i] = -dx/dsdt[i]
    nu = rnorm[:,i]
    tau = [-nu[2];nu[1]]
    octvec[:,i] = boxfmm2d_formoctvec(nu,nu,tau)
    octstr[i] = h*dsdt[i]*cos(5*t)
end

for i = 1:nt
    xt = randn(); yt = randn();
    elnorm = sqrt((xt/a)^2 + (yt/b)^2)
    sc = rand();
    xt = sc*xt/elnorm; yt = sc*yt/elnorm
    targ[:,i] = center + [xt;yt]
end

## make parameters structure

fmmpars = MBHFMM2DParams(lambda,src,targ,ifcharge, ifdipole,
                         ifquad, ifoct, charge, dipstr,
                         dipvec, quadstr, quadvec, octstr,
                         octvec, iprec = 3, ifalltarg = false)

## direct calculation

@printf "PERFORMING DIRECT CALCULATION ...\n"

@time mbhfmm2d_direct!(fmmpars,targ,ifpot,pottarg1,ifgrad,
                 gradtarg1,ifhess,hesstarg1)

## set up fmm

# maxboxes determines size of storage arrays, set it big
# for now
# maxnodes determines max number of sources/targets
# per box

@printf "FORMING FMM ...\n"

@time fmmstor = mbhfmm2d_form(fmmpars,maxnodes=30,maxboxes=100000)

@printf "EVALUATING FMM ...\n"

@time mbhfmm2d_targ!(fmmpars,fmmstor,targ,ifpot,pottarg2,ifgrad,
                 gradtarg2,ifhess,hesstarg2)

println("relative error, pot ",
        norm(pottarg1-pottarg2)/norm(pottarg1))
println("relative error, grad ",
        norm(gradtarg1-gradtarg2)/norm(gradtarg1))
println("relative error, hess ",
        norm(hesstarg1-hesstarg2)/norm(hesstarg1))


text:
###################################################################
#  Encoding related definitions & functions
###################################################################

encoding_mark_prop_spec = Dict(
  :aggregate   =>	AbstractString,
  :timeUnit    =>	AbstractString,
  :bin         =>	Union{Dict, Bool},
  :sort        =>	AbstractString,
  :scale       =>	Dict,
  :axis        =>	Union{Dict, Bool},
  :legend      =>	Union{Dict, Bool},
  :value       =>	Any,
  :field       =>	Symbol,
  :type        =>	AbstractString
)

for channel in [:x, :y, :color, :shape, :size, :text,
                :detail, :order, :path, :row, :column]
  fn = symbol("encoding_" * string(channel))
  path = (:encoding, channel)
  @eval ($fn)(;properties...) = _mkvis($path, encoding_mark_prop_spec, properties)
  @eval ($fn)(val::Symbol;properties...) = ($fn)(;vcat(properties, [(:field, val)])...)
  @eval ($fn)(val;properties...)         = ($fn)(;vcat(properties, [(:value, val)])...)

  for (ts, tn) in [(:quant, "quantitative"), (:temp,   "temporal"),
                   (:ord,   "ordinal")     , (:nominal, "nominal")]
    fn2 = symbol(string(fn) * "_" * string(ts))
    @eval ($fn2)(val;properties...)    = ($fn)(val;vcat(properties, [(:type, $tn)])...)
  end
end

"""
Sets encoding for a channel X with data type Y, with Y equal to `quant` (quantitative),
`temp` (temporal), `ord` (ordinal) or `nominal` (nominal).

```julia
encoding_*X*_*Y*(;properties...)
encoding_*X*_*Y*(val::Symbol;properties...) # binds to field [val]
encoding_*X*_*Y*(val;properties...) # sets to constant value [val]
```

Properties:

- `bin` (`Dict` or `Bool`) : Boolean flag for binning a quantitative field, or a bin property object for binning parameters. Default value: false
- `timeUnit` (`AbstractString`) : Time unit for a temporal field (e.g., year, yearmonth, month, hour). Default value: undefined (None)
- `aggregate` (`AbstractString`) : Aggregation function for the field (e.g., mean, sum, median, min, max, count). Default value: undefined (None)
- `sort` (`AbstractString`) : | Object	Sort order for a particular field. (Default value: "ascending")
  • For quantitative or temporal fields, this can be either "ascending" or , "descending"
  • For quantitative or temporal fields, this can be "ascending", "descending", "none", or a sort field definition object for sorting by an aggregate calculation of a specified sort field.
- `scale`(`Dict`)	: A property object for a scale of a mark property channel.
- `axis`(`Dict` or `Bool`) : Boolean flag for showing an axis (true by default), or a property object for an axis of a position channel (x or y) or a facet channel (row or column).
- `legend`(`Dict` or `Bool`) : Boolean flag for showing a legend (true by default), or a config object for a legend of a non-position mark property channel (color, size, or shape).

"""
encoding_x, encoding_x_quant, encoding_x_temp, encoding_x_ord, encoding_x_nominal,
encoding_y, encoding_y_quant, encoding_y_temp, encoding_y_ord, encoding_y_nominal,
encoding_color, encoding_color_quant, encoding_color_temp,
encoding_color_ord, encoding_color_nominal,
encoding_shape, encoding_shape_quant, encoding_shape_temp,
encoding_shape_ord, encoding_shape_nominal,
encoding_path, encoding_path_quant, encoding_path_temp,
encoding_path_ord, encoding_path_nominal,
encoding_size, encoding_size_quant, encoding_size_temp,
encoding_size_ord, encoding_size_nominal,
encoding_text, encoding_text_quant, encoding_text_temp,
encoding_text_ord, encoding_text_nominal,
encoding_detail, encoding_detail_quant, encoding_detail_temp,
encoding_detail_ord, encoding_detail_nominal,
encoding_order, encoding_order_quant, encoding_order_temp,
encoding_order_ord, encoding_order_nominal,
encoding_row, encoding_row_quant, encoding_row_temp,
encoding_row_ord, encoding_row_nominal,
encoding_column, encoding_column_quant, encoding_column_temp,
encoding_column_ord, encoding_column_nominal

# encoding_x(;properties...) = _mkvis((:encoding, :x), encoding_mark_prop_spec, properties)
#
# encoding_x(val::Symbol;properties...) = encoding_x(;vcat(properties, [(:field, val)])...)
# encoding_x(val;properties...)         = encoding_x(;vcat(properties, [(:value, val)])...)
#
# encoding_x_quant(val;properties...)   = encoding_x(val;vcat(properties, [(:type, "quantitative")])...)
# encoding_x_temp(val;properties...)    = encoding_x(val;vcat(properties, [(:type, "temporal")])...)
# encoding_x_ord(val;properties...)     = encoding_x(val;vcat(properties, [(:type, "ordinal")])...)
# encoding_x_nominal(val;properties...) = encoding_x(val;vcat(properties, [(:type, "nominal")])...)
#
# encoding_y(;properties...) = _mkvis((:encoding, :y), encoding_mark_prop_spec, properties)
#
# encoding_y(val::Symbol;properties...) = encoding_y(;vcat(properties, [(:field, val)])...)
# encoding_y(val;properties...)         = encoding_y(;vcat(properties, [(:value, val)])...)
#
# encoding_y_quant(val;properties...)   = encoding_y(val;vcat(properties, [(:type, "quantitative")])...)
# encoding_y_temp(val;properties...)    = encoding_y(val;vcat(properties, [(:type, "temporal")])...)
# encoding_y_ord(val;properties...)     = encoding_y(val;vcat(properties, [(:type, "ordinal")])...)
# encoding_y_nominal(val;properties...) = encoding_y(val;vcat(properties, [(:type, "nominal")])...)
#
# encoding_color(;properties...) = _mkvis((:encoding, :color), encoding_mark_prop_spec, properties)
#
# encoding_color(val::Symbol;properties...) = encoding_color(;vcat(properties, [(:field, val)])...)
# encoding_color(val;properties...)         = encoding_color(;vcat(properties, [(:value, val)])...)
#
# encoding_color_quant(val;properties...)   = encoding_color(val;vcat(properties, [(:type, "quantitative")])...)
# encoding_color_temp(val;properties...)    = encoding_color(val;vcat(properties, [(:type, "temporal")])...)
# encoding_color_ord(val;properties...)     = encoding_color(val;vcat(properties, [(:type, "ordinal")])...)
# encoding_color_nominal(val;properties...) = encoding_color(val;vcat(properties, [(:type, "nominal")])...)
#
# encoding_shape(;properties...) = _mkvis((:encoding, :shape), encoding_mark_prop_spec, properties)
#
# encoding_shape(val::Symbol;properties...) = encoding_shape(;vcat(properties, [(:field, val)])...)
# encoding_shape(val;properties...)         = encoding_shape(;vcat(properties, [(:value, val)])...)
#
# encoding_shape_quant(val;properties...)   = encoding_shape(val;vcat(properties, [(:type, "quantitative")])...)
# encoding_shape_temp(val;properties...)    = encoding_shape(val;vcat(properties, [(:type, "temporal")])...)
# encoding_shape_ord(val;properties...)     = encoding_shape(val;vcat(properties, [(:type, "ordinal")])...)
# encoding_shape_nominal(val;properties...) = encoding_shape(val;vcat(properties, [(:type, "nominal")])...)
#
# encoding_path(;properties...) = _mkvis((:encoding, :path), encoding_mark_prop_spec, properties)
#
# encoding_path(val::Symbol;properties...) = encoding_path(;vcat(properties, [(:field, val)])...)
# encoding_path(val;properties...)         = encoding_path(;vcat(properties, [(:value, val)])...)
#
# encoding_path_quant(val;properties...)   = encoding_path(val;vcat(properties, [(:type, "quantitative")])...)
# encoding_path_temp(val;properties...)    = encoding_path(val;vcat(properties, [(:type, "temporal")])...)
# encoding_path_ord(val;properties...)     = encoding_path(val;vcat(properties, [(:type, "ordinal")])...)
# encoding_path_nominal(val;properties...) = encoding_path(val;vcat(properties, [(:type, "nominal")])...)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## linalg.jl: Some generic Linear Algebra definitions

# For better performance when input and output are the same array
# See https://github.com/JuliaLang/julia/issues/8415#issuecomment-56608729
function generic_rmul!(X::AbstractArray, s::Number)
    @simd for I in eachindex(X)
        @inbounds X[I] *= s
    end
    X
end

function generic_lmul!(s::Number, X::AbstractArray)
    @simd for I in eachindex(X)
        @inbounds X[I] = s*X[I]
    end
    X
end

function generic_mul!(C::AbstractArray, X::AbstractArray, s::Number)
    if _length(C) != _length(X)
        throw(DimensionMismatch("first array has length $(_length(C)) which does not match the length of the second, $(_length(X))."))
    end
    for (IC, IX) in zip(eachindex(C), eachindex(X))
        @inbounds C[IC] = X[IX]*s
    end
    C
end

function generic_mul!(C::AbstractArray, s::Number, X::AbstractArray)
    if _length(C) != _length(X)
        throw(DimensionMismatch("first array has length $(_length(C)) which does not
match the length of the second, $(_length(X))."))
    end
    for (IC, IX) in zip(eachindex(C), eachindex(X))
        @inbounds C[IC] = s*X[IX]
    end
    C
end

mul!(C::AbstractArray, s::Number, X::AbstractArray) = generic_mul!(C, X, s)
mul!(C::AbstractArray, X::AbstractArray, s::Number) = generic_mul!(C, s, X)

"""
    rmul!(A::AbstractArray, b::Number)

Scale an array `A` by a scalar `b` overwriting `A` in-place.

# Examples
```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rmul!(A, 2)
2×2 Array{Int64,2}:
 2  4
 6  8
```
"""
rmul!(A::AbstractArray, b::Number) = generic_rmul!(A, b)

"""
    lmul!(a::Number, B::AbstractArray)

Scale an array `B` by a scalar `a` overwriting `B` in-place.

# Examples
```jldoctest
julia> B = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> lmul!(2, B)
2×2 Array{Int64,2}:
 2  4
 6  8
```
"""
lmul!(a::Number, B::AbstractArray) = generic_lmul!(a, B)

"""
    cross(x, y)
    ×(x,y)

Compute the cross product of two 3-vectors.

# Examples
```jldoctest
julia> a = [0;1;0]
3-element Array{Int64,1}:
 0
 1
 0

julia> b = [0;0;1]
3-element Array{Int64,1}:
 0
 0
 1

julia> cross(a,b)
3-element Array{Int64,1}:
 1
 0
 0
```
"""
cross(a::AbstractVector, b::AbstractVector) =
    [a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]]

"""
    triu(M)

Upper triangle of a matrix.

# Examples
```jldoctest
julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> triu(a)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0
 0.0  0.0  0.0  1.0
```
"""
triu(M::AbstractMatrix) = triu!(copy(M))

"""
    tril(M)

Lower triangle of a matrix.

# Examples
```jldoctest
julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a)
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0
 1.0  1.0  1.0  0.0
 1.0  1.0  1.0  1.0
```
"""
tril(M::AbstractMatrix) = tril!(copy(M))

"""
    triu(M, k::Integer)

Returns the upper triangle of `M` starting from the `k`th superdiagonal.

# Examples
```jldoctest
julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> triu(a,3)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> triu(a,-3)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
```
"""
triu(M::AbstractMatrix,k::Integer) = triu!(copy(M),k)

"""
    tril(M, k::Integer)

Returns the lower triangle of `M` starting from the `k`th superdiagonal.

# Examples
```jldoctest
julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a,3)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a,-3)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0
```
"""
tril(M::AbstractMatrix,k::Integer) = tril!(copy(M),k)

"""
    triu!(M)

Upper triangle of a matrix, overwriting `M` in the process.
See also [`triu`](@ref).
"""
triu!(M::AbstractMatrix) = triu!(M,0)

"""
    tril!(M)

Lower triangle of a matrix, overwriting `M` in the process.
See also [`tril`](@ref).
"""
tril!(M::AbstractMatrix) = tril!(M,0)

diff(a::AbstractVector) = [ a[i+1] - a[i] for i=1:length(a)-1 ]

"""
    diff(A::AbstractVector)
    diff(A::AbstractMatrix, dim::Integer)

Finite difference operator of matrix or vector `A`. If `A` is a matrix,
specify the dimension over which to operate with the `dim` argument.

# Examples
```jldoctest
julia> a = [2 4; 6 16]
2×2 Array{Int64,2}:
 2   4
 6  16

julia> diff(a,2)
2×1 Array{Int64,2}:
  2
 10

julia> diff(vec(a))
3-element Array{Int64,1}:
  4
 -2
 12
```
"""
function diff(A::AbstractMatrix, dim::Integer)
    if dim == 1
        [A[i+1,j] - A[i,j] for i=1:size(A,1)-1, j=1:size(A,2)]
    elseif dim == 2
        [A[i,j+1] - A[i,j] for i=1:size(A,1), j=1:size(A,2)-1]
    else
        throw(ArgumentError("dimension dim must be 1 or 2, got $dim"))
    end
end

diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to construct a diagonal matrix"))

###########################################################################################
# Inner products and norms

# special cases of vecnorm; note that they don't need to handle isempty(x)
function generic_vecnormMinusInf(x)
    s = start(x)
    (v, s) = next(x, s)
    minabs = norm(v)
    while !done(x, s)
        (v, s) = next(x, s)
        vnorm = norm(v)
        minabs = ifelse(isnan(minabs) | (minabs < vnorm), minabs, vnorm)
    end
    return float(minabs)
end

function generic_vecnormInf(x)
    s = start(x)
    (v, s) = next(x, s)
    maxabs = norm(v)
    while !done(x, s)
        (v, s) = next(x, s)
        vnorm = norm(v)
        maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm)
    end
    return float(maxabs)
end

function generic_vecnorm1(x)
    s = start(x)
    (v, s) = next(x, s)
    av = float(norm(v))
    T = typeof(av)
    sum::promote_type(Float64, T) = av
    while !done(x, s)
        (v, s) = next(x, s)
        sum += norm(v)
    end
    return convert(T, sum)
end

# faster computation of norm(x)^2, avoiding overflow for integers
norm_sqr(x) = norm(x)^2
norm_sqr(x::Number) = abs2(x)
norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x))

function generic_vecnorm2(x)
    maxabs = vecnormInf(x)
    (maxabs == 0 || isinf(maxabs)) && return maxabs
    s = start(x)
    (v, s) = next(x, s)
    T = typeof(maxabs)
    if isfinite(_length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
        sum::promote_type(Float64, T) = norm_sqr(v)
        while !done(x, s)
            (v, s) = next(x, s)
            sum += norm_sqr(v)
        end
        return convert(T, sqrt(sum))
    else
        sum = abs2(norm(v)/maxabs)
        while !done(x, s)
            (v, s) = next(x, s)
            sum += (norm(v)/maxabs)^2
        end
        return convert(T, maxabs*sqrt(sum))
    end
end

# Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p)
# (Not technically a "norm" for p < 1.)
function generic_vecnormp(x, p)
    s = start(x)
    (v, s) = next(x, s)
    if p > 1 || p < -1 # might need to rescale to avoid overflow
        maxabs = p > 1 ? vecnormInf(x) : vecnormMinusInf(x)
        (maxabs == 0 || isinf(maxabs)) && return maxabs
        T = typeof(maxabs)
    else
        T = typeof(float(norm(v)))
    end
    spp::promote_type(Float64, T) = p
    if -1 <= p <= 1 || (isfinite(_length(x)*maxabs^spp) && maxabs^spp != 0) # scaling not necessary
        sum::promote_type(Float64, T) = norm(v)^spp
        while !done(x, s)
            (v, s) = next(x, s)
            sum += norm(v)^spp
        end
        return convert(T, sum^inv(spp))
    else # rescaling
        sum = (norm(v)/maxabs)^spp
        while !done(x, s)
            (v, s) = next(x, s)
            sum += (norm(v)/maxabs)^spp
        end
        return convert(T, maxabs*sum^inv(spp))
    end
end

vecnormMinusInf(x) = generic_vecnormMinusInf(x)
vecnormInf(x) = generic_vecnormInf(x)
vecnorm1(x) = generic_vecnorm1(x)
vecnorm2(x) = generic_vecnorm2(x)
vecnormp(x, p) = generic_vecnormp(x, p)

"""
    vecnorm(A, p::Real=2)

For any iterable container `A` (including arrays of any dimension) of numbers (or any
element type for which `norm` is defined), compute the `p`-norm (defaulting to `p=2`) as if
`A` were a vector of the corresponding length.

The `p`-norm is defined as:
```math
\\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p}
```
with ``a_i`` the entries of ``A`` and ``n`` its length.

`p` can assume any numeric value (even though not all values produce a
mathematically valid vector norm). In particular, `vecnorm(A, Inf)` returns the largest value
in `abs(A)`, whereas `vecnorm(A, -Inf)` returns the smallest. If `A` is a matrix and `p=2`,
then this is equivalent to the Frobenius norm.

# Examples
```jldoctest
julia> vecnorm([1 2 3; 4 5 6; 7 8 9])
16.881943016134134

julia> vecnorm([1 2 3 4 5 6 7 8 9])
16.881943016134134
```
"""
function vecnorm(itr, p::Real=2)
    isempty(itr) && return float(norm(zero(eltype(itr))))
    if p == 2
        return vecnorm2(itr)
    elseif p == 1
        return vecnorm1(itr)
    elseif p == Inf
        return vecnormInf(itr)
    elseif p == 0
        return typeof(float(norm(first(itr))))(count(!iszero, itr))
    elseif p == -Inf
        return vecnormMinusInf(itr)
    else
        vecnormp(itr,p)
    end
end

"""
    vecnorm(x::Number, p::Real=2)

For numbers, return ``\\left( |x|^p \\right) ^{1/p}``.

# Examples
```jldoctest
julia> vecnorm(2, 1)
2

julia> vecnorm(-2, 1)
2

julia> vecnorm(2, 2)
2

julia> vecnorm(-2, 2)
2

julia> vecnorm(2, Inf)
2

julia> vecnorm(-2, Inf)
2
```
"""
@inline vecnorm(x::Number, p::Real=2) = p == 0 ? (x==0 ? zero(abs(x)) : oneunit(abs(x))) : abs(x)

function norm1(A::AbstractMatrix{T}) where T
    m, n = size(A)
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:n
            nrmj::Tsum = 0
            for i = 1:m
                nrmj += norm(A[i,j])
            end
            nrm = max(nrm,nrmj)
        end
    end
    return convert(Tnorm, nrm)
end
function norm2(A::AbstractMatrix{T}) where T
    m,n = size(A)
    if m == 1 || n == 1 return vecnorm2(A) end
    Tnorm = typeof(float(real(zero(T))))
    (m == 0 || n == 0) ? zero(Tnorm) : convert(Tnorm, svdvals(A)[1])
end
function normInf(A::AbstractMatrix{T}) where T
    m,n = size(A)
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for i = 1:m
            nrmi::Tsum = 0
            for j = 1:n
                nrmi += norm(A[i,j])
            end
            nrm = max(nrm,nrmi)
        end
    end
    return convert(Tnorm, nrm)
end

"""
    norm(A::AbstractArray, p::Real=2)

Compute the `p`-norm of a vector or the operator norm of a matrix `A`,
defaulting to the 2-norm.

    norm(A::AbstractVector, p::Real=2)

For vectors, this is equivalent to [`vecnorm`](@ref) and equal to:
```math
\\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p}
```
with ``a_i`` the entries of ``A`` and ``n`` its length.

`p` can assume any numeric value (even though not all values produce a
mathematically valid vector norm). In particular, `norm(A, Inf)` returns the largest value
in `abs(A)`, whereas `norm(A, -Inf)` returns the smallest.

# Examples
```jldoctest
julia> v = [3, -2, 6]
3-element Array{Int64,1}:
  3
 -2
  6

julia> norm(v)
7.0

julia> norm(v, Inf)
6.0
```
"""
norm(x::AbstractVector, p::Real=2) = vecnorm(x, p)

"""
    norm(A::AbstractMatrix, p::Real=2)

For matrices, the matrix norm induced by the vector `p`-norm is used, where valid values of
`p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, `p=2` is currently not
implemented.) Use [`vecnorm`](@ref) to compute the Frobenius norm.

When `p=1`, the matrix norm is the maximum absolute column sum of `A`:
```math
\\|A\\|_1 = \\max_{1 ≤ j ≤ n} \\sum_{i=1}^m | a_{ij} |
```
with ``a_{ij}`` the entries of ``A``, and ``m`` and ``n`` its dimensions.

When `p=2`, the matrix norm is the spectral norm, equal to the largest
singular value of `A`.

When `p=Inf`, the matrix norm is the maximum absolute row sum of `A`:
```math
\\|A\\|_\\infty = \\max_{1 ≤ i ≤ m} \\sum _{j=1}^n | a_{ij} |
```

# Examples
```jldoctest
julia> A = [1 -2 -3; 2 3 -1]
2×3 Array{Int64,2}:
 1  -2  -3
 2   3  -1

julia> norm(A, Inf)
6.0
```
"""
function norm(A::AbstractMatrix, p::Real=2)
    if p == 2
        return norm2(A)
    elseif p == 1
        return norm1(A)
    elseif p == Inf
        return normInf(A)
    else
        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
    end
end

"""
    norm(x::Number, p::Real=2)

For numbers, return ``\\left( |x|^p \\right)^{1/p}``.
This is equivalent to [`vecnorm`](@ref).
"""
@inline norm(x::Number, p::Real=2) = vecnorm(x, p)

"""
    norm(A::Adjoint{<:Any,<:AbstracVector}, q::Real=2)
    norm(A::Transpose{<:Any,<:AbstracVector}, q::Real=2)

For Adjoint/Transpose-wrapped vectors, return the ``q``-norm of `A`, which is
equivalent to the p-norm with value `p = q/(q-1)`. They coincide at `p = q = 2`.

The difference in norm between a vector space and its dual arises to preserve
the relationship between duality and the inner product, and the result is
consistent with the p-norm of `1 × n` matrix.

# Examples
```jldoctest
julia> v = [1; im];

julia> vc = v';

julia> norm(vc, 1)
1.0

julia> norm(v, 1)
2.0

julia> norm(vc, 2)
1.4142135623730951

julia> norm(v, 2)
1.4142135623730951

julia> norm(vc, Inf)
2.0

julia> norm(v, Inf)
1.0
```
"""
norm(v::TransposeAbsVec, q::Real) = q == Inf ? norm(v.parent, 1) : norm(v.parent, q/(q-1))
norm(v::AdjointAbsVec, q::Real) = q == Inf ? norm(conj(v.parent), 1) : norm(conj(v.parent), q/(q-1))
norm(v::AdjointAbsVec) = norm(conj(v.parent))
norm(v::TransposeAbsVec) = norm(v.parent)

function vecdot(x::AbstractArray, y::AbstractArray)
    lx = _length(x)
    if lx != _length(y)
        throw(DimensionMismatch("first array has length $(lx) which does not match the length of the second, $(_length(y))."))
    end
    if lx == 0
        return dot(zero(eltype(x)), zero(eltype(y)))
    end
    s = zero(dot(first(x), first(y)))
    for (Ix, Iy) in zip(eachindex(x), eachindex(y))
        @inbounds  s += dot(x[Ix], y[Iy])
    end
    s
end

"""
    vecdot(x, y)

For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or
any element type for which `dot` is defined), compute the Euclidean dot product (the sum of
`dot(x[i],y[i])`) as if they were vectors.

# Examples
```jldoctest
julia> vecdot(1:5, 2:6)
70

julia> x = fill(2., (5,5));

julia> y = fill(3., (5,5));

julia> vecdot(x, y)
150.0
```
"""
function vecdot(x, y) # arbitrary iterables
    ix = start(x)
    if done(x, ix)
        if !isempty(y)
            throw(DimensionMismatch("x and y are of different lengths!"))
        end
        return dot(zero(eltype(x)), zero(eltype(y)))
    end
    iy = start(y)
    if done(y, iy)
        throw(DimensionMismatch("x and y are of different lengths!"))
    end
    (vx, ix) = next(x, ix)
    (vy, iy) = next(y, iy)
    s = dot(vx, vy)
    while !done(x, ix)
        if done(y, iy)
            throw(DimensionMismatch("x and y are of different lengths!"))
        end
        (vx, ix) = next(x, ix)
        (vy, iy) = next(y, iy)
        s += dot(vx, vy)
    end
    if !done(y, iy)
            throw(DimensionMismatch("x and y are of different lengths!"))
    end
    return s
end

vecdot(x::Number, y::Number) = conj(x) * y

dot(x::Number, y::Number) = vecdot(x, y)

"""
    dot(x, y)
    ⋅(x,y)

Compute the dot product between two vectors. For complex vectors, the first vector is conjugated.
When the vectors have equal lengths, calling `dot` is semantically equivalent to `sum(vx'vy for (vx,vy) in zip(x, y))`.

# Examples
```jldoctest
julia> dot([1; 1], [2; 3])
5

julia> dot([im; im], [1; 1])
0 - 2im
```
"""
function dot(x::AbstractVector, y::AbstractVector)
    if length(x) != length(y)
        throw(DimensionMismatch("dot product arguments have lengths $(length(x)) and $(length(y))"))
    end
    ix = start(x)
    if done(x, ix)
        # we only need to check the first vector, since equal lengths have been asserted
        return zero(eltype(x))'zero(eltype(y))
    end
    @inbounds (vx, ix) = next(x, ix)
    @inbounds (vy, iy) = next(y, start(y))
    s = vx'vy
    while !done(x, ix)
        @inbounds (vx, ix) = next(x, ix)
        @inbounds (vy, iy) = next(y, iy)
        s += vx'vy
    end
    return s
end

# Call optimized BLAS methods for vectors of numbers
dot(x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = vecdot(x, y)

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

"""
    rank(A[, tol::Real])

Compute the rank of a matrix by counting how many singular
values of `A` have magnitude greater than `tol*σ₁` where `σ₁` is
`A`'s largest singular values. By default, the value of `tol` is the smallest
dimension of `A` multiplied by the [`eps`](@ref)
of the [`eltype`](@ref) of `A`.

# Examples
```jldoctest
julia> rank(Matrix(I, 3, 3))
3

julia> rank(diagm(0 => [1, 0, 2]))
2

julia> rank(diagm(0 => [1, 0.001, 2]), 0.1)
2

julia> rank(diagm(0 => [1, 0.001, 2]), 0.00001)
3
```
"""
function rank(A::AbstractMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A))))))
    s = svdvals(A)
    sum(x -> x > tol*s[1], s)
end
rank(x::Number) = x == 0 ? 0 : 1

"""
    trace(M)

Matrix trace. Sums the diagonal elements of `M`.

# Examples
```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> trace(A)
5
```
"""
function trace(A::AbstractMatrix)
    checksquare(A)
    sum(diag(A))
end
trace(x::Number) = x

#kron(a::AbstractVector, b::AbstractVector)
#kron(a::AbstractMatrix{T}, b::AbstractMatrix{S}) where {T,S}

#det(a::AbstractMatrix)

"""
    inv(M)

Matrix inverse. Computes matrix `N` such that
`M * N = I`, where `I` is the identity matrix.
Computed by solving the left-division
`N = M \\ I`.

# Examples
```jldoctest
julia> M = [2 5; 1 3]
2×2 Array{Int64,2}:
 2  5
 1  3

julia> N = inv(M)
2×2 Array{Float64,2}:
  3.0  -5.0
 -1.0   2.0

julia> M*N == N*M == Matrix(I, 2, 2)
true
```
"""
function inv(A::AbstractMatrix{T}) where T
    n = checksquare(A)
    S = typeof(zero(T)/one(T))      # dimensionful
    S0 = typeof(zero(T)/oneunit(T)) # dimensionless
    dest = Matrix{S0}(I, n, n)
    ldiv!(factorize(convert(AbstractMatrix{S}, A)), dest)
end

pinv(v::AbstractVector{T}, tol::Real = real(zero(T))) where {T<:Real} = _vectorpinv(transpose, v, tol)
pinv(v::AbstractVector{T}, tol::Real = real(zero(T))) where {T<:Complex} = _vectorpinv(adjoint, v, tol)
pinv(v::AbstractVector{T}, tol::Real = real(zero(T))) where {T} = _vectorpinv(adjoint, v, tol)
function _vectorpinv(dualfn::Tf, v::AbstractVector{Tv}, tol) where {Tv,Tf}
    res = dualfn(similar(v, typeof(zero(Tv) / (abs2(one(Tv)) + abs2(one(Tv))))))
    den = sum(abs2, v)
    # as tol is the threshold relative to the maximum singular value, for a vector with
    # single singular value σ=√den, σ ≦ tol*σ is equivalent to den=0 ∨ tol≥1
    if iszero(den) || tol >= one(tol)
        fill!(res, zero(eltype(res)))
    else
        res .= dualfn(v) ./ den
    end
    return res
end

# this method is just an optimization: literal negative powers of A are
# already turned by literal_pow into powers of inv(A), but for A^-1 this
# would turn into inv(A)^1 = copy(inv(A)), which makes an extra copy.
@inline Base.literal_pow(::typeof(^), A::AbstractMatrix, ::Val{-1}) = inv(A)

"""
    \\(A, B)

Matrix division using a polyalgorithm. For input matrices `A` and `B`, the result `X` is
such that `A*X == B` when `A` is square. The solver that is used depends upon the structure
of `A`.  If `A` is upper or lower triangular (or diagonal), no factorization of `A` is
required and the system is solved with either forward or backward substitution.
For non-triangular square matrices, an LU factorization is used.

For rectangular `A` the result is the minimum-norm least squares solution computed by a
pivoted QR factorization of `A` and a rank estimate of `A` based on the R factor.

When `A` is sparse, a similar polyalgorithm is used. For indefinite matrices, the `LDLt`
factorization does not use pivoting during the numerical factorization and therefore the
procedure can fail even for invertible matrices.

# Examples
```jldoctest
julia> A = [1 0; 1 -2]; B = [32; -4];

julia> X = A \\ B
2-element Array{Float64,1}:
 32.0
 18.0

julia> A * X == B
true
```
"""
function (\)(A::AbstractMatrix, B::AbstractVecOrMat)
    m, n = size(A)
    if m == n
        if istril(A)
            if istriu(A)
                return Diagonal(A) \ B
            else
                return LowerTriangular(A) \ B
            end
        end
        if istriu(A)
            return UpperTriangular(A) \ B
        end
        return lufact(A) \ B
    end
    return qrfact(A,Val(true)) \ B
end

(\)(a::AbstractVector, b::AbstractArray) = pinv(a) * b
(/)(A::AbstractVecOrMat, B::AbstractVecOrMat) = copy(adjoint(adjoint(B) \ adjoint(A)))
# \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough
# /(x::Number,A::StridedMatrix) = x*inv(A)
/(x::Number, v::AbstractVector) = x*pinv(v)

cond(x::Number) = x == 0 ? Inf : 1.0
cond(x::Number, p) = cond(x)

#Skeel condition numbers
condskeel(A::AbstractMatrix, p::Real=Inf) = norm(abs.(inv(A))*abs.(A), p)

"""
    condskeel(M, [x, p::Real=Inf])

```math
\\kappa_S(M, p) = \\left\\Vert \\left\\vert M \\right\\vert \\left\\vert M^{-1} \\right\\vert \\right\\Vert_p \\\\
\\kappa_S(M, x, p) = \\left\\Vert \\left\\vert M \\right\\vert \\left\\vert M^{-1} \\right\\vert \\left\\vert x \\right\\vert \\right\\Vert_p
```

Skeel condition number ``\\kappa_S`` of the matrix `M`, optionally with respect to the
vector `x`, as computed using the operator `p`-norm. ``\\left\\vert M \\right\\vert``
denotes the matrix of (entry wise) absolute values of ``M``;
``\\left\\vert M \\right\\vert_{ij} = \\left\\vert M_{ij} \\right\\vert``.
Valid values for `p` are `1`, `2` and `Inf` (default).

This quantity is also known in the literature as the Bauer condition number, relative
condition number, or componentwise relative condition number.
"""
condskeel(A::AbstractMatrix, x::AbstractVector, p::Real=Inf) = norm(abs.(inv(A))*(abs.(A)*abs.(x)), p)

issymmetric(A::AbstractMatrix{<:Real}) = ishermitian(A)

"""
    issymmetric(A) -> Bool

Test whether a matrix is symmetric.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> issymmetric(a)
true

julia> b = [1 im; -im 1]
2×2 Array{Complex{Int64},2}:
 1+0im  0+1im
 0-1im  1+0im

julia> issymmetric(b)
false
```
"""
function issymmetric(A::AbstractMatrix)
    indsm, indsn = axes(A)
    if indsm != indsn
        return false
    end
    for i = first(indsn):last(indsn), j = (i):last(indsn)
        if A[i,j] != transpose(A[j,i])
            return false
        end
    end
    return true
end

issymmetric(x::Number) = x == x

"""
    ishermitian(A) -> Bool

Test whether a matrix is Hermitian.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> ishermitian(a)
true

julia> b = [1 im; -im 1]
2×2 Array{Complex{Int64},2}:
 1+0im  0+1im
 0-1im  1+0im

julia> ishermitian(b)
true
```
"""
function ishermitian(A::AbstractMatrix)
    indsm, indsn = axes(A)
    if indsm != indsn
        return false
    end
    for i = indsn, j = i:last(indsn)
        if A[i,j] != adjoint(A[j,i])
            return false
        end
    end
    return true
end

ishermitian(x::Number) = (x == conj(x))

"""
    istriu(A::AbstractMatrix, k::Integer = 0) -> Bool

Test whether `A` is upper triangular starting from the `k`th superdiagonal.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> istriu(a)
false

julia> istriu(a, -1)
true

julia> b = [1 im; 0 -1]
2×2 Array{Complex{Int64},2}:
 1+0im   0+1im
 0+0im  -1+0im

julia> istriu(b)
true

julia> istriu(b, 1)
false
```
"""
function istriu(A::AbstractMatrix, k::Integer = 0)
    m, n = size(A)
    for j in 1:min(n, m + k - 1)
        for i in max(1, j - k + 1):m
            iszero(A[i, j]) || return false
        end
    end
    return true
end
istriu(x::Number) = true

"""
    istril(A::AbstractMatrix, k::Integer = 0) -> Bool

Test whether `A` is lower triangular starting from the `k`th superdiagonal.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> istril(a)
false

julia> istril(a, 1)
true

julia> b = [1 0; -im -1]
2×2 Array{Complex{Int64},2}:
 1+0im   0+0im
 0-1im  -1+0im

julia> istril(b)
true

julia> istril(b, -1)
false
```
"""
function istril(A::AbstractMatrix, k::Integer = 0)
    m, n = size(A)
    for j in max(1, k + 2):n
        for i in 1:min(j - k - 1, m)
            iszero(A[i, j]) || return false
        end
    end
    return true
end
istril(x::Number) = true

"""
    isbanded(A::AbstractMatrix, kl::Integer, ku::Integer) -> Bool

Test whether `A` is banded with lower bandwidth starting from the `kl`th superdiagonal
and upper bandwidth extending through the `ku`th superdiagonal.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> isbanded(a, 0, 0)
false

julia> isbanded(a, -1, 1)
true

julia> b = [1 0; -im -1] # lower bidiagonal
2×2 Array{Complex{Int64},2}:
 1+0im   0+0im
 0-1im  -1+0im

julia> isbanded(b, 0, 0)
false

julia> isbanded(b, -1, 0)
true
```
"""
isbanded(A::AbstractMatrix, kl::Integer, ku::Integer) = istriu(A, kl) && istril(A, ku)

"""
    isdiag(A) -> Bool

Test whether a matrix is diagonal.

# Examples
```jldoctest
julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> isdiag(a)
false

julia> b = [im 0; 0 -im]
2×2 Array{Complex{Int64},2}:
 0+1im  0+0im
 0+0im  0-1im

julia> isdiag(b)
true
```
"""
isdiag(A::AbstractMatrix) = isbanded(A, 0, 0)
isdiag(x::Number) = true

"""
    linreg(x, y)

Perform simple linear regression using Ordinary Least Squares. Returns `a` and `b` such
that `a + b*x` is the closest straight line to the given points `(x, y)`, i.e., such that
the squared error between `y` and `a + b*x` is minimized.

# Examples
```julia
using PyPlot
x = 1.0:12.0
y = [5.5, 6.3, 7.6, 8.8, 10.9, 11.79, 13.48, 15.02, 17.77, 20.81, 22.0, 22.99]
a, b = linreg(x, y)          # Linear regression
plot(x, y, "o")              # Plot (x, y) points
plot(x, a + b*x)             # Plot line determined by linear regression
```

See also:

`\\`, [`cov`](@ref), [`std`](@ref), [`mean`](@ref).

"""
function linreg(x::AbstractVector, y::AbstractVector)
    # Least squares given
    # Y = a + b*X
    # where
    # b = cov(X, Y)/var(X)
    # a = mean(Y) - b*mean(X)
    if size(x) != size(y)
        throw(DimensionMismatch("x has size $(size(x)) and y has size $(size(y)), " *
            "but these must be the same size"))
    end
    mx = mean(x)
    my = mean(y)
    # don't need to worry about the scaling (n vs n - 1) since they cancel in the ratio
    b = Base.covm(x, mx, y, my)/Base.varm(x, mx)
    a = my - b*mx
    return (a, b)
end

"""
    peakflops(n::Integer=2000; parallel::Bool=false)

`peakflops` computes the peak flop rate of the computer by using double precision
[`gemm!`](@ref LinearAlgebra.BLAS.gemm!). By default, if no arguments are specified, it
multiplies a matrix of size `n x n`, where `n = 2000`. If the underlying BLAS is using
multiple threads, higher flop rates are realized. The number of BLAS threads can be set with
[`BLAS.set_num_threads(n)`](@ref).

If the keyword argument `parallel` is set to `true`, `peakflops` is run in parallel on all
the worker processors. The flop rate of the entire parallel computer is returned. When
running in parallel, only 1 BLAS thread is used. The argument `n` still refers to the size
of the problem that is solved on each processor.
"""
function peakflops(n::Integer=2000; parallel::Bool=false)
    a = fill(1.,100,100)
    t = @elapsed a2 = a*a
    a = fill(1.,n,n)
    t = @elapsed a2 = a*a
    @assert a2[1,1] == n
    parallel ? sum(pmap(peakflops, [ n for i in 1:nworkers()])) : (2*Float64(n)^3/t)
end

# BLAS-like in-place y = x*α+y function (see also the version in blas.jl
#                                          for BlasFloat Arrays)
function axpy!(α, x::AbstractArray, y::AbstractArray)
    n = _length(x)
    if n != _length(y)
        throw(DimensionMismatch("x has length $n, but y has length $(_length(y))"))
    end
    for (IY, IX) in zip(eachindex(y), eachindex(x))
        @inbounds y[IY] += x[IX]*α
    end
    y
end

function axpy!(α, x::AbstractArray, rx::AbstractArray{<:Integer}, y::AbstractArray, ry::AbstractArray{<:Integer})
    if _length(rx) != _length(ry)
        throw(DimensionMismatch("rx has length $(_length(rx)), but ry has length $(_length(ry))"))
    elseif !checkindex(Bool, linearindices(x), rx)
        throw(BoundsError(x, rx))
    elseif !checkindex(Bool, linearindices(y), ry)
        throw(BoundsError(y, ry))
    end
    for (IY, IX) in zip(eachindex(ry), eachindex(rx))
        @inbounds y[ry[IY]] += x[rx[IX]]*α
    end
    y
end

function axpby!(α, x::AbstractArray, β, y::AbstractArray)
    if _length(x) != _length(y)
        throw(DimensionMismatch("x has length $(_length(x)), but y has length $(_length(y))"))
    end
    for (IX, IY) in zip(eachindex(x), eachindex(y))
        @inbounds y[IY] = x[IX]*α + y[IY]*β
    end
    y
end

# Elementary reflection similar to LAPACK. The reflector is not Hermitian but
# ensures that tridiagonalization of Hermitian matrices become real. See lawn72
@inline function reflector!(x::AbstractVector)
    n = length(x)
    @inbounds begin
        ξ1 = x[1]
        normu = abs2(ξ1)
        for i = 2:n
            normu += abs2(x[i])
        end
        if normu == zero(normu)
            return zero(ξ1/normu)
        end
        normu = sqrt(normu)
        ν = copysign(normu, real(ξ1))
        ξ1 += ν
        x[1] = -ν
        for i = 2:n
            x[i] /= ξ1
        end
    end
    ξ1/ν
end

# apply reflector from left
@inline function reflectorApply!(x::AbstractVector, τ::Number, A::StridedMatrix)
    m, n = size(A)
    if length(x) != m
        throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m"))
    end
    @inbounds begin
        for j = 1:n
            # dot
            vAj = A[1, j]
            for i = 2:m
                vAj += x[i]'*A[i, j]
            end

            vAj = conj(τ)*vAj

            # ger
            A[1, j] -= vAj
            for i = 2:m
                A[i, j] -= x[i]*vAj
            end
        end
    end
    return A
end

"""
    det(M)

Matrix determinant.

# Examples
```jldoctest
julia> M = [1 0; 2 2]
2×2 Array{Int64,2}:
 1  0
 2  2

julia> det(M)
2.0
```
"""
function det(A::AbstractMatrix{T}) where T
    if istriu(A) || istril(A)
        S = typeof((one(T)*zero(T) + zero(T))/one(T))
        return convert(S, det(UpperTriangular(A)))
    end
    return det(lufact(A))
end
det(x::Number) = x

"""
    logabsdet(M)

Log of absolute value of matrix determinant. Equivalent to
`(log(abs(det(M))), sign(det(M)))`, but may provide increased accuracy and/or speed.

# Examples
```jldoctest
julia> A = [-1. 0.; 0. 1.]
2×2 Array{Float64,2}:
 -1.0  0.0
  0.0  1.0

julia> det(A)
-1.0

julia> logabsdet(A)
(0.0, -1.0)

julia> B = [2. 0.; 0. 1.]
2×2 Array{Float64,2}:
 2.0  0.0
 0.0  1.0

julia> det(B)
2.0

julia> logabsdet(B)
(0.6931471805599453, 1.0)
```
"""
logabsdet(A::AbstractMatrix) = logabsdet(lufact(A))

"""
    logdet(M)

Log of matrix determinant. Equivalent to `log(det(M))`, but may provide
increased accuracy and/or speed.

# Examples
```jldoctest
julia> M = [1 0; 2 2]
2×2 Array{Int64,2}:
 1  0
 2  2

julia> logdet(M)
0.6931471805599453

julia> logdet(Matrix(I, 3, 3))
0.0
```
"""
function logdet(A::AbstractMatrix)
    d,s = logabsdet(A)
    return d + log(s)
end

logdet(A) = log(det(A))

const NumberArray{T<:Number} = AbstractArray{T}

"""
    promote_leaf_eltypes(itr)

For an (possibly nested) iterable object `itr`, promote the types of leaf
elements.  Equivalent to `promote_type(typeof(leaf1), typeof(leaf2), ...)`.
Currently supports only numeric leaf elements.

# Examples
```jldoctest
julia> a = [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]
3-element Array{Any,1}:
  Any[1,2,[3,4]]
 5.0
  Any[0+6im,[7.0,8.0]]

julia> promote_leaf_eltypes(a)
Complex{Float64}
```
"""
promote_leaf_eltypes(x::Union{AbstractArray{T},Tuple{Vararg{T}}}) where {T<:Number} = T
promote_leaf_eltypes(x::Union{AbstractArray{T},Tuple{Vararg{T}}}) where {T<:NumberArray} = eltype(T)
promote_leaf_eltypes(x::T) where {T} = T
promote_leaf_eltypes(x::Union{AbstractArray,Tuple}) = mapreduce(promote_leaf_eltypes, promote_type, Bool, x)

# isapprox: approximate equality of arrays [like isapprox(Number,Number)]
# Supports nested arrays; e.g., for `a = [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]`
# `a ≈ a` is `true`.
function isapprox(x::AbstractArray, y::AbstractArray;
    atol::Real=0,
    rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y),atol),
    nans::Bool=false, norm::Function=vecnorm)
    d = norm(x - y)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(x), norm(y)))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
    end
end

"""
    normalize!(v::AbstractVector, p::Real=2)

Normalize the vector `v` in-place so that its `p`-norm equals unity,
i.e. `norm(v, p) == 1`.
See also [`normalize`](@ref) and [`vecnorm`](@ref).
"""
function normalize!(v::AbstractVector, p::Real=2)
    nrm = norm(v, p)
    __normalize!(v, nrm)
end

@inline function __normalize!(v::AbstractVector, nrm::AbstractFloat)
    # The largest positive floating point number whose inverse is less than infinity
    δ = inv(prevfloat(typemax(nrm)))

    if nrm ≥ δ # Safe to multiply with inverse
        invnrm = inv(nrm)
        rmul!(v, invnrm)

    else # scale elements to avoid overflow
        εδ = eps(one(nrm))/δ
        rmul!(v, εδ)
        rmul!(v, inv(nrm*εδ))
    end

    v
end

"""
    normalize(v::AbstractVector, p::Real=2)

Normalize the vector `v` so that its `p`-norm equals unity,
i.e. `norm(v, p) == vecnorm(v, p) == 1`.
See also [`normalize!`](@ref) and [`vecnorm`](@ref).

# Examples
```jldoctest
julia> a = [1,2,4];

julia> b = normalize(a)
3-element Array{Float64,1}:
 0.2182178902359924
 0.4364357804719848
 0.8728715609439696

julia> norm(b)
1.0

julia> c = normalize(a, 1)
3-element Array{Float64,1}:
 0.14285714285714285
 0.2857142857142857
 0.5714285714285714

julia> norm(c, 1)
1.0
```
"""
function normalize(v::AbstractVector, p::Real = 2)
    nrm = norm(v, p)
    if !isempty(v)
        vv = copy_oftype(v, typeof(v[1]/nrm))
        return __normalize!(vv, nrm)
    else
        T = typeof(zero(eltype(v))/nrm)
        return T[]
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
mutable struct PGEpisode{S<:Real, A<:Real, R<:Real}
    n::Int64
    s::Matrix{S}
    a::Matrix{A}
    r::Vector{R}
    v::Vector{R}
    γ::R
    nstep::Int64
    λ::R
    function PGEpisode{S, A, R}(obssize, actiondim, maxsize;
            γ=0.99, nstep=15, λ=0.97) where {S<:Real, A<:Real, R<:Real}
        new(
            0,
            zeros(S, obssize, maxsize),
            zeros(A, actiondim, maxsize),
            zeros(R, maxsize),
            zeros(R, maxsize),
            γ,
            nstep,
            λ
        )
    end
end

function append!(memory::PGEpisode, s, a, r)
    memory.n += 1
    memory.s[:, memory.n] .= s
    memory.a[:, memory.n] .= a
    memory.r[memory.n] = r
end

function append!(memory::PGEpisode, s, a, r, v)
    append!(memory::PGEpisode, s, a, r)
    memory.v[memory.n] = v
end

discount!(memory::PGEpisode) = discount!(view(memory.r, 1:memory.n), memory.γ)

function nstepr!(memory::PGEpisode)
    nstepr!(view(memory.r, 1:memory.n), view(memory.v, 1:memory.n), memory.γ, memory.nstep)
end

function GAEλ!(memory::PGEpisode, norm=true)
    GAEλ!(view(memory.v, 1:memory.n), view(memory.r, 1:memory.n), memory.γ, memory.λ)
    norm && norm!(view(memory.v, 1:memory.n))
end

function adv!(memory::PGEpisode, norm=false)
    memory.v .= memory.r .- memory.v
    norm && norm!(view(memory.v, 1:memory.n))
end

function batch(memory::PGEpisode, batchsize=0)
    batchsize = (batchsize == 0 || batchsize > memory.n) ? memory.n : batchsize
    inds = sample(1:memory.n, batchsize, replace=false)
    memory.n = 0
    return (
        memory.s[:, inds],
        memory.a[:, inds],
        memory.r[inds],
        memory.v[inds]
    )
end


text:
## a list of reserved names from PyCall.jl
const reserved = Set{String}([
    "while", "if", "for", "try", "return", "break", "continue",
    "function", "macro", "quote", "let", "local", "global", "const",
    "abstract", "typealias", "type", "bitstype", "immutable", "ccall",
    "do", "module", "baremodule", "using", "import", "export", "importall",
    "false", "true", "Tuple", "rmember", "__package__"])

cached_namespaces = Dict{String, Module}()

"""
Import an R package as a julia module.
Sanitizes the imported symbols by default, because otherwise certain symbols cannot be used.

E.g.: `PerformanceAnalytics::charts.Bar` in R becomes `PerformanceAnalytics.charts_Bar` in Julia
```
gg = rimport("ggplot2")
```
"""
function rimport(pkg::String, s::Symbol=:__anonymous__; normalizenames::Bool=true)
    if pkg in keys(cached_namespaces)
        m = cached_namespaces[pkg]
    else
        ns = rcall(:asNamespace, pkg)
        members = rcopy(Vector{String}, rcall(:getNamespaceExports, ns))

        m = Module(s, false)
        id = Expr(:const, Expr(:(=), :__package__, pkg))
        if normalizenames
            exports = [Symbol(replace(x, '.' => '_')) for x in members]
        else
            exports = [Symbol(x) for x in members]
        end

        filtered_indices = filter!(i -> !(string(exports[i]) in reserved),
                                   collect(eachindex(exports)))
        consts = [Expr(:const, Expr(:(=),
                       exports[i],
                       rcall(Symbol("::"), pkg, members[i]))) for i in filtered_indices ]
        Core.eval(m, Expr(:toplevel, id, consts..., Expr(:export, exports...), :(rmember(x) = ($getindex)($ns, x))))
        cached_namespaces[pkg] = m
    end
    m
end
rimport(pkg::Symbol, s::Symbol=:__anonymous__) = rimport(string(pkg), s)

"""
Import an R Package as a Julia module. For example,
```
@rimport ggplot2
```
is equivalent to `ggplot2 = rimport("ggplot2")` with error checking.

You can also use classic Python syntax to make an alias: `@rimport *package-name* as *shorthand*`
```
@rimport ggplot2 as gg
```
which is equivalent to `gg = rimport("ggplot2")`.
"""
macro rimport(x, args...)
    if length(args)==2 && args[1] == :as
        m = Symbol(args[2])
    elseif length(args)==0
        m = Symbol(x)
    else
        throw(ArgumentError("invalid import syntax."))
    end
    pkg = string(x)
    quote
        if !isdefined(@__MODULE__, $(QuoteNode(m)))
            const $(esc(m)) = rimport($pkg)
            nothing
        elseif typeof($(esc(m))) <: Module &&
                    :__package__ in names($(esc(m)), all = true) &&
                    $(esc(m)).__package__ == $pkg
            nothing
        else
            error($pkg * " already exists! Use the syntax `@rimport " * $pkg *" as *shorthand*`.")
            nothing
        end
    end
end

"""
Load all exported functions/objects of an R package to the current module. Almost equivalent to
```
__temp__ = rimport("ggplot2")
using .__temp__
```
"""
macro rlibrary(x)
    m = gensym("RCall")
    quote
        $(esc(m)) = rimport($(QuoteNode(x)))
        Core.eval(@__MODULE__, Expr(:using, Expr(:., :., $(QuoteNode(m)))))
    end
end


text:
struct TrajectoryOptimizationProblem{T}
    obj::Objective{T}
    model::DynamicsModel{T}
    con::ConstraintSet{T}
    x::Vector{Vector{T}}
    u::Vector{Vector{T}}
    w::Vector{Vector{T}}
    λ_dyn::Vector{Vector{T}}
    λ_stage::Vector{Vector{T}}
end

function TrajectoryOptimizationProblem(obj::Objective, model::DynamicsModel, con::ConstraintSet;
    w = [zeros(nw) for nw in model.dim.w])

    x = [zeros(nx) for nx in model.dim.x]
    u = [[zeros(nu) for nu in model.dim.u]..., zeros(0)]

    λ_dyn = [zeros(nx) for nx in model.dim.x[2:end]]
    λ_stage = [zeros(stage.nc) for stage in con.stage for t in stage.idx]

    TrajectoryOptimizationProblem(obj, model, con, x, u, w, λ_dyn, λ_stage)
end

TrajectoryOptimizationProblem(obj::Objective, model::DynamicsModel) = TrajectoryOptimizationProblem(obj, model, ConstraintSet())

struct TrajectoryOptimizationIndices
    obj_hess::Vector{Vector{Int}}
    dyn_con::Vector{Vector{Int}}
    dyn_jac::Vector{Vector{Int}}
    dyn_hess::Vector{Vector{Int}}
    stage_con::Vector{Vector{Int}}
    stage_jac::Vector{Vector{Int}}
    stage_hess::Vector{Vector{Int}}
end

function indices(obj::Objective, dyn::Vector{Dynamics{T}}, stage::StageConstraints{T}, key::Vector{Tuple{Int,Int}},
        nx::Vector{Int}, nu::Vector{Int}) where T
    # Jacobians
    dyn_con = constraint_indices(dyn, shift=0)
    dyn_jac = jacobian_indices(dyn, shift=0)
    stage_con = constraint_indices(stage, shift=num_con(dyn))
    stage_jac = jacobian_indices(stage, shift=num_jac(dyn))

    # Hessian of Lagrangian
    obj_hess = hessian_indices(obj, key, nx, nu)
    dyn_hess = hessian_indices(dyn, key, nx, nu)
    stage_hess = hessian_indices(stage, key, nx, nu)

    return TrajectoryOptimizationIndices(obj_hess, dyn_con, dyn_jac, dyn_hess, stage_con, stage_jac, stage_hess)
end

struct Problem{T} <: MOI.AbstractNLPEvaluator
    trajopt::TrajectoryOptimizationProblem{T}
    num_var::Int
    num_con::Int
    num_jac::Int
    num_hess_lag::Int
    var_bnds::Vector{Vector{T}}
    con_bnds::Vector{Vector{T}}
    idx::TrajectoryOptimizationIndices
    sp_jac
    sp_hess_lag
    hess_lag::Bool
end

function primal_bounds(bnds::Bounds{T}, nz::Int, x_idx::Vector{Vector{Int}}, u_idx::Vector{Vector{Int}}) where T
    zl, zu = -Inf * ones(nz), Inf * ones(nz)
    for bnd in bnds
        for t in bnd.idx
            length(bnd.xl) > 0 && (zl[x_idx[t]] = bnd.xl)
            length(bnd.xu) > 0 && (zu[x_idx[t]] = bnd.xu)
            length(bnd.ul) > 0 && (zl[u_idx[t]] = bnd.ul)
            length(bnd.uu) > 0 && (zu[u_idx[t]] = bnd.uu)
        end
    end
    return zl, zu
end

function constraint_bounds(dyn::Vector{Dynamics{T}}, stage::StageConstraints{T}, nc::Int, idx::TrajectoryOptimizationIndices) where T
    cl, cu = zeros(nc), zeros(nc)
    i = 1
    for con in stage
        for t in con.idx
            if con.type == :inequality
                cu[idx.stage_con[i]] .= Inf
            end
            i += 1
        end
    end
    return cl, cu
end

function Problem(trajopt::TrajectoryOptimizationProblem; eval_hess=false)
    # number of variables
    nz = sum(trajopt.model.dim.x) + sum(trajopt.model.dim.u)

    # number of constraints
    nc_dyn = num_con(trajopt.model.dyn)
    nc_con = num_con(trajopt.con.stage)
    nc = nc_dyn + nc_con

    # number of nonzeros in constraint Jacobian
    nj_dyn = num_jac(trajopt.model.dyn)
    nj_con = num_jac(trajopt.con.stage)
    nj = nj_dyn + nj_con

    # number of nonzeros in Hessian of Lagrangian
    nh = 0

    # primal variable bounds
    zl, zu = primal_bounds(trajopt.con.bounds, nz, trajopt.model.idx.x, trajopt.model.idx.u)

    # constraint Jacobian sparsity
    sp_dyn = sparsity_jacobian(trajopt.model.dyn, trajopt.model.dim.x, trajopt.model.dim.u, row_shift=0)
    sp_con = sparsity_jacobian(trajopt.con.stage, trajopt.model.dim.x, trajopt.model.dim.u, row_shift=nc_dyn)
    sp_jac = collect([sp_dyn..., sp_con...])

    # Hessian of Lagrangian sparsity
    sp_obj_hess = sparsity_hessian(trajopt.obj, trajopt.model.dim.x, trajopt.model.dim.u)
    sp_dyn_hess = sparsity_hessian(trajopt.model.dyn, trajopt.model.dim.x, trajopt.model.dim.u)
    sp_con_hess = sparsity_hessian(trajopt.con.stage, trajopt.model.dim.x, trajopt.model.dim.u)
    sp_hess_lag = [sp_obj_hess..., sp_dyn_hess..., sp_con_hess...]
    sp_hess_lag = !isempty(sp_hess_lag) ? sp_hess_lag : Tuple{Int,Int}[]
    sp_hess_key = sort(unique(sp_hess_lag))
    sp_hess_status = [!isempty(sp_obj_hess), !isempty(sp_dyn_hess), !isempty(sp_con_hess)]

    # indices
    idx = indices(trajopt.obj, trajopt.model.dyn, trajopt.con.stage, sp_hess_key, trajopt.model.dim.x, trajopt.model.dim.u)

    # nonlinear constraint bounds
    cl, cu = constraint_bounds(trajopt.model.dyn, trajopt.con.stage, nc, idx)

    Problem(trajopt, nz, nc, nj, nh, [zl, zu], [cl, cu], idx, sp_jac, sp_hess_key, eval_hess)
end

function trajectory!(x::Vector{Vector{T}}, u::Vector{Vector{T}}, z::Vector{T},
    idx_x::Vector{Vector{Int}}, idx_u::Vector{Vector{Int}}) where T
    for (t, idx) in enumerate(idx_x)
        x[t] .= @views z[idx]
    end
    for (t, idx) in enumerate(idx_u)
        u[t] .= @views z[idx]
    end
end

function duals!(λ_dyn::Vector{Vector{T}}, λ_stage::Vector{Vector{T}}, λ,
    idx_dyn::Vector{Vector{Int}}, idx_stage::Vector{Vector{Int}}) where T
    for (t, idx) in enumerate(idx_dyn)
        λ_dyn[t] .= @views λ[idx]
    end
    for (t, idx) in enumerate(idx_stage)
        λ_stage[t] .= @views λ[idx]
    end
end


text:
function levendist1(s::AbstractString, t::AbstractString)
    ls, lt = length(s), length(t)
    if ls > lt
        s, t = t, s
        ls, lt = lt, ls
    end
    dist = collect(0:ls)
    for (ind2, chr2) in enumerate(t)
        newdist = Vector{Int}(ls+1)
        newdist[1] = ind2
        for (ind1, chr1) in enumerate(s)
            if chr1 == chr2
                newdist[ind1+1] = dist[ind1]
            else
                newdist[ind1+1] = 1 + min(dist[ind1], dist[ind1+1], newdist[end])
            end
        end
        dist = newdist
    end
    return dist[end]
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "dataset_loading" begin
    @test load_station_metadata() isa DataFrame
    @test load_countries_metadata() isa DataFrame
    @test load_inventories() isa DataFrame
end


text:
using AssemblyCalculus: Neuron, empty_synapses, num_synapses
using AssemblyCalculus: NeuralArea, num_neurons, random_firing!, winners!, change_in_assembly
using AssemblyCalculus: Stimulus, rand_stim, zero_stim
using AssemblyCalculus: getneuron, globalindex, select, neuron_attrib_graph, freeze_synapses!
using AssemblyCalculus: MaxIters, Converge
using Graphs
using StatsBase

@testset "Neuron" begin
    findparam(n::Neuron{T, U}) where {T, U} = T
    @testset "Constructors" begin
        neuronf64 = Neuron(Float64)
        neuronf32 = Neuron(Float32)
        @test findparam(neuronf64) == Float64
        @test findparam(neuronf32) == Float32
        neuronf64 = Neuron(NeuralArea(Float64), 0, empty_synapses(Float64))
        T = Float32
        neuronf32 = Neuron(NeuralArea(T), 0, empty_synapses(T))
        @test findparam(neuronf64) == Float64
        @test findparam(neuronf32) == Float32
    end
    @testset "Utilities" begin
        T = Float32
        area1 = NeuralArea(T)
        neuron1 = Neuron(area1, 1, empty_synapses(T))
        push!(area1.neurons, neuron1)
        area2 = NeuralArea(T)
        neuron2 = Neuron(area2, 1, empty_synapses(T))
        push!(area2.neurons, neuron2)
        neuron1[area2, neuron2] = 1.5f0
        neuron2[area1, neuron1] = 1.7f0
        @test neuron1[area2][neuron2] == 1.5f0
        @test neuron2[area1][neuron1] == 1.7f0
        @test get(neuron1, area1, true)
        @test length(get!(neuron1, area2, [])) == 1
        @test length(get!(neuron1, area1, Dict{Neuron{T, NeuralArea{T}}, T}())) == 0
        @test length(collect(keys(neuron1))) == 2
        @test num_synapses(neuron1, area1) == 0
        @test num_synapses(neuron1, area2) == 1
        @test num_synapses(neuron2, area1) == 1
        @test num_synapses(neuron2, area2) == 0
        @test num_synapses(neuron1) == 1
        @test num_synapses(Neuron(Float64)) == 0
    end
end

@testset "NeuralArea" begin
    findparam(na::NeuralArea{T}) where T = T
    @testset "Constructors" begin
        na = NeuralArea(Float64)
        @test findparam(na) == Float64
        @test findparam(NeuralArea(Float32)) == Float32
        @test findparam(NeuralArea(100, 10, 1.0)) == Float64
        @test findparam(NeuralArea(100, 10, 1.0f0)) == Float32
        na = NeuralArea(Float64)
        @test num_neurons(na) == 0
        push!(na.neurons, Neuron(Float64))
        @test num_neurons(na) == 1
        @test_throws UndefVarError push!(na.neurons, Neurons(Float32))
    end
    @testset "Utilities" begin
        T = Float32
        area1 = NeuralArea(T)
        neuron1 = Neuron(area1, 1, empty_synapses(T))
        push!(area1.neurons, neuron1)
        area2 = NeuralArea(T)
        neuron2 = Neuron(area2, 1, empty_synapses(T))
        push!(area2.neurons, neuron2)
        neuron1[area2, neuron2] = 1.5f0
        neuron2[area1, neuron1] = 1.7f0
        @test num_synapses(area1) == 1
        @test num_synapses(area2) == 1
        for n in area1
            @test num_synapses(n) == 1
        end
        @test length(area1.firing) == 0
        area1.assembly_size = 1
        random_firing!(area1)
        @test length(area1.firing) == 1
        area1.assembly_size = 2
        @test_throws ErrorException random_firing!(area1)
        # k-winners
        area1.assembly_size = 3
        @test_throws ErrorException winners!(area1, rand(Float32, 3))
        for i in 2:9
            push!(area1.neurons, Neuron(area1, i, empty_synapses(Float32)))
        end
        currents = zeros(Float32, 10)
        @test_throws ErrorException winners!(area1, currents)
        push!(area1.neurons, Neuron(area1, 10, empty_synapses(Float32)))
        winners!(area1, currents)
        @test length(area1.firing) == 0
        currents[2:4] .= 1.0f0
        winners!(area1, currents)
        @test all(area1.firing .== [2, 3, 4])
        currents = rand(Float32, 10)
        perm = sortperm(currents, rev=true)
        winners!(area1, currents)
        @test all(area1.firing .== perm[1:3])
        # Prev assembly comparison
        area1.firing = [1, 2, 3]
        area1.firing_prev = [1, 2, 3]
        @test change_in_assembly(area1) == 0
        area1.firing = [1, 2, 4]
        @test change_in_assembly(area1) == 1
    end
end

@testset "Stimulus Type" begin
    getparam(ic::Stimulus{T}) where T = T
    @testset "Constructors" begin
        @test getparam(Stimulus()) == Float64
        @test getparam(Stimulus(Float32)) == Float32
        area1 = NeuralArea(4000, 20, 0.01)
        # One edge per neuron
        for neuron in area1
            neuron[area1, neuron] = 1.0
        end
        ic = rand_stim(area1)
        @test ic.area == area1
        @test length(ic.currents) == length(area1)
        @test 0.001 < mean(ic.currents) < 0.01
        ic = rand_stim(area1, p=1. / 400)
        @test ic.area == area1
        @test length(ic.currents) == length(area1)
        @test 0.01 < mean(ic.currents) < 0.1
        @test sum(zero_stim(area1).currents) == 0
        @test_throws ErrorException rand_stim(area1, p=1.0)
    end
end

@testset "Assembly" begin
    getparam(a::Assembly{T}) where T = T
    @testset "Constructors" begin
        @test getparam(Assembly()) == Float64
        @test getparam(Assembly(Float32)) == Float32
        area1 = NeuralArea(100, 10, 0.01)
        ic = zero_stim(area1)
        @test_throws ErrorException Assembly(area1, [ic], Assembly{Float64}[])
        random_firing!(area1)
        a = Assembly(area1, [ic], Assembly{Float64}[])
        @test a.area == area1
        @test length(a.parent_stims) == 1
        @test length(a.parent_assemblies) == 0
        @test length(a.neurons) == area1.assembly_size
    end
end

@testset "Brain" begin
    g = erdos_renyi(10, 1.0, is_directed=true)
    sizes = [3, 3, 4]
    assemblies = [1, 1, 2]
    plasticities = [.01f0, .01f0, .01f0]
    ba = Brain(g, sizes, assemblies, plasticities)
    @testset "Constructors" begin
        @test length(ba) == 3
        @test typeof(ba[2]) == NeuralArea{Float32}
        @test num_synapses(ba) == 90
        @test num_neurons(ba) == 10
        adj = Float32.(rand(10, 10) .< 0.2)
        ba = Brain(adj, sizes, assemblies, plasticities)
        @test all(adj .== adjacency_matrix(ba))
        freeze_synapses!(ba)
        @test all([a.plasticity for a in ba.areas] .== 0)
    end
    @testset "Type Utilities" begin
        g_idx = 5
        @test globalindex(ba, getneuron(ba, g_idx)) == g_idx
        @test select(ba, 1) == ba.areas[1]
        @test_throws ErrorException select(ba, 4)
    end
end

@testset "AtributionGraph" begin
    @testset "Constructors" begin
        T = Float32
        n1, n2, n3 = Neuron(T), Neuron(T), Neuron(T)
        ag = neuron_attrib_graph(T, NeuralArea{T})
        contributors = get!(ag, n1)
        push!(contributors, n2)
        contributors = get!(ag, n1)
        push!(contributors, n3)
        @test all(ag[n1] .== [n2, n3])
        @test all(get(ag, n2, [n1]) .== [n1])

    end
end

@testset "StopCriteria" begin
    areas = [NeuralArea(100, 10, 0.1) for i in 1:3]
    currents = [zeros(100) for i in 1:3]
    random_firing!.(areas)

    @testset "MaxIters" begin
        iters = 9
        mi = MaxIters(max_iters=iters, record=true)
        i = 0
        while ! mi(areas, currents)
            i += 1
        end
        @test i == iters
        @test length(mi.spikes) == 3
        @test all([all(size(s) .== (10, iters + 1)) for s in mi.spikes])
        # Without recording
        mi = MaxIters(max_iters=iters, record=false)
        i = 0
        while ! mi(areas, currents)
            i += 1
        end
        @test i == iters
        @test length(mi.spikes) == 1
        @test all(size(mi.spikes[1]) .== (1, 1))
    end

    @testset "Converge" begin
        iters = 9
        cv = Converge(max_iters=iters, record=true)
        i = 0
        while ! cv(areas, currents)
            i += 1
        end
        @test i == iters
        @test length(cv.spikes) == 3
        @test all([all(size(s) .== (10, iters + 1)) for s in cv.spikes])
        @test all(size(cv.distances) .== (3, iters + 1))
        @test all(cv.distances .== 10)

        # Without recording
        cv = Converge(max_iters=iters, record=false)
        i = 0
        while ! cv(areas, currents)
            i += 1
        end
        @test i == iters
        @test length(cv.spikes) == 1
        @test all(size(cv.spikes[1]) .== (1, 1))

        iters = 20
        tol = 0
        fire_past_convergence = 0
        random_firing!.(areas)
        firing = [a.firing for a in areas]
        cv = Converge(
            tol=tol,
            fire_past_convergence=fire_past_convergence,
            max_iters=iters, record=true)
        i = 0
        while ! cv(areas, currents)
            k = i < 10 ? i + 1 : 10
            for (j, a) in enumerate(areas)
                a.firing_prev = firing[j][1:k]
            end
            i += 1
        end
        @test i == 10 - tol + fire_past_convergence

        iters = 20
        tol = 3
        fire_past_convergence = 0
        random_firing!.(areas)
        firing = [a.firing for a in areas]
        cv = Converge(
            tol=tol,
            fire_past_convergence=fire_past_convergence,
            max_iters=iters, record=true)
        i = 0
        while ! cv(areas, currents)
            k = i < 10 ? i + 1 : 10
            for (j, a) in enumerate(areas)
                a.firing_prev = firing[j][1:k]
            end
            i += 1
        end
        @test i == 10 - tol + fire_past_convergence

        iters = 20
        tol = 3
        fire_past_convergence = 5
        random_firing!.(areas)
        firing = [a.firing for a in areas]
        cv = Converge(
            tol=tol,
            fire_past_convergence=fire_past_convergence,
            max_iters=iters, record=true)
        i = 0
        while ! cv(areas, currents)
            k = i < 10 ? i + 1 : 10
            for (j, a) in enumerate(areas)
                a.firing_prev = firing[j][1:k]
            end
            i += 1
        end
        @test i == 10 - tol + fire_past_convergence
    end
end


text:
@auto_hash_equals struct TransLine
    connecting::Tuple{Bus,Bus}
    impedance::Real # Ohms
    capacity::Real # A

    function TransLine(
        connecting::Tuple{Bus,Bus};
        impedfunc=linear_imped,
        capfunc=volt_cap
    )
        return new(
            connecting,
            impedfunc(connecting[1], connecting[2]),
            capfunc(connecting[1], connecting[2])
        )
    end
end

function TransLine(a::Bus, b::Bus; impedfunc=linear_imped,  capfunc=volt_cap)
    return TransLine((a, b), impedfunc=impedfunc, capfunc=capfunc)
end

function show(io::IO, tl::TransLine)
    println(io, "TransLine(")
    println(io, "\tconnecting: $(tl.connecting),")
    println(io, "\timpedance=$(tl.impedance),")
    println(io, "\tcapacity=$(tl.capacity)")
    print(io, ")")
end

"""
    linear_imped(a::Bus, b::Bus)

Return an impedance value for a transmission line connecting buses `a` and `b`, using a
simple linear approximation. This value is loosely based on the two references below.
Values based on resistances for dc currents at 50C and 350kV lines.

REF 1: Electric Power Research Institute (EPRI).
Transmission Line Reference Book, 345 kV and Above.
Palo Alto, CA: The Institute, 1975.
REF 2: J. Glover, M. Sarma and T. Overbye, Power System Analysis and Design.
Stamford, CT: Cengage Learning, 2012.
"""
function linear_imped(a::Bus, b::Bus)
    respkm = 0.045 # Ohms/km
    d = haversine(a, b) # this is using the
    # straight line distance between nodes. Considering the lines arc, real length
    # is slightly longer. However, this average value is already approximated
    # enough.
    return d * respkm
end

"""
    rand_imped(a::Bus, b::Bus)

Return an impedance value for a transmission line connecting buses `a` and `b`, using a
simple linear approximation with a randomly chosen coefficient. The range we adopted is
based on the values for ACSR cables from the ref: J. Glover, M. Sarma and T. Overbye, Power
System Analysis and Design. Stanford, CT: Cengage Learning, 2012.
"""
function rand_imped(a::Bus, b::Bus)
    respkm = 0.025:0.001:0.250
    d = haversine(a, b)
    return d * rand(respkm)
end

"""
    volt_cap(a::Bus, b::Bus)

Return a current carrying capacity value for a transmission line connecting buses `a` and
`b` based on their operating voltages. This value is loosely based on the reference:
Electric Power Research Institute (EPRI). Transmission Line Reference Book, 345 kV and
Above. Palo Alto, CA: The Institute, 1975.
"""
function volt_cap(a::Bus, b::Bus)
    cappvolt = 14 # A/kV
    volt1 = a.voltage
    volt2 = a.voltage
    bvolt = [v for v in volt1 if v in volt2]
    svolt = vcat(a.voltage, b.voltage)
    if length(bvolt) > 0
        return bvolt[rand(1:length(bvolt))] * cappvolt
    else
        return svolt[rand(1:length(svolt))] * cappvolt
    end
end

end_points(t::TransLine) = t.connecting
impedance(t::TransLine) = t.impedance
capacity(t::TransLine) = t.capacity