text:
struct TstVanillaLoss <: SupervisedLoss end

@testset "Fallback implementations; not prior knowledge" begin
    @test isminimizable(TstVanillaLoss()) == false
    @test isdifferentiable(TstVanillaLoss()) == false
    @test isdifferentiable(TstVanillaLoss(), 0) == false
    @test istwicedifferentiable(TstVanillaLoss()) == false
    @test istwicedifferentiable(TstVanillaLoss(), 0) == false

    @test isstronglyconvex(TstVanillaLoss()) == false
    @test isstrictlyconvex(TstVanillaLoss()) == false
    @test isconvex(TstVanillaLoss()) == false

    @test isnemitski(TstVanillaLoss()) == false
    @test islipschitzcont(TstVanillaLoss()) == false
    @test islocallylipschitzcont(TstVanillaLoss()) == false
    @test isclipable(TstVanillaLoss()) == false
    @test ismarginbased(TstVanillaLoss()) == false
    @test isdistancebased(TstVanillaLoss()) == false
    @test issymmetric(TstVanillaLoss()) == false
end

struct TstStronglyConvexLoss <: SupervisedLoss end
LossFunctions.isstronglyconvex(::TstStronglyConvexLoss) = true

@testset "Fallback implementations; strongly convex" begin
    @test isminimizable(TstStronglyConvexLoss()) == true

    @test isdifferentiable(TstStronglyConvexLoss()) == false
    @test isdifferentiable(TstStronglyConvexLoss(), 0) == false
    @test istwicedifferentiable(TstStronglyConvexLoss()) == false
    @test istwicedifferentiable(TstStronglyConvexLoss(), 0) == false

    @test isstronglyconvex(TstStronglyConvexLoss()) == true
    @test isstrictlyconvex(TstStronglyConvexLoss()) == true
    @test isconvex(TstStronglyConvexLoss()) == true

    @test isnemitski(TstStronglyConvexLoss()) == true
    @test islipschitzcont(TstStronglyConvexLoss()) == false
    @test islocallylipschitzcont(TstStronglyConvexLoss()) == true
    @test isclipable(TstStronglyConvexLoss()) == false
    @test ismarginbased(TstStronglyConvexLoss()) == false
    @test isdistancebased(TstStronglyConvexLoss()) == false
    @test issymmetric(TstStronglyConvexLoss()) == false
end

struct TstTwiceDiffLoss <: SupervisedLoss end
LossFunctions.istwicedifferentiable(::TstTwiceDiffLoss) = true

@testset "Fallback implementations; twice differentiable" begin
    @test isminimizable(TstTwiceDiffLoss()) == false

    @test isdifferentiable(TstTwiceDiffLoss()) == true
    @test isdifferentiable(TstTwiceDiffLoss(), 0) == true
    @test istwicedifferentiable(TstTwiceDiffLoss()) == true
    @test istwicedifferentiable(TstTwiceDiffLoss(), 0) == true

    @test isstronglyconvex(TstTwiceDiffLoss()) == false
    @test isstrictlyconvex(TstTwiceDiffLoss()) == false
    @test isconvex(TstTwiceDiffLoss()) == false

    @test isnemitski(TstTwiceDiffLoss()) == false
    @test islipschitzcont(TstTwiceDiffLoss()) == false
    @test islocallylipschitzcont(TstTwiceDiffLoss()) == false
    @test isclipable(TstTwiceDiffLoss()) == false
    @test ismarginbased(TstTwiceDiffLoss()) == false
    @test isdistancebased(TstTwiceDiffLoss()) == false
    @test issymmetric(TstTwiceDiffLoss()) == false
end

struct TstMarginLoss <: MarginLoss end

@testset "Fallback implementations; margin-based" begin
    @test isminimizable(TstMarginLoss()) == false

    @test isdifferentiable(TstMarginLoss()) == false
    @test isdifferentiable(TstMarginLoss(), 0) == false
    @test istwicedifferentiable(TstMarginLoss()) == false
    @test istwicedifferentiable(TstMarginLoss(), 0) == false

    @test isstronglyconvex(TstMarginLoss()) == false
    @test isstrictlyconvex(TstMarginLoss()) == false
    @test isconvex(TstMarginLoss()) == false

    @test isnemitski(TstMarginLoss()) == true
    @test islipschitzcont(TstMarginLoss()) == false
    @test islocallylipschitzcont(TstMarginLoss()) == false
    @test isclipable(TstMarginLoss()) == false
    @test ismarginbased(TstMarginLoss()) == true
    @test isdistancebased(TstMarginLoss()) == false
    @test issymmetric(TstMarginLoss()) == false

    @test isclasscalibrated(TstMarginLoss()) == false
    @test isfishercons(TstMarginLoss()) == false
    @test isunivfishercons(TstMarginLoss()) == false
end

struct TstDistanceLoss <: DistanceLoss end

@testset "Fallback implementations; distance-based" begin
    @test isminimizable(TstDistanceLoss()) == false

    @test isdifferentiable(TstDistanceLoss()) == false
    @test isdifferentiable(TstDistanceLoss(), 0) == false
    @test istwicedifferentiable(TstDistanceLoss()) == false
    @test istwicedifferentiable(TstDistanceLoss(), 0) == false

    @test isstronglyconvex(TstDistanceLoss()) == false
    @test isstrictlyconvex(TstDistanceLoss()) == false
    @test isconvex(TstDistanceLoss()) == false

    @test isnemitski(TstDistanceLoss()) == false
    @test islipschitzcont(TstDistanceLoss()) == false
    @test islocallylipschitzcont(TstDistanceLoss()) == false
    @test isclipable(TstDistanceLoss()) == true
    @test ismarginbased(TstDistanceLoss()) == false
    @test isdistancebased(TstDistanceLoss()) == true
    @test issymmetric(TstDistanceLoss()) == false
end

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

@testset "LPDistLoss{0.5}" begin
    loss = LPDistLoss(0.5)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == false
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == false

    @test isnemitski(loss) == false
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == false
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "L1DistLoss" begin
    loss = L1DistLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == false
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "L2DistLoss" begin
    loss = L2DistLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true

    @test isstronglyconvex(loss) == true
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "LPDistLoss{3}" begin
    loss = LPDistLoss(3)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true

    @test isstronglyconvex(loss) == true
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "HuberLoss(1)" begin
    loss = HuberLoss(1)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "L1EpsilonInsLoss(1)" begin
    loss = L1EpsilonInsLoss(1)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == false
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "L2EpsilonInsLoss(1)" begin
    loss = L2EpsilonInsLoss(1)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false

    @test isstronglyconvex(loss) == true
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "LogitDistLoss()" begin
    loss = LogitDistLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == false
    @test isdistancebased(loss) == true
    @test issymmetric(loss) == true
end

@testset "QuantileLoss" begin
    l1 = QuantileLoss(.5)
    l2 = QuantileLoss(.7)

    @test issymmetric(l1) == true
    @test issymmetric(l2) == false

    @test isminimizable(l2) == true

    @test isdifferentiable(l2) == false
    @test isdifferentiable(l2, 0) == false
    @test isdifferentiable(l2, 1) == true
    @test istwicedifferentiable(l2) == false
    @test istwicedifferentiable(l2, 0) == false
    @test istwicedifferentiable(l2, 1) == true

    @test isstronglyconvex(l2) == false
    @test isstrictlyconvex(l2) == false
    @test isconvex(l2) == true

    # @test isnemitski(l2) == ?
    @test islipschitzcont(l2) == true
    @test islocallylipschitzcont(l2) == true
    # @test isclipable(l2) == ?
    @test ismarginbased(l2) == false
    @test isdistancebased(l2) == true
end

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

@testset "ZeroOneLoss" begin
    loss = ZeroOneLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == false
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == false

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
end

@testset "PerceptronLoss" begin
    loss = PerceptronLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == false
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == false
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == false
end

@testset "LogitMarginLoss" begin
    loss = LogitMarginLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == false
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

@testset "L1HingeLoss" begin
    loss = L1HingeLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == false
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == false
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == false
end

@testset "L2HingeLoss" begin
    loss = L2HingeLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

@testset "SmoothedL1HingeLoss" begin
    loss = SmoothedL1HingeLoss(2)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, -1) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false
    @test istwicedifferentiable(loss, 2) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
end

@testset "ModifiedHuberLoss" begin
    loss = ModifiedHuberLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == false
    @test istwicedifferentiable(loss, -1) == false
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == false
    @test istwicedifferentiable(loss, 2) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
end

@testset "L2MarginLoss" begin
    loss = L2MarginLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == true
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == true
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

@testset "ExpLoss" begin
    loss = ExpLoss()

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == true
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == false
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == false
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

@testset "SigmoidLoss" begin
    loss = SigmoidLoss()

    @test isminimizable(loss) == false

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == false

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == false
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

@testset "DWDMarginLoss" begin
    loss = DWDMarginLoss(2)

    @test isminimizable(loss) == true

    @test isdifferentiable(loss) == true
    @test isdifferentiable(loss, 0) == true
    @test isdifferentiable(loss, 1) == true
    @test istwicedifferentiable(loss) == true
    @test istwicedifferentiable(loss, 0) == true
    @test istwicedifferentiable(loss, 1) == true

    @test isstronglyconvex(loss) == false
    @test isstrictlyconvex(loss) == false
    @test isconvex(loss) == true

    @test isnemitski(loss) == true
    @test islipschitzcont(loss) == true
    @test islocallylipschitzcont(loss) == true
    @test isclipable(loss) == false
    @test ismarginbased(loss) == true
    @test isdistancebased(loss) == false
    @test issymmetric(loss) == false
    @test isclasscalibrated(loss) == true
    @test isfishercons(loss) == true
    @test isunivfishercons(loss) == true
end

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

function compare_losses(l1, l2, ccal = true)
    @test isminimizable(l1) == isminimizable(l2)

    @test isdifferentiable(l1) == isdifferentiable(l2)
    @test isdifferentiable(l1, 0) == isdifferentiable(l2, 0)
    @test isdifferentiable(l1, 1) == isdifferentiable(l2, 1)
    @test istwicedifferentiable(l1) == istwicedifferentiable(l2)
    @test istwicedifferentiable(l1, 0) == istwicedifferentiable(l2, 0)
    @test istwicedifferentiable(l1, 1) == istwicedifferentiable(l2, 1)
    @test istwicedifferentiable(l1, 2) == istwicedifferentiable(l2, 2)

    @test isstronglyconvex(l1) == isstronglyconvex(l2)
    @test isstrictlyconvex(l1) == isstrictlyconvex(l2)
    @test isconvex(l1) == isconvex(l2)

    @test isnemitski(l1) == isnemitski(l2)
    @test islipschitzcont(l1) == islipschitzcont(l2)
    @test islocallylipschitzcont(l1) == islocallylipschitzcont(l2)
    @test isclipable(l1) == isclipable(l2)
    @test ismarginbased(l1) == ismarginbased(l2)
    @test isdistancebased(l1) == isdistancebased(l2)
    @test issymmetric(l1) == issymmetric(l2)
    @test (ccal && isclasscalibrated(l1)) == isclasscalibrated(l2)
end

compare_losses(PoissonLoss(), 2*PoissonLoss())
compare_losses(PoissonLoss(), 0.5*PoissonLoss())

@testset "Weighted Margin-based" begin
    for loss in margin_losses
        @testset "$loss" begin
            compare_losses(loss, weightedloss(loss,0.2), false)
            compare_losses(loss, weightedloss(loss,0.5), true)
            compare_losses(loss, weightedloss(loss,0.7), false)
        end
    end
