text:
# https://docs.julialang.org/en/latest/manual/parallel-computing/#The-@threads-Macro-1

# july 5, 2018
# threading does not seem to improve performance when applied to the innermost get_deviances loop

# ramp up
include(joinpath(Pkg.dir("DecisionTrees"), "testing", "prepare_a_dtmtable.jl"))

using StatsBase

function get_deviances_threaded(a::DecisionTrees.PoissonDevianceSplit, current_meanl::Float64, current_meanr::Float64, lo, ooo, f, numerator::Array{Float64,1}, denominator::Array{Float64,1}, weight::Array{Float64,1}, elementsInLeftChildBV)
    # for the poisson deviance consider the derivative of the poisson loss (or google, the reacfin paper or the 'axa' master thesis)
	# this is the core function of the modelling process
	# besides copying of the data, the vast majority of time is spent in here!
	# most of the time is spent here, if we can improve the for loop below, that would improve performance greatly!
	# one possibility would be to introduce parallelization here (which is not straightforward, though

	# dr and dl are 'reused' by each iteration
	dr = Threads.Atomic{Float64}(0.0)
	dl = Threads.Atomic{Float64}(0.0)
	# note: inbounds increases efficiency here (about a factor of 2), however if the bounds are violated something nasty might happen (quote: If the subscripts are ever out of bounds, you may suffer crashes or silent corruption.)
	Threads.@threads for count in 1:length(f)
		@inbounds idx = f.parent.refs[count] + ooo
		@inbounds ni = numerator[count]
        @inbounds wi = weight[count]
		@inbounds eli = elementsInLeftChildBV[idx]
        if eli
            if iszero(ni)
                Threads.atomic_add!(dl, -(ni - wi * current_meanl))
            else
                Threads.atomic_add!(dl, -(ni - wi * current_meanl) - ni * log(wi * current_meanl / ni))
            end
        else
            if iszero(ni)
                Threads.atomic_add!(dr, -(ni - wi * current_meanr))
            else
                Threads.atomic_add!(dr, -(ni - wi * current_meanr) - ni * log(wi * current_meanr / ni))
            end
        end
	end
	return dl[], dr[]
end

function get_deviancesALTERNATIVE(a::DecisionTrees.PoissonDevianceSplit, current_meanl::Float64, current_meanr::Float64, lo, ooo, f, numerator::Array{Float64,1}, denominator::Array{Float64,1}, weight::Array{Float64,1}, elementsInLeftChildBV)
    # for the poisson deviance consider the derivative of the poisson loss (or google, the reacfin paper or the 'axa' master thesis)
	# this is the core function of the modelling process
	# besides copying of the data, the vast majority of time is spent in here!
	# most of the time is spent here, if we can improve the for loop below, that would improve performance greatly!
	# one possibility would be to introduce parallelization here (which is not straightforward, though

	# dr and dl are 'reused' by each iteration
	dr = 0.0 # Threads.Atomic{Float64}(0.0)
	dl = 0.0 # Threads.Atomic{Float64}(0.0)
	# note: inbounds increases efficiency here (about a factor of 2), however if the bounds are violated something nasty might happen (quote: If the subscripts are ever out of bounds, you may suffer crashes or silent corruption.)
	for count in 1:length(f)
		@inbounds idx = f.parent.refs[count] + ooo
		@inbounds ni = numerator[count]
        @inbounds wi = weight[count]
		@inbounds eli = elementsInLeftChildBV[idx]
        if eli
            if iszero(ni)
                dl += -(ni - wi * current_meanl)
            else
                dl += -(ni - wi * current_meanl) - ni * log(wi * current_meanl / ni)
            end
        else
            if iszero(ni)
                dr += -(ni - wi * current_meanr)
            else
                dr += -(ni - wi * current_meanr) - ni * log(wi * current_meanr / ni)
            end
        end
	end
	return dl, dr
end

