text:
@testset "Max-min driving force analysis" begin

    model = load_model(model_paths["e_coli_core.json"])

    flux_solution = flux_balance_analysis_dict(
        model,
        GLPK.Optimizer;
        modifications = [add_loopless_constraints()],
    )

    sol = max_min_driving_force(
        model,
        reaction_standard_gibbs_free_energies,
        Tulip.Optimizer;
        flux_solution = flux_solution,
        proton_ids = ["h_c", "h_e"],
        water_ids = ["h2o_c", "h2o_e"],
        concentration_ratios = Dict(
            ("atp_c", "adp_c") => 10.0,
            ("nadh_c", "nad_c") => 0.13,
            ("nadph_c", "nadp_c") => 1.3,
        ),
        concentration_lb = 1e-6,
        concentration_ub = 100e-3,
        ignore_reaction_ids = ["H2Ot"],
        modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
    )

    @test isapprox(sol.mmdf, 1.7661155558545698, atol = TEST_TOLERANCE)

    sols = max_min_driving_force_variability(
        model,
        reaction_standard_gibbs_free_energies,
        Tulip.Optimizer;
        bounds = gamma_bounds(0.9),
        flux_solution = flux_solution,
        proton_ids = ["h_c", "h_e"],
        water_ids = ["h2o_c", "h2o_e"],
        concentration_ratios = Dict{Tuple{String,String},Float64}(
            ("atp_c", "adp_c") => 10.0,
            ("nadh_c", "nad_c") => 0.13,
            ("nadph_c", "nadp_c") => 1.3,
        ),
        constant_concentrations = Dict{String,Float64}(
        # "pi_c" => 10e-3
        ),
        concentration_lb = 1e-6,
        concentration_ub = 100e-3,
        ignore_reaction_ids = ["H2Ot"],
        modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
    )

    pyk_idx = first(indexin(["PYK"], reactions(model)))
    @test isapprox(
        sols[pyk_idx, 1].dg_reactions["PYK"],
        -1.5895040002691128;
        atol = TEST_TOLERANCE,
    )
end


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

using Test

@testset "middle" begin
    @test middle(3) === 3.0
    @test middle(2, 3) === 2.5
    let x = ((realmax(1.0)/4)*3)
        @test middle(x, x) === x
    end
    @test middle(1:8) === 4.5
    @test middle([1:8;]) === 4.5

    # ensure type-correctness
    for T in [Bool,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,Float16,Float32,Float64]
        @test middle(one(T)) === middle(one(T), one(T))
    end
end

@testset "median" begin
    @test median([1.]) === 1.
    @test median([1.,3]) === 2.
    @test median([1.,3,2]) === 2.

    @test median([1,3,2]) === 2.0
    @test median([1,3,2,4]) === 2.5

    @test median([0.0,Inf]) == Inf
    @test median([0.0,-Inf]) == -Inf
    @test median([0.,Inf,-Inf]) == 0.0
    @test median([1.,-1.,Inf,-Inf]) == 0.0
    @test isnan(median([-Inf,Inf]))

    X = [2 3 1 -1; 7 4 5 -4]
    @test all(median(X, 2) .== [1.5, 4.5])
    @test all(median(X, 1) .== [4.5 3.5 3.0 -2.5])
    @test X == [2 3 1 -1; 7 4 5 -4] # issue #17153

    @test_throws ArgumentError median([])
    @test isnan(median([NaN]))
    @test isnan(median([0.0,NaN]))
    @test isnan(median([NaN,0.0]))
    @test isequal(median([NaN 0.0; 1.2 4.5], 2), reshape([NaN; 2.85], 2, 1))

    @test median!([1 2 3 4]) == 2.5
    @test median!([1 2; 3 4]) == 2.5

    @test invoke(median, Tuple{AbstractVector}, 1:10) == median(1:10) == 5.5
end

@testset "mean" begin
    @test_throws ArgumentError mean(())
    @test mean((1,2,3)) === 2.
    @test mean([0]) === 0.
    @test mean([1.]) === 1.
    @test mean([1.,3]) == 2.
    @test mean([1,2,3]) == 2.
    @test mean([0 1 2; 4 5 6], 1) == [2.  3.  4.]
    @test mean([1 2 3; 4 5 6], 1) == [2.5 3.5 4.5]
    @test mean(i->i+1, 0:2) === 2.
    @test mean(isodd, [3]) === 1.
    @test mean(x->3x, (1,1)) === 3.

    @test isnan(mean([NaN]))
    @test isnan(mean([0.0,NaN]))
    @test isnan(mean([NaN,0.0]))

    @test isnan(mean([0.,Inf,-Inf]))
    @test isnan(mean([1.,-1.,Inf,-Inf]))
    @test isnan(mean([-Inf,Inf]))
    @test isequal(mean([NaN 0.0; 1.2 4.5], 2), reshape([NaN; 2.85], 2, 1))
end