end

@testset "Scaled Margin-based" begin
    for loss in margin_losses
        @testset "$loss" begin
            compare_losses(loss, 2*loss)
            compare_losses(loss, 0.5*loss)
        end
    end
end

@testset "Scaled Distance-based" begin
    for loss in distance_losses
        @testset "$loss" begin
            compare_losses(loss, 2*loss)
            compare_losses(loss, 0.5*loss)
        end
    end
end


text:
using CoordinateTransformations, Rotations, TestImages, ImageCore, StaticArrays, OffsetArrays, Interpolations, LinearAlgebra
using Test, ReferenceTests

refambs = detect_ambiguities(CoordinateTransformations, Base, Core)
using ImageTransformations
ambs = detect_ambiguities(ImageTransformations, CoordinateTransformations, Base, Core)
@test isempty(setdiff(ambs, refambs))

function typestring(::Type{T}) where T   # from https://github.com/JuliaImages/ImageCore.jl/pull/133
    buf = IOBuffer()
    show(buf, T)
    String(take!(buf))
end

tests = [
    "autorange.jl",
    "resizing.jl",
    "interpolations.jl",
    "warp.jl",
]

@testset "ImageTransformations" begin
for t in tests
    @testset "$t" begin
        include(t)
    end
end
end

nothing

ext: jl


text:
using JuliaConSubmission
using Test

@testset "AbstractLogicJuliaConSubmission.jl" begin
    # Write your own tests here.
end


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

# Support for @simd for

module SimdLoop

export @simd, simd_outer_range, simd_inner_length, simd_index

# Error thrown from ill-formed uses of @simd
struct SimdError <: Exception
    msg::AbstractString
end

# Parse iteration space expression
#       symbol '=' range
#       symbol 'in' range
function parse_iteration_space(x)
    (isa(x, Expr) && (x.head == :(=) || x.head == :in)) || throw(SimdError("= or in expected"))
    length(x.args) == 2 || throw(SimdError("simd range syntax is wrong"))
    isa(x.args[1], Symbol) || throw(SimdError("simd loop index must be a symbol"))
    x.args # symbol, range
end

# reject invalid control flow statements in @simd loop body
function check_body!(x::Expr)
    if x.head === :break || x.head == :continue
        throw(SimdError("$(x.head) is not allowed inside a @simd loop body"))
    elseif x.head === :macrocall && x.args[1] === Symbol("@goto")
        throw(SimdError("$(x.args[1]) is not allowed inside a @simd loop body"))
    end
    for arg in x.args
        check_body!(arg)
    end
    return true
end
check_body!(x::QuoteNode) = check_body!(x.value)
check_body!(x) = true

# @simd splits a for loop into two loops: an outer scalar loop and
# an inner loop marked with :loopinfo. The simd_... functions define
# the splitting.
# Custom iterators that do not support random access cannot support
# vectorization. In order to be compatible with `@simd` annotated loops,
#they should override `simd_inner_length(v::MyIter, j) = 1`,
#`simd_outer_range(v::MyIter) = v`, and `simd_index(v::MyIter, j, i) = j`.

# Get range for outer loop.
simd_outer_range(r) = 0:0

# Get trip count for inner loop.
@inline simd_inner_length(r, j) = Base.length(r)

# Construct user-level element from original range, outer loop index j, and inner loop index i.
@inline simd_index(r, j, i) = (@inbounds ret = r[i+firstindex(r)]; ret)

# Compile Expr x in context of @simd.
function compile(x, ivdep)
    (isa(x, Expr) && x.head == :for) || throw(SimdError("for loop expected"))
    length(x.args) == 2 || throw(SimdError("1D for loop expected"))
    check_body!(x)

    var,range = parse_iteration_space(x.args[1])
    r = gensym("r") # Range value
    j = gensym("i") # Iteration variable for outer loop
    n = gensym("n") # Trip count for inner loop
    i = gensym("i") # Trip index for inner loop
    quote
        # Evaluate range value once, to enhance type and data flow analysis by optimizers.
        let $r = $range
            for $j in Base.simd_outer_range($r)
                let $n = Base.simd_inner_length($r,$j)
                    if zero($n) < $n
                        # Lower loop in way that seems to work best for LLVM 3.3 vectorizer.
                        let $i = zero($n)
                            while $i < $n
                                local $var = Base.simd_index($r,$j,$i)
                                $(x.args[2])        # Body of loop
                                $i += 1
                                $(Expr(:loopinfo, Symbol("julia.simdloop"), ivdep))  # Mark loop as SIMD loop
                            end
                        end
                    end
                end
            end
        end
        nothing
    end
end

"""
    @simd

Annotate a `for` loop to allow the compiler to take extra liberties to allow loop re-ordering

!!! warning
    This feature is experimental and could change or disappear in future versions of Julia.
    Incorrect use of the `@simd` macro may cause unexpected results.

The object iterated over in a `@simd for` loop should be a one-dimensional range.
By using `@simd`, you are asserting several properties of the loop:

* It is safe to execute iterations in arbitrary or overlapping order, with special consideration for reduction variables.
* Floating-point operations on reduction variables can be reordered, possibly causing different results than without `@simd`.

In many cases, Julia is able to automatically vectorize inner for loops without the use of `@simd`.
Using `@simd` gives the compiler a little extra leeway to make it possible in more situations. In
either case, your inner loop should have the following properties to allow vectorization:

* The loop must be an innermost loop
* The loop body must be straight-line code. Therefore, [`@inbounds`](@ref) is
    currently needed for all array accesses. The compiler can sometimes turn
    short `&&`, `||`, and `?:` expressions into straight-line code if it is safe
    to evaluate all operands unconditionally. Consider using the [`ifelse`](@ref)
    function instead of `?:` in the loop if it is safe to do so.
* Accesses must have a stride pattern and cannot be "gathers" (random-index
    reads) or "scatters" (random-index writes).
* The stride should be unit stride.

!!! note
    The `@simd` does not assert by default that the loop is completely free of loop-carried
    memory dependencies, which is an assumption that can easily be violated in generic code.
    If you are writing non-generic code, you can use `@simd ivdep for ... end` to also assert that:

* There exists no loop-carried memory dependencies
* No iteration ever waits on a previous iteration to make forward progress.
"""
macro simd(forloop)
    esc(compile(forloop, nothing))
end

macro simd(ivdep, forloop)
    if ivdep == :ivdep
        esc(compile(forloop, Symbol("julia.ivdep")))
    else
        throw(SimdError("Only ivdep is valid as the first argument to @simd"))
    end
end

end # module SimdLoop

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
# ComponentSol{T,N,C} = SciMLBase.AbstractODESolution{T,N,C} where C<:AbstractVector{<:ComponentArray}

# Plotting stuff
function SciMLBase.getsyms(sol::SciMLBase.AbstractODESolution{T,N,C}) where {T,N,C<:AbstractVector{<:ComponentArray}}
    if SciMLBase.has_syms(sol.prob.f)
        return sol.prob.f.syms
    else
        return Symbol.(labels(sol.u[1]))
    end
end

max_forks_repo_name: jonniedie/ComponentArrays


text:
# Second version of the new_model.jl
# Parameter iterators
postagv(m)=m[1];
bufferm(m)=m[2];
stackm(m)=m[3];
actm(m)=m[4];
actembed(m)=m[5];
parserv(m)=m[6];

# To initialize parser model from scratch
function initmodel1(o, s; ftype=Float32)
    model = Any[]

    # postag-level initialization, KnetArray initialization
    dpostag = o[:posembed] # 1
    for (k, n, d) in ((:postag, 17, dpostag),)
        push!(model, Any[ initr(d, GPUFEATS=true) for i=1:n ])
    end

    # buffer-lstm initialization # 2
    r_b, wr_b = initbuf(o[:bufembed], o[:bufhidden]) # you may change lstm - gru with kwargs
    push!(model, [r_b, wr_b])

    # stack-lstm initialization # 3
    r_s, wr_s = initbuf(o[:stembed], o[:sthidden])
    push!(model, [r_s, wr_s])

    # action-lstm initialization # 4
    r_a, wr_a = initbuf(o[:actembed], o[:acthidden])
    push!(model, [r_a, wr_a])

    # action embeddings
    p = o[:arctype](s)
    push!(model, initx(o[:actembed], p.nmove)) # 5

    # decision module initialization # 6
    mlpdims = (r_s.hiddenSize + r_b.hiddenSize + r_a.hiddenSize, o[:hidden]..., p.nmove)
    info("mlpdims: $mlpdims")
    decider = Any[]
    for i=2:length(mlpdims)
        push!(decider, initx(mlpdims[i], mlpdims[i-1]))
        push!(decider, initx(mlpdims[i], 1))
    end
    push!(model, decider)
    optims = optimizers(model, o[:optimization])
    return model, optims
end

function mlp(w,x; pdrop=(0,0))
    x = dropout(x,pdrop[1])
    for i=1:2:length(w)-2
        x = relu.(w[i]*x .+ w[i+1])
        x = dropout(x,pdrop[2])
    end
    return w[end-1]*x .+ w[end]
end

###### Buffer-related
# to get the correct index of a matrix
revptr(p::Parser{:ArcHybrid,:R1}) = (p.wptr> p.nword ? -1 : length(p.sentence) - p.wptr + 1)

#### yalls = [wn wn-1.... w3 w2] where we are going to use the reversed version
function scan_buffer3(model, parser, hidden)
    if length(parser.sentence) == 1 # ugly test time!
        t = zeros(Float32, 1078, 1)
        yalls = (gpu()>=0 ? KnetArray(t): Array(t))
        return yalls
    end

    bufmodel, pvecs = bufferm(model), postagv(model)
    rbuf, wbuf = bufmodel[1], bufmodel[2]
    sentence = parser.sentence
    seqe  = length(sentence)
    seqs = parser.wptr
    range_ = seqe:-1:seqs
    xs = Any[]
    for i in range_
        x  = (gpu() >=0 ? KnetArray(sentence.cavec[i]) : Array(sentence.cavec[i]))
        i1 = cat1d(x, pvecs[sentence.postag[i]])
        push!(xs, i1)
    end

    allin = reshape(hcat(xs...), 1078, 1, length(xs))
    yalls, = rnnforw(rbuf, wbuf, allin)
    return reshape(yalls, hidden, length(xs))
end

function cache_bufbatch(model, parsers, hidden)
    # TODO: we may .cache more in Array{Any}

    yalls = Array{Any}(length(parsers))
    B = length(parsers)
    for i in 1:B
        p = parsers[i]
        y1 = scan_buffer3(model, p, hidden)
        yalls[i] = y1
    end
    return yalls
end

function scan_bufbatch(yalls, parsers, hidden)
    #ncols = map(revptr, parsers) # batchsize
    B = length(parsers)
    instate = zeros(Float32, hidden, 1)
    ybuf = (gpu() >=0 ? KnetArray(instate) : Array(instate))

    ybatch = Any[]
    for i in 1:B
        cn = revptr(parsers[i]) # column number
        if cn == -1
            push!(ybatch, ybuf) # zero state for others
        else
            push!(ybatch, yalls[i][:, cn])
        end
    end
    ybatch = cat1d(ybatch...)
    nrow   = div(length(ybatch), B)
    return reshape(ybatch, nrow, B)

end