function check_equality(dataloops, nloops, nsizeNumerator)
# nloops=2#00

    for i = 1:dataloops
        nsize = nsizeNumerator # 15
        denominatoR = rand(nsize)
        numeratoR = rand(nsize)
        weight = rand(nsize)
        trni = sample(1:nsize, floor(Int, .4 * nsize), replace=false)
        sort!(trni)
        psplit = DecisionTrees.PoissonDevianceSplit()

        for thisloop2 = 1:nloops
            strArr = [randstring(1) for i = 1:nsize]
            floatArr = float.(abs.(rand(Int, nsize) .% 99))
            pdStr = DecisionTrees.PooledArray(strArr)
            pdFloat = DecisionTrees.PooledArray(floatArr)

            if rand() < .5
                pd = pdStr
            else
                pd = pdFloat
            end

            fpd = view(pd, trni)
            lo = one(eltype(fpd.parent.refs)) # UInt8(1)
            ooo = one(lo) - lo
            uqrefsmod = sort(unique(fpd.parent.refs)) .+ ooo
            elLBV = BitVector(rand(length(uqrefsmod)) .< rand())
            current_meanl = rand()
            current_meanr = rand()

            resl, resr = DecisionTrees.get_deviances(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
            resl2, resr2 =    get_deviances_threaded(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
            resl3, resr3 = get_deviancesALTERNATIVE(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
            @assert isapprox(resl, resl2)
            @assert isapprox(resl, resl3)
            @assert isapprox(resr, resr2)
            @assert isapprox(resr, resr3)
        end
    end

    return nothing

end

function timeOLD(dataloops, nloops, nsizeNumerator)
    # nloops=2#00

    for i = 1:dataloops
        nsize = nsizeNumerator # 15
        denominatoR = rand(nsize)
        numeratoR = rand(nsize)
        weight = rand(nsize)
        trni = sample(1:nsize, floor(Int, .4 * nsize), replace=false)
        sort!(trni)
        psplit = DecisionTrees.PoissonDevianceSplit()

        for thisloop2 = 1:nloops
            strArr = [randstring(1) for i = 1:nsize]
            floatArr = float.(abs.(rand(Int, nsize) .% 99))
            pdStr = DecisionTrees.PooledArray(strArr)
            pdFloat = DecisionTrees.PooledArray(floatArr)

            if rand() < .5
                pd = pdStr
            else
                pd = pdFloat
            end

            fpd = view(pd, trni)
            lo = one(eltype(fpd.parent.refs)) # UInt8(1)
            ooo = one(lo) - lo
            uqrefsmod = sort(unique(fpd.parent.refs)) .+ ooo
            elLBV = BitVector(rand(length(uqrefsmod)) .< rand())
            current_meanl = rand()
            current_meanr = rand()

            resl, resr = DecisionTrees.get_deviances(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
            # resl2,resr2=    get_deviances_threaded(psplit,current_meanl,current_meanr,lo,ooo,fpd,numeratoR,denominatoR,weight,elLBV)
            # @assert isapprox(resl,resl2)
            # @assert isapprox(resr,resr2)
        end
    end

    return nothing

end

function timeNEW(dataloops, nloops, nsizeNumerator)
    # nloops=2#00

    for i = 1:dataloops
        nsize = nsizeNumerator # 15
        denominatoR = rand(nsize)
        numeratoR = rand(nsize)
        weight = rand(nsize)
        trni = sample(1:nsize, floor(Int, .4 * nsize), replace=false)
        sort!(trni)
        psplit = DecisionTrees.PoissonDevianceSplit()

        for thisloop2 = 1:nloops
            strArr = [randstring(1) for i = 1:nsize]
            floatArr = float.(abs.(rand(Int, nsize) .% 99))
            pdStr = DecisionTrees.PooledArray(strArr)
            pdFloat = DecisionTrees.PooledArray(floatArr)

            if rand() < .5
                pd = pdStr
            else
                pd = pdFloat
            end

            fpd = view(pd, trni)
            lo = one(eltype(fpd.parent.refs)) # UInt8(1)
            ooo = one(lo) - lo
            uqrefsmod = sort(unique(fpd.parent.refs)) .+ ooo
            elLBV = BitVector(rand(length(uqrefsmod)) .< rand())
            current_meanl = rand()
            current_meanr = rand()

            # resl,resr=DecisionTrees.get_deviances(psplit,current_meanl,current_meanr,lo,ooo,fpd,numeratoR,denominatoR,weight,elLBV)
            resl2, resr2 =    get_deviances_threaded(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
            # @assert isapprox(resl,resl2)
            # @assert isapprox(resr,resr2)
        end
    end

    return nothing

end

function timeALT(dataloops, nloops, nsizeNumerator)
        # nloops=2#00

    for i = 1:dataloops
        nsize = nsizeNumerator # 15
        denominatoR = rand(nsize)
        numeratoR = rand(nsize)
        weight = rand(nsize)
        trni = sample(1:nsize, floor(Int, .4 * nsize), replace=false)
        sort!(trni)
        psplit = DecisionTrees.PoissonDevianceSplit()

        for thisloop2 = 1:nloops
            strArr = [randstring(1) for i = 1:nsize]
            floatArr = float.(abs.(rand(Int, nsize) .% 99))
            pdStr = DecisionTrees.PooledArray(strArr)
            pdFloat = DecisionTrees.PooledArray(floatArr)

            if rand() < .5
                pd = pdStr
            else
                pd = pdFloat
            end

            fpd = view(pd, trni)
            lo = one(eltype(fpd.parent.refs)) # UInt8(1)
            ooo = one(lo) - lo
            uqrefsmod = sort(unique(fpd.parent.refs)) .+ ooo
            elLBV = BitVector(rand(length(uqrefsmod)) .< rand())
            current_meanl = rand()
            current_meanr = rand()

                # resl,resr=DecisionTrees.get_deviances(psplit,current_meanl,current_meanr,lo,ooo,fpd,numeratoR,denominatoR,weight,elLBV)
            resl2, resr2 =    get_deviancesALTERNATIVE(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
                # @assert isapprox(resl,resl2)
                # @assert isapprox(resr,resr2)
        end
    end

    return nothing

end

@time check_equality(10, 20, 200000);

@time timeOLD(10, 20, 200000);
@time timeNEW(10, 20, 200000); # 50% more than old (with 1 thread only!)
@time timeALT(10, 20, 200000); # same as old

using BenchmarkTools

@benchmark get_deviances_threaded(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
@benchmark get_deviancesALTERNATIVE($psplit, $current_meanl, $current_meanr, $lo, $ooo, $fpd, $numeratoR, $denominatoR, $weight, $elLBV)
@benchmark DecisionTrees.get_deviances($psplit, $current_meanl, $current_meanr, $lo, $ooo, $fpd, $numeratoR, $denominatoR, $weight, $elLBV)

@code_warntype get_deviancesALTERNATIVE(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
@code_warntype DecisionTrees.get_deviances(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)
@code_warntype get_deviances_threaded(psplit, current_meanl, current_meanr, lo, ooo, fpd, numeratoR, denominatoR, weight, elLBV)


text:
using PyPlot
include("args.jl")

Sw = constant(collect(0:0.001:1))

lambda_brie_3 = Brie(Sw)

ext: jl


text:
using LightGraphs
# using Hargreaves
include("C:/Users/Anton Hinneck/Desktop/Libya/Hargreaves.jl/src/Hargreaves.jl")
using .Hargreaves
# dir = @__DIR__
# cd(dir)

# Amount of vertices
# node_count = 31

# Upper triangular matrix with random values
# weights = Array{Int64, 2}(node_count,node_count)

include("C:/Users/Anton Hinneck/juliaPackages/GitHub/PowerGrids.jl/src/PowerGrids.jl")
using .PowerGrids
dir = @__DIR__
cd(dir)

set_csv_path("C:/Users/Anton Hinneck/Documents/Git/pglib2csv/pglib/2020-08-21.19-54-30-275/csv")
PowerGrids.csv_cases(verbose = true)
#PowerGrids.select_csv_case(48)
PowerGrids.select_csv_case(3)
case = PowerGrids.loadCase() # 118 Bus ieee
g = toGraph(case)

_weights = zeros(length(vertices(g)), length(vertices(g)))
for i in vertices(g)

    for j in vertices(g)

        _weights[j,i] = 0

        if j < i

            _weights[j,i] = 0
            _weights[j,i] = 1#round(Int, rand() * 20 + 1)

        else

            _weights[j,i] = 0

        end
    end
end

# Some random relatively high weight
# weights[2,15] = 30

# Define graph
#g = CompleteGraph(node_count)
cd(dir)
wireplot(g, "wireplot_1", distmx = _weights, wireplot_edge_color_scheme = :purple, static_widths = true, static_colors = true)
wireplot(g, "wireplot_2", distmx = _weights, wireplot_edge_color_scheme = :purple, color_scale = :quadratic, width_scale = :quadratic, max_width = 8)


text:
using Jags, Mamba, Gadfly

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  ## Data
  rats = Dict{String, Any}(
    "Y" =>
      [151 199 246 283 320;
       145 199 249 293 354;
       147 214 263 312 328;
       155 200 237 272 297;
       135 188 230 280 323;
       159 210 252 298 331;
       141 189 231 275 305;
       159 201 248 297 338;
       177 236 285 350 376;
       134 182 220 260 296;
       160 208 261 313 352;
       143 188 220 273 314;
       154 200 244 289 325;
       171 221 270 326 358;
       163 216 242 281 312;
       160 207 248 288 324;
       142 187 234 280 316;
       156 203 243 283 317;
       157 212 259 307 336;
       152 203 246 286 321;
       154 205 253 298 334;
       139 190 225 267 302;
       146 191 229 272 302;
       157 211 250 285 323;
       132 185 237 286 331;
       160 207 257 303 345;
       169 216 261 295 333;
       157 205 248 289 316;
       137 180 219 258 291;
       153 200 244 286 324],
    "x" => [8.0, 15.0, 22.0, 29.0, 36.0]
  )
  rats["N"] = size(rats["Y"], 1)
  rats["T"] = size(rats["Y"], 2)
  rats["x.bar"] = mean(rats["x"])

  ratsmodel = "
  model {
      for (i in 1:N) {
         for (j in 1:T) {
            mu[i,j] <- alpha[i] + beta[i]*(x[j] - x.bar);
            Y[i,j]   ~ dnorm(mu[i,j], tau.c)
         }
         alpha[i] ~ dnorm(alpha.c, tau.alpha);
         beta[i]  ~ dnorm(beta.c, tau.beta);
      }
      alpha.c   ~ dnorm(0, 1.0E-4);
      beta.c    ~ dnorm(0, 1.0E-4);
      tau.c     ~ dgamma(1.0E-3, 1.0E-3);
      tau.alpha ~ dgamma(1.0E-3, 1.0E-3);
      tau.beta  ~ dgamma(1.0E-3, 1.0E-3);
      sigma    <- 1.0/sqrt(tau.c);
      alpha0   <- alpha.c - beta.c*x.bar;
  }
  "

  ## Initial Values
  inits = [
    Dict{String,Any}("alpha" => fill(250, 30), "beta" => fill(6, 30),
    "alpha.c" => 100, "beta.c" => 2, "tau.c" => 1, "tau.alpha" => 1, "tau.beta" => 1),
    Dict{String,Any}("alpha" => fill(150, 30), "beta" => fill(3, 30),
    "alpha.c" => 150, "beta.c" => 2, "tau.c" => 1, "tau.alpha" => 1, "tau.beta" => 1),
    Dict{String,Any}("alpha" => fill(200, 30), "beta" => fill(6, 30),
    "alpha.c" => 200, "beta.c" => 1, "tau.c" => 1, "tau.alpha" => 1, "tau.beta" => 1),
    Dict{String,Any}("alpha" => fill(150, 30), "beta" => fill(3, 30),
    "alpha.c" => 250, "beta.c" => 1, "tau.c" => 1, "tau.alpha" => 1,"tau.beta" => 1)
  ]

  monitors = Dict{String, Bool}(
    "alpha0" => true,
    "beta.c" => true,
    "sigma" =>true
  )

  jagsmodel = Jagsmodel(
    name="rats",
    model=ratsmodel,
    monitor=monitors,
    #ncommands=4, nchains=1,
    adapt=1000, update=10000, thin=10,
    #deviance=true, dic=true, popt=true,
    pdir=ProjDir);

  println("Jagsmodel that will be used:")
  jagsmodel |> display
  println("Input observed data dictionary:")
  rats |> display
  println("\nInput initial values dictionary:")
  inits |> display
  println()

  @time sim = jags(jagsmodel, rats, inits, ProjDir)
  describe(sim)
  println()

  ## Plotting
  p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:svg)
  draw(p, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:pdf)

  # Below will only work on OSX, please adjust for your environment.
  # JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
  @static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
          for i in 1:4
            isfile("$(jagsmodel.name)-summaryplot-$(i).svg") &&
              run(`open -a $(JULIA_SVG_BROWSER) "$(jagsmodel.name)-summaryplot-$(i).svg"`)
          end
        end : println()

end #cd


text:
type Linear{A<:AbstractArray} <: Layer{A, A}
    output :: A
    gradInput :: A
    weight :: A
    bias :: A
    gradWeight :: A
    gradBias :: A
end

# input and output are 1d
Linear(T::Type, batch_size, input_size, output_size) =
    Linear(Array(T, (output_size, batch_size)),
           Array(T, (input_size, batch_size)),
           Array(T, (output_size, input_size)),
           Array(T, (output_size, 1)),
           Array(T, (output_size, input_size)),
           Array(T, (output_size, 1)))

function updateOutput!{A}(self :: Linear{A}, input :: A)
    broadcast!(+, self.output, self.weight * input, self.bias)
    return self.output
end

function updateGradInput!{A}(self :: Linear{A}, input :: A,
                                gradOutput :: A)
    self.gradInput = transpose(self.weight) * gradOutput
    return self.gradInput
end

function accGradParameters!{A}(self :: Linear{A}, input :: A,
                                  gradOutput :: A, scale :: eltype(A))
    self.gradWeight = self.gradWeight + scale * gradOutput * transpose(input)
    self.gradBias = self.gradBias + scale * sum(gradOutput, 2)
end

export Linear


text:
export gemm_batch

doc"""
    gemv(tA::Char, alpha, A::Var, x::Var)

* tA: 'T' (transpose) or 'N' (not transpose)

```math
y = \alpha \times \textrm{tA}(A) \times x
```

```julia
T = Float32
A = Var(rand(T,10,5))
x = Var(rand(T,5))
B = gemv('N', 1, A, x)
```
"""
function gemv(tA::Char, alpha::Number, A::Var, x::Var)
    T = eltype(A)
    ydata = gemv(tA, T(alpha), A.data, x.data)
    Var(ydata, ∇gemv!, (tA,alpha,A,x))
end
gemv(tA, alpha, A::Node, x) = Node(gemv, (tA,alpha,A,x))

function ∇gemv!(y::Var, tA::Char, alpha::Number, A::Var, x::Var)
    T = eltype(A)
    if !isnothing(A.grad)
        gy = reshape(y.grad, length(y.grad), 1)
        xdata = reshape(x.data, length(x.data), 1)
        tA == 'N' ?
        gemm!('N', 'T', T(alpha), gy, xdata, T(1), A.grad) :
        gemm!('N', 'T', T(alpha), xdata, gy, T(1), A.grad)
    end
    if !isnothing(x.grad)
        gemv!(tA=='N' ? 'T' : 'N', T(alpha), A.data, y.grad, T(1), x.grad)
    end
end

doc"""
    gemm(tA::Char, tB::Char, alpha, A::Var, B::Var)
    gemm(A::Var, B::Var, [tA='N'], [tB='N'], [alpha=1])

* tA, tB: 'T' (transpose) or 'N' (not transpose)

```math
C = \alpha \times \textrm{tA}(A) \times \textrm{tB}(B)
```

```julia
T = Float32
A = Var(rand(T,10,5))
B = Var(rand(T,10,7))
C = gemm('T', 'N', 1, A, B)
```
"""
function gemm(tA::Char, tB::Char, alpha::Number, A::Var, B::Var)
    T = eltype(A)
    ydata = gemm(tA, tB, T(alpha), A.data, B.data)
    Var(ydata, ∇gemm!, (tA,tB,alpha,A,B))
end
gemm(tA, tB, alpha, A::Node, B) = Node(gemm, (tA,tB,alpha,A,B))
gemm(tA, tB, alpha, A, B::Node) = Node(gemm, (tA,tB,alpha,A,B))

function ∇gemm!(C::Var, tA::Char, tB::Char, alpha::Number, A::Var, B::Var)
    isa(A.grad,Nothing) || ∇gemm_A!(C.grad, tA, tB, alpha, A.grad, B.data)
    isa(B.grad,Nothing) || ∇gemm_B!(C.grad, tA, tB, alpha, A.data, B.grad)
end

function ∇gemm_A!(gC, tA::Char, tB::Char, alpha, gA, B)
    T = eltype(gC)
    if tA == 'N'
        gemm!('N', tB=='N' ? 'T' : 'N', T(alpha), gC, B, T(1), gA)
    else
        gemm!(tB, 'T', T(alpha), B, gC, T(1), gA)
    end
end

function ∇gemm_B!(gC, tA::Char, tB::Char, alpha, A, gB)
    T = eltype(gC)
    if tB == 'N'
        gemm!(tA=='N' ? 'T' : 'N', 'N', T(alpha), A, gC, T(1), gB)
    else
        gemm!('T', tA, T(alpha), gC, A, T(1), gB)
    end
end


text:
using EngEconomics, Roots, Plots, ColorSchemes
plotly()

uwEngColorScheme = ColorScheme([parse(Colorant, "#D0B4E7"), parse(Colorant, "#57058B"), parse(Colorant, "#8100B4"), parse(Colorant, "#000000")])

# Given
PT = 100
PA = 150
lifeT = 5
lifeA = 5
AT = 50
AA = 62
FT = 20
FA = 30

# Find
# (a.) Construct a break even graph showing the present worth of each
#      alternative as a function of interest rates between 6 percent and 20
#      percent.
i = collect(0.06:0.01:0.2)
NPWT = -PT .+ AT * seriesPresentAmountFactor.(i, lifeT) .+ FT * presentWorthFactor.(i, lifeT)
NPWA = -PA .+ AA * seriesPresentAmountFactor.(i, lifeA) .+ FA * presentWorthFactor.(i, lifeA)

plot(i, NPWT, label="Model T", xlabel="Interest Sensitivity (6%-20%)",
	ylabel="Present Worth (\$1000)",
	title="Sensitivity of Net Present Worth of Models T and A to Interest Rate",
	palette=cgrad(uwEngColorScheme))
plot!(i, NPWA, label="Model A")

# (b.) Which is the preferred choice at 8 percent interest?
i8 = findall(x->x == 0.08, i)
@assert length(i8) == 1
i8 = i8[1]
NPWT_i8 = NPWT[i8]
NPWA_i8 = NPWA[i8]

if NPWT_i8 > NPWA_i8
	println("Model T is preferred because its NPW at i=8% is greater than that of Model A")
else
	println("Model A is preferred because its NPW at i=8% is greater than that of Model T")
end

# (c.) Which is the preferred choice at 16 percent interest?
i16 = findall(x->x == 0.16, i)
@assert length(i16) == 1
i16 = i16[1]
NPWT_i16 = NPWT[i16]
NPWA_i16 = NPWA[i16]

if NPWT_i16 > NPWA_i16
	println("Model T is preferred because its NPW at i=16% is greater than that of Model A")
else
	println("Model A is preferred because its NPW at i=16% is greater than that of Model T")
end

# (d.) What is the break-even interest rate?
NPWT_func(interest) = -PT .+ AT * seriesPresentAmountFactor.(interest, lifeT) .+ FT * presentWorthFactor.(interest, lifeT)
NPWA_func(interest) = -PA .+ AA * seriesPresentAmountFactor.(interest, lifeA) .+ FA * presentWorthFactor.(interest, lifeA)
f(interest) = NPWT_func(interest) - NPWA_func(interest)
iBreakEven = find_zero(f, 0.11)

println("The break even interest rate is roughly ", iBreakEven)


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4a.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    GU.read("ospy/data4/sites.shp") do shp
        GU.read("ospy/data4/aster.img") do img
            shplayer = GU.fetchlayer(shp, 0)
            featuredefn = GU.getlayerdefn(shplayer)
            id_index = GU.getfieldindex(featuredefn, "ID")

            # get georeference info
            transform = GU.getgeotransform(img)
            xOrigin = transform[1]
            yOrigin = transform[4]
            pixelWidth = transform[2]
            pixelHeight = transform[6]

            # loop through the features in the shapefile
            for feature in shplayer
                geom = GU.getgeom(feature)
                x = GU.getx(geom, 0)
                y = GU.gety(geom, 0)
                # compute pixel offset
                xOffset = round(Int, (x - xOrigin) / pixelWidth)
                yOffset = round(Int, (y - yOrigin) / pixelHeight)
                # create a string to print out
                s = string(GU.fetchfield(feature, id_index)) * " "
                for j in 1:GU.nraster(img)
                    data = GU.fetch(img, j, 1, 1, xOffset, yOffset)
                    s = s * string(data[1,1]) * " "
                end
                # print out the data string
                println(s)
            end
        end
    end
end

# Expected Output:
# 1 46 30 47
# 2 85 77 76
# 3 64 51 66
# 4 71 60 73
# 5 59 40 68
# 6 97 95 80
# 7 51 33 52
# 8 79 69 62
# 9 78 74 63
# 10 90 81 67
# 11 82 81 67
# 12 100 88 68
# 13 90 80 65
# 14 76 65 81
# 15 83 75 61
# 16 73 61 63
# 17 81 71 68
# 18 69 57 60
# 19 98 90 75
# 20 91 79 76
# 21 86 77 64
# 22 78 73 63
# 23 82 66 71
# 24 64 46 82
# 25 69 55 68
# 26 98 88 78
# 27 62 43 79
# 28 68 57 62
# 29 66 51 74
# 30 60 39 83
# 31 80 62 80
# 32 86 73 86
# 33 97 86 84
# 34 47 29 49
# 35 91 74 90
# 36 67 53 63
# 37 71 59 63
# 38 54 23 16
# 39 53 22 17
# 40 76 66 61
# 41 0 0 0
# 42 0 0 0


text:
using DashBootstrapComponents

pagination = html_div([
    html_div("Collapse long pagination objects using ellipsis"),
    dbc_pagination(max_value = 10, fully_expanded = false),
    html_div("If space won't be saved, it won't be collapsed"),
    dbc_pagination(max_value = 5, fully_expanded = false),
])

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
"""
    savitsky_golay(x::Vector, windowSize::Integer, polyOrder::Integer; deriv::Integer=0)

Polynomial smoothing with the Savitsky Golay filters.

Copied from https://github.com/blakejohnson/Qlab.jl/blob/master/src/SavitskyGolay.jl

#### Sources

Theory: http://www.ece.rutgers.edu/~orfanidi/intro2sp/orfanidis-i2sp.pdf
Python example: http://wiki.scipy.org/Cookbook/SavitzkyGolay
"""
function savitsky_golay(x::Vector, windowSize::Integer, polyOrder::Integer; deriv::Integer=0)

    #Some error checking
    @assert isodd(windowSize) "Window size must be an odd integer."
    @assert polyOrder < windowSize "Polynomial order must me less than window size."

    halfWindow = round(Int, (windowSize-1) / 2)

    #Setup the S matrix of basis vectors.
    S = zeros(windowSize, polyOrder+1)
    for ct = 0:polyOrder
        S[:,ct+1] = collect(-halfWindow:halfWindow).^(ct)
    end

    #Compute the filter coefficients for all orders
    #From the scipy code it seems pinv(S) and taking rows should be enough
    G = S*pinv(S'*S)

    #Slice out the derivative order we want
    filterCoeffs = G[:,deriv+1] * factorial(deriv);

    #Pad the signal with the endpoints and convolve with filter
    paddedX = [x[1]*ones(halfWindow); x; x[end]*ones(halfWindow)]
    y = conv(filterCoeffs[end:-1:1], paddedX)

    #Return the valid midsection
    return y[2*halfWindow+1:end-2*halfWindow]

end

"""
    trend(::SSR, channel::Int)

Determine the trend for a steady state response signal.
"""
function trend(a::SSR, chan::Int; secs = 0.5, order = 2, kwargs...)

    ls = size(a.data, 1)

    windowSize = round(Int, 8192 * secs)
    if iseven(windowSize)
        windowSize += 1
    end

    savitsky_golay(vec(a.data[:, chan]), windowSize, order)
end

"""
    detrend(::SSR)

Remove the trend from a steady state response.
"""
function detrend(a::SSR; kwargs...)

    for c in 1:size(a.data, 2)
        a.data[:, c] = a.data[:, c] - trend(a, c; kwargs...)
    end

    return a
end


text:
function optimize_linear_constraints(f, A, f_g, f_h, yinit=0)
    """
    Convert constrained optimization problem min f(x) s.t. Ax = b
        to unconstrained problem min g(y) = f(B*y + x0), where x0 is a solution to Ax = b
        and B is a basis for the kernel of A.

    Defaults to using Newton's line search algorithm.

    f: constrained objective
    A: constraint matrix
    f_g: gradient of f
    f_h: hessian of f
    yinit: init y with a vector of all yinit values
    """

    m, n = size(A)
    b = zeros(m)
    x0 = A \ b
    B = nullspace(A)

    y0 = ones(size(B)[2]) .* yinit

    g(y) = f(B*y + x0)
    g_g(y) = B' * f_g(B*y + x0)
    g_h(y) = B' * f_h(B*y + x0) * B

    yvals = line_search(g, y0, g_g, g_h, "newton", 5000)
    xvals = [B*y + x0 for y in yvals]
    lambdas = [A' \ f_g(x) for x in xvals]

    return yvals, xvals, lambdas
end


text:
# project to interval if we are not on the interview
# TODO: need to work out how to set piecewise domain

/(c::Number,f::Fun{S}) where {S<:ContinuousSpace} = Fun(map(f->c/f,components(f)),PiecewiseSpace)
^(f::Fun{S},c::Integer) where {S<:ContinuousSpace} = Fun(map(f->f^c,components(f)),PiecewiseSpace)
^(f::Fun{S},c::Number) where {S<:ContinuousSpace} = Fun(map(f->f^c,components(f)),PiecewiseSpace)

^(f::Fun{<:PolynomialSpace},k::Integer) = intpow(f,k)

#TODO: implement
^(f::Fun{Jacobi},k::Integer) = intpow(f,k)
^(f::Fun{Jacobi},k::Real) = Fun(f,Chebyshev)^k

function log(f::Fun{<:PolynomialSpace{<:IntervalOrSegment}})
    g = log(setdomain(f, ChebyshevInterval()))
    setdomain(g, domain(f))
end

# ODE gives the first order ODE a special function op satisfies,
# RHS is the right hand side
# growth says what to use to choose a good point to impose an initial condition
for (op,ODE,RHS,growth) in ((:(exp),"D-f'","0",:(real)),
                            (:(asinh),"sqrt(f^2+1)*D","f'",:(real)),
                            (:(acosh),"sqrt(f^2-1)*D","f'",:(real)),
                            (:(atanh),"(1-f^2)*D","f'",:(real)),
                            (:(erfcx),"D-2f*f'","-2f'/sqrt(π)",:(real)),
                            (:(dawson),"D+2f*f'","f'",:(real)))
    L,R = Meta.parse(ODE),Meta.parse(RHS)
    @eval $op(f::Fun{<:ContinuousSpace}) = Fun(map(f->$op(f),components(f)),PiecewiseSpace)
end

for OP in (:abs,:sign,:log,:angle)
    @eval $OP(f::Fun{<:ContinuousSpace{<:Any,<:Real},<:Real}) =
            Fun(map($OP,components(f)),PiecewiseSpace)
end

sin(f::Fun{S,T}) where {S<:Union{Ultraspherical,Chebyshev},T<:Real} = imag(exp(im*f))
cos(f::Fun{S,T}) where {S<:Union{Ultraspherical,Chebyshev},T<:Real} = real(exp(im*f))

## Second order functions with parameter ν

for (op,ODE,RHS,growth) in ((:(hankelh1),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D+(f^2-ν^2)*f'^3","0",:(imag)),
                            (:(hankelh2),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D+(f^2-ν^2)*f'^3","0",:(imag)),
                            (:(besselj),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D+(f^2-ν^2)*f'^3","0",:(imag)),
                            (:(bessely),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D+(f^2-ν^2)*f'^3","0",:(imag)),
                            (:(besseli),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D-(f^2+ν^2)*f'^3","0",:(real)),
                            (:(besselk),"f^2*f'*D^2+(f*f'^2-f^2*f'')*D-(f^2+ν^2)*f'^3","0",:(real)),
                            (:(besselkx),"f^2*f'*D^2+((-2f^2+f)*f'^2-f^2*f'')*D-(f+ν^2)*f'^3","0",:(real)),
                            (:(hankelh1x),"f^2*f'*D^2+((2im*f^2+f)*f'^2-f^2*f'')*D+(im*f-ν^2)*f'^3","0",:(imag)),
                            (:(hankelh2x),"f^2*f'*D^2+((-2im*f^2+f)*f'^2-f^2*f'')*D+(-im*f-ν^2)*f'^3","0",:(imag)))
    L,R = Meta.parse(ODE),Meta.parse(RHS)
    @eval begin
        function $op(ν,fin::Fun{S,T}) where {S<:Union{Ultraspherical,Chebyshev},T}
            f=setcanonicaldomain(fin)

            g=chop($growth(f),eps(T))
            xmin = isempty(g.coefficients) ? leftendpoint(domain(g)) : argmin(g)
            xmax = isempty(g.coefficients) ? rightendpoint(domain(g)) : argmax(g)
            opfxmin,opfxmax = $op(ν,f(xmin)),$op(ν,f(xmax))
            opmax = maximum(abs,(opfxmin,opfxmax))
            while opmax≤10eps(T) || abs(f(xmin)-f(xmax))≤10eps(T)
                xmin,xmax = rand(domain(f)),rand(domain(f))
                opfxmin,opfxmax = $op(ν,f(xmin)),$op(ν,f(xmax))
                opmax = maximum(abs,(opfxmin,opfxmax))
            end
            D=Derivative(space(f))
            B=[Evaluation(space(f),xmin),Evaluation(space(f),xmax)]
            u=\([B;eval($L)],[opfxmin;opfxmax;eval($R)];tolerance=eps(T)*opmax)

            setdomain(u,domain(fin))
        end
    end
end

#TODO ≤,≥

## Piecewise Space

# Return the locations of jump discontinuities
#
# Non Piecewise Spaces are assumed to have no jumps.
function jumplocations(f::Fun)
    eltype(domain(f))[]
end

# Return the locations of jump discontinuities
function jumplocations(f::Fun{S}) where{S<:Union{PiecewiseSpace,ContinuousSpace}}
    d = domain(f)

    if ncomponents(d) < 2
      return eltype(domain(f))[]
    end

    dtol=10eps(eltype(d))
    ftol=10eps(cfstype(f))

    dc = components(d)
    fc = components(f)

    isjump = isapprox.(leftendpoint.(dc[2:end]), rightendpoint.(dc[1:end-1]), rtol=dtol) .&
           .!isapprox.(first.(fc[2:end]), last.(fc[1:end-1]), rtol=ftol)

    locs = rightendpoint.(dc[1:end-1])
    locs[isjump]
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(pset::PointSet{Dim,T}) where {Dim,T}
  X = coordinates(pset, 1:nelements(pset))

  seriestype --> :scatter
  seriescolor --> :black
  legend --> false

  if Dim == 1
    @series begin
      X[1,:], fill(0, nelements(pset))
    end
  elseif Dim == 2
    aspect_ratio --> :equal
    @series begin
      X[1,:], X[2,:]
    end
  elseif Dim == 3
    aspect_ratio --> :equal
    @series begin
      X[1,:], X[2,:], X[3,:]
    end
  else
    @error "cannot plot in more than 3 dimensions"
  end
end

@recipe function f(pset::PointSet{Dim,T}, data::AbstractVector) where {Dim,T}
  X = coordinates(pset, 1:nelements(pset))

  seriestype --> :scatter

  if Dim == 1
    X[1,:], data
  elseif Dim == 2
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:]
  elseif Dim == 3
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:], X[3,:]
  else
    @error "cannot plot in more than 3 dimensions"
  end
end


text:
"""
2D Cartesian grid functionality
"""
module Cartesian

export CartesianGrid, nb_nodes, origin, coordinates, cartesianindices, nb_neighbors, setneighbors!, dirichletrhs, xrange, yrange

"""
Represents a 2D cartesian grid with uniform spacing
"""
struct CartesianGrid
  nx::Int # Number of points in the X direction
  ny::Int # Number of points in the Y direction
  h::Float64 # Spacing
end

cartesianindices(g::CartesianGrid) = CartesianIndices((g.nx, g.ny))
Base.LinearIndices(g::CartesianGrid) = LinearIndices((g.nx, g.ny))

"""
Origin of the mesh
"""
origin(g::CartesianGrid) = (g.h*(g.nx-1)/2, g.h*(g.ny-1)/2)

coordrange(g::CartesianGrid, dim) = -origin(g)[dim]:g.h:origin(g)[dim]
xrange(g::CartesianGrid) = coordrange(g,1)
yrange(g::CartesianGrid) = coordrange(g,2)

"""
Number of nodes in the mesh
"""
nb_nodes(g::CartesianGrid) = g.nx*g.ny

struct LinearCoordinates{T<:Number} <: AbstractArray{T,2}
  grid::CartesianGrid
  origin::NTuple{2,T}
  LinearCoordinates{T}(g::CartesianGrid) where T = new(g, origin(g))
end

coordinates(g::CartesianGrid) = LinearCoordinates{Float64}(g)

Base.IndexStyle(::Type{LinearCoordinates{T}}) where {T} = IndexCartesian()

@inline function Base.getindex(coords::LinearCoordinates{T}, i::Int, j::Int) where {T}
  tocartesian = cartesianindices(coords.grid)
  return (tocartesian[i][j]-1)*coords.grid.h - coords.origin[j]
end

Base.size(coords::LinearCoordinates) = (nb_nodes(coords.grid), 2)

@inline function laplace2D_stencil(I::CartesianIndex)
  return CartesianIndex.((0,0,-1, 0,1),
                         (0,1, 0,-1,0)) .+ Ref(I)
end

function nb_neighbors(grid::CartesianGrid, gid::Integer, f::Function)
  return count(i -> i ∈ cartesianindices(grid) && f(grid, i), laplace2D_stencil(cartesianindices(grid)[gid]))
end

function setneighbors!(array, startidx, grid::CartesianGrid, gid::Integer, f::Function)
  cartinds = cartesianindices(grid)
  lininds = LinearIndices(grid)
  i = startidx
  for I in laplace2D_stencil(cartinds[gid])
    if I ∈ cartinds && f(grid,I)
      array[i] = lininds[I]
      i += 1
    end
  end
  return i - startidx
end

function dirichletrhs(grid::CartesianGrid, gid::Integer, isdirichlet::Function, dirichletval::Function)
  cartinds = cartesianindices(grid)
  lininds = LinearIndices(grid)
  result = 0.0
  for I in laplace2D_stencil(cartinds[gid])
    if I ∈ cartinds && isdirichlet(grid,I)
      result += dirichletval(grid, I)
    end
  end
  return result
end

end


text:
## menus in JGUI

using JGUI

w = window(title="menus")
b = button(w, "button"); push!(w,b)

help = action(w, "help", () -> println("hi"))
file = action(w, "file", () -> println("file"))

m = menubar(w)
sm = addMenu(m, "File")
push!(sm, help)
push!(sm, separator(sm))
push!(sm, file) ## push!

sm = addMenu(m, "Hey")
push!(sm, help)
push!(sm, file) ## push!

## radio button and check button values
if !JGUI.isgtk()
    sm = addMenu(m, "buttons")
    ## radio button
    rb = radiogroup(sm, ["one", "two", "three"])
    widget = push!(sm, rb)
    push!(sm,separator(sm))
    ## check button
    cb = checkbox(sm, true, "label")
    push!(sm, cb)
    push!(sm, separator(sm))

    push!(sm, action(w, "values", () -> println((rb[:value], cb[:value]))))
end

## popup menus
pm = menu(b)
b[:context] = "ping"
push!(pm, action(w, "foo", () -> begin println(b[:context]); b[:context] = "foo" end))
push!(pm, action(w, "bar", () -> begin println(b[:context]); b[:context] = "bar" end))

raise(w)


text:
# predefined adaptors for working with types from the Julia standard library

adapt_structure(to, xs::Union{Tuple,NamedTuple}) = map(x->adapt(to,x), xs)

## Broadcast

import Base.Broadcast: Broadcasted, Extruded

adapt_structure(to, bc::Broadcasted{Style}) where Style =
  Broadcasted{Style}(bc.f, map(arg->adapt(to, arg), bc.args), bc.axes)

adapt_structure(to, ex::Extruded) =
    Extruded(adapt(to, ex.x), ex.keeps, ex.defaults)


text:
function point_cloud(R::AbstractRoute; remove_station=true)
    zone_idx, _, _ = split_by_zone(R, remove_station=remove_station)
    S = stops(R, remove_station=remove_station)

    lat = latitude.(S)
    long = longitude.(S)

    points = [Singleton([x, y]) for (x, y) in zip(long, lat)]

    out = [points[idx] for idx in zone_idx]

    return out
end

function path(R::AbstractRoute; remove_station=true)
    P = point_cloud(R, remove_station=remove_station)

    num_zones = length(P)
    out = [Vector{LineSegment{Float64, Vector{Float64}}}() for _ in 1:num_zones]

    for (k, points) in enumerate(P) # loop over zones
        m = length(points)
        X = Vector{LineSegment{Float64, Vector{Float64}}}()
        for i in 1:m-1 # loop over points for the current zone
            L = LineSegment(element(points[i]), element(points[i+1]))
            push!(out[k], L)
        end

        # add last segment
        if k < num_zones
            push!(out[k], LineSegment(element(points[m]), element(P[k+1][1])))
        end
    end
    return out
end

#=
FIXME add plot recipe

function plot_sequence(R; kwargs...)
    fig = plot()
    plot_sequence!(fig, R; kwargs...)
end

function plot_sequence!(fig, R; path_color=:black, path_style=:solid, path_width=1.0)
    vertices = point_cloud(R, remove_station=true);
    edges = path(R, remove_station=true)

    colores = Vector{String}()
    i = 1
    for (k, v) in Colors.color_names
        if mod(i, 10) == 0
            push!(colores, k)
        end
        i += 1
    end

    for i in 1:length(edges)
        plot!(fig, edges[i], marker=:none, seriestype=:path, lw=path_width, alpha=1., c=path_color, ls=path_style)
        plot!(fig, vertices[i], alpha=1, color=colores[i])

        #for (k, vi) in enumerate(vertices[i])
        #    s = vi.element
        #    annotate!(s[1], s[2], "$k", :black)
        #end
    end
    fig
end
=#


text:
set_data(d::Dict) = false

# macro to simplify concept declaration

const _supported_concept_symbols = Set(
    :data)

function _concept_set(d::Symbol, concepts::Symbol...)
    stmts = Expr[]
    for c in concepts
        if !(c in _supported_concept_symbols)
            error("Invalid concept name: $c")
        end

        fun = symbol(string("set_", string(c)))
        stmt = :( $(fun)(::$(d)) = true )
        push!(stmts, stmt)
    end
    Expr(:block, stmts...)
end

macro concept_set(d, concepts...)
    esc(_concept_set(d, concepts...))
end

# macro to check interface requirements

function _concept_check(d::Symbol, concepts::Symbol...)
    stmts = Expr[]
    for c in concepts
        if !(c in _supported_concept_symbols)
            error("Invalid concept name: $c")
        end

        fun = symbol(string("set_", string(c)))
        msg = "Require: $(c)."
        stmt = :( $(fun)($d) ? nothing : throw(ArgumentError($msg)))
        push!(stmts, stmt)
    end
    Expr(:block, stmts...)
end

macro concept_check(d, concepts...)
    esc(_concept_check(d, concepts...))
end


text:
include("../imports.jl")

@testset ExtendedTestSet "projective" begin
    keypoints = Keypoints([SVector(2., 2)], (32, 32))
    tfm = Project(Translation(20, 10))

    @test_nowarn apply(tfm, keypoints)
    @test boundsranges(apply(tfm, keypoints).bounds) == (21:52, 11:42)

    @test_nowarn apply(tfm, Polygon(keypoints))
end


text:
"""
    struct Attribute
        name::String
    end

Specify an attribute to be returned by the query.

Example :

    BioMart.Attribute("external_gene_name")
"""
struct Attribute
    name::String
end

function Base.push!(query::Query, f::Attribute)
    xa = new_child(query.dataset, "Attribute")
    set_attributes(xa, name = f.name)
end

struct Attributes
    attributes::Vector{Attribute}
end

"""
    Attributes(args...)

Specify a list of attributes to be returned by the query.

Example :

    BioMart.Attributes(
        "external_gene_name",
        "strand"
    )
"""
function Attributes(args...)
    Attributes([Attribute(arg) for arg in args])
end

Base.push!(query::Query, f::Attributes) = push!.(Ref(query), f.attributes)


text:
using Dolfin

import Dolfin:dot, Function, det
import PyCall: PyDict

mesh = UnitSquareMesh(32, 32)

# Define function spaces and mixed (product) space
DRT = FunctionSpace(mesh, "DRT", 2)
CG  = FunctionSpace(mesh, "CG", 3)
W = DRT * CG

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source functions
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5.0*x[0])")

# Define variational form
a = (dot(sigma, tau) + dot(grad(u), tau) + dot(sigma, grad(v)))*dx
L = - f*v*dx - g*v*ds

bc = DirichletBC(W[:sub](1), 0.0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")

# Compute solution
w = Function(W)
solve(a==L, w, bc)
(sigma, u) = w["split"]()

plot(sigma, interactive=true)


text:
@testset "Core" begin
        A = [NaN 1 2 3;
             NaN 6 5 4]
        @test ImageContrastAdjustment.minfinite(A) == 1
        @test ImageContrastAdjustment.maxfinite(A) == 6
        A = rand(10:20, 5, 5)
        @test ImageContrastAdjustment.minfinite(A) == minimum(A)
        @test ImageContrastAdjustment.maxfinite(A) == maximum(A)
        A = reinterpret(N0f8, rand(0x00:0xff, 5, 5))
        @test ImageContrastAdjustment.minfinite(A) == minimum(A)
        @test ImageContrastAdjustment.maxfinite(A) == maximum(A)
        A = rand(Float32,3,5,5)
        img = colorview(RGB, A)
        dc = ImageContrastAdjustment.minfinite(img)-RGB{Float32}(minimum(A, dims=(2,3))...)
        @test abs(dc) < 1e-6
        dc = ImageContrastAdjustment.maxfinite(img)-RGB{Float32}(maximum(A, dims=(2,3))...)
        @test abs(dc) < 1e-6
        @test ImageContrastAdjustment.minfinite(x->x^2,[NaN,10,2]) == 4
        @test ImageContrastAdjustment.maxfinite(x->x^2,[NaN,10,2]) == 100
end


text:
# ---
# title: 1419. Minimum Number of Frogs Croaking
# id: problem1419
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/minimum-number-of-frogs-croaking/description/>
# hidden: true
# ---
#
# Given the string `croakOfFrogs`, which represents a combination of the string
# "croak" from different frogs, that is, multiple frogs can croak at the same
# time, so multiple "croak" are mixed.  _Return the minimum number of_ different
# _frogs to finish all the croak in the given string._
#
# A valid "croak" means a frog is printing 5 letters 'c', 'r', 'o', 'a', 'k'
# **sequentially**. The frogs have to print all five letters to finish a croak.
# If the given string is not a combination of valid "croak" return -1.
#
#
#
# **Example 1:**
#
#
#
#     Input: croakOfFrogs = "croakcroak"
#     Output: 1
#     Explanation: One frog yelling "croak **"** twice.
#
#
# **Example 2:**
#
#
#
#     Input: croakOfFrogs = "crcoakroak"
#     Output: 2
#     Explanation: The minimum number of frogs is two.
#     The first frog could yell " **cr** c **oak** roak".
#     The second frog could yell later "cr **c** oak **roak** ".
#
#
# **Example 3:**
#
#
#
#     Input: croakOfFrogs = "croakcrook"
#     Output: -1
#     Explanation: The given string is an invalid combination of "croak **"** from different frogs.
#
#
# **Example 4:**
#
#
#
#     Input: croakOfFrogs = "croakcroa"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= croakOfFrogs.length <= 10^5`
#   * All characters in the string are: `'c'`, `'r'`, `'o'`, `'a'` or `'k'`.
#
#
## @lc code=start
using LeetCode

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


text:
using JuLIP: neighbourlist,
             @D,
             JVec

import JuLIP.Potentials: evaluate,
                         evaluate_d,
                         evaluate_dd,
                         evaluate_d!

using NeighbourLists: nbodies,
                      maptosites!,
                      maptosites_d!,
                      max_neigs,
                      sites

using KahanSummation

function _get_loop_ex(N)
   # generate the expression for the multiple for loops, e.g. for 4B:
   # for i_1 = 1:(nR-2), i_2 = (i_1+1):(nR-1), i_3 = (i_2+1):nR
   str_loop = "for i_1 = 1:(nR-$(N-2))"
   for n = 2:N-1
      str_loop *= ", i_$n = (i_$(n-1)+1):(nR-$(N-1-n))"
   end
   str_loop *= "\n end"
   return Meta.parse(str_loop)
end

function _get_Jvec_ex(N)
   # inside these N-1 loops we collect the loop indices into an SVector, e.g.
   # J = @SVector [i_1, i_2, i_3]
   str = "J = @SVector Int[i_1"
   for n = 2:N-1
      str *= ", i_$n"
   end
   return Meta.parse(str * "]")
end

skipunorderedsimplices(V::NBodyFunction) =
      skipunorderedsimplices(descriptor(V))

skipunorderedsimplices(D::NBSiteDescriptor) = false

skipunorderedsimplices(D::NBClusterDescriptor) = true

function unordered(i, j, J)
   for jj in J
      i > j[jj] && return true
   end
   return false
end

@generated function eval_site_nbody!( ::Val{N},
                                      ii::Int, jj::AbstractVector{TI},
                                      Rs::AbstractVector{JVec{T}},
                                      rcut::T,
                                      skipunordered::Bool,
                                      reducefun,
                                      out,
                                      temp ) where {N, T, TI <: Integer}
   code = Expr[]
   # initialise the output
   push!(code, :( nR = length(Rs)  ))
   push!(code, :( skip = false ))

   # generate the multi-for-loop
   ex_loop = _get_loop_ex(N)

   # inside the loop
   # ---------------
   code_inner = Expr[]
   # collect the indices into a vector
   push!(code_inner,      _get_Jvec_ex(N) )
   push!(code_inner, :(   if skipunordered; skip = unordered(ii, jj, J); end   ))
   # now call `V` with the simplex-corner vectors and "add" this to the site energy
   push!(code_inner, :(   if !skip; out = reducefun(out, Rs, ii, J, temp); end ))

   # put code_inner into the loop expression
   ex_loop.args[2] = Expr(:block, code_inner...)

   # now append the loop to the main code
   push!(code, ex_loop)

   quote
      @inbounds $(Expr(:block, code...))
      return out
   end
end

function site_energies(V::NBodyFunction{N, DT}, at::Atoms{T}
                  ) where {N, T, DT <: NBSiteDescriptor}
   Es = zeros(T, length(at))
   for (i, j, r, R) in sites(at, cutoff(V))
      Es[i] = eval_site_nbody!(Val(N), i, j, R, cutoff(V),
                               skipunorderedsimplices(V),
                               ((out, R, ii, J, temp) -> out + evaluate(V, R, ii, J)),
                               zero(T), nothing)
   end
   return Es
end

# this is probably already in JuLIP??? if not, it should be moved to JuLIP??
energy(V::NBodyFunction, at::Atoms) = sum_kbn(site_energies(V, at))

# this appears to be a nice generic implementation of forces with a
# temporary array => move this to JuLIP!
function forces(V::NBodyFunction{N}, at::Atoms{T}) where {N, T}
   nlist = neighbourlist(at, cutoff(V))
   maxneigs = max_neigs(nlist)
   F = zeros(JVec{T}, length(at))
   dVsite = zeros(JVec{T}, maxneigs)
   for (i, j, r, R) in sites(nlist)
      fill!(dVsite, zero(JVec{T}))
      eval_site_nbody!(
            Val(N), i, j, R, cutoff(V), skipunorderedsimplices(V),
            (out, R, ii, J, temp) -> evaluate_d!(out, V, R, ii, J),
            dVsite, nothing )   # dVsite == out, nothing == temp
      # write site energy gradient into forces
      for n = 1:length(j)
         F[j[n]] -= dVsite[n]
         F[i] += dVsite[n]
      end
   end
   return F
end

function virial(V::NBodyFunction{N}, at::Atoms{T}) where {N, T}
   nlist = neighbourlist(at, cutoff(V))
   maxneigs = max_neigs(nlist)
   S = @SMatrix zeros(3,3)
   dVsite = zeros(JVec{T}, maxneigs)
   for (i, j, r, R) in sites(nlist)
      dVsite .*= 0.0
      # eval_site_d!(dVsite, V, R)
      eval_site_nbody!(
            Val(N), i, j, R, cutoff(V), skipunorderedsimplices(V),
            (out, R, ii, J, temp) -> evaluate_d!(out, V, R, ii, J),
            dVsite, nothing )
      S += JuLIP.Potentials.site_virial(dVsite, R)
   end
   return S
end

# ========================= assembly support for LSQ system ====================

# For assembling the LSQ system efficiently we need a way to evaluate all basis
# functions of the same body-order at the same time. Otherwise we would be
# re-computing the invariants many many times, which is very expensive.
# To achieve this we just wrap all basis functions of a body-order into
# a new type `NBBasis` which evaluates to a long vector
#
# at the moment, it seems we need to hard-code this to the Polys
# sub-module, but it would be good if this can be fixed, so we keep this
# "interface" here.

function energy(B::AbstractVector{TB}, at::Atoms{T}
                ) where {TB <: NBodyFunction{N}, T} where {N}
   # TODO: assert that all B[j] have the same invariants
   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   E = zeros(T, length(B))
   for (i, j, r, R) in sites(nlist)
      # evaluate all the site energies at the same time
      # for each simplex, write the nB energies into temp
      # then add them to E, which is just passed through all the
      # various loops, so no need to update it here again
      eval_site_nbody!(Val(N), i, j, R, rcut, skipunorderedsimplices(B[1]),
                       (out, R, ii, J, temp) -> evaluate_many!(out, B, R, ii, J),
                       E, nothing)
   end
   return E
end

function forces(B::AbstractVector{TB}, at::Atoms{T}
              ) where {TB <: NBodyFunction{N}, T} where {N}
   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   maxneigs = max_neigs(nlist)
   nB = length(B)
   # forces
   F =      [ zeros(JVec{T}, length(at)) for n = 1:nB ]
   # site gradient
   dVsite = [ zeros(JVec{T}, maxneigs)   for n = 1:nB ]

   for (i, j, r, R) in sites(nlist)
      # clear dVsite
      for n = 1:nB; fill!(dVsite[n], zero(JVec{T})); end
      # fill dVsite
      eval_site_nbody!(Val(N), i, j, R, rcut, skipunorderedsimplices(B[1]),
                       (out, R, ii, J, temp) -> evaluate_many_d!(out, B, R, ii, J),
                       dVsite, nothing)
      # write it into the force vectors
      for ib = 1:nB, n = 1:length(j)
         F[ib][j[n]] -= dVsite[ib][n]
         F[ib][i] += dVsite[ib][n]
      end
   end
   return F
end

function virial(B::AbstractVector{TB}, at::Atoms{T}
                ) where {TB <: NBodyFunction{N}, T} where {N}
   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   maxneigs = max_neigs(nlist)
   nB = length(B)
   # virials (main output)
   S = fill((@SMatrix zeros(3,3)), nB)
   # site gradient
   dVsite = [ zeros(JVec{T}, maxneigs)   for n = 1:nB ]

   for (i, j, r, R) in sites(nlist)
      # clear dVsite
      for n = 1:nB; dVsite[n] .*= 0.0; end
      # fill dVsite
      eval_site_nbody!(Val(N), i, j, R, rcut, skipunorderedsimplices(B[1]),
                       (out, R, ii, J, temp) -> evaluate_many_d!(out, B, R, ii, J),
                       dVsite, nothing)
      # update the virials
      for iB = 1:nB
         S[iB] += JuLIP.Potentials.site_virial(dVsite[iB], R)
      end
   end
   return S
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -5.483106975000851 <= q <= 6.193029789058127)

                     add_NL_constraint(m, :(softplus(0.7547513043737739 + -0.48113681936364117*softplus(-0.10262045764995253 + -0.12743000480166566*softplus(0.5778005478095731 + 0.8852363638482865*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*softplus(0.07154777431849713 + -0.5862593464659436*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.9820460002458025*softplus(0.1263157800600383 + 0.1285425867505281*softplus(0.5778005478095731 + 0.8852363638482865*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*softplus(0.07154777431849713 + -0.5862593464659436*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) + softplus(-0.6861294460251406 + 0.36167652828680286*softplus(-0.10262045764995253 + -0.12743000480166566*softplus(0.5778005478095731 + 0.8852363638482865*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*softplus(0.07154777431849713 + -0.5862593464659436*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.5440699588803759*softplus(0.1263157800600383 + 0.1285425867505281*softplus(0.5778005478095731 + 0.8852363638482865*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*softplus(0.07154777431849713 + -0.5862593464659436*softplus(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*softplus(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
@axpy
@ascal
@dot

BLAS Level 1 functions.
"""

macro axpy(a, x, y)
    expr = quote
        let u = $(esc(a))
            for i in eachindex($(esc(x)))
                @inbounds $(esc(y))[i] += u * $(esc(x))[i]
            end
        end
    end
    expr
end

macro scal(a, x)
    expr = quote
        let u = $(esc(a))
            for i in eachindex($(esc(x)))
                @inbounds $(esc(x))[i] *= u
            end
        end
    end
    expr
end

macro dot(x, y)
    expr = quote
        s = 0
        for i in eachindex($(esc(x)))
            @inbounds s += $(esc(x))[i] * $(esc(y))[i]
        end
        s
    end
    expr
end

"""
itime(t)

Get interval time from a given cumulative time vector t.
The first element is t[1]

Retuen value:
dt: interval time vector
maxt: maximum interval time
"""

function itime(t::AbstractVector{Tv}) where Tv
    dt = similar(t)
    prev = Tv(0)
    maxt = Tv(0)
    @inbounds for i = eachindex(t)
        dt[i] = t[i] - prev
        prev = t[i]
        if dt[i] > maxt
            maxt = dt[i]
        end
    end
    return dt, maxt
end


text:
module EveInvent

import Requests
import HttpServer
import URIParser
import JSON
import LightXML
import DataFrames

include("crest.jl")
include("crestauth.jl")
include("apiv2.jl")
include("evec.jl")
include("dataimp.jl")

global config
global services, regions, itemCategories, itemTypes, itemNames
global itemPrices = nothing
global marketRegion, marketGroups, marketPrices, industrySystems, industryFacilities, systemFacilities
global shipTypes, blueprintTypes, blueprints, inventableProducts, productionProducts, shipInventable
global decryptors, skills

global charSkills, charBlueprints, charBPOTypes

global crestAuth = nothing
global apiAuth

global dataDir = joinpath(dirname(dirname(@__FILE__)), "deps")
global cfgDir = joinpath(dataDir, "cfg")
global cacheDir = joinpath(dirname(dirname(@__FILE__)), "cache")

const priceTimeout = 6.0

id_from_url(url::AbstractString) = parse(Int, rsplit(url, '/'; limit = 2, keep = false)[end])

function access_path(o, p::Array, default = nothing)
	for i in p
		if !(isa(o, Associative) && haskey(o, i) || isa(o, AbstractArray) && 1 <= i <= length(o))
			return default
		end
		o = o[i]
	end
	return o
end

function server_version(serv)
	strVer = match(r"[0-9|.?]+", serv["serverVersion"])
	map(v->parse(Int, v), split(strVer.match, '.'))
end

function get_services()
	oldServ = read_from_cache(urlCrest, Inf64)
	newServ = get_crest(urlCrest, crestAuth, 0.0)
	@assert newServ != nothing
	if oldServ == nothing || server_version(oldServ)[1:2] != server_version(newServ)[1:2]
		info("Server version differs from cached data, clearing cache")
		clear_cache(".json")
		write_to_cache(urlCrest, newServ)
	end
	return newServ
end

get_service(path::Array) = access_path(services, path)["href"]

make_row_dict(colNames, vals) = Dict(filter(nv->!DataFrames.isna(nv[2]), zip(colNames, vals)))

function make_array(df::DataFrames.DataFrame)
	colNames = map(string, names(df))
	return [make_row_dict(colNames, vals) for vals in zip(DataFrames.columns(df)...)]
end

function make_index(key::Symbol, df::DataFrames.DataFrame)
	colNames = map(string, names(df))
	keyIndex = indexin([key], names(df))[1]
	keyType = eltype(df[key])
	return Dict(vals[keyIndex]::keyType => make_row_dict(colNames, vals) for vals in zip(DataFrames.columns(df)...))
end

function make_index(kvFunc::Function, collection, target)
	for e in collection
		k,v = kvFunc(e)
		target[k] = v
	end
	return target
end

function get_regions()
	jsonRegions = get_crest(get_service(["regions"]), crestAuth)
	make_index(jsonRegions, Dict{Int, Any}()) do r
		id_from_url(r["href"]), r
	end
end

get_region(id::Int) = get_crest(regions[id]["href"], crestAuth)

get_item_types() = make_index(t->(id_from_url(t["href"]), t), get_crest(get_service(["itemTypes"]), crestAuth), Dict{Int, Any}())
get_item_names() = Dict(t["name"]=>id for (id, t) in itemTypes)

global nextFakeId = 9999000
function add_fake_item(name::AbstractString)
	global nextFakeId
	assert(!haskey(itemTypes, nextFakeId))
	itemTypes[nextFakeId] = Dict("name"=>name, "id"=>nextFakeId)
	nextFakeId += 1
	return nextFakeId - 1
end

fake_item(id::Int) = !haskey(itemTypes[id], "href")

function get_item_orders(typeID::Int, sellOrders::Bool)
	if fake_item(typeID)
		return [Dict("price"=>0.0, "id"=>typeID, "buy"=>!sellOrders, "volume"=>0)]
	end
	get_crest(marketRegion[sellOrders? "marketSellOrders" : "marketBuyOrders"]["href"] * "?type=" * itemTypes[typeID]["href"], crestAuth, priceTimeout)
end

function resolve_item_prices(types)
	global itemPrices
	if itemPrices == nothing
		itemPrices = Dict{Int, Any}()
	end
	missing = filter(t->!haskey(itemPrices, t), types)
	for typeID in missing
		sell = get_item_orders(typeID, true)
		sellPrice = reduce(min, Inf64, map(o->o["price"], sell))
		if sellPrice == Inf64
			info("Warning: Cannot find price for typeID $typeID - $(itemTypes[typeID]["name"])")
		end
		# buy = get_item_orders(typeID, false)
		# buyPrice = reduce(max, 0.0, map(o->o["price"], buy))
		itemPrices[typeID] = Dict(#="buy"=>buyPrice,=# "sell"=>sellPrice)
	end
end

function resolve_item_prices_evec(types)
	global itemPrices
	itemPricesFile = joinpath(cacheDir, "itemPrices.json")
	if itemPrices == nothing
		itemPrices = json_read(itemPricesFile)
		if itemPrices == nothing
			itemPrices = Dict{Int, Any}()
		end
	end
	timeNow = time()
	missing = filter(t->!haskey(itemPrices, t) || itemPrices[t]["timeStamp"]<=timeNow, types)
	if isempty(missing)
		return
	end
	timeStamp = timeNow + priceTimeout * 60 * 60
	res = get_evec("marketstat", query=Dict{Any, Any}("typeid"=>join(missing, ','), "regionlimit"=>config["marketRegionID"]); timeoutHours=0.0)
	for typeID in missing
		buy = xml_find(res, "marketstat/id:$typeID/buy/max")[1, 1]
		sell = xml_find(res, "marketstat/id:$typeID/sell/min")[1, 1]
		itemPrices[typeID] = Dict("buy"=>buy, "sell"=>sell, "timeStamp"=>timeStamp)
	end
	json_write(itemPricesFile, itemPrices)
end

get_item_categories() = get_crest(get_service(["itemCategories"]), crestAuth)
get_item_category(catName::AbstractString) = get_crest(filter(c->c["name"]==catName, itemCategories)[1]["href"], crestAuth)

function get_category_types(category)
	types = Dict{Int, Any}()
	if category == nothing
		return types
	end
	if haskey(category, "types")
		make_index(category["types"], types) do t
		    id_from_url(t["href"]), t["name"]
		end
	end
	if haskey(category, "groups")
		for group in category["groups"]
			subtypes = get_category_types(get_crest(group["href"], crestAuth))
			merge!(types, subtypes)
		end
	end
	return types
end

get_category_types(catName::AbstractString) = get_category_types(get_item_category(catName))

get_market_groups() = get_crest(get_service(["marketGroups"]), crestAuth)
get_market_group(group::AbstractString) = filter(g->g["name"]==group, marketGroups)[1]

function get_group_types(group::AbstractString)
	typeGroup = get_market_group(group)
	jTypes = get_crest(typeGroup["types"]["href"], crestAuth)
	make_index(jTypes, Dict{Int, AbstractString}()) do t
		t["type"]["id"], t["type"]["name"]
	end
end

function get_market_prices()
	priceArray = get_crest(get_service(["marketPrices"]), crestAuth, priceTimeout)
	make_index(priceArray, Dict{Int, Any}()) do p
		p["type"]["id"], p
	end
end

function get_industry_facilities()
	facilities = get_crest(get_service(["industry", "facilities"]), crestAuth)
	make_index(facilities, Dict{Int, Any}()) do f
		f["facilityID"], f
	end
end

function get_industry_systems()
	systems = get_crest(get_service(["industry", "systems"]), crestAuth)
	make_index(systems, Dict{Int, Any}()) do s
		s["solarSystem"]["id"], s
	end
end

function get_system_facilities()
	facilities = Dict{Int, Array{Any}}()
	for (facId, fac) in industryFacilities
		facs = get!(facilities, fac["solarSystem"]["id"], Array(Any, 0))
		push!(facs, fac)
	end
	return facilities
end

function get_or_import(fileName::AbstractString, importFunc::Function)
	data = json_read(fileName)
	if data == nothing
		info("Importing $fileName from static data...")
		importFunc(fileName)
		data = json_read(fileName)
	end
	return data
end

function get_decryptors_crest()
	typeGroup = get_market_group("Decryptors")
	jTypes = get_crest(typeGroup["types"]["href"], crestAuth)
	decr = []
	for t in jTypes
		jd = get_crest(t["type"]["href"], crestAuth)
		attr = Dict(a["attribute"]["name"]=>a["value"] for a in jd["dogma"]["attributes"])
		d = Dict("name"=>jd["name"], "id"=>jd["id"], "attributes"=>attr)
		push!(decr, d)
	end
	return decr
end

function get_decryptors()
	decr = get_decryptors_crest()
	noDecrName = "No Decryptor"
	noDecrId = add_fake_item(noDecrName)
	itemTypes[noDecrId] = Dict("name"=>noDecrName, "id"=>noDecrId)
	noDecr = Dict("name"=>noDecrName, "id"=>noDecrId, "attributes"=>Dict("inventionMEModifier"=>0.0, "inventionMaxRunModifier"=>0.0, "inventionTEModifier"=>0.0, "inventionPropabilityMultiplier"=>1.0))
	push!(decr, noDecr)
	return decr
end

get_blueprints() = get_or_import(joinpath(cacheDir, "blueprints.json"), import_blueprints)

get_crest_authorization() = request_access(get_service(["authEndpoint"]), config["appInfo"])
get_api_authorization() = Dict("accessToken"=>crestAuth["access_token"])

function get_config()
	config = json_read(joinpath(cfgDir,"config.json"))
	config["appInfo"] = json_read(joinpath(cfgDir, "appinfo.json"))
	set_api_user_agent(config["appInfo"]["userAgent"])
	return config
end

get_char_skills() = xml_find(get_api("Char/Skills"; auth=apiAuth), "result/:skills/*/*:")
get_char_blueprints() = xml_find(get_api("Char/Blueprints"; auth=apiAuth), "result/:blueprints/*/*:")
get_skills() = xml_find(get_api("Eve/SkillTree"), "result/:skillGroups/*/:skills/*/*:")

function process_market_groups()
	global shipTypes = get_category_types("Ship")
	global blueprintTypes = get_category_types("Blueprint")
end

get_products(blueprint, activity::AbstractString) = access_path(blueprint, ["activities", activity, "products"])

function get_inventable_blueprints(groupTypes::Dict{Int})
	groupPrints = Set{Int}()
	for (b,bv) in blueprints
		invProducts = get_products(bv, "invention")
		if invProducts == nothing
			continue
		end
		for prod in invProducts
			t2BlueprintID = prod["typeID"]
			t2Blueprint = get(blueprints, t2BlueprintID, nothing)
			t2Blueprint == nothing && continue
			t2Products = get_products(t2Blueprint, "manufacturing")
			t2Products == nothing && continue
			t2Product = t2Products[1]["typeID"]
			if haskey(groupTypes, t2Product) && haskey(blueprintTypes, b) # need to check if we're looking at a blueprint because of the t3 inventable research items
				push!(groupPrints, t2Product)
			end
		end
	end
	return groupPrints
end

function process_products(f::Function, blueprint, activity::AbstractString)
	prods = access_path(blueprint, ["activities", activity, "products"])
	if prods != nothing
		for p in prods
			f(p)
		end
	end
end

function process_blueprints()
	global inventableProducts = Dict{Int, Any}()
	global productionProducts = Dict{Int, Any}()
	for (bid, bp) in blueprints
		process_products(bp, "invention") do prod
			inventableProducts[prod["typeID"]] = bp
		end
		process_products(bp, "manufacturing") do prod
			productionProducts[prod["typeID"]] = bp
		end
	end
	global shipInventable = get_inventable_blueprints(shipTypes)

	global charBPOTypes = Dict{Int, Array}()
	for (itemID, item) in charBlueprints
		if item["runs"] >= 0
			continue
		end
		itemArr = get!(charBPOTypes, item["typeID"], Array(Dict, 0))
		push!(itemArr, item)
	end
end

system_id(systemName::AbstractString) = first(filter((id,s)->s["solarSystem"]["name"]==systemName, industrySystems))[1]

function process_industry()
	config["productionSystemID"] = system_id(config["productionSystem"])
	config["inventionSystemID"] = system_id(config["inventionSystem"])
	config["marketRegionID"] = first(filter(f->f["solarSystem"]["id"]==config["productionSystemID"], values(industryFacilities)))["region"]["id"]
	global marketRegion = get_region(config["marketRegionID"])
end

function production_amount(runs::Float64, baseQuantity::Float64, ME::Int, facilityMultiplier::Float64 = 1.0)
	materialModifier = facilityMultiplier * (1 - ME / 100.0)
	return max(runs, ceil(round(runs * baseQuantity * materialModifier, 2)))
end

function production_time(runs::Float64, baseProductionTime::Float64, TE::Int, skills::Array{Int}, facilityMultiplier::Float64 = 1.0)
	timeModifier = facilityMultiplier * (1 - TE / 100.0)
	skillModifier = prod(map(s->(1 - s / 100.0), skills))
	return baseProductionTime * timeModifier * skillModifier * runs
end

adjusted_price(itemID::Int) = marketPrices[itemID]["adjustedPrice"]

function base_job_cost(productID::Int)
	if haskey(blueprints, productID)
		materials = blueprints[productID]["activities"]["manufacturing"]["materials"]
	else
		materials = productionProducts[productID]["activities"]["manufacturing"]["materials"]
	end
	cost = 0.0
	for m in materials
		cost += m["quantity"] * adjusted_price(m["typeID"])
	end
	return cost
end

function production_install_cost(productID::Int, systemCostIndex::Float64, facilityTaxRate::Float64)
	baseCost = base_job_cost(productID)
	return baseCost * systemCostIndex * (1.0 + facilityTaxRate)
end

function invention_copy_install_cost(runs::Float64, blueprintID::Int, systemCostIndex::Float64, facilityTaxRate::Float64)
	baseCost = base_job_cost(blueprintID)
	return baseCost * 0.02 * systemCostIndex * (1.0 + facilityTaxRate) * runs
end

function is_encryption(skillID::Int)
	skillName = skills[skillID]["typeName"]
	return !isempty(search(skillName, "Encryption Methods"))
end

function invention_chance(baseChance::Float64, skills::Array, decryptorModifier::Float64)
	skillModifier = 1.0
	for s in skills
		skillID = s["typeID"]
		if haskey(charSkills, skillID)
			level = charSkills[skillID]["level"]
		else
			level = 0
			info("Missing character skill " * itemTypes[skillID]["name"])
		end
		skillModifier += is_encryption(skillID)? level/40.0 : level/30.0
	end
	return baseChance * skillModifier * decryptorModifier
end

function get_cost_index(system, activityName::AbstractString)
	filter(c->c["activityName"]==activityName, system["systemCostIndices"])[1]["costIndex"]
end

get_tax_rate(system) = access_path(systemFacilities, [system["solarSystem"]["id"], 1, "tax"], 0.1)

function blueprint_materials(runs::Float64, blueprintItem, systemID::Int)
	blueprint = blueprints[blueprintItem["typeID"]]
	ME = blueprintItem["materialEfficiency"]
	productID = blueprint["activities"]["manufacturing"]["products"][1]["typeID"]
	baseMaterials = blueprint["activities"]["manufacturing"]["materials"]
	materials = Dict(m["typeID"]::Int => production_amount(runs, Float64(m["quantity"]), ME) for m in baseMaterials)
	system = industrySystems[systemID]
	installCost = production_install_cost(productID, get_cost_index(system, "Manufacturing"), get_tax_rate(system))
	return materials, installCost
end

function material_cost(materials)
	resolve_item_prices(keys(materials))
	cost = 0.0
	for (matID, matAmount) in materials
		cost += itemPrices[matID]["sell"] * matAmount
	end
	return cost
end

function invention_result(targetID::Int, decryptor, systemID::Int=config["inventionSystemID"])
	baseRuns = 1.0
	baseME = 2.0
	basePE = 4.0
	baseProbability = 1.0
	t2BP = productionProducts[targetID]
	t2BPID = t2BP["blueprintTypeID"]
	t1BP = inventableProducts[t2BPID]
	invention = t1BP["activities"]["invention"]
	product = filter(p->p["typeID"]==t2BPID, invention["products"])[1]
	baseProbability = invention_chance(baseProbability * product["probability"], invention["skills"], decryptor["attributes"]["inventionPropabilityMultiplier"])
	baseME += decryptor["attributes"]["inventionMEModifier"]
	basePE += decryptor["attributes"]["inventionTEModifier"]
	baseRuns += decryptor["attributes"]["inventionMaxRunModifier"]
	resultItem = Dict{AbstractString, Any}()
	resultItem["typeID"] = t2BPID
	resultItem["materialEfficiency"] = round(Int, baseME)
	resultItem["timeEfficiency"] = round(Int, basePE)
	resultItem["probability"] = baseProbability
	resultItem["runs"] = baseRuns
	resultItem["quantity"] = product["quantity"]
	materials = Dict{Int, Float64}()
	for m in invention["materials"]
		materials[m["typeID"]] = m["quantity"]
	end
	materials[decryptor["id"]] = 1
	system = industrySystems[systemID]
	taxRate = get_tax_rate(system)
	copyCost = invention_copy_install_cost(1.0, t1BP["blueprintTypeID"], get_cost_index(system, "Copying"), taxRate)
	installCost = invention_copy_install_cost(1.0, targetID, get_cost_index(system, "Invention"), taxRate)
	materialCost = material_cost(materials)
	attemptsPerCopy = 1 / (baseProbability * baseRuns * product["quantity"])
	return resultItem, (materialCost + installCost + copyCost) * attemptsPerCopy
end

function add_materials!(target, source, multiplier)
	for (matID, matAmount) in source
		target[matID] = get(target, matID, 0.0) + matAmount * multiplier
	end
	return target
end

function get_char_bpo(productID::Int, systemID::Int)
	if !haskey(productionProducts, productID)
		return nothing
	end
	bp = productionProducts[productID]
	bpID = bp["blueprintTypeID"]
	if !haskey(charBPOTypes, bpID)
		return nothing
	end
	localBPOs = filter(charBPOTypes[bpID]) do item
		access_path(industryFacilities, [item["locationID"], "solarSystem", "id"]) == systemID
	end
	if isempty(localBPOs)
		localBPOs = charBPOTypes[bpID]
	end
	return localBPOs[1]
end

function break_materials(runs::Float64, blueprintItem, systemID::Int=config["productionSystemID"])
	materials, installCost = blueprint_materials(runs, blueprintItem, systemID)
	intermediates = Dict{Int, Float64}()
	while true
		bpItem = nothing
		local matID, matAmount
		for (matID, matAmount) in materials
			bpItem = get_char_bpo(matID, systemID)
			if bpItem != nothing
				break
			end
		end
		if bpItem == nothing
			break
		end
		intermediates[matID] = get(intermediates, matID, 0.0) + matAmount
		mats, cost = blueprint_materials(matAmount, bpItem, systemID)
		delete!(materials, matID)
		add_materials!(materials, mats, 1.0)
		installCost += cost
	end
	return materials, intermediates, installCost
end

function optimal_invention(targetID::Int)
	resolve_item_prices([targetID])
	price = itemPrices[targetID]["sell"]
	optimal = (-Inf64, nothing)
	for decr in decryptors
		resultItem, inventCost = invention_result(targetID, decr)
		materials, intermediates, produceInstall = break_materials(1.0, resultItem)
		totalCost = inventCost + material_cost(materials) + produceInstall
		profit = price - totalCost
		if optimal[1] < profit
			optimal = (profit, decr)
		end
	end
	return optimal
end

function optimal_invention(checkBPOs::Bool = true)
	result = []
	for ship in shipInventable
		t2BP = productionProducts[ship]
		if checkBPOs
			t1BP = inventableProducts[t2BP["blueprintTypeID"]]
			if !haskey(charBPOTypes, t1BP["blueprintTypeID"])
				continue
			end
		end
		profit, decr = optimal_invention(ship)
		if profit > 0.0
			isFrigate = any(t2BP["activities"]["manufacturing"]["skills"]) do s
							itemTypes[s["typeID"]]["name"] == "Advanced Small Ship Construction"
						end
			number = isFrigate? 10 : 1
			push!(result, Dict("typeID"=>ship, "decryptor"=>decr, "profit"=>profit*number, "number"=>number))
		end
	end
	sort!(result; lt=(x,y)->x["profit"] < y["profit"], rev=true)
	return result
end

fmt_float(f::Float64) = @sprintf("%.2f", f)
fmt_amount(f::Float64) = @sprintf("%.1f", f)

function fmt_item(typeID::Int; markup::Bool = false)
	name = itemTypes[typeID]["name"]
	if markup
		name = "<loc><url=showinfo:$typeID>$name</url>"
	end
	return name
end

function print_materials(materials, intermediates, installCost::Float64, inventCost::Float64, number::Int=1; io::IO = STDOUT, markup::Bool = false)
	resolve_item_prices(keys(materials))
	println(io, "Materials needed:")
	totalPrice = 0.0
	for (matID, matAmount) in materials
		matAmount *= number
		price = itemPrices[matID]["sell"] * matAmount
		println(io, fmt_item(matID; markup=markup) * " => $(fmt_amount(matAmount)) ($(fmt_float(price)) isk)")
		totalPrice += price
	end
	installCost *= number
	println(io, "Materials $(fmt_float(totalPrice)) + Install $(fmt_float(installCost)) + Invent $(fmt_float(inventCost)) = $(fmt_float(totalPrice+installCost+inventCost)) isk")

	if !isempty(intermediates)
		println(io, "\nIntermediate products:")
		for (matID, matAmount) in intermediates
			matAmount *= number
			println(io, fmt_item(matID; markup=markup) * " => $(fmt_amount(matAmount))")
		end
	end
end

function print_opt_row(i::Int, o; io::IO = STDOUT, markup::Bool = false)
	name = fmt_item(o["typeID"], markup=markup)
	number = o["number"]
	if haskey(o, "profit")
		profit = "=> " * fmt_float(o["profit"])
	else
		profit = ""
	end
	decryptor = fmt_item(o["decryptor"]["id"], markup=markup)
	println(io, "$i. $name x $number $profit using $decryptor")
end

function print_optimal(opt::Array; io::IO = STDOUT, markup::Bool = false)
	for i = 1:length(opt)
		print_opt_row(i, opt[i]; io=io, markup=markup)
	end
end

function print_plan(plan, filename::AbstractString = "")
	io = STDOUT
	if !isempty(filename)
		io = open(filename, "w+")
	end
	markup = io != STDOUT
	println(io, "\nMaterials to manufacture")
	materials = Dict{Int, Float64}()
	intermediates = Dict{Int, Float64}()
	installCost = 0.0
	inventionCost = 0.0
	totalSell = 0.0
	resolve_item_prices([entry["typeID"] for entry in plan])
	for i = 1:length(plan)
		entry = plan[i]
		print_opt_row(i, entry; io=io, markup=markup)
		productID = entry["typeID"]
		decr = entry["decryptor"]
		number = entry["number"]
		invItem, invCost = invention_result(productID, decr)
		mats, intermeds, instCost = break_materials(1.0, invItem)
		inventionCost += invCost * number
		installCost += instCost * number
		totalSell += itemPrices[productID]["sell"] * number
		add_materials!(materials, mats, number)
		add_materials!(intermediates, intermeds, number)
	end
	materialCost = material_cost(materials)
	totalProfit = totalSell - (inventionCost + installCost + materialCost)
	println(io, "Expected profit $(fmt_float(totalProfit))\n")
	print_materials(materials, intermediates, installCost, inventionCost; io=io, markup=markup)
	if io != STDOUT
		seekstart(io)
		clipboard(readstring(io)) # send all output to the system clipboard so it can be readily pasted into Eve
		close(io)
	end
end

function find_optimal(checkBPOs::Bool = true)
	opt = optimal_invention(checkBPOs)
	print_optimal(opt)
	print("\nEnter rows to invent & manufacture: ")
	rows = map(n->parse(Int, n), split(strip(readline()), [' ', ',']; keep = false))
	plan = [opt[rows[i]] for i = 1:length(rows)]
	print_plan(plan)
	print_plan(plan, joinpath(cacheDir, "result.txt"))
end

function add_to_plan(product, decryptor, number, plan = Array(Any, 0))
	if isa(product, AbstractString)
		product = itemNames[product]
	end
	if isa(decryptor, AbstractString)
		decryptor = filter(d->contains(d["name"], decryptor), decryptors)[1]
	end
	push!(plan, Dict{Any, Any}("typeID"=>product, "decryptor"=>decryptor, "number"=>number))
end

function init()
	global config = get_config()
	global services = get_services()
	global crestAuth = get_crest_authorization()
	global apiAuth = get_api_authorization()

	global charSkills = make_index(:typeID, get_char_skills())
	global charBlueprints = make_index(:itemID, get_char_blueprints())
	global skills = make_index(:typeID, get_skills())

	global regions = get_regions()
	global itemTypes = get_item_types()
	global itemNames = get_item_names()
	global itemCategories = get_item_categories()
	global marketGroups = get_market_groups()
	global marketPrices = get_market_prices()
	global industryFacilities = get_industry_facilities()
	global industrySystems = get_industry_systems()
	global systemFacilities = get_system_facilities()
	global blueprints = get_blueprints()
	global decryptors = get_decryptors()

	process_market_groups()
	process_industry()
	process_blueprints()
end

init()
find_optimal()

end