@testset "var & std" begin
    # edge case: empty vector
    # iterable; this has to throw for type stability
    @test_throws ArgumentError var(())
    @test_throws ArgumentError var((); corrected=false)
    @test_throws ArgumentError var((); mean=2)
    @test_throws ArgumentError var((); mean=2, corrected=false)
    # reduction
    @test isnan(var(Int[]))
    @test isnan(var(Int[]; corrected=false))
    @test isnan(var(Int[]; mean=2))
    @test isnan(var(Int[]; mean=2, corrected=false))
    # reduction across dimensions
    @test isequal(var(Int[], 1), [NaN])
    @test isequal(var(Int[], 1; corrected=false), [NaN])
    @test isequal(var(Int[], 1; mean=[2]), [NaN])
    @test isequal(var(Int[], 1; mean=[2], corrected=false), [NaN])

    # edge case: one-element vector
    # iterable
    @test isnan(@inferred(var((1,))))
    @test var((1,); corrected=false) === 0.0
    @test var((1,); mean=2) === Inf
    @test var((1,); mean=2, corrected=false) === 1.0
    # reduction
    @test isnan(@inferred(var([1])))
    @test var([1]; corrected=false) === 0.0
    @test var([1]; mean=2) === Inf
    @test var([1]; mean=2, corrected=false) === 1.0
    # reduction across dimensions
    @test isequal(@inferred(var([1], 1)), [NaN])
    @test var([1], 1; corrected=false) ≈ [0.0]
    @test var([1], 1; mean=[2]) ≈ [Inf]
    @test var([1], 1; mean=[2], corrected=false) ≈ [1.0]

    @test var(1:8) == 6.
    @test varm(1:8,1) == varm(collect(1:8),1)
    @test isnan(varm(1:1,1))
    @test isnan(var(1:1))
    @test isnan(var(1:-1))

    @test @inferred(var(1.0:8.0)) == 6.
    @test varm(1.0:8.0,1.0) == varm(collect(1.0:8.0),1)
    @test isnan(varm(1.0:1.0,1.0))
    @test isnan(var(1.0:1.0))
    @test isnan(var(1.0:-1.0))

    @test @inferred(var(1.0f0:8.0f0)) === 6.f0
    @test varm(1.0f0:8.0f0,1.0f0) == varm(collect(1.0f0:8.0f0),1)
    @test isnan(varm(1.0f0:1.0f0,1.0f0))
    @test isnan(var(1.0f0:1.0f0))
    @test isnan(var(1.0f0:-1.0f0))

    @test varm([1,2,3], 2) ≈ 1.
    @test var([1,2,3]) ≈ 1.
    @test var([1,2,3]; corrected=false) ≈ 2.0/3
    @test var([1,2,3]; mean=0) ≈ 7.
    @test var([1,2,3]; mean=0, corrected=false) ≈ 14.0/3

    @test varm((1,2,3), 2) ≈ 1.
    @test var((1,2,3)) ≈ 1.
    @test var((1,2,3); corrected=false) ≈ 2.0/3
    @test var((1,2,3); mean=0) ≈ 7.
    @test var((1,2,3); mean=0, corrected=false) ≈ 14.0/3
    @test_throws ArgumentError var((1,2,3); mean=())

    @test var([1 2 3 4 5; 6 7 8 9 10], 2) ≈ adjoint([2.5 2.5])
    @test var([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) ≈ adjoint([2.0 2.0])

    @test stdm([1,2,3], 2) ≈ 1.
    @test std([1,2,3]) ≈ 1.
    @test std([1,2,3]; corrected=false) ≈ sqrt(2.0/3)
    @test std([1,2,3]; mean=0) ≈ sqrt(7.0)
    @test std([1,2,3]; mean=0, corrected=false) ≈ sqrt(14.0/3)

    @test stdm((1,2,3), 2) ≈ 1.
    @test std((1,2,3)) ≈ 1.
    @test std((1,2,3); corrected=false) ≈ sqrt(2.0/3)
    @test std((1,2,3); mean=0) ≈ sqrt(7.0)
    @test std((1,2,3); mean=0, corrected=false) ≈ sqrt(14.0/3)

    @test std([1 2 3 4 5; 6 7 8 9 10], 2) ≈ sqrt.(adjoint([2.5 2.5]))
    @test std([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) ≈ sqrt.(adjoint([2.0 2.0]))

    let A = ComplexF64[exp(i*im) for i in 1:10^4]
        @test varm(A, 0.) ≈ sum(map(abs2, A)) / (length(A) - 1)
        @test varm(A, mean(A)) ≈ var(A)
    end
end

function safe_cov(x, y, zm::Bool, cr::Bool)
    n = length(x)
    if !zm
        x = x .- mean(x)
        y = y .- mean(y)
    end
    dot(vec(x), vec(y)) / (n - Int(cr))
end
X = adjoint([1. 2. 3. 4. 5.; 5. 4. 6. 2. 1.])
Y = adjoint([6. 1. 5. 3. 2.; 2. 7. 8. 4. 3.])

@testset "covariance" begin
    for vd in [1, 2], zm in [true, false], cr in [true, false]
        # println("vd = $vd: zm = $zm, cr = $cr")
        if vd == 1
            k = size(X, 2)
            Cxx = zeros(k, k)
            Cxy = zeros(k, k)
            for i = 1:k, j = 1:k
                Cxx[i,j] = safe_cov(X[:,i], X[:,j], zm, cr)
                Cxy[i,j] = safe_cov(X[:,i], Y[:,j], zm, cr)
            end
            x1 = vec(X[:,1])
            y1 = vec(Y[:,1])
        else
            k = size(X, 1)
            Cxx = zeros(k, k)
            Cxy = zeros(k, k)
            for i = 1:k, j = 1:k
                Cxx[i,j] = safe_cov(X[i,:], X[j,:], zm, cr)
                Cxy[i,j] = safe_cov(X[i,:], Y[j,:], zm, cr)
            end
            x1 = vec(X[1,:])
            y1 = vec(Y[1,:])
        end

        c = zm ? Base.covm(x1, 0, corrected=cr) :
                 cov(x1, corrected=cr)
        @test isa(c, Float64)
        @test c ≈ Cxx[1,1]
        @inferred cov(x1, corrected=cr)

        @test cov(X) == Base.covm(X, mean(X, 1))
        C = zm ? Base.covm(X, 0, vd, corrected=cr) :
                 cov(X, vd, corrected=cr)
        @test size(C) == (k, k)
        @test C ≈ Cxx
        @inferred cov(X, vd, corrected=cr)

        @test cov(x1, y1) == Base.covm(x1, mean(x1), y1, mean(y1))
        c = zm ? Base.covm(x1, 0, y1, 0, corrected=cr) :
                 cov(x1, y1, corrected=cr)
        @test isa(c, Float64)
        @test c ≈ Cxy[1,1]
        @inferred cov(x1, y1, corrected=cr)

        if vd == 1
            @test cov(x1, Y) == Base.covm(x1, mean(x1), Y, mean(Y, 1))
        end
        C = zm ? Base.covm(x1, 0, Y, 0, vd, corrected=cr) :
                 cov(x1, Y, vd, corrected=cr)
        @test size(C) == (1, k)
        @test vec(C) ≈ Cxy[1,:]
        @inferred cov(x1, Y, vd, corrected=cr)

        if vd == 1
            @test cov(X, y1) == Base.covm(X, mean(X, 1), y1, mean(y1))
        end
        C = zm ? Base.covm(X, 0, y1, 0, vd, corrected=cr) :
                 cov(X, y1, vd, corrected=cr)
        @test size(C) == (k, 1)
        @test vec(C) ≈ Cxy[:,1]
        @inferred cov(X, y1, vd, corrected=cr)

        @test cov(X, Y) == Base.covm(X, mean(X, 1), Y, mean(Y, 1))
        C = zm ? Base.covm(X, 0, Y, 0, vd, corrected=cr) :
                 cov(X, Y, vd, corrected=cr)
        @test size(C) == (k, k)
        @test C ≈ Cxy
        @inferred cov(X, Y, vd, corrected=cr)
    end
end

function safe_cor(x, y, zm::Bool)
    if !zm
        x = x .- mean(x)
        y = y .- mean(y)
    end
    x = vec(x)
    y = vec(y)
    dot(x, y) / (sqrt(dot(x, x)) * sqrt(dot(y, y)))
end
@testset "correlation" begin
    for vd in [1, 2], zm in [true, false]
        # println("vd = $vd: zm = $zm")
        if vd == 1
            k = size(X, 2)
            Cxx = zeros(k, k)
            Cxy = zeros(k, k)
            for i = 1:k, j = 1:k
                Cxx[i,j] = safe_cor(X[:,i], X[:,j], zm)
                Cxy[i,j] = safe_cor(X[:,i], Y[:,j], zm)
            end
            x1 = vec(X[:,1])
            y1 = vec(Y[:,1])
        else
            k = size(X, 1)
            Cxx = zeros(k, k)
            Cxy = zeros(k, k)
            for i = 1:k, j = 1:k
                Cxx[i,j] = safe_cor(X[i,:], X[j,:], zm)
                Cxy[i,j] = safe_cor(X[i,:], Y[j,:], zm)
            end
            x1 = vec(X[1,:])
            y1 = vec(Y[1,:])
        end

        c = zm ? Base.corm(x1, 0) : cor(x1)
        @test isa(c, Float64)
        @test c ≈ Cxx[1,1]
        @inferred cor(x1)

        @test cor(X) == Base.corm(X, mean(X, 1))
        C = zm ? Base.corm(X, 0, vd) : cor(X, vd)
        @test size(C) == (k, k)
        @test C ≈ Cxx
        @inferred cor(X, vd)

        @test cor(x1, y1) == Base.corm(x1, mean(x1), y1, mean(y1))
        c = zm ? Base.corm(x1, 0, y1, 0) : cor(x1, y1)
        @test isa(c, Float64)
        @test c ≈ Cxy[1,1]
        @inferred cor(x1, y1)

        if vd == 1
            @test cor(x1, Y) == Base.corm(x1, mean(x1), Y, mean(Y, 1))
        end
        C = zm ? Base.corm(x1, 0, Y, 0, vd) : cor(x1, Y, vd)
        @test size(C) == (1, k)
        @test vec(C) ≈ Cxy[1,:]
        @inferred cor(x1, Y, vd)

        if vd == 1
            @test cor(X, y1) == Base.corm(X, mean(X, 1), y1, mean(y1))
        end
        C = zm ? Base.corm(X, 0, y1, 0, vd) : cor(X, y1, vd)
        @test size(C) == (k, 1)
        @test vec(C) ≈ Cxy[:,1]
        @inferred cor(X, y1, vd)

        @test cor(X, Y) == Base.corm(X, mean(X, 1), Y, mean(Y, 1))
        C = zm ? Base.corm(X, 0, Y, 0, vd) : cor(X, Y, vd)
        @test size(C) == (k, k)
        @test C ≈ Cxy
        @inferred cor(X, Y, vd)
    end

    @test cor(repmat(1:17, 1, 17))[2] <= 1.0
    @test cor(1:17, 1:17) <= 1.0
    @test cor(1:17, 18:34) <= 1.0
    let tmp = linspace(1, 85, 100)
        tmp2 = collect(tmp)
        @test cor(tmp, tmp) <= 1.0
        @test cor(tmp, tmp2) <= 1.0
    end
end

@testset "quantile" begin
    @test quantile([1,2,3,4],0.5) == 2.5
    @test quantile([1,2,3,4],[0.5]) == [2.5]
    @test quantile([1., 3],[.25,.5,.75])[2] == median([1., 3])
    @test quantile(100.0:-1.0:0.0, 0.0:0.1:1.0) == collect(0.0:10.0:100.0)
    @test quantile(0.0:100.0, 0.0:0.1:1.0, sorted=true) == collect(0.0:10.0:100.0)
    @test quantile(100f0:-1f0:0.0, 0.0:0.1:1.0) == collect(0f0:10f0:100f0)
    @test quantile([Inf,Inf],0.5) == Inf
    @test quantile([-Inf,1],0.5) == -Inf
    @test quantile([0,1],1e-18) == 1e-18
    @test quantile([1, 2, 3, 4],[]) == []
    @test quantile([1, 2, 3, 4], (0.5,)) == (2.5,)
    @test quantile([4, 9, 1, 5, 7, 8, 2, 3, 5, 17, 11], (0.1, 0.2, 0.4, 0.9)) == (2.0, 3.0, 5.0, 11.0)
    @test quantile([1, 2, 3, 4], ()) == ()
end

# StatsBase issue 164
let y = [0.40003674665581906, 0.4085630862624367, 0.41662034698690303, 0.41662034698690303, 0.42189053966652057, 0.42189053966652057, 0.42553514344518345, 0.43985732442991354]
    @test issorted(quantile(y, linspace(0.01, 0.99, 17)))
end

@testset "variance of complex arrays (#13309)" begin
    z = rand(ComplexF64, 10)
    @test var(z) ≈ invoke(var, Tuple{Any}, z) ≈ cov(z) ≈ var(z,1)[1] ≈ sum(abs2, z .- mean(z))/9
    @test isa(var(z), Float64)
    @test isa(invoke(var, Tuple{Any}, z), Float64)
    @test isa(cov(z), Float64)
    @test isa(var(z,1), Vector{Float64})
    @test varm(z, 0.0) ≈ invoke(varm, Tuple{Any,Float64}, z, 0.0) ≈ sum(abs2, z)/9
    @test isa(varm(z, 0.0), Float64)
    @test isa(invoke(varm, Tuple{Any,Float64}, z, 0.0), Float64)
    @test cor(z) === 1.0
    v = varm([1.0+2.0im], 0; corrected = false)
    @test v ≈ 5
    @test isa(v, Float64)
end

@testset "cov and cor of complex arrays (issue #21093)" begin
    x = [2.7 - 3.3im, 0.9 + 5.4im, 0.1 + 0.2im, -1.7 - 5.8im, 1.1 + 1.9im]
    y = [-1.7 - 1.6im, -0.2 + 6.5im, 0.8 - 10.0im, 9.1 - 3.4im, 2.7 - 5.5im]
    @test cov(x, y) ≈ 4.8365 - 12.119im
    @test cov(y, x) ≈ 4.8365 + 12.119im
    @test cov(x, reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1)
    @test cov(reshape(x, :, 1), y) ≈ reshape([4.8365 - 12.119im], 1, 1)
    @test cov(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1)
    @test cov([x y]) ≈ [21.779 4.8365-12.119im;
                        4.8365+12.119im 54.548]
    @test cor(x, y) ≈ 0.14032104449218274 - 0.35160772008699703im
    @test cor(y, x) ≈ 0.14032104449218274 + 0.35160772008699703im
    @test cor(x, reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
    @test cor(reshape(x, :, 1), y) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
    @test cor(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
    @test cor([x y]) ≈ [1.0                                          0.14032104449218274-0.35160772008699703im
                        0.14032104449218274+0.35160772008699703im  1.0]
end

@testset "Issue #17153 and PR #17154" begin
    a = rand(10,10)
    b = deepcopy(a)
    x = median(a, 1)

    @test b == a
    x = median(a, 2)
    @test b == a
    x = mean(a, 1)
    @test b == a
    x = mean(a, 2)
    @test b == a
    x = var(a, 1)
    @test b == a
    x = var(a, 2)
    @test b == a
    x = std(a, 1)
    @test b == a
    x = std(a, 2)
    @test b == a
end

# dimensional correctness
isdefined(Main, :TestHelpers) || @eval Main include("TestHelpers.jl")
using Main.TestHelpers: Furlong
@testset "Unitful elements" begin
    r = Furlong(1):Furlong(1):Furlong(2)
    a = collect(r)
    @test sum(r) == sum(a) == Furlong(3)
    @test cumsum(r) == Furlong.([1,3])
    @test mean(r) == mean(a) == median(a) == median(r) == Furlong(1.5)
    @test var(r) == var(a) == Furlong{2}(0.5)
    @test std(r) == std(a) == Furlong{1}(sqrt(0.5))

    # Issue #21786
    A = [Furlong{1}(rand(-5:5)) for i in 1:2, j in 1:2]
    @test mean(mean(A, 1), 2)[1] === mean(A)
    @test var(A, 1)[1] === var(A[:, 1])
    @test_broken std(A, 1)[1] === std(A[:, 1])
end

# Issue #22901
@testset "var and quantile of Any arrays" begin
    x = Any[1, 2, 4, 10]
    y = Any[1, 2, 4, 10//1]
    @test var(x) === 16.25
    @test var(y) === 65//4
    @test std(x) === sqrt(16.25)
    @test quantile(x, 0.5)  === 3.0
    @test quantile(x, 1//2) === 3//1
end

@testset "Promotion in covzm. Issue #8080" begin
    A = [1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
    @test Base.covzm(A) - mean(A, 1)'*mean(A, 1)*size(A, 1)/(size(A, 1) - 1) ≈ cov(A)
    A = [1//1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
    @test (A'A - size(A, 1)*Base.mean(A, 1)'*Base.mean(A, 1))/4 == cov(A)
end

@testset "Mean along dimension of empty array" begin
    a0  = zeros(0)
    a00 = zeros(0, 0)
    a01 = zeros(0, 1)
    a10 = zeros(1, 0)
    @test isequal(mean(a0, 1)      , fill(NaN, 1))
    @test isequal(mean(a00, (1, 2)), fill(NaN, 1, 1))
    @test isequal(mean(a01, 1)     , fill(NaN, 1, 1))
    @test isequal(mean(a10, 2)     , fill(NaN, 1, 1))
end

@testset "cov/var/std of Vector{Vector}" begin
    x = [[2,4,6],[4,6,8]]
    @test var(x) ≈ vec(var([x[1] x[2]], 2))
    @test std(x) ≈ vec(std([x[1] x[2]], 2))
    @test cov(x) ≈ cov([x[1] x[2]], 2)
end

ext: jl

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
"""
$(SIGNATURES)

Direct inline-style links are properly processed by Julia's Markdown processor but not:

* `[link title][some reference]` and later `[some reference]: http://www.reddit.com`
* `[link title]` and later `[link title]: https://www.mozilla.org`
* `[link title](https://www.google.com "Google's Homepage")` (we don't either)
"""
function find_and_fix_md_links(hs::String)::String
    # 1. find all occurences of things that look like links
    m_link_refs = collect(eachmatch(ESC_LINK_PAT, hs))

    # recuperate the appropriate id which has a chance to match def_names
    ref_names = [
        # no second bracket or empty second bracket ?
        # >> true then the id is in the first bracket   A --> [id] or [id][]
        # >> false then the id is in the second bracket B --> [...][id]
        ifelse(isnothing(ref.captures[3]) ||
               isempty(ref.captures[3]),
                         ref.captures[2], # A. first bracket
                         ref.captures[3]) # B. second bracket
                #
                for ref in m_link_refs]

    isempty(ref_names) || (ref_names = strip.(ref_names))

    # reconstruct the text
    h = IOBuffer()
    head = 1
    i = 0
    for (m, refn) in zip(m_link_refs, ref_names)
        # write what's before
        (head < m.offset) && write(h, subs(hs, head, prevind(hs, m.offset)))
        #
        def = get(PAGE_LINK_DEFS, refn) do
            ""
        end
        if isempty(def)
            # no def found --> just leave it as it was
            write(h, m.match)
        else
            if !isnothing(m.captures[1])
                # CASE: ![alt][id] --> image
                write(h, html_img(def, refn))
            else
                # It's a link
                if isnothing(m.captures[3]) || isempty(m.captures[3])
                    # CASE: [id] or [id][] the id is also the link text
                    write(h, html_ahref(def, refn))
                else
                    # It's got a second, non-empty bracket
                    # CASE: [name][id]
                    name = m.captures[2]
                    write(h, html_ahref(def, name))
                end
            end
        end
        # move the head after the match
        head = nextind(hs, m.offset + lastindex(m.match) - 1)
    end
    strlen = lastindex(hs)
    (head ≤ strlen) && write(h, subs(hs, head, strlen))

    return String(take!(h))
end

"""
$(SIGNATURES)

For a project website, for instance `username.github.io/project/` all relative
paths should eventually be pre-pended with `/project/`. This would happen just
before you publish the website (see `optimize` or `publish`).
"""
function fix_links(pg::String)::String
    prepath = globvar(:prepath)
    isempty(prepath) && return pg
    pp = strip(prepath, '/')
    ss = SubstitutionString("\\1=\"/$(pp)/")
    # replace things that look like href="/..." with href="/$prepath/..."
    return replace(pg, r"(src|href|formaction|url)\s*?=\s*?\"\/" => ss)
end


text:
using StanSample, Test

import CompatHelperLocal as CHL
CHL.@check()

if haskey(ENV, "JULIA_CMDSTAN_HOME") || haskey(ENV, "CMDSTAN")

  TestDir = @__DIR__
  tmpdir = mktempdir()
  println()

  test_bernoulli = [
    "test_keyedarray/test_bernoulli_keyedarray_01.jl",
    "test_keyedarray/test_keyedarray.jl",
  ]

  @testset "Bernoulli array tests" begin
      include(joinpath(TestDir, "test_bernoulli/test_bernoulli_keyedarray_01.jl"))

      if success(rc)

        sdf = read_summary(sm)
        @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05

        (samples, parameters) = read_samples(sm, :array;
          return_parameters=true)
        @test size(samples) == (1000, 1, 6)
        @test length(parameters) == 1

        (samples, parameters) = read_samples(sm, :array;
          return_parameters=true, include_internals=true)
        @test size(samples) == (1000, 8, 6)
        @test length(parameters) == 8

        samples = read_samples(sm, :array;
          include_internals=true)
        @test size(samples) == (1000, 8, 6)

        samples = read_samples(sm, :array)
        @test size(samples) == (1000, 1, 6)
      end

      include(joinpath(TestDir, "test_bernoulli/test_bernoulli_array_02.jl"))
      if success(rc)
        sdf = read_summary(sm)
        @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05

        (samples, parameters) = read_samples(sm, :array;
          return_parameters=true)
        @test size(samples) == (250, 1, 4)
        @test length(parameters) == 1

        samples = read_samples(sm, :array;
          include_internals=true)
        @test size(samples) == (250, 8, 4)
      end
  end

  basic_run_tests = [
    "test_bernoulli/test_bernoulli_array_01.jl",
    "test_basic_runs/test_bernoulli_dict.jl",
    "test_basic_runs/test_bernoulli_array_dict_1.jl",
    "test_basic_runs/test_bernoulli_array_dict_2.jl",
    "test_basic_runs/test_parse_interpolate.jl"
  ]

  @testset "Bernoulli cpp tests" begin
      println("\nTesting bernoulli_cpp.jl")
      include(joinpath(TestDir, "test_bernoulli", "bernoulli_cpp.jl"))
  end

  @testset "Bernoulli basic run tests" begin
    for test in basic_run_tests
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
      @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05
    end
    println()
  end

  generate_quantities_tests = [
    "test_generate_quantities/test_generate_quantities.jl"
  ]

  @testset "Generate_quantities tests" begin
    for test in generate_quantities_tests
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
    end
    println()
  end

  test_tables_interface = [
    "test-tables-interface/ex-00.jl",
    "test-tables-interface/ex-01.jl",
    "test-tables-interface/ex-02.jl",
    "test-tables-interface/ex-03.jl",
    "test-tables-interface/ex-04.jl",
    "test-tables-interface/ex-05.jl"
  ]

  @testset "Tables.jl interface" begin
    for test in test_tables_interface
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
    end
    println()
  end

  test_dimensionaldata = [
    "test_dimensionaldata/test_dimensionaldata.jl",
  ]

  @testset "DimensionalData interface" begin
    for test in test_dimensionaldata
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
    end
    println()
  end

  test_keywords = [
    "test_keywords/test_bernoulli_keyedarray_01.jl",
    "test_keywords/test_bernoulli_keyedarray_02.jl",
    "test_keywords/test_bernoulli_keyedarray_03.jl",
  ]

  @testset "Seed and num_chains keywords" begin
    for test in test_keywords
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
    end
    println()
  end

  include(joinpath(TestDir, ))
  println()

else
  println("\nCMDSTAN and JULIA_CMDSTAN_HOME not set. Skipping tests")
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
using AutoPreallocation
using BenchmarkTools
using Test

# Its not worth the effort of writing these tests to work both on 64bit and 32bit.
# Pointer size basically gets counted in a ton of things we measure, and distangling
# number of pointer allocations from content allocations is too much work.
sizeof(1) == 8 || error("These tests cn only be run on 64-bit systems.")

const files = (
    "integration_tests.jl",
    "record_types.jl",
    "readme_examples.jl"
)
@testset "AutoPreallocation.jl" begin
    @testset "$file" for file in files
        include(file)
    end
end


text:
"""
将一个图形进行三角切分
"""
module Triangulated

using Printf
using ..Basics

export split_square, split_hexagon
export save_triangulated, load_triangulated

export TYPE_LADJS

const TYPE_LADJS{P} = Vector{Tuple{Union{Missing, P}, Union{Missing, P}, Union{Missing, P}}}

"""
将一个正方形切分成4个直角三角形
"""
function split_single_square(squ::Basics.AbstractRectangle{:AXISSQUARE})
    rtver = squ.center
    top = RtTriangle(rtver, squ.vertex[4], squ.vertex[1])
    left = RtTriangle(rtver, squ.vertex[1], squ.vertex[2])
    btm = RtTriangle(rtver, squ.vertex[2], squ.vertex[3])
    right = RtTriangle(rtver, squ.vertex[3], squ.vertex[4])
    return [top, left, btm, right]
end

const RtTriOrMissing = Union{Basics.AbstractTriangle{:RT}, Missing}
const EqTriOrMissing = Union{Basics.AbstractTriangle{:EQ}, Missing}

"""
将一个大正方形切分成很多的直角三角形
"""
function split_square(squ::Basics.AbstractRectangle{:AXISSQUARE}, nps::Int64)
    lsqus = Matrix{Basics.AbstractRectangle{:AXISSQUARE}}(undef, (nps, nps))
    width = squ.width / nps
    statx, staty = squ.vertex[2].x, squ.vertex[2].y
    #从左下角开始，先水平的顺序进行标号
    # [idx2]
    # nps
    # ...
    # 1   x  x
    # 0   x  x
    #     0  1  ...  nps [idx1]
    for idx12 in CartesianIndices(lsqus)
        idx1, idx2 = Tuple(idx12)
        xct = (idx1 - 0.5) * width + statx
        yct = (idx2 - 0.5) * width + staty
        lsqus[idx1, idx2] = Square(Point2D(xct, yct), width)
    end
    #将每个小正方切分
    ltris = Array{Basics.AbstractTriangle{:RT}, 3}(undef, (nps, nps, 4))
    for idx12 in CartesianIndices(lsqus)
        idx1, idx2 = Tuple(idx12)
        ltris[idx1, idx2, :] = split_single_square(lsqus[idx1, idx2])
    end
    #找到每个小三角形挨着的小三角形
    #这个时候的顺序应该是和小三角形的边的顺序一致的
    ladjs = Array{Tuple{RtTriOrMissing, RtTriOrMissing, RtTriOrMissing}, 3}(
        undef, (nps, nps, 4)
    )
    for idx12 in CartesianIndices(lsqus)
        idx1, idx2 = Tuple(idx12)
        #上面的那个小三角，第一条边是右侧，第二条是上边
        adjs2 = idx2 < nps ? ltris[idx1, idx2+1, 3] : missing
        adjs = (ltris[idx1, idx2, 4], adjs2, ltris[idx1, idx2, 2])
        ladjs[idx1, idx2, 1] = adjs
        #左面的小三角，第一条边是上边，第二条是左面
        adjs2 = idx1 > 1 ? ltris[idx1-1, idx2, 4] : missing
        adjs = (ltris[idx1, idx2, 1], adjs2, ltris[idx1, idx2, 3])
        ladjs[idx1, idx2, 2] = adjs
        #下面的小三角，第一条边是左面，第二条是下面
        adjs2 = idx2 > 1 ? ltris[idx1, idx2-1, 1] : missing
        adjs = (ltris[idx1, idx2, 2], adjs2, ltris[idx1, idx2, 4])
        ladjs[idx1, idx2, 3] = adjs
        #右面的小三角，第一条边是下面，第二条是右边
        adjs2 = idx1 < nps ? ltris[idx1+1, idx2, 2] : missing
        adjs = (ltris[idx1, idx2, 3], adjs2, ltris[idx1, idx2, 1])
        ladjs[idx1, idx2, 4] = adjs
    end
    #将小三角和它对应的近邻都reshape
    ltris = reshape(ltris, nps*nps*4)
    ladjs = reshape(ladjs, nps*nps*4)
    return ltris, ladjs
end

"""
切分一个水平的等腰梯形，按照底边的差进行切分，所以角度必须是60度
的等腰梯形
"""
function split_trapezoid(lseg::Segment, sseg::Segment)
    dlen = lseg.length - sseg.length
    part = sseg.length / dlen
    @assert isapprox(part, round(part); atol=1e-6)
    part = Int64(round(part))
    #两个边上平均分得的所有点
    lpt = Vector{Point2D}(undef, part+2)
    spt = Vector{Point2D}(undef, part+1)
    #
    spt[1] = sseg.pt1
    lpt[1] = lseg.pt1
    #
    for idx in 1:1:part
        coef = idx / part
        coefb = 1 - coef
        spt[idx+1] = middle_point(sseg.pt1, sseg.pt2, sc1=coefb, sc2=coef)
        coef = idx / (part+1)
        coefb = 1 - coef
        lpt[idx+1] = middle_point(lseg.pt1, lseg.pt2, sc1=coefb, sc2=coef)
    end
    lpt[end] = lseg.pt2
    #
    #从长的开始计，第二个点到倒数第二个点
    eqtris = Vector{Basics.AbstractTriangle{:EQ}}(undef, 1+2*part)
    eqtris[1] = EqTriangle(lpt[1], spt[1], lpt[2])
    for idx in 1:1:part
        eqtris[2*idx] = EqTriangle(lpt[idx+1], spt[idx+1], spt[idx])
        eqtris[2*idx+1] = EqTriangle(lpt[idx+1], spt[idx+1], lpt[idx+2])
    end
    return eqtris
end

"""
切分一个水平的六边型
"""
function split_hexagon(hexa::Basics.AbstractHexagon{:EQ}, nps::Int64)
    #六边形分成两个等腰梯形
    topedge = hexa.edges[2]
    midedge = Segment(hexa.vertex[1], hexa.vertex[4])
    btmedge = hexa.edges[5]
    #
    #所有的小三角形
    eqtris = Vector{Basics.AbstractTriangle{:EQ}}(undef, 6*nps^2)
    #
    #先处理上面半个梯形
    #注意这里由于Hexagon是逆时针对顶点进行的排序
    #所以梯形和三角形增长的方向是从右往左
    llines = Vector{Segment}(undef, nps+1)
    llines[1] = topedge
    for idx in 1:1:nps
        coef = idx / nps
        coefb = 1 - coef
        llines[1+idx] = Segment(
            middle_point(topedge.pt1, midedge.pt1, sc1=coefb, sc2=coef),
            middle_point(topedge.pt2, midedge.pt2, sc1=coefb, sc2=coef)
        )
    end
    #llines中一共有nps+1个线
    #注意由于梯形的点是逆时针的，所以这边三角形的增长也是从右向左的
    sptr = 1::Int64
    for idx in 1:1:nps
        ets = split_trapezoid(llines[idx+1], llines[idx])
        @assert length(ets) == 2*nps + 1 + 2*(idx-1)
        eqtris[sptr:sptr+length(ets)-1] = ets
        sptr += length(ets)
    end
    @assert sptr == (length(eqtris) // 2)+1
    #再处理下面半个梯形的
    #这里需要注意下半个梯形的点的顺序问题
    llines = Vector{Segment}(undef, nps+1)
    llines[1] = midedge
    for idx in 1:1:nps
        coef = idx / nps
        coefb = 1 - coef
        llines[1+idx] = Segment(
            middle_point(midedge.pt1, btmedge.pt2, sc1=coefb, sc2=coef),
            middle_point(midedge.pt2, btmedge.pt1, sc1=coefb, sc2=coef)
        )
    end
    #llines中共有nps+1条线
    for idx in 1:1:nps
        ets = split_trapezoid(llines[idx], llines[idx+1])
        @assert length(ets) == 4*nps - 1 - 2*(idx-1)
        eqtris[sptr:sptr+length(ets)-1] = ets
        sptr += length(ets)
    end
    @assert sptr == length(eqtris)+1
    #开始处理相邻关系
    ladjs = Vector{Tuple{EqTriOrMissing, EqTriOrMissing, EqTriOrMissing}}(
        undef, length(eqtris)
    )
    # Warning: 这个ladjs中的相邻的三角形的顺序一定和Eqtriangle的边的顺序是一样的
    #先是nps个上半层，其中三角形的数量是2*nps+1+2*idx
    #idx从0到nps-1
    #先处理最上面的一层的相邻
    #最上面一层就是从[1,2*nps+1]这个范围内的所有
    for idx in 1:1:(2*nps+1)
        #每个三角都有三条边，于是有三个相邻
        #这些是底边在下面的相邻的有下一个，和下一行正对的
        #上半部分中，底边在下面的三角形的点是按照右-左-下的顺序来的
        #底边在上面的三角形是左-上-右
        if idx % 2 == 1
            right = idx == 1 ? missing : eqtris[idx-1]
            left = idx == 2*nps+1 ? missing : eqtris[idx+1]
            ladjs[idx] = (right, left, eqtris[idx+2*nps+2])
        else
            ladjs[idx] = (eqtris[idx+1], missing, eqtris[idx-1])
        end
    end
    #之后，还有nps-1个上半部分的梯形
    sptr = 2*nps + 1
    for npsidx in 2:1:nps
        trinum = 2*nps+1+2*(npsidx-1)
        #如果到了最后一个，下面的一行的三角形的数量和这一行是一样多的
        #计算的时候不需要多添加一个1
        offset = npsidx == nps ? 0 : 1
        for idx in 1:1:trinum
            tot_idx = sptr + idx
            #上半部分中，底边在下面的三角形的点是按照右-左-下的顺序来的
            #底边在上面的三角形是左-上-右
            if idx % 2 == 1
                right = idx == 1 ? missing : eqtris[tot_idx-1]
                left = idx == trinum ? missing : eqtris[tot_idx+1]
                ladjs[tot_idx] = (right, left, eqtris[tot_idx+trinum+offset])
            else
                ladjs[tot_idx] =
                (eqtris[tot_idx+1], eqtris[tot_idx-trinum+1], eqtris[tot_idx-1])
            end
        end
        sptr += trinum
    end
    #下半部分从后往前开始处理
    #最下面的一个梯形
    for idx in 1:1:2*nps+1
        #每个三角都有三条边，于是有三个相邻
        #这些是底边在上面的，相邻的有下一个，上面对的，还有上一个
        #下半部分中，底边在上面的三角形是按照右-左-上的顺序来的
        #底边在下面的三角形是按照左-下-右的顺序来的
        tot_idx = length(eqtris) - (2*nps+1) + idx
        if idx % 2 == 1
            rightone = idx == 1 ? missing : eqtris[tot_idx-1]
            leftone = idx == 2*nps+1 ? missing : eqtris[tot_idx+1]
            ladjs[tot_idx] = (rightone, leftone, eqtris[tot_idx-2-2*nps])
        else
            ladjs[tot_idx] = (eqtris[tot_idx+1], missing, eqtris[tot_idx-1])
        end
    end
    sptr = length(eqtris) - (2*nps+1)
    for npsidx in 2:1:nps
        offset = npsidx != nps ? 1 : 0
        trinum = 2*nps+1+2*(npsidx-1)
        for idx in 1:1:trinum
            tot_idx = sptr - trinum + idx
            #下半部分中，底边在上面的三角形是按照右-左-上的顺序来的
            if idx % 2 == 1
                rightone = idx == 1 ? missing : eqtris[tot_idx-1]
                leftone = idx == trinum ? missing : eqtris[tot_idx+1]
                ladjs[tot_idx] = (rightone, leftone, eqtris[tot_idx-offset-trinum])
            #下半部分中，底边在下面的三角形是按照左-下-右的顺序来的
            else
                ladjs[tot_idx] =
                (eqtris[tot_idx+1], eqtris[tot_idx+trinum-1], eqtris[tot_idx-1])
            end
        end
        sptr -= trinum
    end
    return eqtris, ladjs
end

"""
平行于坐标轴的图形的字符串
"""
macro axis_polygon_string(poly)
    type2name = Dict(
        Basics.AbstractRectangle{:AXISSQUARE} => "Square",
        Basics.AbstractHexagon{:EQ} => "EqHexagon"
    )
    return quote
        #宏里面不要加return，代换回去以后直接就把函数返回了
        par1 = @sprintf "%s\n" $(esc(type2name))[typeof($(esc(poly)))]
        par2 = @sprintf "%.12f,%.12f\n" $(esc(poly)).center.x $(esc(poly)).center.y
        par3 = @sprintf "%.12f\n" $(esc(poly)).edges[1].length
        par1 * par2 * par3
    end
end

"""
三角型的字符串
"""
macro triangle_string(tri)
    type2name = Dict(
        Basics.AbstractTriangle{:COM} => "Triangle",
        Basics.AbstractTriangle{:EQ} => "EqTriangle",
        Basics.AbstractTriangle{:RT} => "RtTriangle"
    )
    return esc(quote
        par1 = @sprintf "%s:" $(type2name)[typeof($tri)]
        pt1 = @sprintf "%.12f,%.12f:" ($tri).vertex[1].x $(tri).vertex[1].y
        pt2 = @sprintf "%.12f,%.12f:" ($tri).vertex[2].x ($tri).vertex[2].y
        pt3 = @sprintf "%.12f,%.12f\n" ($tri).vertex[3].x ($tri).vertex[3].y
        par1*pt1*pt2*pt3
    end)
end

"""
保存结果到文件, 会分别保存到三个文件里面
"""
function save_triangulated(
    poly::AbstractPolygon,
    ltris::Vector{T},
    ladjs::Vector{Tuple{Union{Missing, T}, Union{Missing, T}, Union{Missing, T}}},
    prefix
) where T <: Basics.AbstractTriangle
    @assert isa(poly, Union{
        Basics.AbstractRectangle{:AXISSQUARE},
        Basics.AbstractHexagon{:EQ}}
        )
    #保存描述图形的文件
    plystr = @axis_polygon_string poly
    outf = open(prefix*".ply", "w")
    write(outf, plystr)
    close(outf)
    #保存描述三角形的文件
    outf = open(prefix*".tri", "w")
    tri2idx = Dict{
        Union{Missing, AbstractPolygon},
        Int64
    }(missing => -1)
    for (idx, tri) in enumerate(ltris)
        tristr = @triangle_string tri
        tri2idx[tri] = idx
        write(outf, tristr)
    end
    close(outf)
    #保存描述相邻三角型的文件
    outf2 = open(prefix*".adj", "w")
    for idx in 1:1:length(ladjs)
        adjs = ladjs[idx]
        aidx1, aidx2, aidx3 = tri2idx[adjs[1]], tri2idx[adjs[2]], tri2idx[adjs[3]]
        write(outf2, @sprintf "%d:%d:%d\n" aidx1 aidx2 aidx3)
    end
    close(outf2)
end

"""
读取配置
"""
function load_triangulated(prefix)
    #从名称到类型
    name2type = Dict(
        "Triangle" => Basics.AbstractTriangle{:COM},
        "EqTriangle" => Basics.AbstractTriangle{:EQ},
        "RtTriangle" => Basics.AbstractTriangle{:RT}
    )
    #读取多边形
    inf = open(prefix*".ply", "r")
    plystr = readlines(inf)
    cntstr = split(plystr[2], ',')
    poly = eval(Symbol(plystr[1]))(
        Point2D(parse(Float64, cntstr[1]), parse(Float64, cntstr[2])),
        parse(Float64, plystr[3])
    )
    close(inf)
    #读取切分的三角形
    inf = open(prefix*".tri", "r")
    ltristr = readlines(inf)
    typestr = split(ltristr[1], ':')[1]
    ltris = Vector{name2type[typestr]}(undef, length(ltristr))
    for (idx, lstr) in enumerate(ltristr)
        context = split(lstr, ':')
        pt1 = split(context[2], ',')
        pt2 = split(context[3], ',')
        pt3 = split(context[4], ',')
        ltris[idx] = eval(Symbol(context[1]))(
            Point2D(parse(Float64, pt1[1]), parse(Float64, pt1[2])),
            Point2D(parse(Float64, pt2[1]), parse(Float64, pt2[2])),
            Point2D(parse(Float64, pt3[1]), parse(Float64, pt3[2]))
        )
    end
    close(inf)
    #读取近邻的配置
    inf = open(prefix*".adj", "r")
    ladjstr = readlines(inf)
    typeormis = Union{Missing, name2type[typestr]}
    ladjs = Vector{Tuple{typeormis, typeormis, typeormis}}(
        undef, length(ladjstr)
    )
    for (idx, lstr) in enumerate(ladjstr)
        context = split(lstr, ':')
        adj1 = context[1] == "-1" ? missing : ltris[parse(Int64, context[1])]
        adj2 = context[2] == "-1" ? missing : ltris[parse(Int64, context[2])]
        adj3 = context[3] == "-1" ? missing : ltris[parse(Int64, context[3])]
        ladjs[idx] = (adj1, adj2, adj3)
    end
    close(inf)
    return poly, ltris, ladjs
end

end # end module


text:
module PriceDist

using Dates

include("eventIO.jl")
include("stateResult.jl")
include("forward.jl")

export MktEvent, IODataCache, fillBinaryEvent!
export StateResult, processPairs
export walkForward

end # module

max_stars_repo_name: Andy-P/PriceDist

max_issues_repo_name: Andy-P/PriceDist

max_forks_repo_name: Andy-P/PriceDist


text:
using Zygote, OptimKit
using ..ITensorAutoHOOT
using ..ITensorNetworks
using ..ITensorAutoHOOT: batch_tensor_contraction
using ..ITensorNetworks:
  PEPS,
  inner_network,
  inner_networks,
  flatten,
  insert_projectors,
  split_network,
  rayleigh_quotient,
  Models,
  tree_w_projectors

function loss_grad_wrap(
  peps::PEPS,
  Hs::Array{Models.LineMPO},
  ::typeof(insert_projectors);
  cutoff=1e-15,
  maxdim=100,
)
  Hs_row = [H for H in Hs if H.coord[2] isa Colon]
  Hs_column = [H for H in Hs if H.coord[1] isa Colon]
  init_call = true
  cache = NetworkCache()
  function loss(peps::PEPS)
    tn_split_row, tn_split_column, projectors_row, projectors_column = insert_projectors(
      peps, cutoff, maxdim
    )
    peps_bra = addtags(linkinds, peps, "bra")
    peps_ket = addtags(linkinds, peps, "ket")
    peps_bra_rot = addtags(linkinds, peps, "brarot")
    peps_ket_rot = addtags(linkinds, peps, "ketrot")
    sites = commoninds(peps_bra, peps_ket)
    peps_bra_split = split_network(peps_bra)
    peps_ket_split = split_network(peps_ket)
    peps_ket_split_ham = prime(sites, peps_ket_split)
    peps_bra_split_rot = split_network(peps_bra_rot, true)
    peps_ket_split_rot = split_network(peps_ket_rot, true)
    peps_ket_split_rot_ham = prime(sites, peps_ket_split_rot)
    # generate network
    network_list_row = inner_networks(
      peps_bra_split,
      peps_ket_split,
      peps_ket_split_ham,
      projectors_row,
      Hs_row,
      tree_w_projectors,
    )
    network_list_column = inner_networks(
      peps_bra_split_rot,
      peps_ket_split_rot,
      peps_ket_split_rot_ham,
      projectors_column,
      Hs_column,
      tree_w_projectors,
    )
    network_inner = inner_network(
      peps_bra_split, peps_ket_split, projectors_row[1], tree_w_projectors
    )
    network_list = vcat(network_list_row, network_list_column, [network_inner])
    variables = flatten([
      peps_bra_split,
      peps_ket_split,
      peps_ket_split_ham,
      peps_bra_split_rot,
      peps_ket_split_rot,
      peps_ket_split_rot_ham,
    ])
    if init_call == true
      cache = NetworkCache(network_list)
      init_call = false
    end
    inners = batch_tensor_contraction(network_list, cache, variables...)
    return rayleigh_quotient(inners)
  end
  loss_w_grad(peps::PEPS) = loss(peps), gradient(loss, peps)[1]
  return loss_w_grad
end

function gradient_descent(
  peps::PEPS,
  Hs::Array,
  ::typeof(insert_projectors);
  stepsize::Float64,
  num_sweeps::Int,
  cutoff=1e-15,
  maxdim=100,
)
  loss_w_grad = loss_grad_wrap(peps, Hs, insert_projectors; cutoff=cutoff, maxdim=maxdim)
  return gradient_descent(peps, loss_w_grad; stepsize=stepsize, num_sweeps=num_sweeps)
end

function gradient_descent(
  peps::PEPS,
  Hs::Array,
  ::typeof(insert_projectors),
  ::typeof(backtracking_linesearch);
  beta::Float64,
  num_sweeps::Int,
  cutoff=1e-15,
  maxdim=100,
)
  loss_w_grad = loss_grad_wrap(peps, Hs, insert_projectors; cutoff=cutoff, maxdim=maxdim)
  return gradient_descent(
    peps, loss_w_grad, backtracking_linesearch; beta=beta, num_sweeps=num_sweeps
  )
end

function OptimKit.optimize(
  peps::PEPS,
  Hs::Array,
  ::typeof(insert_projectors);
  num_sweeps::Int,
  method="GD",
  cutoff=1e-15,
  maxdim=100,
)
  loss_w_grad = loss_grad_wrap(peps, Hs, insert_projectors; cutoff=cutoff, maxdim=maxdim)
  return optimize(peps, loss_w_grad; num_sweeps=num_sweeps, method=method)
end

function gd_error_tracker(
  peps::PEPS, Hs::Vector; stepsize::Float64, num_sweeps::Int, cutoff=1e-15, maxdim=100
)
  ITensors.set_warn_order(40)
  loss_w_grad = loss_grad_wrap(peps, Hs)
  loss_w_grad_approx = loss_grad_wrap(
    peps, Hs, insert_projectors; cutoff=cutoff, maxdim=maxdim
  )
  return gd_error_tracker(
    peps,
    loss_w_grad,
    loss_w_grad_approx;
    stepsize=stepsize,
    num_sweeps=num_sweeps,
    cutoff=cutoff,
    maxdim=maxdim,
  )
end


text:
## Dependencies
using Pkg
Pkg.activate("/ibox/afeinbe2/jabante1/CpelNano/CpelNano.jl/")     # <== Modify with your environment path
using CpelNano

## Arguments
out_dir = ARGS[1]       # output directory
out_prefix = ARGS[2]    # output prefix
fa_file = ARGS[3]       # fasta file
mod_file = ARGS[4]      # model file

## Function
function gen_seqs_sims(out_dir,out_prefix,fa_file,mod_file)

    # Set variables
    count = 1
    avg_cov = 0.0

    # Load fasta record chromosome 22
    CpelNano.print_log("Loading fasta record ...")
    fa_record = CpelNano.get_chr_fa_rec("chr22",fa_file)
    chrstart = fa_record.sequence[1]
    chrend = fa_record.sequence[length(fa_record.sequence)]

    # Get chromosome 22 partition (use default partition by nanopolish)
    CpelNano.print_log("Partitioning chromosome ...")
    config = CpelNano.CpelNanoConfig()
    chr_part = CpelNano.get_chr_part("chr22",config,fa_file)

    # Get models from chromosome 22
    CpelNano.print_log("Read in models ...")
    model_regs = CpelNano.read_reg_chr(mod_file,"chr22",fa_file,fa_record)

    # Loop until we reach average 25X coverage
    CpelNano.print_log("Starting read generation...")
    while(avg_cov < 25.0)

        # Generate read_length and read
        read_len = CpelNano.gen_read_len(1)[1]
        genom_pos,new_read = CpelNano.gen_read("chr22",fa_file,read_len)

        # Filter based on read length
        read_len>100 || continue

        # Discard read if it contains N's at all
        n_match = map(x->getfield(x,:offset),eachmatch(r"N",new_read))
        length(n_match)/read_len<0.05 || continue

        # Get indices of estimation regions
        startidx = CpelNano.linearsearch_readstart(genom_pos,chr_part)
        testidx = CpelNano.linearsearch_readend(genom_pos+read_len-1,startidx,chr_part)

        # Get Ising realization and CpG positions. Skip if no CpG sites
        x_vec,cpg_pos = CpelNano.get_methyl_vector(genom_pos,genom_pos+read_len-1,startidx,testidx,model_regs)
        length(cpg_pos)>0 || continue

        # Methylate read and leave if not successful
        new_methyl_read = CpelNano.meth_seq(new_read,x_vec)
        length(new_methyl_read)>0 || continue

        # Print methylated sequence to a fasta file
        io = open(string(out_dir,out_prefix,"_",count,"_reads.fa"),"w")
        println(io,">Read $(count) @ $(genom_pos):$(genom_pos+read_len-1)")
        println(io,new_methyl_read)
        close(io)

        # Output methylation lattice
        io = open(string(out_dir,out_prefix,"_",count,"_x_vec.txt"),"w")
        @inbounds for i in 1:length(x_vec)
            print(io,cpg_pos[i])
            print(io,'\t')
            println(io,x_vec[i])
        end
        close(io)

	    # Calculate new average coverage (substracting N's in chr22)
        avg_cov += read_len / (chrend - chrstart + 1 - 11658691)

        # Increase count counter
        count += 1

        # Print progress every 500 reads
	    mod(count,100)==0 &&  CpelNano.print_log("Read ID: $(count). Avg Cov: $(avg_cov)")

    end

    CpelNano.print_log("Finished read generation. Total number reads: $(count)")

    # Return nothing
    return nothing

end

# Function call
gen_seqs_sims(out_dir,out_prefix,fa_file,mod_file)


text:
export proj_reps, cover_group, check_proj_coeff
"""
    proj_reps(g, coeff, p; R, tol)

Calculate the projective representation of `g` with coefficients `coeff`.

Inputs:
-------
g     : Finite group object.
coeff : Coefficients represented by an integer matrix.
p     : Modulus of the the integer.
R     : Real or not, default = `false`.
tol   : Tolerance, default = `1e-7`.
"""
function proj_reps(
    g::AbstractFiniteGroup,
    coeff::AbstractMatrix{<:Integer},
    p::Integer;
    R::Bool=false,
    tol::Real=1e-7
)
    cg = cover_group(g, coeff, p)
    reps = irreps(cg, R=R)
    check = [check_proj_coeff(g, rep, coeff, p, tol=tol) for rep in reps]
    [reps[i][1:order(g)] for i=1:length(reps) if check[i]]
end
#-------------------------------------------------------------------------------
"""
    cover_group(g::AbstractFiniteGroup, coeff::AbstractMatrix{<:Integer}, p::Integer)

Construct the cover group of `g`, given the projective coefficient `coeff`.
The projective coefficient is unitary and the phase is module p.

Inputs:
-------
g     : Finite group.
coeff : D(gᵢ) * D(gⱼ) = exp(i2π/p * coeff[i,j]) * D(gᵢ⋅gⱼ).
p     : Module of coefficient.
"""
function cover_group(
    g::AbstractFiniteGroup,
    coeff::AbstractMatrix{<:Integer},
    p::Integer
)
    n = order(g)
    mt = Matrix{Int64}(undef, n * p, n * p)
    for i = 1:n, j = 1:n
        k = g[i, j]
        c = coeff[i, j]
        for a = 0:p-1, b = 0:p-1
            row = i + n * a
            col = j + n * b
            val = k + n * mod(a + b + c, p)
            mt[row, col] = val
        end
    end
    FiniteGroup(mt, "Cover Group for $(name(g))")
end
#-------------------------------------------------------------------------------
"""
Check whether a linear representation of the cover group of `g` is indeed the
projective representation obeying the projective coefficients.
"""
function check_proj_coeff(
    g::AbstractFiniteGroup,
    r::AbstractVector{<:AbstractMatrix},
    coeff::AbstractMatrix{<:Integer},
    p::Integer;
    tol::Real=1e-7
)
    n = order(g)
    for i = 1:n, j = 1:n
        k = g[i,j]
        err = r[i] * r[j] - exp(1im * 2π/p * coeff[i, j]) * r[k]
        norm(err) > tol && (return false)
    end
    true
end

#-------------------------------------------------------------------------------
# Chiral projective representation.
#-------------------------------------------------------------------------------
export chiral_proj_reps
"""
    chiral_proj_reps(g, χ; R=true)

Representation of `g` with a chiral operator S, which is represented by
Diagonal([1, …, 1, -1, …, -1]).
"""
function chiral_proj_reps(
    g::AbstractFiniteGroup,
    χ::AbstractVector{<:Integer};
    R::Bool=true
)
    n = order(g)
    bg, coeff = double_group(g, χ)
    reps = proj_reps(bg, coeff, 2, R=R)
    out = []
    for i = 1:length(reps)
        rep = reps[i]
        S = rep[n+1]
        e, v = eigen(Symmetric(S))
        dim = length(e)÷2
        @assert norm(e[1:dim] .+ 1) < 1e-7 "Chiral operator not unitary, eigvals = $e."
        @assert norm(e[dim+1:2dim] .- 1) < 1e-7 "Chiral operator not unitary, eigvals = $e."
        U = hcat(svd(v[:,dim+1:2dim]).U, svd(v[:, 1:dim]).U)
        @assert size(U) == (2dim, 2dim) "Dimension net correct, got $(size(U)), expect $((2dim, 2dim))"
        Ui = U'
        push!(out, [Ui * rep[i] * U for i = 1:n])
    end
    out
end
#-------------------------------------------------------------------------------
"""
    double_group(g::AbstractFiniteGroup, χ::AbstractVector{<:Integer})

Construct a double group with element g⊕S⋅g, satisfying:
    gᵢS⋅gⱼ = S⋅gᵢgⱼ,
The projective coefficients is:
    D(gᵢ)D(S⋅gⱼ) = χᵢD(S⋅gᵢgⱼ)

The function takes Finitegroup `g` and coefficent `χ`, and return the doubled
group together with the projective coefficients.
"""
function double_group(g::AbstractFiniteGroup, χ::AbstractVector{<:Integer})
    @assert length(χ) == length(g.cls) "Length of χ should be equal to classes of g, while length(χ) = $(length(χ))."
    n = order(g)
    coeff = begin
        ct = zeros(Int64, 2n, 2n)
        for i = 1:n
            p = χ[inclass(g, i)]
            if p == -1
                ct[i, n+1:2n] .= 1
                ct[n+i, n+1:2n] .= 1
            elseif p != 1
                error("p should be ±1, get $p.")
            end
        end
        ct
    end
    mt = Matrix{Int64}(undef, 2n, 2n)
    for i=1:n, j=1:n
        k = g[i,j]
        mt[i, j] = k
        mt[n+i, j] = n+k
        mt[i, n+j] = n+k
        mt[n+i, n+j] = k
    end
    FiniteGroup(mt), coeff
end


text:
using Base.Test
using DensityEstimationML
using Distributions

using ForwardDiff
@testset "cdf consistent with pdf $dataset" for dataset in [GenerateDatasets.Likas2()]

    pdf_l2(x) = pdf(dataset, x)
    cdf_l2(x) = cdf(dataset, x)
    for x in 10*rand(1_000_000)-2.0
        abs(x-(3+√2))<10.0^-3 && continue
        abs(x-(2))<10.0^-3 && continue
        abs(x-0)<10.0^-3 && continue

        num = Calculus.gradient(cdf_l2, x)
        ana =  pdf_l2(x)
        diff = abs(num-ana)
        if diff>10.0^-3.0
            error("At x=$x, \t analytical pdf=$ana \t numerical pdf=$num")
        end
    end
end

@testset "Sample valid $dataset" for dataset in [GenerateDatasets.Likas1(), GenerateDatasets.Likas2(), GenerateDatasets.MagdonIsmailAndAtiya()]
    samples = original_sample(dataset)
    @test all(pdf.(dataset, samples) .> 0.0)
end


text:
module Taps
#import Base.getproperty, Base.setproperty!
abstract type Coords{T, N} <: AbstractArray{T, N} end
abstract type Model end
abstract type Database end
# abstract type Finder end

export *
export Paths, get_properties, getproperty, setproperty!, get_distances, get_momentum, get_kinetics
export Coords, Cartesian, Atomic, Fouriered
export Model, Gaussian
export ImageData

mutable struct Paths
    coords::Coords
    model::Model
    prefix::String
    directory::String
    imgdata::Database
    tag::Dict
    cache::Dict
end

function Paths(;coords="Cartesian", coords_kwargs=Dict(), model="MullerBrown", model_kwargs=Dict(),
       prefix="paths", directory=".", imgdata="ImgData", imgdata_kwargs=Dict(),
       tag=Dict(), cache=Dict())
    coords = getfield(Main, Symbol(coords))(;coords_kwargs...)
    model = getfield(Main, Symbol(model))(;model_kwargs...)
    imgdata = getfield(Main, Symbol(imgdata))(;imgdata_kwargs...)
    Paths(coords, model, prefix, directory, imgdata, tag, cache)
end

function Base.getproperty(paths::Paths, name::Symbol, x)
    if length(String(name)) > 7 && String(name)[end-6:end] == "_kwargs"
        name = Symbol(String(name)[1:end-7])
        field = getfield(paths, name)
        x = x == nothing ? Dict([(f, nothing) for f in fieldnames(field)]) : x
        ans = Dict()
        for (key, value) in pairs(x)
            ans[key] = getproperty(field, key, value)
        end
        return ans
    elseif parentmodule(fieldtype(typeof(paths), name)) ∉ [Base, Core]
        typename = typeof(getfield(paths, name))
        return "$typename"
    else
        return getfield(paths, name)
    end
end

function Base.setproperty!(obj::T, name::Symbol, x) where {T<:Union{Paths, Coords, Database, Model}}
    if length(String(name)) > 7 && String(name)[end-6:end] == "_kwargs"
        name = Symbol(String(name)[1:end-7])
        field = getfield(obj, name)
        for (key, value) in pairs(x)
            setproperty!(field, key, value)
        end
    elseif parentmodule(fieldtype(typeof(obj), name)) ∉ [Base, Core]
        # Change the class
        if typeof(x) == String
            newfield = getfield(Main, Symbol(x))()
            setfield!(obj, name, newfield)
        # Complex update
        elseif typeof(x) <: Dict
            # :coords => {kwargs...} case. Update kwargs.
            if get(x, :name, nothing) == nothing
                field = getfield(obj, name)
                for (key, value) in pairs(x)
                    setproperty!(field, name, x)
                end
            # :coords => {:name=> "Sine", :kwargs => {...}} case. Reconstruct coords
            else
                newname = x[:name]
                new_kwargs = x[:kwargs]
                newfield = getfield(Main, Symbol(newname))(;new_kwargs...)
                setfield!(obj, name, newfield)
            end
        end
    else
        setfield!(obj, name, x)
    end
end

@inline get_distances(paths::Paths, args...; kwargs...) = get_distances(paths.model, paths, args...; kwargs...)
@inline get_momentums(paths::Paths, args...; kwargs...) = get_momentums(paths.model, paths, args...; kwargs...)
@inline get_kinetic_energies(paths::Paths, args...; kwargs...) = get_kinetic_energies(paths.model, paths, args...; kwargs...)
@inline get_kenetic_energy_graident(paths::Paths, args...; kwargs...) = get_kenetic_energy_graident(paths.model, paths, args...; kwargs...)
@inline get_velocities(paths::Paths, args...; kwargs...) = get_velocities(paths.model, paths, args...; kwargs...)
@inline get_accelerations(paths::Paths, args...; kwargs...) = get_accelerations(paths.model, paths, args...; kwargs...)
const get_accelerations = get_acceleration
@inline get_masses(paths::Paths, args...; kwargs...) = get_masses(paths.model, paths, args...; kwargs...)
@inline get_effective_mass(paths::Paths, args...; kwargs...) = get_effective_mass(paths.model, paths, args...; kwargs...)
""" Directly calls the :meth:`get_properties` in ``paths.model``"""
@inline get_properties(paths::Paths, args...; kwargs...) = get_properties(paths.model, paths, args...; kwargs...)
""" Calculate potential """
@inline get_potential(paths::Paths, args...; kwargs...) = get_potential(paths.model, paths, args...; kwargs...)
""" Calculate potential( energy) """
const get_potential_energy = get_potential
""" Calculate potentials, individual energy of each atoms"""
@inline get_potentials(paths::Paths, args...; kwargs...) = get_potentials(paths.model, paths, args...; kwargs...)
""" Equivalanet to Calculate potentials"""
const get_potential_energies = get_potentials
""" Calculate potential gradient(s)"""
@inline get_gradients(paths::Paths, args...; kwargs...) = get_gradients(paths.model, paths, args...; kwargs...)
""" Calculate - potential gradient"""
const get_forces = get_gradients
""" Calculate potential gradient"""
const get_gradient = get_gradients

""" Calculate Hessian of a potential"""
@inline get_hessian(paths::Paths, args...; kwargs...) = get_hessian(paths.model, paths, args...; kwargs...);
""" Calculate kinetic + potential energy"""
function get_total_energy(paths::Paths, args...; kwargs...)
    V = get_potential_energy(paths.model, paths, args...; kwargs...)
    T = get_kinetic_energies(paths.model, paths, args...; kwargs...)
    return V + T
end

""" Calculate covariance. It only applies when potential is guassian"""
@inline get_covariance(paths::Paths, args...; kwargs...) = get_covariance(paths.model, paths, args...; kwargs...)

""" Get index of highest potential energy simplified"""
function get_higest_energy_idx(paths::Paths)
    E = get_potential_energy(path)
    return argmax(E)
end

""" Get index of lowest of covariance simplified"""
function get_lowest_confident_idx(paths::Paths)
    cov = get_covariance(path)
    return argmax(diag(cov))
end

"""" list of int; Get rowid of data"""
@inline get_image_data(paths::Paths, args...; kwargs...) = get_image_data(paths.model::Model, paths, args...; kwargs...)

include("./utils/utils.jl")
include("./coordinates/transformations.jl")
include("./coordinates/coordinates.jl")
include("./database/database.jl")
include("./models/models.jl")

end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule psurface_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("psurface")
JLLWrappers.@generate_main_file("psurface", UUID("d6cb37aa-0f65-5fd6-8d86-c55d2c9ef130"))
end  # module psurface_jll


text:
#system & io tools
using DelimitedFiles

# WARNING: nothing here works yet!!
# Everything you see is a work in progress

export ReadConf

#sreaddlm
#swritedlm

#print_shortest
#print_coolest

#ReadConf
############################################
#Read parameters from configure file
#
#input: file=path_to_configure_file
#list of variables to search eg. "teff"
#output: list of values read from conf-file
############################################
function ReadConf(file, vars...)
    f_stream=open(file)
    conf_values=readlines(f_stream)
    close(f_stream)
    var_arr=Any[]
    nvars=length(vars)

    for var in vars
        for line in conf_values
            #Reading all non-comment lines
            if !(ismatch(r"^\s*(?:#|$)",line))
                 m=search(line,var*"=")
                 if m != (0:-1);
                     var_s=strip(line[m[end]+1:end])
                     m2=search(var_s,".")
                     if m2 == (0:-1)
                         var_f=int64(var_s)
                     else
                         var_f=float64(var_s)
                     end
#                     println("found ",var," = ",var_f)
                     if nvars == 1
                         return var_f
                     else
                         push!(var_arr,var_f)
                     end
                 end
            end
        end
    end
    return var_arr
end

#function print_coolest(io::IO, x::FloatingPoint)
function print_coolest(io::IO, x::AbstractFloat)
    if isnan(x); return write(io, "NaN"); end
    if x < 0 write(io,'-') end
    if isinf(x); return write(io, "Inf"); end

    #smallest change
    #<100 use f
    #>100 use E
    #integer?

end
#print_coolest(io::IO, x::Union(FloatingPoint, Integer)) = print_coolest(io, float(x))

function swritedlm(io::IO, a::AbstractVecOrMat, dlm)
    pb = PipeBuffer()
    nr = size(a,1)
    nc = size(a,2)
    for i = 1:nr
        for j = 1:nc
            writedlm_cell(pb, a[i,j], dlm)
            j == nc ? write(pb,'\n') : print(pb,dlm)
        end
        (nb_available(pb) > (16*1024)) && write(io, takebuf_array(pb))
    end
    write(io, takebuf_array(pb))
    nothing
end
#writedlm(io::IO, a::AbstractArray{T,0}, dlm) = writedlm(io, reshape(a,1), dlm) where T
##writedlm{T}(io::IO, a::AbstractArray{T,0}, dlm) = writedlm(io, reshape(a,1), dlm)

#function fwritedlm(io::IO, a::AbstractVecOrMat; header::ASCIIString="")
function fwritedlm(io::IO, a::AbstractVecOrMat; header::String="")
#    pb = PibeBuffer()
    nr = size(a,1)
    nc = size(a,2)

    hnts = split(header, " ")
    if lenght(types) != nc
        hnts = split(header, ",")
        if length(types) != nc
            error("Could not parse header, use space or comma as delimiter")
        end
    end

    headers = Array(String, nc)
    types = Array(String, nc)
    for i = 1:nc
        hnt = split(hnts[i], "%")
        headers[i] = hnt[1]
        types[i] = length(hnt) == 2 ? hnt[2] : ""
    end

    println(headers)
    println(types)

#    write(io, takebuf_array(pb))
#    nothing
end

max_stars_repo_name: thjsal/toolbox

max_issues_repo_name: thjsal/toolbox

max_forks_repo_name: thjsal/toolbox


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file contains functions and structs for models

mutable struct JouliaModel
	parameters
	solver
	f_build_model
	results

	function JouliaModel(args...)
		self = new()
		self.parameters = Dict{Symbol, Dict}()
		self.f_build_model = create_model(args...)
		self.results = Dict{Symbol, DataFrame}()

		return self
	end
end

"""
	create_model(pp::PowerPlants,
		res::RenewableEnergySource,
		nodes::Nodes
		)

Create the model with 'pp', 'res', and 'nodes'.
"""
function create_model(pp::PowerPlants,
		res::RenewableEnergySource,
		nodes::Nodes
		)

	N = nodes.id
	P = pp.id
	RES = res.id

	mc = pp.mc
	exchange = nodes.exchange
	load = nodes.load
	input_gen = pp.capacity
	input_res = res.infeed

	function build_model(T::Array{Int, 1})

		m = Model()

	    @variables m begin
	        G[P,T] >= 0
	        G_RES[RES,T] >= 0
	        LOST_LOAD[T] >= 0
	        LOST_GENERATION[T] >= 0
	    end

		println("Set model objective")
		@objective(m, Min,
			sum(mc[p][t] * G[p,t] for p in P, t in T)
			+ 1000 * sum(LOST_LOAD[t] + LOST_GENERATION[t] for t in T)
		)

		println("Building constraints:")

	    prog = Progress(length(T), 0.2, "MarketClearing...   ", 50)

	    @constraintref MarketClearing[T]
	    for t=T
	    MarketClearing[t] = @constraint(m,

	        sum(G[p,t] for p in P)
	        + sum(G_RES[r,t] for r in RES)
	        - LOST_GENERATION[t]

	        ==

	        sum(load[n][t] for n in N)
	        - sum(exchange[n][t] for n in N)
	        - LOST_LOAD[t] );

	    next!(prog)
	    end
	    JuMP.registercon(m, :MarketClearing, MarketClearing)

	    @constraintref GenerationRestriction[P,T]
	    prog = Progress(length(P)*length(T), 0.2, "Max generation...   ", 50)

	    for p=P, t=T
	    GenerationRestriction[p,t] = @constraint(m,

	        G[p,t] <= input_gen[p][t] );

	    next!(prog)
	    end

	    @constraintref ResRestriction[RES,T]
	    prog = Progress(length(RES)*length(T), 0.2, "Renewable infeed... ", 50)

	    for r=RES, t=T
	    ResRestriction[r,t] = @constraint(m,

	        G_RES[r,t] <= input_res[r][t] );

	    next!(prog)

	    end

		return m
	end # end of build model

	return build_model
end # end function `create_model`

"""
	create_model(pp::PowerPlants,
		res::RenewableEnergySource,
		storages::Storages,
		nodes::Nodes
		)

Create the model with 'pp', 'res', 'storages', and 'nodes'.
"""
function create_model(
	pp::PowerPlants,
	res::RenewableEnergySource,
	storages::Storages,
	nodes::Nodes
	)

	N = nodes.id
	P = pp.id
	S = storages.id
	RES = res.id

	mc = pp.mc
	exchange = nodes.exchange
	load = nodes.load
	input_gen = pp.capacity
	input_res = res.infeed
	sto_capacity = storages.storage
	sto_power = storages.power
	eff = storages.efficiency

	function build_model(T::Array{Int, 1})

		m = Model()

	    @variables m begin
	        G[P,T] >= 0
	        G_RES[RES,T] >= 0
	        G_STOR[S,T] >= 0
	        W_STOR[S,T] >= 0
	        STOR_LEVEL[S,T] >= 0
	        LOST_LOAD[T] >= 0
	        LOST_GENERATION[T] >= 0
	    end

	    for s in S
	        JuMP.fix(STOR_LEVEL[s,T[1]], 0)
	        JuMP.fix(G_STOR[s,T[1]], 0)
	        JuMP.fix(W_STOR[s,T[1]], 0)
	    end

		println("Set model objective")
		@objective(m, Min,
			sum(mc[p][t] * G[p,t] for p in P, t in T)
			+ 1000 * sum(LOST_LOAD[t] + LOST_GENERATION[t] for t in T)
		)

		println("Building constraints:")

	    prog = Progress(length(T), 0.2, "MarketClearing...   ", 50)

	    @constraintref MarketClearing[T]
	    for t=T
	    MarketClearing[t] = @constraint(m,

	        sum(G[p,t] for p in P)
	        + sum(G_RES[r,t] for r in RES)
	        + sum(G_STOR[s,t] for s in S)
	        - LOST_GENERATION[t]

	        ==

	        sum(load[n][t] for n in N)
	        - sum(exchange[n][t] for n in N)
	        + sum(W_STOR[s,t] for s in S)
	        - LOST_LOAD[t] );

	    next!(prog)
	    end
	    JuMP.registercon(m, :MarketClearing, MarketClearing)

	    @constraintref GenerationRestriction[P,T]
	    prog = Progress(length(P)*length(T), 0.2, "Max generation...   ", 50)

	    for p=P, t=T
	    GenerationRestriction[p,t] = @constraint(m,

	        G[p,t] <= input_gen[p][t] );

	    next!(prog)
	    end

	    @constraintref ResRestriction[RES,T]
	    prog = Progress(length(RES)*length(T), 0.2, "Renewable infeed... ", 50)

	    for r=RES, t=T
	    ResRestriction[r,t] = @constraint(m,

	        G_RES[r,t] <= input_res[r][t] );

	    next!(prog)

	    end

		@constraintref StorageCapGeneration[S,T]
	    @constraintref StorageCapWithdraw[S,T]
	    @constraintref StorageCapLevel[S,T]
	    @constraintref StorageBalance[S,T]
	    prog = Progress(length(S)*length(T), 0.2, "Storages...         ", 50)

	    for s=S, t=T
	    StorageCapGeneration[s,t] = @constraint(m,

	        G_STOR[s,t] <= sto_power[s] );

	    StorageCapWithdraw[s,t] = @constraint(m,

	        W_STOR[s,t] <= sto_power[s] );

	    StorageCapLevel[s,t] = @constraint(m,

	        STOR_LEVEL[s,t] <= sto_capacity[s] );

	    if t != T[end]
	    StorageBalance[s,t] = @constraint(m,

	        STOR_LEVEL[s,t+1]
	        ==
	        STOR_LEVEL[s,t] + (eff[s] * W_STOR[s,t]) - G_STOR[s,t] );

	    elseif t == T[end]
	    StorageBalance[s,t] = @constraint(m,

	        STOR_LEVEL[s,T[1]]
	        ==
	        STOR_LEVEL[s,t] + (eff[s] * W_STOR[s,t]) - G_STOR[s,t] );

	    end #end of if statement
	    next!(prog)
	    end

		return m
	end # end of build model

	return build_model
end # end function `create_model`

"""
	create_model(pp::PowerPlants,
		res::RenewableEnergySource,
		storages::Storages,
		nodes::Nodes,
		lines::Lines
		)

Create the model with 'pp', 'res', 'storages', 'nodes', and 'lines'.
"""
function create_model(pp::PowerPlants,
		res::RenewableEnergySource,
		storages::Storages,
		nodes::Nodes,
		lines::Lines
		)

	N = nodes.id
	P = pp.id
	S = storages.id
	L = lines.id
	RES = res.id

	mc = pp.mc
	map_n_p = Dict(map(n-> n=> [p for p in P if Symbol(pp.node[p]) == n], N))
	map_n_s = Dict(map(n-> n=> [s for s in S if Symbol(storages.node[s]) == n], N))
	map_n_res = Dict(map(n-> n=> [r for r in RES if Symbol(res.node[r]) == n], N))
	h,b = calc_h_b(lines, nodes)
	p_max = lines.pmax
	exchange = nodes.exchange
	load = nodes.load
	input_gen = pp.capacity
	input_res = res.infeed
	sto_capacity = storages.storage
	sto_power = storages.power
	eff = storages.efficiency

	function build_model(T::Array{Int, 1})

		m = Model()

	    @variables m begin
	        G[P,T] >= 0
	        G_RES[RES,T] >= 0
	        G_STOR[S,T] >= 0
	        W_STOR[S,T] >= 0
	        STOR_LEVEL[S,T] >= 0
	        LINEFLOW[L,T]
	        DELTA[N,T]
	        LOST_LOAD[N,T] >= 0
	        LOST_GENERATION[N,T] >= 0
	    end

	    for t in T JuMP.fix(DELTA[Symbol(158),t], 0) end

	    for s in S
	        JuMP.fix(STOR_LEVEL[s,T[1]], 0)
	        JuMP.fix(G_STOR[s,T[1]], 0)
	        JuMP.fix(W_STOR[s,T[1]], 0)
	    end

		println("Set model objective")
		@objective(m, Min,
			sum(mc[p][t] * G[p,t] for p in P, t in T)
			+ 1000 * sum(LOST_LOAD[n,t] + LOST_GENERATION[n,t] for n in N, t in T)
		)

		println("Building constraints:")

	    prog = Progress(length(N)*length(T), 0.2, "MarketClearing...   ", 50)

	    @constraintref MarketClearing[N,T]
	    for n=N, t=T
	    MarketClearing[n,t] = @constraint(m,

	        sum(G[p,t] for p in map_n_p[n])
	        + sum(G_RES[r,t] for r in map_n_res[n])
	        + sum(G_STOR[s,t] for s in map_n_s[n])
	        + 500 * sum(b[n,nn] * DELTA[nn,t] for nn in N)
	        - LOST_GENERATION[n,t]

	        ==

	        load[n][t]
	        - exchange[n][t]
	        + sum(W_STOR[s,t] for s in map_n_s[n])
	        - LOST_LOAD[n,t] );

	    next!(prog)
	    end
	    JuMP.registercon(m, :MarketClearing, MarketClearing)

	    @constraintref GenerationRestriction[P,T]
	    prog = Progress(length(P)*length(T), 0.2, "Max generation...   ", 50)

	    for p=P, t=T
	    GenerationRestriction[p,t] = @constraint(m,

	        G[p,t] <= input_gen[p][t] );

	    next!(prog)
	    end

	    @constraintref ResRestriction[RES,T]
	    prog = Progress(length(RES)*length(T), 0.2, "Renewable infeed... ", 50)

	    for r=RES, t=T
	    ResRestriction[r,t] = @constraint(m,

	        G_RES[r,t] <= input_res[r][t] );

	    next!(prog)

	    end

		@constraintref StorageCapGeneration[S,T]
	    @constraintref StorageCapWithdraw[S,T]
	    @constraintref StorageCapLevel[S,T]
	    @constraintref StorageBalance[S,T]
	    prog = Progress(length(S)*length(T), 0.2, "Storages...         ", 50)

	    for s=S, t=T
	    StorageCapGeneration[s,t] = @constraint(m,

	        G_STOR[s,t] <= sto_power[s] );

	    StorageCapWithdraw[s,t] = @constraint(m,

	        W_STOR[s,t] <= sto_power[s] );

	    StorageCapLevel[s,t] = @constraint(m,

	        STOR_LEVEL[s,t] <= sto_capacity[s] );

	    if t != T[end]
	    StorageBalance[s,t] = @constraint(m,

	        STOR_LEVEL[s,t+1]
	        ==
	        STOR_LEVEL[s,t] + (eff[s] * W_STOR[s,t]) - G_STOR[s,t] );

	    elseif t == T[end]
	    StorageBalance[s,t] = @constraint(m,

	        STOR_LEVEL[s,T[1]]
	        ==
	        STOR_LEVEL[s,t] + (eff[s] * W_STOR[s,t]) - G_STOR[s,t] );

	    end #end of if statement
	    next!(prog)
	    end

	    @constraintref LinecapPositive[L,T]
	    @constraintref LinecapNegative[L,T]
	    @constraintref Lineflow[L,T]

	    prog = Progress(length(L)*length(T), 0.2, "Lineflow...         ", 50)

	    for l=L, t=T

	    LinecapPositive[l,t] = @constraint(m,

	        500 * LINEFLOW[l,t] <= p_max[l] );

	    LinecapNegative[l,t] = @constraint(m,

	        500 * LINEFLOW[l,t] >= - p_max[l] );

	    Lineflow[l,t] = @constraint(m,

	        LINEFLOW[l,t] == sum(h[l,n] * DELTA[n,t] for n in N) );

	    next!(prog)
	    end

		return m
	end # end of build model

	return build_model
end # end function `create_model`

"""
    run_model(em::JouliaModel, T::UnitRange{Int}; solver=error("Set a solver!"))

Run the model.
"""
function run_model(em::JouliaModel, T::UnitRange{Int}; solver=error("Set a solver!"))

	total = @elapsed begin
	elap_build = @elapsed m = em.f_build_model(collect(T))
	@info "Building model took $(round(elap_build, digits=1)) s"

	setsolver(m, solver)
	elap_solve = @elapsed status = solve(m)

	@info "Optimization took $(round(elap_solve, digits=1)) s to solve"

	vars = get_variable_names(m)

	for v in vars
		if haskey(em.results, v)
			df = jumpvar_to_df(getindex(m, v), dim_names=[:id,:hour])
			append!(em.results[v], df)
		else
			em.results[v] = jumpvar_to_df(getindex(m, v), dim_names=[:id,:hour])
		end
	end

	price = jumpvar_to_df(getindex(m, :MarketClearing), dual=true)
	haskey(em.results, :Price) ?
		append!(em.results[:Price], price) : em.results[:Price] = price

end # end of total elapsed

	return (Modelstat=status, Buildtime=elap_build, Solvetime=elap_solve,
		Totaltime=total, Starthour=T[1], Endhour=T[end])
end

"""
    run_model(em::JouliaModel, arr_T::Array{UnitRange{Int}, 1};
		solver=error("Set a solver!"))

Run the model.
"""
function run_model(em::JouliaModel, arr_T::Array{UnitRange{Int}, 1}
	;solver=Error("Set a solver!"))

	stat_arr = []
	for T in arr_T
		@info "Starting interval from $(T[1]) to $(T[end])"
		status = run_model(em, T, solver=solver)
		push!(stat_arr, status)
		status.Modelstat != :Optimal && @warn("Instance from $(T[1]) until $(T[end]) was not optimal!")
	end

	return stat_arr
end


text:
# startup script for solving an equation

push!(LOAD_PATH, "/users/creanj/julialib_fork/PUMI.jl")
using PumiInterface # pumi interface
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
include("../equation/Equation.jl")  # equation types
include("../rk4/rk4.jl")  # timestepping
include("./euler/euler.jl")  # solver functions
include("./euler/ic.jl")  # initial conditions functions
include("./euler/adaptfuncs.jl")  # mesh adapt function
#=
# timestepping parameters
delta_t = 0.5
t_max = 1.00
=#

# create operator
sbp = TriSBP{Float64}()  # create linear sbp operator

# create mesh
dmg_name = ".null"
#smb_name = "tri2l.smb"
smb_name = "tri18l.smb"
#smb_name = "adapt_big.smb"
mesh = PumiMesh2(dmg_name, smb_name, 1; dofpernode=4)  #create linear mesh with 1 dof per node

# create euler equation
eqn = EulerEquation(sbp)
println("eqn.bigQT_xi = \n", eqn.bigQT_xi)
println("eqn.bigQT_eta = \n", eqn.bigQT_eta)
println("sbp.QT_xi' = \n", sbp.Q[:,:,1].')
println("sbp.QT_eta' = \n", sbp.Q[:,:,2].')

# create vectors to hold solution at current, previous timestep
u0 = zeros(mesh.numDof)  # solution at previous timestep
u = zeros(mesh.numDof) # solution at current timestep

# populate u0 with initial condition
ICLinear(mesh, sbp, eqn, u0) # change this
println("u0 = \n", u0)
saveSolutionToMesh(mesh, u0)

cfunc3 = cfunction(func3, Void, (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}))

createAnisoFunc(mesh.m_ptr, cfunc3, u)
writeVtkFiles("output_pre", mesh.m_ptr)

  runAnisoAdapt(mesh.m_ptr)

writeVtkFiles("output_post", mesh.m_ptr)
reinitPumiMesh2(mesh)
u = zeros(mesh.numDof)
retrieveSolutionFromMesh(mesh, u);
println("u = \n", u)

# test code, please ignore
#=
# test getF1
for i=1:16
  u0[i] = i
end

f1 = zeros(12)
getF1(mesh, sbp, eqn, u0, 1, f1)
f2 = zeros(12)
getF2(mesh, sbp, eqn, u0, 1, f2)

# test assembleU
#assembleU(f1, 1, u)

vec = [f1[1], f1[5], f1[9]]
assembleU(vec, 1, 1, u)

edgenum_local = getBoundaryEdgeLocalNum(mesh, 1)
println("edgenum_local = ", edgenum_local)
=#

#=
function evalEuler(t, x)
# this function is called by time stepping algorithm
# t is the current time
# x is the solution value at the previous timestep
# u = output, the function value at the current timestep
# u is declared outside this function to avoid reallocating memory

u[:] = 0.0  # zero out u before starting
evalVolumeIntegrals(mesh, sbp, eqn, u, x)
evalBoundaryIntegrals(mesh, sbp, eqn, u, x)
applyMassMatrixInverse(mesh, sbp, eqn, u, x)

return u

end  # end evalEuler

# call timestepper
rk4(evalEuler, delta_t, u0, t_max)
=#


text: println("Hello OIST")

max_stars_repo_name: oist/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
#__precompile__()

#=
A [Julia](http://julialang.org/) library to get remote data
via [Requests.jl](https://github.com/JuliaWeb/Requests.jl).

Inspired by [Pandas-DataReader](https://github.com/pydata/pandas-datareader).

get_components_yahoo
get_data_famafrench
get_data_google
get_data_yahoo
get_data_yahoo_actions
get_quote_google
get_quote_yahoo
DataReader
Options

=#

module DataReaders

    export get, DataReader, DataSymbol
    export DataFrame, TimeArray
    export D_DATAREADERS

    import Base: convert
    import Requests: get, get_streaming
    import DataFrames: readtable, DataFrame, rename!
    import DataStructures: OrderedDict
    import RequestsCache: create_query, execute
    using TimeSeriesIO: TimeArray

    abstract type AbstractDataReader end
    struct DataReaderGoogleDaily <: AbstractDataReader
        retry_count::Int
        pause::AbstractFloat
        timeout::AbstractFloat
        session
    end
    struct DataReaderGoogleQuotes <: AbstractDataReader
        retry_count::Int
        pause::AbstractFloat
        timeout::AbstractFloat
        session
    end
    struct DataReaderYahooDaily <: AbstractDataReader
        retry_count::Int
        pause::AbstractFloat
        timeout::AbstractFloat
        session
    end

    abstract type AbstractDataReaderResponse end

    D_DATAREADERS = OrderedDict(
        "google" => DataReaderGoogleDaily,
        "google-daily" => DataReaderGoogleDaily,
        "google-quotes" => DataReaderGoogleQuotes,
        "yahoo" => DataReaderYahooDaily,
        "yahoo-daily" => DataReaderYahooDaily
    )

    function DataReader(s_source::String; retry_count=3, pause=0.1, timeout=30, session=nothing)
        s_source = lowercase(s_source)
        if haskey(D_DATAREADERS, s_source)
            dr = D_DATAREADERS[s_source]
            dr(retry_count, pause, timeout, session)
        else
            allowed_data_source = keys(D_DATAREADERS)
            s_allowed_data_source = "[" * join(map(s->"\""*s*"\"", allowed_data_source), ", ") * "]"
            error("'$s_source' is not an allowed data source. It must be in $s_allowed_data_source.")
        end
    end

    DEFAULT_DT_END = now(Dates.UTC)
    DEFAULT_DT_START = DEFAULT_DT_END - Dates.Day(7)

    struct DataSymbol
        s::AbstractString
    end
    Base.hash(symb::DataSymbol, h::UInt) = hash(symb.s, hash(:DataSymbol, h))
    Base.:(==)(symb1::DataSymbol, symb2::DataSymbol) = isequal(symb1.s, symb2.s)

    convert(::Type{DataSymbol}, s::String) = DataSymbol(s)

    include("google/daily.jl")
    include("google/quotes.jl")
    include("yahoo/daily.jl")

    function get_several_symbols_to_ordereddict(dr::AbstractDataReader, symbols::Vector{DataSymbol}, args...; kwargs...)
        d = OrderedDict{DataSymbol,AbstractDataReaderResponse}()
        for symb in symbols
            data = get(dr, symb, args...; kwargs...)
            d[symb] = data
        end
        return d
    end

    function get_response(url::AbstractString, query::Dict{AbstractString,AbstractString}, dr::AbstractDataReader)
        if dr.session==nothing
            #r = get_streaming(url; query = query, timeout = dr.timeout)
            r = get(url; query = query, timeout = dr.timeout)
        else
            #error("cached session not supported")
            #r = get(dr.session, url; query = query)
            prepared_query = create_query("get", url; query = query, timeout = dr.timeout)
            r = execute(prepared_query; session=dr.session)
        end
        #println(r)
        #println(typeof(r))
        #println(r.response)
        #if r.response.status / 100 != 2
        if r.status / 100 != 2
            error("Error $(r.status) - not able to process request on the URL: $url with $query")
        end
        return r
    end

    function DataFrame(multi_symbol_response::OrderedDict{DataSymbol,AbstractDataReaderResponse})
        d = OrderedDict{DataSymbol,DataFrame}()
        for (symb, response) in multi_symbol_response
            d[symb] = DataFrame(response)
        end
        return d
    end

end


text:
using Documenter

makedocs(
    debug = true,
    doctestfilters = [r"Ptr{0x[0-9]+}"],
    sitename = "Documenter example",
    pages = ["index.md"],
)


text:
using ArrayViews
using Base.Test

avparent = rand(7, 8)
bvparent = rand(8, 5)

_ab = avparent * bvparent

@test _ab == view(avparent, :, :) * bvparent
@test _ab == avparent * view(bvparent, :, :)
@test _ab == view(avparent, :, :) * view(bvparent, :, :)

for j = 1:size(bvparent,2)
    @test_approx_eq _ab[:,j] view(avparent,:,:) * view(bvparent,:,j)
end

@test avparent[:, 2:2:7] * bvparent[1:3, 1:2:5] == view(avparent, :, 2:2:7) * view(bvparent, 1:3, 1:2:5)


text:
### Activation functions.

function activation_function_selector(activation::ASCIIString)
	if activation == "exponential"
		return activation, exponential_activation!, exponential_backward!
	elseif activation == "leaky_rectified_linear"
		return activation, leaky_rectified_linear_activation!, leaky_rectified_linear_backward!
	elseif activation == "rectified_linear"
		return activation, rectified_linear_activation!, rectified_linear_backward!
	elseif activation == "sigmoid"
		return activation, sigmoid_activation!, sigmoid_backward!
	elseif activation == "softmax"
		return activation, softmax_activation!, softmax_backward!
	elseif activation == "softplus"
		return activation, softplus_activation!, softplus_backward!
	elseif activation == "tanh"
		return activation, tanh_activation!, tanh_backward!
	# Else default to linear activations.
	else
		return "linear", linear_activation!, linear_backward!
	end
end

### Exponential activations.

function exponential_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			ACT[i] = exp(NET[i])
		end
	end
end

function exponential_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			DE_DNET[o] = DELTAS_ABOVE[o] * ACT[o]
		end
	end
end

### Linear activations.

function linear_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			ACT[i] = NET[i]
		end
	end
end

function linear_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			DE_DNET[o] = DELTAS_ABOVE[o]
		end
	end
end

### Leaky rectified linear activations.

function leaky_rectified_linear_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			if NET[i] > 0.0
				ACT[i] = NET[i]
			else
				ACT[i] = 0.01 * NET[i]
			end
		end
	end
end

function leaky_rectified_linear_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			if NET[o] > 0.0
				DE_DNET[o] = DELTAS_ABOVE[o]
			else
				DE_DNET[o] = DELTAS_ABOVE[o] * 0.01
			end
		end
	end
end

### Rectified linear activations.

function rectified_linear_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			if NET[i] > 0.0
				ACT[i] = NET[i]
			else
				ACT[i] = 0.0
			end
		end
	end
end

function rectified_linear_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			if NET[o] > 0.0
				DE_DNET[o] = DELTAS_ABOVE[o]
			else
				DE_DNET[o] = 0.0
			end
		end
	end
end

### Sigmoid activations.

function sigmoid_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			ACT[i] = 1.0 / (1.0 + exp(-NET[i]))
		end
	end
end

function sigmoid_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			DE_DNET[o] = DELTAS_ABOVE[o] * ACT[o] * (1.0 - ACT[o])
		end
	end
end

### Softmax activations.

function softmax_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for j = 1:size(NET, 2)
			# Get maximum net value (for numerical stability).
			maxnet::T = -Inf
			for i = 1:length(NET)
				if NET[i] > maxnet
					maxnet = NET[i]
				end
			end
			# Get sum of exponentials.
			expsum::T = 0.0
			for i = 1:length(NET)
				ACT[i] = exp(NET[i] - maxnet)
				expsum += ACT[i]
			end
			# Set activations.
			for i = 1:length(NET)
				ACT[i] ./= expsum
			end
		end
	end
end

function softmax_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		no = length(NET)
		for o = 1:no
			DE_DNET[o] = 0.0
			for oo = 1:no
				if o == oo
					DE_DNET[o] += DELTAS_ABOVE[oo] * ACT[o] * (1.0 - ACT[o])
				else
					DE_DNET[o] -= DELTAS_ABOVE[oo] * ACT[o] * ACT[oo]
				end
			end
		end
	end
end

### Softplus activations.

function softplus_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			ACT[i] = log(1.0 + exp(NET[i]))
		end
	end
end

function softplus_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			expx = exp(NET[o])
			DE_DNET[o] = DELTAS_ABOVE[o] * expx / (1.0 + expx)
		end
	end
end

### Tanh activations.

function tanh_activation!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T})
	@inbounds begin
		for i = 1:length(NET)
			ACT[i] = tanh(NET[i])
		end
	end
end

function tanh_backward!{T<:FloatingPoint}(NET::Vector{T}, ACT::Vector{T}, DELTAS_ABOVE::Vector{T}, DE_DNET::Vector{T})
	@inbounds begin
		for o = 1:length(NET)
			DE_DNET[o] = DELTAS_ABOVE[o] * (1.0 - ACT[o] * ACT[o])
		end
	end
end


text:
using Coordinates
using Test

function check_coordinate(coord::Coordinate{N, ElType}, axes::Tuple) where {N, ElType}
    @test all(CartesianIndices(coord)) do I
        (@inferred coord[I])::ElType == getindex.(axes, Tuple(I))
    end
    @test size(coord) == length.(axes)
    @test coordinateaxes(coord) == axes
end
check_coordinate(coord::Coordinate{N}, ax) where {N} = check_coordinate(coord, ntuple(i -> ax, Val(N)))

@testset "Constructors" begin
    @test (@inferred Coordinate{2, Tuple{Int, Float64}, Tuple{UnitRange{Int}, Vector{Float64}}}((1:2, [2.0,3.0]))).axes == (1:2, [2.0,3.0])
    @test (@inferred Coordinate{2, Tuple{Int, Char}, Tuple{UnitRange{Int}, String}}((1:2, "abc"))).axes == (1:2, "abc")
    @test_throws ArgumentError Coordinate{2, Tuple{Float64, Float64}, Tuple{UnitRange{Int}, Vector{Float64}}}((1:2, [2.0,3.0]))
    # from axes
    for axes in ((1:3,),
                (1:3, 2:4),
                (1:3, 2:4, 5:8),
                (1:3, 2:2.0:4, "abcd"),
                (1:3, [1,2,3.0], "abcd"))
        N = length(axes)
        ElType = Tuple{map(eltype, axes)...}
        C = (@inferred Coordinate(axes))::Coordinate{N, ElType, typeof(axes)}
        C = (@inferred Coordinate(axes...))::Coordinate{N, ElType, typeof(axes)}
        check_coordinate(C, axes)
    end
    # from axis
    for axis in (1:3, [2,4,4.0], "abcd")
        for N in 1:3
            axes = ntuple(i -> axis, Val(N))
            ElType = Tuple{map(eltype, axes)...}
            C = (@inferred Coordinate{N}(axis))::Coordinate{N, ElType, typeof(axes)}
            check_coordinate(C, axis)
        end
    end
end

@testset "Operations" begin
    @testset "subarray" begin
        x = Coordinate(1:4, [2,3,4])
        sub = Coordinate(2:3, [2,3])
        vec = Coordinate(1:4, 3)
        rowvec = Coordinate(3, [3,4])
        @test (@inferred x[2:3, 1:2])::typeof(sub) == sub
        @test (@inferred x[CartesianIndices((2:3, 1:2))])::typeof(sub) == sub
        @test (@inferred x[:, 2])::typeof(vec) == vec
        @test (@inferred x[3, 2:3])::typeof(rowvec) == rowvec
    end
    @testset "transpose" begin
        x = Coordinate(1:2, [1,2,3])
        xᵀ = Coordinate([1,2,3], 1:2)
        @test (@inferred adjoint(x))::typeof(xᵀ) == xᵀ
        @test (@inferred transpose(x))::typeof(xᵀ) == xᵀ
    end
end