####### Stack-related

function scan_stack(model, parser, hstack)
    stmodel, pvecs = stackm(model), postagv(model)
    rst, wst = stmodel[1], stmodel[2]
    sentence = parser.sentence
    seqe = parser.sptr
    instate = zeros(Float32, hstack, 1)
    yst = hout = cout = (gpu() >= 0 ? KnetArray(instate) : Array(instate))
    range_ = 1:seqe # words in stack p.stack[1:p.sptr]

    if length(range_) < 1
        return yst
    end
    inputs = Any[]
    for i in range_
        indx  = parser.stack[i]
        #word = sentence.word[indx]; print(" $word") # dbg line
        x = (gpu() >= 0 ? KnetArray(sentence.cavec[indx]) : Array(sentence.cavec[indx]))
        input = cat1d(x, pvecs[sentence.postag[indx]])
        push!(inputs, input)
        #yst, hout, cout, _ = rnnforw(rst, wst, input, hout, cout, hy=true, cy=true)
    end
    in2rnn = reshape(hcat(inputs...), 1078, 1, length(inputs))
    _, hstck, = rnnforw(rst, wst, in2rnn, hy=true)
    return reshape(hstck, hstack, 1)
end

function scan_stackbatch(model, parsers, hstack)
    yalls = Any[]
    for p in parsers
        yst = scan_stack(model, p, hstack)
        push!(yalls, yst)
    end
    yalls = cat1d(yalls...)
    ncols = length(parsers)
    nrows = div(length(yalls), ncols)
    return reshape(yalls, nrows, ncols)
end

###### Action related

# compute batched actions
function scan_action(model, actions, state)
    acemb, actmodel = actembed(model), actm(model)
    yout0, hout0, cout0 = state[1], state[2], state[3]
    input = acemb[:, actions]
    yout, hout, cout = rnnforw(actmodel[1], actmodel[2], input, hout0, cout0, hy=true, cy=true)
    return (yout, hout, cout)
end

function oracleloss(allmodel, sentences, arctype, hiddens; losses=nothing, pdrop=(0.0, 0.0))
    hidstack, hidact, hidbuf = hiddens
    parsers     = map(arctype, sentences)
    mcosts      = Array{Cost}(parsers[1].nmove)
    parserdone  = falses(length(parsers))
    cachedbufs  = cache_bufbatch(allmodel, parsers, hidbuf)
    mlpmodel = parserv(allmodel)

    B = length(parsers) # bathcsize
    t0 = zeros(Float32, hidact, B, 1)
    h0 = c0 = (gpu()>=0 ? KnetArray(t0) : Array(t0))
    y0 = reshape(h0, hidact, B)
    totalloss = 0.0

    # action lstm-init
    nstate = (y0, h0, c0);
    inacts = map(t->1, 1:B)

    while !all(parserdone)
        # batched score
        yact, hact, cact = scan_action(allmodel, inacts, nstate)
        nstate = (yact, hact, cact);
        ybufs  = scan_bufbatch(cachedbufs, parsers, hidbuf);
        ysts     = scan_stackbatch(allmodel, parsers, hidstack)
        encoded  = vcat(yact, ybufs, ysts)
        scores   = mlp(mlpmodel, encoded, pdrop=pdrop)
        logprobs = logp(scores, 1)

        inacts = Int[] # to get the next step's actions
        # iterative loss-val
        for (i, p) in enumerate(parsers)
            if parserdone[i]
                push!(inacts, 1); # not changing batchsize
                continue
            end
            movecosts(p, p.sentence.head, p.sentence.deprel, mcosts)
            goldmove = indmin(mcosts)

            if mcosts[goldmove] == typemax(Cost)
                parserdone[i] = true
                p.sentence.parse = p
                push!(inacts, 1); # not changing batchsize
            else
                totalloss -= logprobs[goldmove, i]
                move!(p, goldmove);
                push!(inacts, goldmove);
                if losses != nothing
                    loss1 = -getval(logprobs)[goldmove,i]
                    losses[1] += loss1
                    losses[2] += 1
                    if losses[2] < 1000
                        losses[3] = losses[1]/losses[2]
                    else
                        losses[3] = 0.999 * losses[3] + 0.001 * loss1
                    end
                end
            end
        end
    end
    return totalloss / length(sentences)
end

oraclegrad = grad(oracleloss)

# Implement test related

# endofparse to use in test
endofparse(p)=(p.sptr == 1 && p.wptr > p.nword)

function oracletest(allmodel, corpus, arctype, hiddens, batchsize; pdrop=(0.0, 0.0))
    sentbatches = minibatch(corpus, batchsize)
    hidstack, hidact, hidbuf = hiddens
    mlpmodel = parserv(allmodel)

    for sentences in sentbatches
        parsers = map(arctype, sentences)
        B = length(parsers) # bathcsize
        mcosts  = Array{Cost}(parsers[1].nmove)
        parserdone = map(endofparse, parsers)
        cachedbufs  = cache_bufbatch(allmodel, parsers, hidbuf)
        #parserdone = falses(length(parsers))
        t0 = zeros(Float32, hidact, B, 1)
        h0 = c0 = (gpu()>=0 ? KnetArray(t0) : Array(t0))
        y0 = reshape(h0, hidact, B)
        nstate = (y0, h0, c0);
        inacts = map(t->1, 1:B)
        yact, hact, cact = scan_action(allmodel, inacts, nstate)

        if all(parserdone) # shows how careless they prepared the dataset!
            for p in parsers
                p.sentence.parse = p
            end
        end

        while !all(parserdone)
            yact, hact, cact = scan_action(allmodel, inacts, nstate)
            nstate = (yact, hact, cact);
            ybufs    = scan_bufbatch(cachedbufs, parsers, hidbuf)
            ysts     = scan_stackbatch(allmodel, parsers, hidstack)
            encoded  = vcat(yact, ybufs, ysts)
            scores   = mlp(mlpmodel, encoded, pdrop=pdrop)
            logprobs = Array(logp(scores, 1))

            inacts = Int[] # to get the next step's actions
            for (i, p) in enumerate(parsers)
                if parserdone[i]
                    if p.sentence.parse == nothing # to avoid from the first time
                        p.sentence.parse = p
                    end
                    push!(inacts, 1) # not to break batchsize
                    continue
                end
                isorted = sortperm(logprobs[:,i], rev=true) # ith col ith inst.
                for m in isorted # best&val action
                    if moveok(p, m)
                        move!(p, m)
                        push!(inacts, m) # not to break batchsize
                        break
                    end
                end
                if endofparse(p)
                    parserdone[i] = true
                    p.sentence.parse = p
                end
            end
        end
    end
    return las(corpus)
end


text:
using Flux
using Random

function get_batch(dataset, max_sequence::Integer, batch_size::Integer; set::Symbol = :train, max_regions=nothing)
    data = shuffle(dataset[set])

    # Positional data
    word_size = size(data[1])[1]
    region_length = max_sequence
    num_ins = size(data[1])[3]

    # Chop songs up along sequence length
    for j in 1:length(data)
        num_regions = Integer(ceil(size(data[j])[2] / max_sequence))
        if max_regions != nothing
            num_regions = min(num_regions, max_regions)
        end

        data[j] = [data[j][:,max(1, min(i, size(data[j])[2] - max_sequence)):i+min(max_sequence-1, size(data[j])[2]-i-1),:]
                    for i in 1:max_sequence:num_regions*max_sequence]
    end

    # Batch songs
    num_batches = Integer(ceil(length(data) / batch_size))
    data = [data[i:i+min(batch_size-1, length(data)-i)]
                for i in 1:batch_size:batch_size*num_batches]

    return data
end

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
using MAT, Plots

@views function generate_grf2D(lx=100.0, ly=100.0, sf=1.0, cl=(10.0, 8.0), k_m=100.0, nh=10000, nx=64, ny=64;
                               cov_typ="expon", do_reset=true, do_viz=true, do_save=false)
    # Derived numerics
    dx, dy  = lx/nx, ly/ny  # numerical grid step size
    me = 0
    co1, co2 = 0, 0
    # Array allocation
    Yf      = @zeros(nx, ny)
    # Visu init
    if do_viz
        X, Y = -lx/2:dx:lx/2, -ly/2:dy:ly/2
    end

    if cov_typ=="expon"
        # Generate the 2D exponential covariance function
        wtime_it = grf2D_expon!(Yf, sf, cl, nh, nx, ny, dx, dy; me, co1, co2, do_reset)
    elseif cov_typ=="gauss"
        # Generate the 2D Gaussian covariance function
        wtime_it = grf2D_gauss!(Yf, sf, cl[1]/2.0, nh, k_m, nx, ny, dx, dy; me, co1, co2, do_reset)
    else
        error("Undefined covariance function.")
    end

    # Performance
    A_eff    = 2/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    T_eff    = A_eff/wtime_it                   # Effective memory throughput [GB/s]
    if (me==0) @printf("T_eff = %1.2f GB/s \n", round(T_eff, sigdigits=3)) end

    # Visualisation
    if do_viz
        display(heatmap(X, Y, Array(Yf)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:hot, title="2D RandomField"))
    end
    if do_save  file = matopen("grf2D_$(cov_typ).mat", "w"); write(file, "grf2D", Array(Yf)); close(file)  end

    return Yf
end


text:
# Just make sure they run for now.
module TestHeatEquationExample
    mktempdir() do dir
        cd(dir) do
            include(joinpath(@__DIR__, "..","docs", "src", "examples", "heat_equation.jl"))
        end
    end
end


text:
module PrettyNumbers

using Printf

################################################################################
#                                   Includes
################################################################################

include("./constants.jl")
include("./misc.jl")
include("./print.jl")

include("./backends/text/print.jl")
include("./backends/text/render.jl")

end # module


text:
"""
Parsing routines for BioNetGen .net files. Note, this only handles a subset
of the BioNetGen .net file format. In particular, any functions in the file
that use non-Julia expressions (i.e. like conditional logic) or time-dependent
features, will definitely not work.
"""

"""
Seek to a line matching `start_string`.
"""
function seek_to_block(lines, idx, start_string)
    while lines[idx] != start_string
        idx += 1
        (idx > length(lines)) && error("Block: ", start_string, " was never found.")
    end
    idx += 1
end

# for swapping out BioNetGen symbols
const REPLACEMENT_DICT = Dict( :ln => :log )
function recursive_replace!(expr::Any, replace_requests::Dict{Symbol,Symbol}=REPLACEMENT_DICT)
    if typeof(expr) == Symbol
        haskey(replace_requests,expr) && return replace_requests[expr]
    elseif typeof(expr) == Expr
        for i = 1:length(expr.args)
            expr.args[i] = recursive_replace!(expr.args[i], replace_requests)
        end
    end
    return expr
end

const PARAM_BLOCK_START = "begin parameters"
const PARAM_BLOCK_END = "end parameters"
function parse_params(ft::BNGNetwork, lines, idx)

    idx = seek_to_block(lines, idx, PARAM_BLOCK_START)

    # parse params
    pvals = []
    ptoids = OrderedDict{Symbol,Int}()
    psyms = OrderedSet{Symbol}()
    while lines[idx] != PARAM_BLOCK_END
        vals = split(lines[idx])
        pidx = parse(Int,vals[1])
        psym = Symbol(vals[2])
        ptoids[psym] = pidx

        # value could be an expression
        pval = Meta.parse(vals[3])
        (typeof(pval) <: Expr) && recursive_replace!(pval)

        # BNG Constant is a number, so save sym
        (vals[end] == "Constant") && push!(psyms, psym)

        push!(pvals,pval)
        idx += 1
        (idx > length(lines)) && error("Block: ", PARAM_BLOCK_END, " was never found.")
    end

    ptoids,pvals,psyms,idx
end

const CHARS_TO_STRIP = ",~!()."
stripchars = (s) -> replace(s, Regex("[$CHARS_TO_STRIP]") => "")
hasstripchars = (s) -> occursin(Regex("[$CHARS_TO_STRIP]"), s)
const SPECIES_BLOCK_START = "begin species"
const SPECIES_BLOCK_END = "end species"
const MAX_SYM_LEN = 8
function parse_species(ft::BNGNetwork, lines, idx, ptoids, pvals)

    idx = seek_to_block(lines, idx, SPECIES_BLOCK_START)

    # parse species
    shortsymstoids = Dict{Symbol,Int}()
    shortsymstosyms = Dict{Symbol,Symbol}()
    u0exprs   = Vector{Any}()
    while lines[idx] != SPECIES_BLOCK_END
        vals = split(lines[idx])
        sidx = parse(Int,vals[1])
        sym = Symbol(vals[2])
        shortsym = length(vals[2]) < MAX_SYM_LEN ? Symbol(stripchars(vals[2])) : Symbol(string("S$sidx"))
        shortsymstoids[shortsym] = sidx
        shortsymstosyms[shortsym] = sym
        push!(u0exprs,  Meta.parse(vals[3]))
        idx += 1
        (idx > length(lines)) && error("Block: ", SPECIES_BLOCK_END, " was never found.")
    end

    u0exprs,shortsymstoids,shortsymstosyms,idx
end

const REACTIONS_BLOCK_START = "begin reactions"
const REACTIONS_BLOCK_END = "end reactions"
function parse_reactions!(rxiobuf, ft, lines, idx, ptoids, pvals, psyms, symstoids, idstosyms)

    # map from species id to shortsym
    idtosymstr = id -> (id==0) ? ∅ : string(idstosyms[id])

    idx = seek_to_block(lines, idx, REACTIONS_BLOCK_START)
    write(rxiobuf, "begin\n")
    cntdict = Dict{Int,Int}()
    while lines[idx] != REACTIONS_BLOCK_END
        vals        = split(lines[idx])
        reactantids = (parse(Int,rid) for rid in split(vals[2],","))
        productids  = (parse(Int,pid) for pid in split(vals[3],","))
        rateexpr    = Meta.parse(vals[4])
        (typeof(rateexpr) <: Expr) && recursive_replace!(rateexpr)
        reactstr    = join((idtosymstr(rid) for rid in reactantids), " + ")
        productstr  = join((idtosymstr(pid) for pid in productids), " + ")

        # correct for higher-order rate rescalings by BioNetGen
        empty!(cntdict)
        foreach(rid -> haskey(cntdict,rid) ? (cntdict[rid] += 1) : (cntdict[rid]=1), reactantids)
        scalefactor = prod(factorial, values(cntdict))
        pstr = string(scalefactor, "*", string(rateexpr))

        # create string for this reaction
        write(rxiobuf, "  ", pstr, ", ", reactstr, " --> ", productstr, "\n")

        idx += 1
        (idx > length(lines)) && error("Block: ", REACTIONS_BLOCK_END, " was never found.")
    end
    write(rxiobuf, "end")

    idx
end

const GROUPS_BLOCK_START = "begin groups"
const GROUPS_BLOCK_END = "end groups"
function parse_groups(ft::BNGNetwork, lines, idx, idstoshortsyms, rn)

    idx = seek_to_block(lines, idx, GROUPS_BLOCK_START)
    namestoids = Dict{Symbol,Vector{Int}}()
    while lines[idx] != GROUPS_BLOCK_END
        vals = split(lines[idx])
        name = Symbol(vals[2])

        # map from BioNetGen id to reaction_network id
        ids = [rn.syms_to_ints[idstoshortsyms[parse(Int,val)]] for val in split(vals[3],",")]

        namestoids[name] = ids
        idx += 1
        (idx > length(lines)) && error("Block: ", GROUPS_BLOCK_END, " was never found.")
    end

    namestoids,idx
end

function exprs_to_nums(ptoids, pvals, u0exprs)
    p    = zeros(Float64, length(ptoids))
    pmod = Module()
    for (psym,pid) in ptoids
        p[pid] = Base.eval(pmod, :($psym = $(pvals[pid])))
    end

    u0 = zeros(Float64, length(u0exprs))
    for (i,u0expr) in enumerate(u0exprs)
        u0[i] = Base.eval(pmod, :($u0expr))
    end

    p,u0
end

# for parsing a subset of the BioNetGen .net file format
function loadrxnetwork(ft::BNGNetwork, networkname, rxfilename; kwargs...)

    file  = open(rxfilename, "r");
    lines = readlines(file)
    idx   = 1
    print("Parsing parameters...")
    ptoids,pvals,psyms,idx = parse_params(ft, lines, idx)
    println("done")

    print("Parsing species...")
    u0exprs,shortsymstoids,shortsymstosyms,idx = parse_species(ft, lines, idx, ptoids, pvals)
    println("done")

    # map from species id to short sym
    idstoshortsyms = Vector{Symbol}(undef,length(shortsymstoids))
    for (k,v) in shortsymstoids
        idstoshortsyms[v] = k
    end

    print("Parsing reactions...")
    rxiobuf = IOBuffer()
    write(rxiobuf, string("@min_reaction_network ", networkname, " "))
    idx = parse_reactions!(rxiobuf, ft, lines, idx, ptoids, pvals, psyms, shortsymstoids, idstoshortsyms)
    println("done")

    # finish the min_reaction_network string
    foreach(psym -> write(rxiobuf, " ", string(psym)), keys(ptoids))
    write(rxiobuf, "\n")
    rxstrs = String(take!(rxiobuf))

    # build the DiffEqBiological representation of the network
    print("Building network...")
    rn = eval(Meta.parse(rxstrs))
    println("done")

    print("Parsing groups...")
    groupstoids,idx = parse_groups(ft, lines, idx, idstoshortsyms, rn)
    println("done")
    close(file)

    # get numeric values for parameters and u₀
    p,u0 = exprs_to_nums(ptoids, pvals, u0exprs)

    # reorder species to reaction_network ordering
    u₀ = similar(u0)
    for i in eachindex(u0)
        u₀[rn.syms_to_ints[idstoshortsyms[i]]] = u0[i]
    end

    ParsedReactionNetwork(rn, u₀; p = p,
                                  paramexprs = pvals,
                                  symstonames = shortsymstosyms,
                                  groupstoids = groupstoids,
                                  rnstr = rxstrs)
end

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
# Downsampling a PointCloud using a VoxelGrid filter
# http://www.pointclouds.org/documentation/tutorials/voxel_grid.php

using PCLCommon
using PCLIO
using PCLFilters
using PCLVisualization

pcd_file = Pkg.dir("PCLIO", "test", "data", "table_scene_lms400.pcd")
cloud = PointCloud{PointXYZI}(pcd_file)
cloud_filtered = PointCloud{PointXYZI}()

sor = VoxelGrid{PointXYZI}()
setInputCloud(sor, cloud)
setLeafSize(sor, 0.01, 0.01, 0.01)
filter(sor, cloud_filtered)

println("PointCloud before filtering: $(length(cloud)) data points")
println("PointCloud before filtering: $(length(cloud_filtered)) data points")

if isdefined(:vis) && vis
    viewer = PCLVisualizer("pcl visualizer")
    red_handler = PointCloudColorHandlerCustom(cloud, 255, 0, 0)
    green_handler = PointCloudColorHandlerCustom(cloud, 0, 255, 0)
    addPointCloud(viewer, cloud, red_handler, id="cloud")
    addPointCloud(viewer, cloud_filtered, green_handler, id="cloud_filtered")
    spin(viewer)
    close(viewer)
end


text:
using Test
using Cleaner: CleanTable, names, polish_names!, polish_names, generate_polished_names
using DataFrames: DataFrame

@testset "name polishing is working" begin
    testDF = DataFrame(
        "  _aName with_loTsOfProblems" => [1, 2, 3],
        "  _aName with_loTsOfProblems1" => [4, 5, 6],
        "  _aName with_loTsOfProblems2" => [7, 8, 9],
    )
    testCT = CleanTable(testDF)

    @test names(polish_names!(testCT; style=:snake_case)) == Vector{Symbol}([
        :a_name_with_lo_ts_of_problems,
        :a_name_with_lo_ts_of_problems1,
        :a_name_with_lo_ts_of_problems2,
    ])

    @test names(polish_names!(testCT; style=:camelCase)) == Vector{Symbol}([
        :aNameWithLoTsOfProblems, :aNameWithLoTsOfProblems1, :aNameWithLoTsOfProblems2
    ])

    @test polish_names(testDF) isa CleanTable

    @test generate_polished_names(
        [
            "  _aName with_loTsOfProblems",
            "  _aName with_loTsOfProblems",
            "  _aName with_loTsOfProblems_1",
        ];
        style=:snake_case,
    ) == Vector{Symbol}([
        :a_name_with_lo_ts_of_problems,
        :a_name_with_lo_ts_of_problems_1,
        :a_name_with_lo_ts_of_problems_1_1,
    ])

    @test generate_polished_names(
        [
            "  _aName with_loTsOfProblems",
            "  _aName with_loTsOfProblems",
            "  _aName with_loTsOfProblems_1",
        ];
        style=:camelCase,
    ) == Vector{Symbol}([
        :aNameWithLoTsOfProblems, :aNameWithLoTsOfProblems_1, :aNameWithLoTsOfProblems1
    ])

    let err = nothing
        try
            generate_polished_names(["something"]; style=:any)
        catch err
        end
        @test err isa Exception
        @test sprint(showerror, err) ==
            "Invalid style selected. Options are :snake_case, :camelCase"
    end
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
module KernelSpectralClustering
# core files
include("KSCdata.jl")
include("KSCselect.jl")
include("KSCmodel.jl")
include("KSCevaluation.jl")
# algorithm implementation files
include("KSCnet.jl")
include("KSCbignet.jl")

export runKSCnetwork, rerunKSCnetwork
export kscnet, kscbignet
export coverage, coverageApprox, modularity, modularityApprox
export convertDSVNet, createMirrored

function runKSCnetwork(networkfile::String, delimiter::Char, mink::Int, maxk::Int; sizeMB=4000::Int, undirected=true::Bool, eval=false::Bool)
  filedir = dirname(networkfile)
  filename = splitext(basename(networkfile))[1]
  if !undirected
    createMirrored(networkfile, delimiter, "$(filedir)/$(filename)_mirror.txt")
    networkfile = "$(filedir)/$(filename)_mirror.txt"
  end
  if sizeMB > 0
    # kscbignet algo
    qTest, baf, k, firstNode, lastNode, trainSubset, numFiles, fileWeighted = @time kscbignet(networkfile, filedir, filename, delimiter, sizeMB; fraction=0.15, method="furs", decomp="eigs", minK=mink, maxK=maxk, convertF=true)
    println("$(numFiles) files with training subset size $(length(trainSubset)) and BAF value $(baf[k-mink+1]) for $(k) clusters")
    if eval
      println("Evaluation step")
      metricsfile = open("$(filedir)/$(filename)/metricsKSC.txt", "w")
      # Coverage
      tic()
      coverageM, totalNodes = coverageApprox(numFiles, trainSubset, firstNode, lastNode)
      timeCoverage = toq()
      println("Coverage $(coverageM) finished, time elapsed $(timeCoverage)s")
      write(metricsfile, "Coverage: value $(coverageM) for ratio $(length(trainSubset)) / $(totalNodes), time $(timeCoverage)s\n")
      # Modularity
      tic()
      modularityM = modularityApprox(numFiles, k, fileWeighted)
      timeModularity = toq()
      println("Modularity $(modularityM) finished, time elapsed $(timeModularity)s")
      write(metricsfile, "Modularity: value $(modularityM) for $(k) communities, time $(timeModularity)s\n")
      close(metricsfile)
    end
    return qTest, baf, k, firstNode, lastNode, trainSubset, numFiles, fileWeighted
  else
    # kscnet algo
    result = @time kscnet()
    return result
  end
end

function rerunKSCnetwork(networkfile::String, delimiter::Char, mink::Int, maxk::Int; sizeMB=4000::Int, undirected=true::Bool, eval=false::Bool)
  filedir = dirname(networkfile)
  filename = splitext(basename(networkfile))[1]
  if !undirected
    networkfile = "$(filedir)$(filename)_mirror.txt"
  end
  if sizeMB > 0
    # kscbignet algo
    qTest, baf, k, firstNode, lastNode, trainSubset, numFiles, fileWeighted = @time kscbignet(networkfile, filedir, filename, delimiter, sizeMB; fraction=0.15, method="furs", decomp="eigs", minK=mink, maxK=maxk, convertF=false)
    println("KSCbignet: $(numFiles) files with training subset size $(length(trainSubset)) and BAF value $(baf[k-mink+1]) for $(k) clusters")
    if eval
      println("Evaluation step")
      metricsfile = open("$(filedir)/$(filename)/metricsKSC.txt", "w")
      # Coverage
      tic()
      coverageM, totalNodes = coverageApprox(numFiles, trainSubset, firstNode, lastNode)
      timeCoverage = toq()
      println("Coverage $(coverageM) finished ratio $(length(trainSubset)) / $(totalNodes), time elapsed $(timeCoverage)s")
      write(metricsfile, "Coverage: value $(coverageM) for ratio $(length(trainSubset)) / $(totalNodes), time $(timeCoverage)s\n")
      # Modularity
      tic()
      modularityM = modularityApprox(numFiles, k, fileWeighted)
      timeModularity = toq()
      println("Modularity $(modularityM) finished, time elapsed $(timeModularity)s")
      write(metricsfile, "Modularity: value $(modularityM) for $(k) communities, time $(timeModularity)s\n")
      close(metricsfile)
    end
    return qTest, baf, k, firstNode, lastNode, trainSubset, numFiles, fileWeighted
  else
    # kscnet algo
    result = @time kscnet()
    return result
  end
end
end # module


text:
# Autogenerated wrapper script for libsass_jll for aarch64-linux-musl-cxx11
export libsass_so

JLLWrappers.@generate_wrapper_header("libsass")
JLLWrappers.@declare_library_product(libsass_so, "libsass.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsass_so,
        "lib/libsass.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function initjumpmodel(options)
    @unpack solver, threads, showsolverlog = options
    if solver == :cplex
    # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.6.0/ilog.odms.studio.help/CPLEX/ReleaseNotes/topics/releasenotes126/newParameterNames.html
    # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.8.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/introListAlpha.html
    # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.8.0/ilog.odms.cplex.help/CPLEX/UsrMan/topics/cont_optim/barrier/19_tuning_title_synopsis.html
        m = Model(solver=CplexSolver(CPXPARAM_LPMethod=4, CPXPARAM_Threads=threads, CPXPARAM_Read_Scale=1,
                    CPXPARAM_Preprocessing_Dual=-1, CPXPARAM_SolutionType=2,
                    CPXPARAM_ScreenOutput=Int(showsolverlog)))
                    # CPXPARAM_Emphasis_Numerical=1,
                    # CPXPARAM_Barrier_ConvergeTol=1e-8 (default), CPXPARAM_Barrier_Display=2,
                    # CPXPARAM_Barrier_Crossover=-1 or 2, CPXPARAM_Emphasis_Memory=1,
                    # Tom: CPXPARAM_Barrier_ColNonzeros=200, CPXPARAM_Preprocessing_Fill=0 (doesn't help)

    elseif solver == :gurobi
        m = Model(solver=GurobiSolver(Method=2, Threads=threads, BarConvTol=1.5e-8, Crossover=0))
                # BarConvTol=1e-8, NumericFocus=2 (0-3), ScaleFlag=3 (-1-3), Crossover=0
    elseif solver == :glpk
        m = Model(solver=GLPKSolverLP(method=:InteriorPoint, msg_lev=GLPK.MSG_ON))
    elseif solver == :clp
        m = Model(solver=ClpSolver(SolveType=4, LogLevel=4))
    else
        println("No solver installed called $solver.")
    end
    return m
end

function makevariables(m, sets)
    @unpack REGION, FUEL, TECH, CLASS, STORAGECLASS, HOUR, techtype = sets

    storagetechs = [k for k in TECH if techtype[k] == :storage]

    @variables m begin
        Systemcost[r in REGION]                                                         # M€/year
        CO2emissions[r in REGION]                                                       # kton CO2/year
        FuelUse[r in REGION, f in FUEL] >= 0                                            # GWh fuel/year
        Electricity[r in REGION, k in TECH, c in CLASS[k], h in HOUR] >= 0              # GWh elec/period
        AnnualGeneration[r in REGION, k in TECH] >= 0                                   # GWh elec/year
        Charging[r in REGION, k in storagetechs, h in HOUR] >= 0       # GWh elec/period (electricity used to charge)
        StorageLevel[r in REGION, k in storagetechs, c in STORAGECLASS[k], h in HOUR] >= 0  # TWh elec (in storage)
        Transmission[r1 in REGION, r2 in REGION, h in HOUR] >= 0                        # GWh elec/period
        TransmissionCapacity[r1 in REGION, r2 in REGION] >= 0                           # GW elec
        Capacity[r in REGION, k in TECH, c in CLASS[k]] >= 0                            # GW elec
        SolarCapacity[r in REGION, k in [:pv, :csp], pv in CLASS[:pv], csp in CLASS[:csp]] >= 0     # GW elec
    end #variables

    return Vars(Systemcost, CO2emissions, FuelUse, Electricity, AnnualGeneration, Charging, StorageLevel,
                    Transmission, TransmissionCapacity, Capacity, SolarCapacity)
end

function setbounds(sets, params, vars, options)
    @unpack REGION, TECH, CLASS, techtype = sets
    @unpack Capacity, TransmissionCapacity = vars
    @unpack classlimits, hydrocapacity, transmissionislands, demand = params
    @unpack hydroinvestmentsallowed, nuclearallowed, transmissionallowed, disabletechs = options
    for r in REGION, k in TECH
        if techtype[k] == :vre || k == :csp
            for c in CLASS[k]
                setupperbound(Capacity[r,k,c], classlimits[r,k,c])
            end
        end
    end
    for r in REGION
        setlowerbound(Capacity[r,:hydro,:x0], hydrocapacity[r,:x0])
        for c in CLASS[:hydro]
            if hydroinvestmentsallowed
                setupperbound(Capacity[r,:hydro,c], hydrocapacity[r,c])
            else
                setupperbound(Capacity[r,:hydro,c], c == :x0 ? hydrocapacity[r,c] : 0.0)
            end
        end
    end
    if !nuclearallowed
        for r in REGION
            setupperbound(Capacity[r,:nuclear,:_], 0.0)
        end
    end
    for r1 in REGION, r2 in REGION
        if transmissionallowed == :none || (transmissionallowed == :islands && !transmissionislands[r1,r2])
            setupperbound(TransmissionCapacity[r1,r2], 0.0)
        end
    end
    for r in REGION, k in disabletechs, c in CLASS[k]
        setupperbound(Capacity[r,k,c], 0.0)
    end
    # for r in REGION
    #   setupperbound(AnnualGeneration[r,:csp], 0.50 * sum(demand[r,h] for h in HOUR) * hoursperperiod)
    # end
end

function makeconstraints(m, sets, params, vars, hourinfo, options)
    @unpack REGION, FUEL, TECH, CLASS, STORAGECLASS, HOUR, techtype, techfuel, reservoirclass = sets
    @unpack cf, transmissionlosses, demand, cfhydroinflow, efficiency, rampingrate, dischargetime, initialstoragelevel,
            minflow_existinghydro, emissionsCO2, fuelcost, variablecost, smalltransmissionpenalty, investcost, crf, fixedcost,
            transmissioninvestcost, transmissionfixedcost, hydroeleccost, solarcombinedarea,
            pv_density, csp_density, cspsolarmultiple = params
    @unpack Systemcost, CO2emissions, FuelUse, Electricity, AnnualGeneration, Charging, StorageLevel,
            Transmission, TransmissionCapacity, Capacity, SolarCapacity = vars
    @unpack hoursperperiod = hourinfo
    @unpack carbontax, carboncap, rampingconstraints, maxbioenergy, globalnuclearlimit = options

    storagetechs = [k for k in TECH if techtype[k] == :storage]

    @constraints m begin
        ElecCapacity[r in REGION, k in TECH, c in CLASS[k], h in HOUR],
            Electricity[r,k,c,h] <= Capacity[r,k,c] * (k == :csp ? 1 : cf[r,k,c,h]) * hoursperperiod

        ElecDemand[r in REGION, h in HOUR],
            sum(Electricity[r,k,c,h] for k in TECH, c in CLASS[k]) - sum(Charging[r,k,h] for k in TECH if techtype[k] == :storage) +
                + sum((1-transmissionlosses[r2,r])*Transmission[r2,r,h] - Transmission[r,r2,h] for r2 in REGION) >=
                    demand[r,h] * hoursperperiod

        # <= instead of == to avoid need of slack variable to deal with spillage during spring floods, etc
        StorageBalance[r in REGION, k in storagetechs, sc in STORAGECLASS[k], h in HOUR],
            (StorageLevel[r,k,sc,h] - StorageLevel[r,k,sc, (h>1) ? h-1 : length(HOUR)]) / 1 <=  # unit: energy diff per period (TWh/period)
                0.001 * Charging[r,k,h] +
                + (k == :hydro ? 0.001 * hoursperperiod * sum(cfhydroinflow[r,c,h] * Capacity[r,:hydro,c] for c in reservoirclass[sc])
                                : 0.0) +
                + (k == :csp ? 0.001 * hoursperperiod * sum(cf[r,:csp,c,h] * Capacity[r,:csp,c] for c in reservoirclass[sc])
                                : 0.0) +
                - 0.001 * sum(Electricity[r,k,c,h]/efficiency[k] for c in reservoirclass[sc])

        MaxStorageCapacity[r in REGION, k in storagetechs, sc in STORAGECLASS[k], h in HOUR],
            StorageLevel[r,k,sc,h] <= sum(Capacity[r,k,c] * dischargetime[r,k,c] for c in reservoirclass[sc]) / 1000

        InitialStorageLevel[r in REGION, k in storagetechs, sc in STORAGECLASS[k]],
            StorageLevel[r,k,sc,1] ==
                initialstoragelevel * sum(Capacity[r,k,c] * dischargetime[r,k,c] for c in reservoirclass[sc]) / 1000

        # consider turbine efficiencies later
        # POSSIBLE BUG (why only class x0 here?)
        MinHydroFlow[r in REGION, h in HOUR],
            Electricity[r,:hydro,:x0,h] >= minflow_existinghydro * hoursperperiod * cfhydroinflow[r,:x0,h] * Capacity[r,:hydro,:x0]

        # We'll add pumped hydro later
        NoCharging[r in REGION, h in HOUR, k in [:hydro, :csp]],
            Charging[r,k,h] == 0

        ChargingNeedsBattery[r in REGION, h in HOUR],
            Charging[r,:battery,h] <= Capacity[r,:battery, :_] * hoursperperiod

        MaxTransmissionCapacity[r1 in REGION, r2 in REGION, h in HOUR],
            Transmission[r1,r2,h] <= TransmissionCapacity[r1,r2] * hoursperperiod

        TwoWayStreet[r1 in REGION, r2 in REGION],
            TransmissionCapacity[r1,r2] == TransmissionCapacity[r2,r1]

        NoTransmission[r1 in REGION, r2 in REGION; transmissioninvestcost[r1,r2] == 0],
            TransmissionCapacity[r1,r2] == 0

        Calculate_AnnualGeneration[r in REGION, k in TECH],
            AnnualGeneration[r,k] == sum(Electricity[r,k,c,h] for c in CLASS[k], h in HOUR)

        Calculate_FuelUse[r in REGION, f in FUEL; f != :_],
            FuelUse[r,f] == sum(AnnualGeneration[r,k]/efficiency[k] for k in TECH if techfuel[k]==f)

        BioLimit[r in REGION],
            sum(AnnualGeneration[r,k] for k in [:bioGT, :bioCCGT]) <= maxbioenergy * sum(demand[r,h] for h in HOUR) * hoursperperiod

        # This does not quite make the variable bound redundant, because e.g. some pixels in PV class 1 are class "0" for CSP,
        # and are therefore unaffected by this constraint.
        SolarOverlap[r in REGION, pv in CLASS[:pv], csp in CLASS[:csp]],
            SolarCapacity[r,:pv,pv,csp]/pv_density + SolarCapacity[r,:csp,pv,csp]/csp_density/cspsolarmultiple <=
                solarcombinedarea[r,pv,csp]

        # There could be a minor bug here depending on the numerical class limits. If there are pixels in e.g. PV class a1
        # which do not belong to any CSP class, then we may miscount PV capacity. But this does not happen with our parameters.
        PVcapacity[r in REGION, pv in CLASS[:pv]],
            Capacity[r,:pv,pv] == sum(SolarCapacity[r,:pv,pv,csp] for csp in CLASS[:csp])

        CSPcapacity[r in REGION, csp in CLASS[:csp]],
            Capacity[r,:csp,csp] == sum(SolarCapacity[r,:csp,pv,csp] for pv in CLASS[:pv])

        TotalCO2[r in REGION],
            CO2emissions[r] == sum(FuelUse[r,f] * emissionsCO2[f] for f in FUEL)

        GlobalCO2Cap,
            sum(CO2emissions[r] for r in REGION) <= carboncap * sum(demand) * hoursperperiod

        # Storage costs included in ordinary Capacity costs
        # Transmission costs halved since they are counted twice (allocate half the cost to sending and receiving regions)
        # make this regional later
        Totalcosts[r in REGION],
            Systemcost[r] ==
                0.001 * sum(FuelUse[r,f] * fuelcost[f] for f in FUEL) +
                + 0.001 * carbontax * CO2emissions[r] +
                + 0.001 * sum(Electricity[r,k,c,h] * variablecost[k] for k in TECH, c in CLASS[k], h in HOUR) +
                + 0.001 * sum(Electricity[r,:hydro,c,h] * hydroeleccost[r,c] for c in CLASS[:hydro], h in HOUR) +
                + 0.001 * sum(Transmission[r,r2,h] * smalltransmissionpenalty for r2 in REGION, h in HOUR) +
                + sum(Capacity[r,k,c] * (investcost[k,c] * crf[k] + fixedcost[k]) for k in TECH, c in CLASS[k]) +
                + 0.5 * sum(TransmissionCapacity[r,r2] *
                            (transmissioninvestcost[r,r2] * crf[:transmission] + transmissionfixedcost[r,r2]) for r2 in REGION)
        # =#
    end #constraints

    if rampingconstraints
        @constraints m begin
            RampingDown[r in REGION, k in TECH, c in CLASS[k], h in HOUR; rampingrate[k] < 1],
                Electricity[r,k,c,h] - Electricity[r,k,c, (h>1) ? h-1 : length(HOUR)] >=
                    -rampingrate[k] * Capacity[r,k,c] * cf[r,k,c,h] * hoursperperiod

            RampingUp[r in REGION, k in TECH, c in CLASS[k], h in HOUR; rampingrate[k] < 1],
                Electricity[r,k,c,h] - Electricity[r,k,c, (h>1) ? h-1 : length(HOUR)] <=
                    rampingrate[k] * Capacity[r,k,c] * cf[r,k,c,h] * hoursperperiod
        end
    else
        RampingDown = RampingUp = nothing
    end

    if globalnuclearlimit != Inf
        @constraints m begin
            NuclearCapacityLimit,
                sum(Capacity[r,:nuclear,:_] for r in REGION) <= globalnuclearlimit
        end
    end

    return Constraints(ElecCapacity, ElecDemand, RampingDown, RampingUp, StorageBalance, MaxStorageCapacity, InitialStorageLevel,
                MaxTransmissionCapacity, TwoWayStreet, NoTransmission, NoCharging, ChargingNeedsBattery,
                Calculate_AnnualGeneration, Calculate_FuelUse, TotalCO2, Totalcosts)
end

function makeobjective(m, sets, vars)
    @unpack REGION = sets
    @unpack Systemcost = vars

    @objective m Min begin
        sum(Systemcost[r] for r in REGION)
    end
end

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
using NPZ, JLD

function main()
	in_file = ARGS[end]
	if isfile(in_file)
		if endswith(in_file, ".jld")
			model = load(in_file)

			w1 = model["w1"]
			w2 = model["w2"]

			npzwrite(replace(in_file, ".jld" => "_w1") * ".npy", w1)
			npzwrite(replace(in_file, ".jld" => "_w2") * ".npy", w2)
		else
			println("$in_file needs to be a JLD file, have suffix .jld")
		end
	else
		println("$in_file is not a file")
	end
end

main()


text:
#73
na = NamedArray([1, 2], ([1, missing],), ("A",))
let one = na[Name(1)], two = na[Name(missing)]
    @test one == 1
    @test two == 2
end

#39
include("init-namedarrays.jl")
v = n[1, :]
@test sin.(v).array == sin.(v.array)
@test namesanddim(sin.(v)) == namesanddim(v)


text:
using FixedPointNumbers
using LRUCache
using FileIO
using Lazy: @>, @>>
using IntervalSets

import FileIO: save
import DSP: resample
import LibSndFile
import SampledSignals: samplerate
import SampledSignals
import Distributions: nsamples
import Base: show, length, start, done, next, getindex, size,
  setindex!, vcat, similar, convert, .*, .+, *, minimum, maximum
importall IntervalSets # I just need .., but there's a syntax parsing bug

export sound, playable, duration, nchannels, nsamples, save, samplerate, length,
  samples, vcat, leftright, similar, left, right, resample,
  audiofn, limit, .., ends, data

immutable Sound{R,T,N} <: AbstractArray{T,N}
  data::Array{T,N}
  function Sound{R,T,N}(a::Array{T,N}) where {R,T,N}
    if T <: Integer
      error("Cannot use integer arrays for a sound. ",
            "Use FixedPointNumbers instead.")
    end

    if N == 1
      new(a)
    elseif N == 2
      new(a)
    else
      error("Unexpected dimension count $N for sound array, should be 1",
            " (for mono) or 2 (for mono or stereo).")
    end
  end
end
convert(::Type{Sound{R,T,N}},x) where {R,T,N} =
  Sound{R,T,N}(convert(Array{T,N},x))
function convert(::Type{Sound{R,T,N}},x::Sound{R,S,N}) where {R,T,S,N}
  Sound{R,T,N}(convert(Array{T,N},x.data))
end
function convert{R,Q,T,S}(::Type{Sound{R,T}},x::Sound{Q,S})
  error("Cannot convert a sound with sampling rate $(Q*Hz) to a sound with ",
        "sampling rate $(R*Hz). Use `resample` to change the sampling rate.")
end

function Base.Array(x::Sound{R,T,N}) where {R,T,N}
  x.data
end

"""
    samplerate([sound])

Report the sampling rate of the sound or of any object
that can be turned into a sound.

The sampling rate of an object determines how many samples per second are used
to represent the sound. Objects that can be converted to sounds are assumed to
be at the sampling rate of the current hardware settings as defined by
[`setup_sound`](@ref).
"""
samplerate(x::Vector) = samplerate()
samplerate(x::Matrix) = samplerate()
samplerate{R}(x::Sound{R}) = R*Hz
rtype{R}(x::Sound{R}) = R

length(x::Sound) = length(x.data)

asstereo{R,T}(x::Sound{R,T,1}) = hcat(x.data,x.data)
asstereo{R,T}(x::Sound{R,T,2}) = size(x,2) == 1 ? hcat(x.data,x.data) : x.data
asmono{R,T}(x::Sound{R,T,1}) = x.data
asmono{R,T}(x::Sound{R,T,2}) = squeeze(x.data,2)

vcat(xs::Sound...) = error("Sample rates differ, fix with `resample`.")
function vcat(xs::Sound{R}...) where R
  T = promote_type(map(eltype,xs)...)
  vcat((map(el -> convert(T,el),x) for x in xs)...)
end
vcat(xs::Sound{R,T,1}...) where {R,T} =
    Sound{R,T,1}(vcat(map(x -> x.data,xs)...))
function vcat(xs::Sound{R,T}...) where {R,T}
  if any(x -> nchannels(x) == 2,xs)
    Sound{R,T,2}(vcat(map(asstereo,xs)...))
  else
    Sound{R,T,1}(vcat(map(asmono,xs)...))
  end
end

function *(x::Number,y::Sound)
  z = similar(y)
  z .= x .* y
end

function *(y::Sound,x::Number)
  z = similar(y)
  z .= x .* y
end

"""
    duration(x)

Get the duration of the given sound in seconds.
"""
duration(x::Sound{R}) where R = uconvert(s,nsamples(x) / (R*Hz))

"""
    nchannels(sound)

Return the number of channels (1 for mono, 2 for stereo) in this sound.
"""
nchannels(x::Sound) = size(x.data,2)

"""
    nsamples(sound::Sound)

Returns the number of samples in the sound.
"""
nsamples(x::Sound) = size(x.data,1)
size(x::Sound) = size(x.data)
Base.IndexStyle(::Type{Sound}) = IndexLinear()

# adapted from:
# https://github.com/JuliaAudio/SampledSignals.jl/blob/0a31806c3f7d382c9aa6db901a83e1edbfac62df/src/SampleBuf.jl#L109-L139
rounded_time(x,R) = round(ustrip(inseconds(x,R)),ceil(Int,log(10,R)))*s
function show(io::IO, x::Sound{R}) where R
  seconds = rounded_time(duration(x),R)
  typ = if eltype(x) == Q0f15
    "16 bit PCM"
  elseif eltype(x) <: AbstractFloat
    "$(sizeof(eltype(x))*8) bit floating-point"
  else
    eltype(x)
  end

  channel = size(x.data,2) == 1 ? "mono" : "stereo"

  println(io, "$seconds $typ $channel sound")
  print(io, "Sampled at $(R*Hz)")
  nsamples(x) > 0 && showchannels(io, x)
end
show(io::IO, ::MIME"text/plain", x::Sound) = show(io,x)

const ticks = ['_','▁','▂','▃','▄','▅','▆','▇']
function showchannels(io::IO, x::Sound, widthchars=80)
  # number of samples per block
  blockwidth = round(Int, nsamples(x)/widthchars, RoundUp)
  nblocks = round(Int, nsamples(x)/blockwidth, RoundUp)
  blocks = Array{Char}(nblocks, nchannels(x))
  for blk in 1:nblocks
    i = (blk-1)*blockwidth + 1
    n = min(blockwidth, nsamples(x)-i+1)
    peaks = sqrt.(mean(float(x[(1:n)+i-1,:]).^2,1))
    # clamp to -60dB, 0dB
    peaks = clamp.(20log10.(peaks), -60.0, 0.0)
    idxs = trunc.(Int, (peaks+60)/60 * (length(ticks)-1)) + 1
    blocks[blk, :] = ticks[idxs]
  end
  for ch in 1:nchannels(x)
    println(io)
    print(io, convert(String, blocks[:, ch]))
  end
end

@inline function getindex(x::Sound,i::Int)
  @boundscheck checkbounds(x.data,i)
  @inbounds return x.data[i]
end

@inline function setindex!{R,T,S}(x::Sound{R,T},v::S,i::Int)
  @boundscheck checkbounds(x.data,i)
  @inbounds return x.data[i] = convert(T,v)
end

@inline function getindex(x::Sound,i::Int,j::Int)
  @boundscheck checkbounds(x.data,i,j)
  @inbounds return x.data[i,j]
end

@inline function setindex!{R,T,S}(x::Sound{R,T},v::S,i::Int,j::Int)
  @boundscheck checkbounds(x.data,i,j)
  @inbounds return x.data[i,j] = convert(T,v)
end

"""
    left(sound::Sound)

Extract the left channel a stereo sound or stream.
"""
function left(sound::Sound{R,T,N}) where {R,T,N}
  if size(sound.data,2) == 1
    error("Expected stereo sound.")
  else
    sound[:,1]
  end
end

"""
    right(sound::Sound)

Extract the right channel of a stereo sound or stream.
"""
function right(sound::Sound{R,T,N}) where {R,T,N}
  if size(sound.data,2) == 1
    error("Expected stereo sound.")
  else
    sound[:,2]
  end
end

immutable EndSecs
end
const ends = EndSecs()

immutable ClosedIntervalEnd{N}
  from::Quantity{N}
end
minimum(x::ClosedIntervalEnd) = x.from

..(x::Time,::EndSecs) = ClosedIntervalEnd(x)
..(x::SampleQuant,::EndSecs) = ClosedIntervalEnd(x)
..(x,::EndSecs) = error("Unexpected quantity $x in interval.")

function checktime(time)
  if ustrip(time) < 0
    throw(BoundsError("Unexpected negative time."))
  end
end

@inline @Base.propagate_inbounds function getindex(x::Sound,ixs,js::Symbol)
  if js == :left
    getindex(x,ixs,1)
  elseif js == :right
    getindex(x,ixs,2)
  else
    throw(BoundsError(x,js))
  end
end
@inline @Base.propagate_inbounds function setindex!(x::Sound,vals,ixs,js::Symbol)
  if js == :left
    setindex!(x,vals,ixs,1)
  elseif js == :right
    setindex!(x,vals,ixs,2)
  else
    throw(BoundsError(x,js))
  end
end

########################################
# getindex
const Index = Union{Integer,Range,AbstractVector,Colon}
@inline function getindex(x::Sound{R,T},
                          ixs::ClosedIntervalEnd,js::I) where {R,T,I <: Index}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  @boundscheck checkbounds(x.data,from,js)
  @inbounds return Sound{R,T,2}(x.data[from:end,js])
end

@inline function
    getindex(x::Sound{R,T,N},
             ixs::ClosedInterval{TM},js::I) where {R,T,I <: Index,N,TM <: Quantity}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  to = insamples(maximum(ixs),R*Hz)-1
  @boundscheck checkbounds(x.data,from:to,js)
  @inbounds result = x.data[from:to,js]
  return Sound{R,T,ndims(result)}(result)
end

@inline function getindex(x::Sound{R,T},ixs::ClosedIntervalEnd) where {R,T}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  if size(x,2) == 1
    @boundscheck checkbounds(x.data,from)
    @inbounds return Sound{R,T,1}(x.data[from:end])
  else
    @boundscheck checkbounds(x.data,from,:)
    @inbounds return Sound{R,T,2}(x.data[from:end,:])
  end
end

@inline function getindex(x::Sound{R,T},
                          ixs::ClosedInterval{TM}) where {R,T,TM <: Quantity}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  to = insamples(maximum(ixs),R*Hz)-1
  if size(x,2) == 1
    @boundscheck checkbounds(x.data,from)
    @boundscheck checkbounds(x.data,to)
    @inbounds return Sound{R,T,1}(x.data[from:to])
  else
    @boundscheck checkbounds(x.data,from,:)
    @boundscheck checkbounds(x.data,to,:)
    @inbounds return Sound{R,T,2}(x.data[from:to,:])
  end
end

########################################
# setindex

@inline function setindex!{R,T,I}(
  x::Sound{R,T},vals::AbstractArray,ixs::ClosedIntervalEnd,js::I)
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  @boundscheck checkbounds(x.data,from,js)
  @inbounds x.data[from:end,js] = vals
  vals
end

@inline function setindex!(
  x::Sound{R,T},vals::AbstractArray,
  ixs::ClosedInterval{TM},js::I) where {R,T,I,TM <: Quantity}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  to = insamples(maximum(ixs),R*Hz)-1
  @boundscheck checkbounds(x.data,from,js)
  @boundscheck checkbounds(x.data,to,js)
  @inbounds x.data[from:to,js] = vals
  vals
end

@inline function setindex!{R,T}(
  x::Sound{R,T},vals::AbstractArray,ixs::ClosedIntervalEnd)
    @boundscheck checktime(minimum(ixs))
    from = max(1,insamples(minimum(ixs),R*Hz))
    if size(x,2) == 1
      @boundscheck checkbounds(x.data,from)
      @inbounds x.data[from:end] = vals
    else
      @boundscheck checkbounds(x.data,from,:)
      @inbounds x.data[from:end,:] = vals
    end
    vals

end

@inline function setindex!(
    x::Sound{R,T},vals::AbstractArray,
    ixs::ClosedInterval{TM}) where {R,T,TM <: Quantity}
  @boundscheck checktime(minimum(ixs))
  from = max(1,insamples(minimum(ixs),R*Hz))
  to = insamples(maximum(ixs),R*Hz)-1
  if size(x,2) == 1
    @boundscheck checkbounds(x.data,from:to)
    @inbounds x.data[from:to] = vals
  else
    @boundscheck checkbounds(x.data,from:to,:)
    @inbounds x.data[from:to,:] = vals
  end
  vals
end

function similar(x::Sound{R,T,N},::Type{S},dims::NTuple{M,Int}) where {R,T,S,N,M}
  if M ∉ [1,2] || (M == 2 && dims[2] ∉ [1,2])
    similar(x.data,S,dims)
  else
    Sound{R,S,M}(similar(x.data,S,dims))
  end
end

save(file::Union{AbstractString,IO},sound::Sound) = save(file,assampled(sound))
assampled(x::Sound{R}) where R = SampledSignals.SampleBuf(x.data,float(R))

"""
    resample(x::Sound,samplerate)

Returns a new sound representing the sound `x` at the given sampling rate.

You will loose all frequencies in the sound above samplerate/2. Resampling
occurs automatically when you call [`sound`](@ref)--which is called inside
[`play`](@ref)--anytime the sampling rate of the sound and the current audio
playback settings (determined by [`setup_sound`](@ref)) are not the same.

To avoid automatic resampling you can either create sounds at the appropriate
sampling rate, as determined by [`samplerate`](@ref) (recommended), or change
the sampling rate initialized during [`setup_sound`](@ref) (not recommended).
"""
function resample(x::Sound{R,T,N},new_sample_rate) where {R,T,N}
  new_rate = floor(Int,ustrip(inHz(new_sample_rate)))
  if new_rate < R
    warn("The function `resample` reduced the sample rate, high freqeuncy",
         " information above $(new_rate/2) Hz will be lost.")
  end
  if size(x,2) == 2
    Sound{new_rate,T,N}(hcat(T.(resample(x.data[:,1],new_rate // R)),
                             T.(resample(x.data[:,2],new_rate // R))))
  else
    Sound{new_rate,T,1}(T.(resample(x.data[:],new_rate // R)))
  end
end

function duration(x::Array{Float64};sample_rate_Hz=samplerate())
  uconvert(s,nsamples(x) / inHz(sample_rate_Hz))
end

const sound_cache = LRU{Tuple{UInt,Int,Function},Sound}(256)
clear_sound_cache() = empty!(sound_cache)
function with_cache(fn,usecache,x,sr)
  if usecache
    get!(fn,sound_cache,(object_id(x),ustrip(inHz(Int,sr)),fn))
  else
    fn()
  end
end

"""
    sound(x::Array,[cache=TimedSound.usecache()];[sample_rate=samplerate()])

Creates a sound object from an arbitrary array.

Assumes 1 is the loudest and -1 the softest. The array should be 1d for mono
signals, or an array of size (N,2) for stereo sounds.

When cache is set to true, sound will cache its results thus avoiding repeatedly
creating a new sound for the same object.

!!! note "Called Implicitly"

    This function is normally called implicitly in a call to `play(x)`, so it
    need not normally be called directly.

"""
function sound(x::Array{T,N},cache=usecache();
               sample_rate=samplerate()) where {T <: Number,N}
  if N ∉ [1,2]
    error("Array must have 1 or 2 dimensions to be converted to a sound.")
  end

  with_cache(cache,x,sample_rate) do
    R = ustrip(inHz(sample_rate))
    Sound{R,T,N}(x)
  end
end

function sound(x::SampledSignals.SampleBuf,cache=usecache();
               sample_rate=samplerate(x)*Hz)
  if ustrip(sample_rate) != samplerate(x)
    error("Unexpected sample rate $sample_rate. Use `playable` or `resample`",
          " to change the sampling rate.")
  end

  with_cache(cache,x,sample_rate) do
    R,T = ustrip(inHz(Int,samplerate(x)*Hz)),eltype(x.data)
    Sound{R,T,ndims(x)}(x.data)
  end
end

function sound(x::Sound{R},cache=usecache();sample_rate=R*Hz) where R
  if ustrip(sample_rate) != R
    error("Unexpected sample rate $sample_rate. Use `playable` or `resample`",
          " to change the sampling rate.")
  end
  x
end

"""
    playable(x,[cache=usecache()],[sample_rate=samplerate()])

Prepare a sound or stream to be played.

A call to `playable` will ensure the sound is in the format required by
[`play`](@ref).  This automatically calls [`sound`](@ref) on `x` if it does not
appear to already be a sound or a stream.

!!! note "Called Implicitly"

    This need not be called explicitly, as play will call it for you, if need be.
"""

function playable(x,cache=usecache(),sample_rate=samplerate())
  with_cache(cache,x,sample_rate) do
    playable(sound(x,sample_rate=sample_rate),false,sample_rate)
  end
end

function playable{R,T,N}(x::Sound{R,T,N},cache=usecache(),sample_rate=samplerate())
  with_cache(cache,x,sample_rate) do
    bounded = clamp.(x.data,typemin(Q0f15),typemax(Q0f15))
    T2 = Q0f15
    playable(Sound{R,T2,N}(Q0f15.(bounded)),false,sample_rate)
  end
end

function playable{R}(x::Sound{R,Q0f15},cache=usecache(),sample_rate=samplerate())
  if size(x.data,2) == 2
    if R == ustrip(sample_rate)
      x
    else
      with_cache(cache,x,sample_rate) do
        warn("Reampling sound")
        resample(Sound{R,Q0f15,2}(x.data),sample_rate)
      end
    end
  else
    with_cache(cache,x,sample_rate) do
      data = hcat(x.data,x.data)
      if R == ustrip(sample_rate)
        Sound{R,Q0f15,2}(data)
      else
        warn("Reampling sound")
        resample(Sound{R,Q0f15,2}(data),sample_rate)
      end
    end
  end
end

"""
    sound(file,[cache=usecache()];[sample_rate=samplerate(file)])

Load a specified file as a sound.
"""
sound(file::File,cache=usecache();keys...) = sound(load(file),cache;keys...)
sound(file::String,cache=usecache();keys...) = sound(load(file),cache;keys...)
function sound(stream::IOStream,cache=false;keys...)
  if cache
    error("Cannot cache a sound from an IOStream. Please use the filename.")
  end
  sound(load(stream),cache;keys...)
end

"""
    leftright(left,right)

Create a stereo sound from two vectors or two monaural sounds.

For vectors, one can specify a sample_rate other than the default,
if desired.
"""
function leftright(x::Sound{R,T},y::Sound{R,T},sample_rate=R*Hz) where {R,T}
  if sample_rate != R*Hz
    error("Unexpected sampling rate $sample_rate. ",
          "Use `playable` or `resample` to change the sampling rate.")
  end
  if size(x.data,2) == size(y.data,2) == 1
    mix(hcat(x,silence(nsamples(x)*samples)),
        hcat(silence(nsamples(y)*samples),y))
  else
    error("Expected two monaural sounds.")
  end
end

limit(sound::Sound{R},len::Time) where R = limit(sound,insamples(len,R*Hz))
limit(sound::Sound{R,T,1},len::Int) where {R,T} = sound[1:len]
limit(sound::Sound{R,T,2},len::Int) where {R,T} = sound[1:len,:]

"""
    audiofn(fn,x)

Apply `fn` to x for both sounds and streams.

For a sound this is the same as calling `fn(x)`.
"""
audiofn(fn::Function,x::Sound{R}) where R = sound(fn(x),sample_rate=R*Hz)

function soundop(op,xs::Union{Sound{R},Array}...) where R
  len = maximum(map(x -> size(x,1),xs))
  channels = maximum(map(x -> size(x,2),xs))
  y = similar(xs[1],(len,channels))

  for i in 1:size(y,1)
    used = false
    for j in 1:length(xs)
      if i <= size(xs[j],1)
        if !used
          used = true
          @inbounds y[i,:] .= xs[j][i,:]
        else
          @inbounds y[i,:] .= op(y[i,:],xs[j][i,:])
        end
      end
    end
  end

  y
end


text:
export zoom_generic

function zoom_generic(h :: AbstractNLPModel,
                  	  t₀ :: Float64,
	                  t₁ :: Float64;
					  direction :: String = "Nwt",
	                  tol :: Float64 = 1e-7,
	                  maxiter :: Int = 100,
	                  verbose :: Bool = false,
					  kwargs...)
    if obj(h, [t₀])[1] < obj(h, [t₁])[1]
    tl = t₀; th = t₁
    else
    tl = t₁; th = t₀
    end

    hl = obj(h, [tl])[1]; dhl = grad(h, [tl])[1]
    hh = obj(h, [th])[1]; dhh = grad(h, [th])[1]
    h₀ = obj(h, [t₀])[1]; dh₀ = grad(h, [t₀])[1]

    γ = 0.8; t = t₁; tp = t₀; tqnp = t₀; hk = 0; hkm1 = 0; gkm1 = 0; hp = 0
    i = 0

    hk = obj(h, [t])[1]; gk = grad(h, [t])[1]
    hkm1 = obj(h, [tqnp])[1]; gkm1 = grad(h, [tqnp])[1]

    verbose &&
        @printf(" iter        tₗ        tₕ         t        hₗ        hₕ")
    verbose &&
        @printf("         hk         gk\n")
    verbose &&
        @printf(" %7.2e %7.2e  %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e\n",
                i, tl, th, t, hl, hh, hk, gk)

    while ((i < maxiter) && (abs(gk) > tol)) || i == 0

    if (hk > h₀) || (hl <= hk)
      th = t; tlast = th; hlast = hh
      dhlast = dhh; hh = hk; dhh = gk
    else
      if (abs(gk) < tol)
        topt = t; iter = i
        return (topt, iter)
      elseif gk * (th - tl) >= 0
        th = tl; hh = hl; dhh = dhl
      end
      tl = t; tlast = tl; hlast = hl
      dhlast = dhl; hl = hk; dhl = gk
    end
    if direction == "Nwt"
        kₖ = hess(h, [t])[1]
        dN = -gk/kₖ #direction de Newton
    elseif direction == "Sec"
        s = t - tqnp; y = gk - gkm1
        dN =-gk * s / y
    elseif direction == "SecA"
        s = t - tqnp; y = gk - gkm1
        Γ = 3 * (gk + gkm1) * s - 6 * (hk - hkm1)
        if y * s + Γ < eps(Float64) * (s^2)
          yt = y
        else
          yt = y + Γ / s
        end
        dN = -gk * s / yt
    elseif direction == "Cub"
        s = t - tqnp; y = gk - gkm1; α = -s
        z = gk + gkm1 + 3 * (hk - hkm1) / α
        discr = z^2 - gk * gkm1; denom = gk + gkm1 + 2 * z
        if (discr > 0.0) && (abs(denom) > eps(Float64))
          w = sqrt(discr)
          dN = -s * (gk + z + sign(α) * w) / (denom)
        else
          dN=-gk*s/y
        end
    end

    if ((tp - t) * dN > 0) & (dN/(tp - t) < γ)
      tplus = t + dN
      hplus = obj(h, [tplus])[1]; gplus = grad(h, [tplus])[1]
      verbose && println("N")
    else
      tplus = (t + tp) / 2
      hplus = obj(h, [tplus])[1]; gplus = grad(h, [tplus])[1]
      verbose && println("B")
    end

    if t > tp
      if gplus < 0.0
        tp = t; tqnp = t; t = tplus
      else
        tqnp = t; t = tplus
      end
    else
      if gplus > 0.0
        tp = t; tqnp = t; t = tplus
      else
        tqnp = t; t = tplus
      end
    end

    #mise à jour des valeurs
    hkm1 = hk; gkm1 = gk; hk = hplus; gk = gplus
    i = i + 1
    #println("on continue d'itérer")
    verbose &&
        @printf(" %7.2e %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e %7.2e\n",
                i, tl, th, t, hl, hh, hk, gk)

    end
    topt = t; iter = i

    return (topt, iter)
end #function


text:
using LinearAlgebra
using Optim
using GFX

import StatsFuns: @irrational
@irrational log2π  1.8378770664093454836 log(big(2.)*π)

function KLDivergence(μ0, σ0, μ1, σ1)
    "KL of 2 univariate Normals"
#     return log(σ_q / σ_p) + (σ_p + (μ_p - μ_q)^2)/ (2*σ_q) - 0.5
    return 0.5*( (σ0/σ1)^2 + (μ0 - μ1)^2/σ1^2 +1 + 2*log(σ1/σ0))
end

function minEFE(current_state,
                goal_state,
                model_params;
                Δt=1.0,
                num_iters=100,
                plan_horizon=2)

    # Objective function
    G(π) = EFE(π, current_state, goal_state, model_params, Δt=Δt, time_horizon=plan_horizon)

    # Minimize
    results = optimize(G, rand(plan_horizon,), Optim.Options(iterations=num_iters), LBFGS(); autodiff=:forward)

    return Optim.minimizer(results)
end

function EFE(action, state_k, goal_state, params; Δt = 1., time_horizon=1)
    "Expected Free Energy"

    # Unpack goal state
    μ_star = goal_state[1]
    Σ_star = goal_state[2]

    # Unpack model parameters
    θ, η, τ, ξ = params

    # Process noise
    Σ_z = inv(τ)*[Δt^3/3 Δt^2/2; Δt^2/2 Δt]

    # Utility matrices
    S = [1. Δt; 0. 1.]
    s = [0., Δt]

    # Model matrices
    A = S + s*θ'
    B = s*η
    C = [1., 0.]'

    # Unpack prior state
    μ_tmin1 = mn(state_k)
    Σ_tmin1 = cv(state_k)

    cumEFE = 0
    for t in 1:time_horizon

        # State transition
        μ_t = A*μ_tmin1 + B*action[t]
        Σ_t = A*Σ_tmin1*A' + Σ_z

        # Predicted observation
        y_t = C*μ_t

        # Block covariance matrix of joint
        Σ_11 = Σ_t
        Σ_21 = C*Σ_t
        Σ_12 = Σ_t*C'
        Σ_22 = C*Σ_t*C' .+ inv(ξ)

        # Calculate conditional entropy
        Σ_cond = Σ_22 - Σ_21*inv(Σ_11)*Σ_12
        ambiguity = 0.5(log2π + log(Σ_cond[1]) - 1)

        # Risk as KL between marginal and goal prior
        risk = 0.5*(log(det(Σ_star)/det(Σ_22)) + 1 + (y_t - μ_star)'*inv(Σ_star)*(y_t - μ_star) + tr(inv(Σ_star)*Σ_22))

        # Add to cumulative EFE
        cumEFE += risk + ambiguity

        # Update previous state
        μ_tmin1 = μ_t
        Σ_tmin1 = Σ_t

    end
    return cumEFE
end;

max_stars_repo_name: wmkouw/actinf-oscillator

max_issues_repo_name: wmkouw/actinf-oscillator

max_forks_repo_name: wmkouw/actinf-oscillator


text:
function solve_euler(timestep::Float64, n::Int64)
    euler_result = Vector{Float64}(undef, n)

    # Setting the initial condition
    euler_result[1] = 1;
    for i = 2:length(euler_result)
        euler_result[i] = euler_result[i-1] - 3.0*euler_result[i-1]*timestep
    end
    return euler_result
end

function check_result(euler_result::Vector{Float64}, threshold::Float64,
                      timestep::Float64)
    is_approx = true

    for i = 1:length(euler_result)
        time = (i - 1)*timestep
        solution = exp(-3*time);
        if (abs(euler_result[i] - solution) > threshold)
            println(euler_result[i], solution)
            is_approx = false
        end
    end

    return is_approx
end

function main()
    timestep = 0.01
    n = 100
    threshold = 0.01

    euler_result = solve_euler(timestep,n)
    is_approx = check_result(euler_result, threshold, timestep)

    println(is_approx)
end

main()

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive