text:
# 1D diffusion problem

# TODO: Add more complex tests.

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test,OrdinaryDiffEq

# Tests
@testset "Test 00: Dt(u(t,x)) ~ Dxx(u(t,x))" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ cos(x),
           u(t,0) ~ exp(-t),
           u(t,Float64(pi)) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = prob.space[2]
    t = sol.t

    # Plot and save results
    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    #     scatter!(x, u_exact(x, t[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test00.png")

    # Test against exact solution
    for i in 1:size(t,1)
        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
    end
end

@testset "Test 01: Dt(u(t,x)) ~ D*Dxx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x D
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dxx''~x

    D = 1.1

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ D*Dxx(u(t,x))
    bcs = [u(0,x) ~ -x*(x-1)*sin(x),
           u(t,0) ~ 0.0,
           u(t,1) ~ 0.0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x,D],[u])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    # Plot and save results
    # x = prob.space[2]
    # t = sol.t

    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test01.png")

    # Test
    n = size(sol,1)
    t_f = size(sol,3)
    @test sol[:,1,t_f] ≈ zeros(n) atol = 0.001;
end

@testset "Test 02: Dt(u(t,x)) ~ Dx(D(t,x))*Dx(u(t,x))+D(t,x)*Dxx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..) D(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions

    eq  = [ Dt(u(t,x)) ~ Dx(D(t,x))*Dx(u(t,x))+D(t,x)*Dxx(u(t,x)),
            D(t,x) ~ 0.999 + 0.001 * t * x  ]

    bcs = [u(0,x) ~ -x*(x-1)*sin(x),
           u(t,0) ~ 0.0,
           u(t,1) ~ 0.0,
           D(0,x) ~ 0.999,
           D(t,0) ~ 0.999,
           D(t,1) ~ 0.999 + 0.001 * t ]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u,D])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    # Plot and save results
    # x = prob.space[2]
    # t = sol.t

    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test02.png")

    # Test
    n = size(sol,1)
    t_f = size(sol,3)
    @test sol[:,1,t_f] ≈ zeros(n) atol=0.01;
end

@testset "Test 03: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           Dx(u(t,0)) ~ exp(-t),
           Dx(u(t,Float64(pi))) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = prob.space[2]
    t = sol.t

    # Plot and save results
    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    #     scatter!(x, u_exact(x, t[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test03.png")

    # Test against exact solution
    for i in 1:size(t,1)
        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
    end
end

@testset "Test 04: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann + Dirichlet BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           u(t,0) ~ 0.0,
           Dx(u(t,Float64(pi))) ~ -exp(-t)]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,Float64(pi))]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = prob.space[2]
    t = sol.t

    # Plot and save results
    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    #     scatter!(x, u_exact(x, t[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test04.png")

    # Test against exact solution
    for i in 1:size(t,1)
        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
    end
end

@testset "Test 05: Dt(u(t,x)) ~ Dxx(u(t,x)), Robin BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (sin(-1.0) + 3cos(-1.0)),
           4u(t,1.0) + Dx(u(t,1.0)) ~ exp(-t) * (4sin(1.0) + cos(1.0))]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(-1.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.01
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = prob.space[2]
    t = sol.t

    # Plot and save results
    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    #     scatter!(x, u_exact(x, t[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test05.png")

    # Test against exact solution
    for i in 1:size(t,1)
        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
        @test u_approx ≈ u_exact(x, t[i]) atol=0.05
    end
end

@testset "Test 06: Dt(u(t,x)) ~ Dxx(u(t,x)), time-dependent Robin BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    bcs = [u(0,x) ~ sin(x),
           t^2 * u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (t^2 * sin(-1.0) + 3cos(-1.0)),
           4u(t,1.0) + t * Dx(u(t,1.0)) ~ exp(-t) * (4sin(1.0) + t * cos(1.0))]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(-1.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.01
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)
    x = prob.space[2]
    t = sol.t

    # Plot and save results
    # using Plots
    # plot()
    # for i in 1:4
    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
    #     scatter!(x, u_exact(x, t[i]))
    # end
    # savefig("MOL_1D_Linear_Diffusion_Test06.png")

    # Test against exact solution
    for i in 1:size(t,1)
        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
        @test u_approx ≈ u_exact(x, t[i]) atol=0.05
    end
end

@testset "Test errors" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..) v(..)
    @derivatives Dt'~t
    @derivatives Dx'~x
    @derivatives Dxx''~x

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0), x ∈ IntervalDomain(0.0,1.0)]

    # Method of lines discretization
    dx = 0.1
    order = 2
    discretization = MOLFiniteDifference(dx,order)

    # Missing boundary condition
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Boundary condition not at t=0
    bcs = [u(1,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           Dx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Boundary condition not at an edge of the domain
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           Dx(u(t,0.5)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Second-order derivative in BC
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dxx(u(t,0)) ~ 0.0,
           Dx(u(t,0.5)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Wrong format for Robin BCs
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           u(t,1) * Dx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           u(t,1) + Dxx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           u(t,1) + 2Dxx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           Dx(u(t,1)) + u(t,1) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           u(t,1) / 2 + Dx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
           Dx(u(t,0)) ~ 0.0,
           u(t,1) + Dx(u(t,1)) / 2 ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Mismatching arguments
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
    Dx(u(t,0)) ~ 0.0,
    u(t,0) + Dx(u(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

    # Mismatching variables
    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
    Dx(u(t,0)) ~ 0.0,
    u(t,1) + Dx(v(t,1)) ~ 0.0]
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
    @test_throws BoundaryConditionError discretize(pdesys,discretization)

end


text:
using Base.Test
using QuantumOptics

b = SpinBasis(1//2)

psi1 = spinup(b)
psi2 = spindown(b)

rho = tensor(psi1, dagger(psi1))
sigma = tensor(psi2, dagger(psi2))

@test tracedistance(rho, sigma) == 1.
@test tracedistance(rho, rho) == 0.
@test tracedistance(sigma, sigma) == 0.

ext: jl


text:
function maxThroughputPrediction(this::PredictionCenter)

	meanAbsError = {};
	meanRelError = {};
	errorHeader = {};
	extra = {};
	title = ""
	legends={}
	Xdata={}
	Ydata={}
	Xlabel=""
	Ylabel=""

	if this.testMode == 1 #PredictionCenter.TEST_MODE_DATASET
		range = (1:15000)'';
		maxFlushRate = 300; # DY: This is hard-coded... should we do something about this?

		#%cfFlushRateApprox_conf = struct('io_conf', this.testConfig.io_conf, 'workloadName', this.workloadName);
		cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf,"TPCC");
		myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testConfig.transactionMixture, maxFlushRate, cfFlushRateApprox_conf);
		modelP = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.transactionCount);

		testMaxThroughputIdx, testMaxThroughput = findMaxThroughput(this.testConfig.TPSUngrouped);
		trainMaxThroughputIdx, trainMaxThroughput = findMaxThroughput(this.trainConfig.TPSUngrouped);

		actualThr = testMaxThroughput;
		if isempty(actualThr)
			actualThr = Max(vec(this.testConfig.TPSUngrouped));
		end

		idx=[]
		if isdefined(:trainMaxThroughputIdx) && !isempty(trainMaxThroughputIdx)
			idx=(1:trainMaxThroughputIdx)';
		else
			idx=(1:size(this.trainConfig.averageCpuUsage,1))';
		end

		#%CPU-based throughput with classification
		cpuC = barzanLinInvoke(modelP, range*this.testConfig.transactionMixture);

		tempLen=length(cpuC)
		cpuCLThroughput=[]
		while tempLen>=1
			if cpuC[tempLen]>88 && cpuC[tempLen]<90
				cpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(cpuC)
		cpuCUThroughput=[]
		while tempLen>=1
			if cpuC[tempLen]>98 && cpuC[tempLen]<100
				cpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		#%CPU-based without classification
		cpuTModel = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.TPS);
		cpuT = barzanLinInvoke(cpuTModel, range);
		tempLen=length(cpuT)
		cpuTLThroughput=[]
		while tempLen>=1
			if cpuT[tempLen]>88 && cpuT[tempLen]<90
				cpuTLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(cpuT)
		cpuTUThroughput=[]
		while tempLen>=1
			if cpuT[tempLen]>98 && cpuT[tempLen]<100
				cpuTUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		myModelP = barzanLinSolve(this.trainConfig.averageCpuUsage[vec(idx),:], this.trainConfig.transactionCount[vec(idx),:]);
		myCpuC = barzanLinInvoke(myModelP, range*this.testConfig.transactionMixture);
		tempLen=length(myCpuC)
		myCpuCLThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>88 && myCpuC[tempLen]<90
				myCpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCUThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>98 && myCpuC[tempLen]<100
				myCpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCLThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>44 && myCpuC[tempLen]<45
				myCpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCUThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>59 && myCpuC[tempLen]<50
				myCpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		#%Our IO-based throughput
		cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf,"TPCC");
		myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testConfig.transactionMixture, maxFlushRate, cfFlushRateApprox_conf);

		#%Lock-based throughput
		#%getConcurrencyLebel_conf = struct('lock_conf', this.testConfig.lock_conf, 'workloadName', this.workloadName);
		getConcurrencyLebel_conf = typeOfConcurrencyLevel( this.lockConf, "TPCC");
		concurrencyThroughput = findClosestValue("getConcurrencyLevel", (1:10000)''*this.testConfig.transactionMixture, 160, getConcurrencyLebel_conf);

		#%Linear IO-based throughput
		modelFlushRate = barzanLinSolve(this.trainConfig.pagesFlushed, this.trainConfig.transactionCount);
		linFlushRate = barzanLinInvoke(modelFlushRate, range*this.testConfig.transactionMixture);
		linFlushRateThroughput = 0
		for i=length(linFlushRate):-1:1
			if linFlushRate[i]<maxFlushRate
				linFlushRateThroughput = i
				break
			end
		end

		if isempty(linFlushRateThroughput)
			linFlushRateThroughput=0;
		end

		myMaxThroughput1 = min(myCpuCLThroughput, myFlushRateThroughput, concurrencyThroughput);
		if myMaxThroughput1 == myCpuCLThroughput
			PredReasonIdx1 = 1
		elseif myMaxThroughput1 == myFlushRateThroughput
			PredReasonIdx1 = 2
		elseif myMaxThroughput1 == concurrencyThroughput
			PredReasonIdx1 = 3
		end
		myMaxThroughput2 = min(myCpuCUThroughput, myFlushRateThroughput, concurrencyThroughput);
		if myMaxThroughput2 == myCpuCUThroughput
			PredReasonIdx2 = 1
		elseif myMaxThroughput2 == myFlushRateThroughput
			PredReasonIdx2 = 2
		elseif myMaxThroughput2 == concurrencyThroughput
			PredReasonIdx2 = 3
		end

		Xdata = {[1:size(this.testConfig.TPS, 1)]''};
		Ydata = {this.testConfig.TPS};
		legends = {"Original Signal"};

		num_row = size(Xdata[1], 1);
		num_col = size(Xdata[1], 2);

		if !isempty(actualThr)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'');
			push!(Ydata,repmat([actualThr], num_row, num_col));
			push!(legends,"Actual Max Throughput");
		end
		if !isempty(cpuCLThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([cpuCLThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on adjusted LR for CPU+classification");
			push!(errorHeader,"Max Throughput based on adjusted LR for CPU+classification");
			push!(meanAbsError,mae(cpuCLThroughput, actualThr));
			push!(meanRelError,mre(cpuCLThroughput, actualThr));
		end
		if !isempty(cpuCUThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([cpuCUThroughput], num_row, num_col));
			push!(legends,"Max Throughput on LR for CPU+classification");
			push!(errorHeader,"Max Throughput on LR for CPU+classification");
			push!(meanAbsError,mae(cpuCUThroughput, actualThr));
			push!(meanRelError,mre(cpuCUThroughput, actualThr));
		end
		if !isempty(cpuTLThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([cpuTLThroughput], num_row, num_col));
			push!(legends,"Max Throughput on adjusted LR for CPU");
			push!(errorHeader,"Max Throughput on adjusted LR for CPU");
			push!(meanAbsError,mae(cpuTLThroughput, actualThr));
			push!(meanRelError,mre(cpuTLThroughput, actualThr));
		end
		if !isempty(cpuTUThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([cpuTUThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on LR for CPU");
			push!(errorHeader,"Max Throughput based on LR for CPU");
			push!(meanAbsError,mae(cpuTUThroughput, actualThr));
			push!(meanRelError,mre(cpuTUThroughput, actualThr));
		end
		if !isempty(myFlushRateThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([myFlushRateThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on our flush rate model");
			push!(errorHeader,"Max Throughput based on our flush rate model");
			push!(meanAbsError,mae(myFlushRateThroughput, actualThr));
			push!(meanRelError,mre(myFlushRateThroughput, actualThr));
		end
		if !isempty(linFlushRateThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([linFlushRateThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on LR for flush rate");
			push!(errorHeader,"Max Throughput based on LR for flush rate");
			push!(meanAbsError,mae(linFlushRateThroughput, actualThr));
			push!(meanRelError,mre(linFlushRateThroughput, actualThr));
		end
		if !isempty(concurrencyThroughput)
			push!(Xdata,[1:size(this.testConfig.TPS, 1)]'')
			push!(Ydata,repmat([concurrencyThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on our contention model");
			push!(errorHeader,"Max Throughput based on our contention model");
			push!(meanAbsError,mae(concurrencyThroughput, actualThr));
			push!(meanRelError,mre(concurrencyThroughput, actualThr));
		end

                title = "Max Throughput Prediction";
                Ylabel = "TPS";
                Xlabel = "Time";
	elseif this.testMode == 0 #PredictionCenter.TEST_MODE_MIXTURE_TPS
                #% testTPS = [this.testMinTPS:(this.testMaxTPS-this.testMinTPS)/(this.NUM_TPS_VALUES-1):this.testMaxTPS];
                #% testTransactionCount(1,:) = this.testMixture * this.testMinTPS;
                #% testTransactionCount(2,:) = this.testMixture * ((this.testMinTPS + this.testMaxTPS) / 2);
                #% testTransactionCount(3,:) = this.testMixture * this.testMaxTPS;

		testTPS = this.testSampleTPS;
		testTransactionCount = this.testSampleTransactionCount;

		range = (1:15000)'';
		maxFlushRate = 6000;

		#%cfFlushRateApprox_conf = struct('io_conf', this.testConfig.io_conf, 'workloadName', this.workloadName);
		cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf,"TPCC");

		myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testMixture, maxFlushRate, cfFlushRateApprox_conf);
		modelP = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.transactionCount);

		#%[testMaxThroughputIdx testMaxThroughput] = findMaxThroughput(testTPS);
		trainMaxThroughputIdx, trainMaxThroughput = findMaxThroughput(this.trainConfig.TPSUngrouped);

		actualThr = [];

		idx=[]
		if isdefined(:trainMaxThroughputIdx) && !isempty(trainMaxThroughputIdx)
			idx=(1:trainMaxThroughputIdx)';
		else
			idx=(1:size(this.trainConfig.averageCpuUsage,1))';
		end

		#%CPU-based throughput with classification
		cpuC = barzanLinInvoke(modelP, range*this.testMixture);

		tempLen=length(cpuC)
		cpuCLThroughput=[]
		while tempLen>=1
			if cpuC[tempLen]>88 && cpuC[tempLen]<90
				cpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(cpuC)
		cpuCUThroughput=[]
		while tempLen>=1
			if cpuC[tempLen]>98 && cpuC[tempLen]<100
				cpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		#%CPU-based without classification
		cpuTModel = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.TPS);
		cpuT = barzanLinInvoke(cpuTModel, range);

		tempLen=length(cpuT)
		cpuTLThroughput=[]
		while tempLen>=1
			if cpuT[tempLen]>88 && cpuT[tempLen]<90
				cpuTLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(cpuT)
		cpuTUThroughput=[]
		while tempLen>=1
			if cpuT[tempLen]>98 && cpuT[tempLen]<100
				cpuTUThroughput = tempLen
				break
			end
			tempLen -= 1
		end
		myModelP = barzanLinSolve(this.trainConfig.averageCpuUsage[vec(idx),:], this.trainConfig.transactionCount[vec(idx),:]);
		myCpuC = barzanLinInvoke(myModelP, range*this.testMixture);

		tempLen=length(myCpuC)
		myCpuCLThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>88 && myCpuC[tempLen]<90
				myCpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCUThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>98 && myCpuC[tempLen]<100
				myCpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCLThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>44 && myCpuC[tempLen]<45
				myCpuCLThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		tempLen=length(myCpuC)
		myCpuCUThroughput=[]
		while tempLen>=1
			if myCpuC[tempLen]>59 && myCpuC[tempLen]<50
				myCpuCUThroughput = tempLen
				break
			end
			tempLen -= 1
		end

		#Our IO-based throughput
		#% cfFlushRateApprox_conf = struct('io_conf', this.testConfig.io_conf, 'workloadName', this.workloadName);
		cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf, "TPCC");
		myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testMixture, maxFlushRate, cfFlushRateApprox_conf);

		#%Lock-based throughput
		#%getConcurrencyLebel_conf = struct('lock_conf', this.testConfig.lock_conf, 'workloadName', this.workloadName);
		getConcurrencyLebel_conf = typeOfConcurrencyLevel(this.lockConf, "TPCC");
		concurrencyThroughput = findClosestValue("getConcurrencyLevel", (1:10000)''*this.testMixture, 160, getConcurrencyLebel_conf);

		#%Linear IO-based throughput
		modelFlushRate = barzanLinSolve(this.trainConfig.pagesFlushed, this.trainConfig.transactionCount);
		linFlushRate = barzanLinInvoke(modelFlushRate, range*this.testMixture);
		for i=length(linFlushRate):-1:1
			if linFlushRate[i]<maxFlushRate
				linFlushRateThroughput = i
				break
			end
		end
		if isempty(linFlushRateThroughput)
			linFlushRateThroughput=0
		end

		myMaxThroughput1 = min(myCpuCLThroughput, myFlushRateThroughput, concurrencyThroughput);
		if myMaxThroughput1 == myCpuCLThroughput
			PredReasonIdx1 = 1
		elseif myMaxThroughput1 == myFlushRateThroughput
			PredReasonIdx1 = 2
		elseif myMaxThroughput1 == concurrencyThroughput
			PredReasonIdx1 = 3
		end
		myMaxThroughput2 = min(myCpuCUThroughput, myFlushRateThroughput, concurrencyThroughput);
		if myMaxThroughput2 == myCpuCUThroughput
			PredReasonIdx2 = 1
		elseif myMaxThroughput2 == myFlushRateThroughput
			PredReasonIdx2 = 2
		elseif myMaxThroughput2 == concurrencyThroughput
			PredReasonIdx2 = 3
		end

		Xdata = {[1:size(testTPS, 1)]''};
		Ydata = {testTPS};
		legends = {"Signal Generated From User Input"};

		num_row = size(Xdata[1], 1);
		num_col = size(Xdata[1], 2);

		if !isempty(actualThr)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([actualThr], num_row, num_col));
			push!(legends,"Actual Max Throughput");
		end
		if !isempty(cpuCLThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([cpuCLThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on adjusted LR for CPU+classification");
		end
		if !isempty(cpuCUThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([cpuCUThroughput], num_row, num_col));
			push!(legends,"Max Throughput on LR for CPU+classification");
		end
		if !isempty(cpuTLThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([cpuTLThroughput], num_row, num_col));
			push!(legends,"Max Throughput on adjusted LR for CPU");
		end
		if !isempty(cpuTUThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([cpuTUThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on LR for CPU");
		end
		if !isempty(myFlushRateThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([myFlushRateThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on our flush rate model");
		end
		if !isempty(linFlushRateThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([linFlushRateThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on LR for flush rate");
		end
		if !isempty(concurrencyThroughput)
			push!(Xdata,[1:size(testTPS, 1)]'')
			push!(Ydata,repmat([concurrencyThroughput], num_row, num_col));
			push!(legends,"Max Throughput based on our contention model");
		end

		title = "Max Throughput Prediction";
		Ylabel = "TPS";
		Xlabel = "Time";
	end
	return title,legends,Xdata,Ydata,Xlabel,Ylabel,meanAbsError,meanRelError,errorHeader,extra
end # end function

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
using Test
using NQCBase
using NQCModels
using LinearAlgebra
using FiniteDiff

function finite_difference_gradient(model::NQCModels.AdiabaticModels.AdiabaticModel, R)
    f(x) = potential(model, x)
    FiniteDiff.finite_difference_gradient(f, R)
end

function finite_difference_gradient(model::NQCModels.DiabaticModels.DiabaticModel, R)
    f(x, i, j) = potential(model, x)[i,j]
    grad = [Hermitian(zeros(nstates(model), nstates(model))) for i in CartesianIndices(R)]
    for k in eachindex(R)
        for i=1:nstates(model)
            for j=1:nstates(model)
                grad[k].data[i,j] = FiniteDiff.finite_difference_gradient(x->f(x,i,j), R)[k]
            end
        end
    end
    grad
end

function test_model(model::NQCModels.Model, atoms; rtol=1e-5)
    R = rand(ndofs(model), atoms)
    D = derivative(model, R)
    finite_diff = finite_difference_gradient(model, R)
    return isapprox(finite_diff, D, rtol=rtol)
end

function test_model(model::NQCModels.FrictionModels.AdiabaticFrictionModel, atoms)
    R = rand(ndofs(model), atoms)
    D = derivative(model, R)
    friction(model, R)
    return finite_difference_gradient(model, R) ≈ D
end

@testset "Potential abstraction" begin
    struct TestModel <: NQCModels.Model end

    NQCModels.ndofs(::TestModel) = 3
    @test_throws MethodError potential(TestModel(), rand(3,1))

    NQCModels.ndofs(::TestModel) = 1
    NQCModels.potential(::TestModel, ::Real) = 1
    NQCModels.potential(::TestModel, ::AbstractVector) = 2
    @test potential(TestModel(), rand(1,1)) == 1
    @test potential(TestModel(), rand(1,2)) == 2
end

@testset "Plot" begin
    using Plots
    plot(-10:0.1:10, Harmonic())
    plot(-10:0.1:10, DoubleWell())
end

@testset "DiatomicHarmonic" begin
    model = DiatomicHarmonic()
    @test test_model(model, 2)

    R = [0 0; 0 0; 1 0]
    @test potential(model, R) ≈ 0
    R = [sqrt(3) 0; sqrt(3) 0; sqrt(3) 0]
    @test potential(model, R) ≈ 2
end

@testset "AdiabaticModels" begin
    @test test_model(Harmonic(), 10)
    @test test_model(Free(), 10)
    @test test_model(AveragedPotential((Harmonic(), Harmonic()), zeros(1, 10)), 10)
    @test test_model(BosonBath(OhmicSpectralDensity(2.5, 0.1), 10), 10)
    @test test_model(DarlingHollowayElbow(), 2)
end

@testset "DiabaticModels" begin
    @test test_model(DoubleWell(), 1)
    @test test_model(TullyModelOne(), 1)
    @test test_model(TullyModelTwo(), 1)
    @test test_model(TullyModelThree(), 1)
    @test test_model(Scattering1D(), 1)
    @test test_model(ThreeStateMorse(), 1)
    @test test_model(SpinBoson(DebyeSpectralDensity(0.25, 0.5), 10, 1.0, 1.0), 10)
    @test test_model(OuyangModelOne(), 1)
    @test test_model(GatesHollowayElbow(), 2)
    @test test_model(MiaoSubotnik(), 1)
    @test test_model(AnanthModelOne(), 1)
    @test test_model(AnanthModelTwo(), 1)
end

@testset "FrictionModels" begin
    @test test_model(CompositeFrictionModel(Free(2), ConstantFriction(2, 1)), 3)
    @test test_model(CompositeFrictionModel(Free(3), RandomFriction(3)), 3)
end

@testset "JuLIP" begin
    using JuLIP: JuLIP
    at = JuLIP.bulk(:Si, cubic=true)
    deleteat!(at, 1)
    JuLIP.set_calculator!(at, JuLIP.StillingerWeber())
    model = AdiabaticModels.JuLIPModel(at)
    @test test_model(model, length(at))
end

@testset "ASE" begin
    using PyCall

    ase = pyimport("ase")

    h2 = ase.Atoms("H2", [(0, 0, 0), (0, 0, 0.74)])
    h2.center(vacuum=2.5)

    @testset "EMT" begin
        emt = pyimport("ase.calculators.emt")
        h2.calc = emt.EMT()
        model = AdiabaticASEModel(h2)
        @test test_model(model, 2)
    end

    @testset "GPAW" begin
        gpaw = pyimport("gpaw")
        h2.calc = gpaw.GPAW(xc="PBE", mode=gpaw.PW(300), txt="h2.txt")
        model = AdiabaticASEModel(h2)
        @test test_model(model, 2, rtol=1e-3)
    end

end


text:
#Additional constructors for the PlanetEnvironment with different number of inputs
PlanetEnvironment() = PlanetEnvironment(CentralGravity)
PlanetEnvironment(gravfunc::Function) = PlanetEnvironment(gravfunc, SphericalAltitude)
PlanetEnvironment(gravfunc::Function, planetfunc::Function) =  PlanetEnvironment(gravfunc, planetfunc, false, x-> nothing)
function PlanetEnvironment(gravfunc::Function, planetfunc::Function, atmosphereflag::Bool, atmospherefunc::Function)
    return PlanetEnvironment(gravfunc, planetfunc, atmosphereflag, atmospherefunc, false, x->nothing)
end

#Define Constructor for Planet and Moon Types
function Planet(Name::Type{<:abstractCelestialBody})
    return Planet(Name, PlanetEnvironment())
end

function NaturalSatellite(Name::Type{<:abstractCelestialBody})
    return NaturalSatellite(Name, PlanetEnvironment())
end


text:
type ScreenPower
    current::Nullable{Bool}
end
ScreenPower() = ScreenPower(Nullable{Bool}())

function setOn(power::ScreenPower)
    if !get(power.current, false)
        run(`tvservice -p`)
        sleep(.5)
        run(`xset dpms force on`)
        power.current = true
    end
end

function setOff(power::ScreenPower)
    if get(power.current, true)
        run(`tvservice -o`)
        power.current = false
    end
end

@enum(Direction, IN, OUT)

immutable Gpio
    pin
    gpio_dir

    function Gpio(pin, direction::Direction=IN)
        write("/sys/class/gpio/export", dec(pin))
        g = new(pin, "/sys/class/gpio/gpio$pin/")
        retry(() -> write(joinpath(g.gpio_dir, "direction"),
                          lowercase(string(direction))))
        return g
    end
end

function getState(gpio::Gpio)
    parse(Int, strip(readstring(joinpath(gpio.gpio_dir, "value"))))
end

function main(motionSensorGpioPin::Int, minimumScreenOnTime::Dates.Second)
    run(`xset s off`)
    run(`xset dpms 0 0 0`)
    gpio = Gpio(motionSensorGpioPin)
    power = ScreenPower()

    info("watching pin ", gpio.pin)
    lastMotionTime = now()
    while true
        motion = getState(gpio) == 1
        if motion
            lastMotionTime = now()
        end

        if now() - lastMotionTime < minimumScreenOnTime
            setOn(power)
        else
            setOff(power)
        end
        sleep(1)
    end
end
main(parse(Int, ARGS[1]), Dates.Second(parse(Float64, ARGS[2])))

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# saveTable("test.txt", "<tab:testJulia1>", [[1 2]; [3 4]])
# saveTable("test.txt", "<tab:testJulia2>", rand(2, 2))

using DelimitedFiles

function saveTable(outfile, tag, outmatrix)
    open(outfile, "a") do file
        write(file, "$(tag)\n")
        writedlm(file, outmatrix, "\t")
    end
end

max_stars_repo_name: kylebarron/tablefill

max_issues_repo_name: kylebarron/tablefill

max_forks_repo_name: kylebarron/tablefill


text:
solve(x, offset) = sum(x[x .== circshift(x, offset)])
d = parse.(Int, split(readline(open("input.txt")), ""))
println(solve(d, 1))
println(solve(d, Int(length(d) / 2)))


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

JLLWrappers.@generate_main_file_header("LLVMOpenMP")
JLLWrappers.@generate_main_file("LLVMOpenMP", UUID("1d63c593-3942-5779-bab2-d838dc0a180e"))
end  # module LLVMOpenMP_jll


text:
#-------------IFS.jl-----------------------------------------------------------#
#
# Purpose: quick test of iterated function system fractals
#
#------------------------------------------------------------------------------#

using Images

mutable struct Pixel
    color::Vector{Float64}
    val::Float64
end

# Functions for fractal flames
function linear(point::Array{Float64, 1}, c::Float64, color::Array{Float64, 1})
    ret_point = zeros(5)
    ret_point[1] = point[1] * c
    ret_point[2] = point[2] * c
    ret_point[3] = color[1]
    ret_point[4] = color[2]
    ret_point[5] = color[3]

    return ret_point
end

function heart(point::Array{Float64, 1}, color::Array{Float64, 1})
    ret_point = zeros(5)
    r = sqrt(point[1]^2 + point[2]^2)
    theta = atan2(point[2],point[1])
    ret_point[1] = r*(sin(r*theta))
    ret_point[2] = -r*(cos(r*theta))
    ret_point[3] = color[1]
    ret_point[4] = color[2]
    ret_point[5] = color[3]

    return ret_point
end

function sinusoidal(point::Array{Float64, 1}, color::Array{Float64, 1})
    ret_point = zeros(5)
    ret_point[1] = sin(point[1])
    ret_point[2] = cos(point[2])
    ret_point[3] = color[1]
    ret_point[4] = color[2]
    ret_point[5] = color[3]

    return ret_point
end

function swirl(point::Array{Float64, 1}, color::Array{Float64, 1})
    ret_point = zeros(5)
    r = sqrt(point[1]^2 + point[2]^2)
    ret_point[1] = point[1]*sin(r^2) - point[2]*cos(r^2)
    ret_point[2] = point[1]*cos(r^2) + point[2]*sin(r^2)
    ret_point[3] = color[1]
    ret_point[4] = color[2]
    ret_point[5] = color[3]

    return ret_point
end

function sierpensky_flame(point::Array{Float64, 1}, chosen_function::Int64)
    if (chosen_function == 1)
        point[1:2] *= 0.5
        point[3] = 0.5
        point[4] = 0.5
        point[5] = 0.0
    elseif(chosen_function == 2)
        point[1:2] = [(point[1] - 1)/2, (point[2])/2]
        point[3] = 0.0
        point[4] = 0.5
        point[5] = 0.5
    else
        point[1:2] = [(point[1])/2, (point[2]-1)/2]
        point[3] = 0.5
        point[4] = 0.0
        point[5] = 0.5
    end

    return point
end

function barnsley_flame(point::Array{Float64, 1}, chosen_function::Int64)
    if (chosen_function == 1)
        point[1:2] = [0 0; 0 0.16]*point[1:2]
        point[3] = 1.0
        point[4] = 1.0
        point[5] = 1.0
    elseif(chosen_function == 2)
        point[1:2] = [0.85 0.04; -0.04 0.85]*point[1:2] + [0, 1.60]
        point[3] = 1.0
        point[4] = 0.0
        point[5] = 0.0
    elseif(chosen_function == 3)
        point[1:2] = [0.2 -0.26; 0.23 0.22]*point[1:2] + [0, 1.60]
        point[3] = 0.0
        point[4] = 1.0
        point[5] = 0.0
    else
        point[1:2] = [-0.15 0.28; 0.26 0.24]*point[1:2] + [0, 0.44]
        point[3] = 0.0
        point[4] = 0.0
        point[5] = 1.0
    end

    return point

end

function write_image(p::Array{Pixel,2})
    a = Array{RGB{N0f8},2}(size(p,1), size(p,2))
    for i = 1:length(pixels)
        a[i] = RGB{N0f8}(pixels[i].color[1], pixels[i].color[2],
                         pixels[i].color[3])
    end
    save("fractal.png", a)
end

function barnsley_fern(n::Int64)
    point = [0.0, 0.0]
    f = open("barnsley.dat", "w")
    for i = 1:n
        rnd = rand()
        if (rnd <= 0.01)
            point = [0 0; 0 0.16]*point
        elseif(rnd > 0.01 && rnd <= 0.86)
            point = [0.85 0.04; -0.04 0.85]*point + [0, 1.60]
        elseif(rnd > 0.86 && rnd <= 0.93)
            point = [0.2 -0.26; 0.23 0.22]*point + [0, 1.60]
        else
            point = [-0.15 0.28; 0.26 0.24]*point + [0, 0.44]
        end

        println(f, point[1], '\t', point[2])
    end
    close(f)
end

function sierpensky(n::Int64)
    element_x = rand()*2 - 1
    element_y = rand()*2 - 1
    point = [element_x element_y]

    f = open("sierpensky.dat", "w")

    for i = 1:n
        if (i > 20)
            println(f, point[1], '\t', point[2])
        end
        rnd = rand()

        if (rnd <= 0.33)
            point = 0.5*point
        elseif(rnd > 0.33 && rnd <= 0.66)
            point[1] = (point[1] - 1)/2
            point[2] = (point[2])/2
        else
            point[1] = (point[1])/2
            point[2] = (point[2]-1)/2
        end
    end

    close(f)
end

function flame(n::Int64, resx::Int64, resy::Int64,
               scalex::Float64, scaley::Float64, gamma::Float64)
    element_x = rand()*scalex*0.5 - 1
    element_y = rand()*scaley*0.5 - 1
    point = [element_x, element_y, 0.0, 0.0, 0.0]

    pixels = [Pixel([0.0, 0.0, 0.0], 0.0) for i = 1:resx, j = 1:resy]

    pixel_file = open("pixels.dat","w")

    color = zeros(3)
    for i = 1:n
        rnd = rand()
        r = sqrt(point[1]^2 + point[2]^2)
        if (abs(r) < 0.0001)
            r = 0.0001
        end

#=
        if (rnd <= 0.01)
            point = barnsley_flame(point, 1)
        elseif(rnd > 0.01 && rnd <= 0.86)
            point = barnsley_flame(point, 2)
        elseif(rnd > 0.86 && rnd <= 0.93)
            point = barnsley_flame(point, 3)
        else
            point = barnsley_flame(point, 4)
        end
=#

        if (rnd <= 0.33)
            #point = sierpensky_flame(point, 1)
            point = sinusoidal(point, [0.0, 0.0, 1.0])
        elseif (rnd > 0.33 && rnd <= 0.66)
            #point = sierpensky_flame(point, 2)
            point = linear(point, 0.5, [1.0, 0.0, 0.0])
        else
            #point = sierpensky_flame(point, 3)
            point = swirl(point, [1.0, 1.0, 1.0])
        end

        # Final function attempts
        final_point = heart(point, [1.0, 0.0, 1.0])

        # blending the colors
        if (i > 20)
            xval = Int(floor((final_point[1]+scalex*0.5) / (scalex / resx)))
            yval = Int(floor((final_point[2]+scaley*0.5) / (scaley / resy)))
            if (xval < resx && yval < resy && xval > 0 && yval > 0)
                pixels[yval,xval].val += 1
                pixels[yval,xval].color = (color+point[3:5])*0.5
                color = pixels[yval,xval].color
                pixels[yval,xval].color = (color+final_point[3:5])*0.5
            end
        end

    end

    # find maximum in a dumb way
    max = 0;
    for pixel in pixels
        if (pixel.val > 0)
            pixel.val = log10(pixel.val)
        end
        if pixel.val > max
            max = pixel.val
        end
    end

    for i = 1:resy-1
        index_string = ""
        for j = 1:resx-1
            index_string *= string(log10(pixels[i,j].val), '\t')
        end
        println(pixel_file, index_string)
    end
    close(pixel_file)

    # Setting log plot and changing brightness
    for pixel in pixels
        pixel.val /= max
        if (pixel.val > 0.0000001)
            pixel.color = (pixel.color * pixel.val).^(1/gamma)
            for c in pixel.color
                if c > 1
                    c = 1
                end
            end
        end
    end

    return pixels

end

pixels = flame(100000000, 1000, 1000, 4., 4., 1.0)
#pixels = flame(1000000, 1000, 1000, 2., 2., 2.2)
write_image(pixels)
#sierpensky(1000000)
#barnsley_fern(100000)

max_stars_repo_name: shiffman/simuleios

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
import Distributions: cdf, Normal

N(x) = cdf(Normal(), x)

function d1(S, K, τ, r, σ, q)
    return (log(S / K) + (r - q + σ^2 / 2) * τ) / (σ * √(τ))
end

function d2(S, K, τ, r, σ, q)
    return d1(S, K, τ, r, σ, q) - σ * √(τ)
end

"""
    eurocall(;S=1.,K=1.,τ=1,r,σ,q=0.)

Calculate the Black-Scholes implied option price for a european call, where:

- `S` is the current asset price
- `K` is the strike or exercise price
- `τ` is the time remaining to maturity (can be typed with \\tau[tab])
- `r` is the continuously compounded risk free rate
- `σ` is the (implied) volatility (can be typed with \\sigma[tab])
- `q` is the continuously paid dividend rate

Rates should be input as rates (not percentages), e.g.: `0.05` instead of `5` for a rate of five percent.

!!! Experimental: this function is well-tested, but the derivatives functionality (API) may change in a future version of ActuaryUtilities.

# Extended Help

This is the same as the formulation presented in the [dividend extension of the BS model in Wikipedia](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model#Black%E2%80%93Scholes_equation).

## Other general comments:

- Swap/OIS curves are generally better sources for `r` than government debt (e.g. US Treasury) due to the collateralized nature of swap instruments.
- (Implied) volatility is characterized by a curve that is a function of the strike price (among other things), so take care when using
- Yields.jl can assist with converting rates to continuously compounded if you need to perform conversions.

"""
function eurocall(; S = 1.0, K = 1.0, τ = 1, r, σ, q = 0.0)
    d₁ = d1(S, K, τ, r, σ, q)
    d₂ = d2(S, K, τ, r, σ, q)
    return (N(d₁) * S * exp(τ * (r - q)) - N(d₂) * K) * exp(-r * τ)
end

"""
    europut(;S=1.,K=1.,τ=1,r,σ,q=0.)

Calculate the Black-Scholes implied option price for a european call, where:

- `S` is the current asset price
- `K` is the strike or exercise price
- `τ` is the time remaining to maturity (can be typed with \\tau[tab])
- `r` is the continuously compounded risk free rate
- `σ` is the (implied) volatility (can be typed with \\sigma[tab])
- `q` is the continuously paid dividend rate

Rates should be input as rates (not percentages), e.g.: `0.05` instead of `5` for a rate of five percent.

!!! Experimental: this function is well-tested, but the derivatives functionality (API) may change in a future version of ActuaryUtilities.

# Extended Help

This is the same as the formulation presented in the [dividend extension of the BS model in Wikipedia](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model#Black%E2%80%93Scholes_equation).

## Other general comments:

- Swap/OIS curves are generally better sources for `r` than government debt (e.g. US Treasury) due to the collateralized nature of swap instruments.
- (Implied) volatility is characterized by a curve that is a function of the strike price (among other things), so take care when using
- Yields.jl can assist with converting rates to continuously compounded if you need to perform conversions.

"""
function europut(; S = 1.0, K = 1.0, τ = 1, r, σ, q = 0.0)
    d₁ = d1(S, K, τ, r, σ, q)
    d₂ = d2(S, K, τ, r, σ, q)
    return (N(-d₂) * K - N(-d₁) * S * exp(τ * (r - q))) * exp(-r * τ)
end

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
using BinaryBuilder

name = "ODEInterface"
version = v"0.0.1"

function get_autogenerated_build_script()
  return raw"""
    cd $WORKSPACE/srcdir/ODEInterfac*/src
    if [[ $target == *mingw* ]] ; then

    # dopri
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5.o ./dopri5.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5_i32.o ./dopri5.f
    gfortran -shared -o ./dopri5.dll ./dopri5.o
    gfortran -shared -o ./dopri5_i32.dll ./dopri5_i32.o
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853.o ./dop853.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853_i32.o ./dop853.f
    gfortran -shared -o ./dop853.dll ./dop853.o
    gfortran -shared -o ./dop853_i32.dll ./dop853_i32.o

    # odex
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex.o ./odex.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex_i32.o ./odex.f
    gfortran -shared -o ./odex.dll ./odex.o
    gfortran -shared -o ./odex_i32.dll ./odex_i32.o

    # lapack
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack.o ./dc_lapack.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack_i32.o ./dc_lapack.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack.o ./lapack.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack_i32.o ./lapack.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc.o ./lapackc.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc_i32.o ./lapackc.f

    # radau
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5.o ./radau5.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5_i32.o ./radau5.f
    gfortran -shared -o ./radau5.dll ./radau5.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -o ./radau5_i32.dll ./radau5_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau.o ./radau.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau_i32.o ./radau.f
    gfortran -shared -o ./radau.dll ./radau.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -o ./radau_i32.dll ./radau_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # seulex
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex.o ./seulex.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex_i32.o ./seulex.f
    gfortran -shared -o ./seulex.dll ./seulex.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -o ./seulex_i32.dll ./seulex_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # rodas
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas.o ./rodas.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas_i32.o ./rodas.f
    gfortran -shared -o ./rodas.dll ./rodas.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -o ./rodas_i32.dll ./rodas_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # slatec
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec.o ./slatec.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec_i32.o ./slatec.f

    # ddeabm
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm.o ./ddeabm.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm_i32.o ./ddeabm.f
    gfortran -shared -o ./ddeabm.dll ./ddeabm.o ./slatec.o
    gfortran -shared -o ./ddeabm_i32.dll ./ddeabm_i32.o ./slatec_i32.o

    # ddebdf
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf.o ./ddebdf.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf_i32.o ./ddebdf.f
    gfortran -shared -o ./ddebdf.dll ./ddebdf.o ./slatec.o
    gfortran -shared -o ./ddebdf_i32.dll ./ddebdf_i32.o ./slatec_i32.o

    # bvpsol
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol.o ./bvpsol.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol_i32.o ./bvpsol.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol.o ./linalg_bvpsol.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol_i32.o ./linalg_bvpsol.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst.o ./zibconst.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst_i32.o ./zibconst.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol.o ./ma28_bvpsol.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol_i32.o ./ma28_bvpsol.f
    gfortran -shared -o ./bvpsol.dll ./bvpsol.o ./linalg_bvpsol.o ./zibconst.o ./ma28_bvpsol.o
    gfortran -shared -o ./bvpsol_i32.dll ./bvpsol_i32.o ./linalg_bvpsol_i32.o ./zibconst_i32.o ./ma28_bvpsol_i32.o

    # colnew
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew.o ./colnew.f
    gfortran -c -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew_i32.o ./colnew.f
    gfortran -shared -o ./colnew.dll ./colnew.o
    gfortran -shared -o ./colnew_i32.dll ./colnew_i32.o

    # bvpm2
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvp_la-2.o ./bvp_la-2.f
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -std=f2008 -o ./bvp_m-2.o ./bvp_m-2.f90
    gfortran -c -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -Wall -Wextra -Wimplicit-interface -std=f2008ts -o ./bvp_m_proxy.o ./bvp_m_proxy.f90
    gfortran -shared -o ./bvp_m_proxy.dll ./bvp_m_proxy.o ./bvp_m-2.o ./bvp_la-2.o
    cp *.dll $libdir

    elif [[ $target == *apple* ]] ; then

    # dopri
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5.o ./dopri5.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5_i32.o ./dopri5.f
    gfortran -shared -fPIC -o ./dopri5.dylib ./dopri5.o
    gfortran -shared -fPIC -o ./dopri5_i32.dylib ./dopri5_i32.o
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853.o ./dop853.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853_i32.o ./dop853.f
    gfortran -shared -fPIC -o ./dop853.dylib ./dop853.o
    gfortran -shared -fPIC -o ./dop853_i32.dylib ./dop853_i32.o

    # odex
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex.o ./odex.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex_i32.o ./odex.f
    gfortran -shared -fPIC -o ./odex.dylib ./odex.o
    gfortran -shared -fPIC -o ./odex_i32.dylib ./odex_i32.o

    # lapack
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack.o ./dc_lapack.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack_i32.o ./dc_lapack.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack.o ./lapack.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack_i32.o ./lapack.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc.o ./lapackc.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc_i32.o ./lapackc.f

    # radau
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5.o ./radau5.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5_i32.o ./radau5.f
    gfortran -shared -fPIC -o ./radau5.dylib ./radau5.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./radau5_i32.dylib ./radau5_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau.o ./radau.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau_i32.o ./radau.f
    gfortran -shared -fPIC -o ./radau.dylib ./radau.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./radau_i32.dylib ./radau_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # seulex
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex.o ./seulex.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex_i32.o ./seulex.f
    gfortran -shared -fPIC -o ./seulex.dylib ./seulex.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./seulex_i32.dylib ./seulex_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # rodas
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas.o ./rodas.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas_i32.o ./rodas.f
    gfortran -shared -fPIC -o ./rodas.dylib ./rodas.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./rodas_i32.dylib ./rodas_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # slatec
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec.o ./slatec.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec_i32.o ./slatec.f

    # ddeabm
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm.o ./ddeabm.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm_i32.o ./ddeabm.f
    gfortran -shared -fPIC -o ./ddeabm.dylib ./ddeabm.o ./slatec.o
    gfortran -shared -fPIC -o ./ddeabm_i32.dylib ./ddeabm_i32.o ./slatec_i32.o

    # ddebdf
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf.o ./ddebdf.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf_i32.o ./ddebdf.f
    gfortran -shared -fPIC -o ./ddebdf.dylib ./ddebdf.o ./slatec.o
    gfortran -shared -fPIC -o ./ddebdf_i32.dylib ./ddebdf_i32.o ./slatec_i32.o

    # bvpsol
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol.o ./bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol_i32.o ./bvpsol.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol.o ./linalg_bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol_i32.o ./linalg_bvpsol.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst.o ./zibconst.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst_i32.o ./zibconst.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol.o ./ma28_bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol_i32.o ./ma28_bvpsol.f
    gfortran -shared -fPIC -o ./bvpsol.dylib ./bvpsol.o ./linalg_bvpsol.o ./zibconst.o ./ma28_bvpsol.o
    gfortran -shared -fPIC -o ./bvpsol_i32.dylib ./bvpsol_i32.o ./linalg_bvpsol_i32.o ./zibconst_i32.o ./ma28_bvpsol_i32.o

    # colnew
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew.o ./colnew.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew_i32.o ./colnew.f
    gfortran -shared -fPIC -o ./colnew.dylib ./colnew.o
    gfortran -shared -fPIC -o ./colnew_i32.dylib ./colnew_i32.o

    # bvpm2
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvp_la-2.o ./bvp_la-2.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -std=f2008 -o ./bvp_m-2.o ./bvp_m-2.f90
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -Wall -Wextra -Wimplicit-interface -std=f2008ts -o ./bvp_m_proxy.o ./bvp_m_proxy.f90
    gfortran -shared -fPIC -o ./bvp_m_proxy.dylib ./bvp_m_proxy.o ./bvp_m-2.o ./bvp_la-2.o
    cp *.dylib $libdir

    else

    # dopri
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5.o ./dopri5.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dopri5_i32.o ./dopri5.f
    gfortran -shared -fPIC -o ./dopri5.so ./dopri5.o
    gfortran -shared -fPIC -o ./dopri5_i32.so ./dopri5_i32.o
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853.o ./dop853.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dop853_i32.o ./dop853.f
    gfortran -shared -fPIC -o ./dop853.so ./dop853.o
    gfortran -shared -fPIC -o ./dop853_i32.so ./dop853_i32.o

    # odex
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex.o ./odex.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./odex_i32.o ./odex.f
    gfortran -shared -fPIC -o ./odex.so ./odex.o
    gfortran -shared -fPIC -o ./odex_i32.so ./odex_i32.o

    # lapack
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack.o ./dc_lapack.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./dc_lapack_i32.o ./dc_lapack.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack.o ./lapack.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapack_i32.o ./lapack.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc.o ./lapackc.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./lapackc_i32.o ./lapackc.f

    # radau
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5.o ./radau5.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau5_i32.o ./radau5.f
    gfortran -shared -fPIC -o ./radau5.so ./radau5.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./radau5_i32.so ./radau5_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau.o ./radau.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./radau_i32.o ./radau.f
    gfortran -shared -fPIC -o ./radau.so ./radau.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./radau_i32.so ./radau_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # seulex
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex.o ./seulex.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./seulex_i32.o ./seulex.f
    gfortran -shared -fPIC -o ./seulex.so ./seulex.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./seulex_i32.so ./seulex_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # rodas
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas.o ./rodas.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./rodas_i32.o ./rodas.f
    gfortran -shared -fPIC -o ./rodas.so ./rodas.o ./dc_lapack.o ./lapack.o ./lapackc.o
    gfortran -shared -fPIC -o ./rodas_i32.so ./rodas_i32.o ./dc_lapack_i32.o ./lapack_i32.o ./lapackc_i32.o

    # slatec
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec.o ./slatec.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./slatec_i32.o ./slatec.f

    # ddeabm
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm.o ./ddeabm.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddeabm_i32.o ./ddeabm.f
    gfortran -shared -fPIC -o ./ddeabm.so ./ddeabm.o ./slatec.o
    gfortran -shared -fPIC -o ./ddeabm_i32.so ./ddeabm_i32.o ./slatec_i32.o

    # ddebdf
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf.o ./ddebdf.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ddebdf_i32.o ./ddebdf.f
    gfortran -shared -fPIC -o ./ddebdf.so ./ddebdf.o ./slatec.o
    gfortran -shared -fPIC -o ./ddebdf_i32.so ./ddebdf_i32.o ./slatec_i32.o

    # bvpsol
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol.o ./bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvpsol_i32.o ./bvpsol.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol.o ./linalg_bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./linalg_bvpsol_i32.o ./linalg_bvpsol.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst.o ./zibconst.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./zibconst_i32.o ./zibconst.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol.o ./ma28_bvpsol.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./ma28_bvpsol_i32.o ./ma28_bvpsol.f
    gfortran -shared -fPIC -o ./bvpsol.so ./bvpsol.o ./linalg_bvpsol.o ./zibconst.o ./ma28_bvpsol.o
    gfortran -shared -fPIC -o ./bvpsol_i32.so ./bvpsol_i32.o ./linalg_bvpsol_i32.o ./zibconst_i32.o ./ma28_bvpsol_i32.o

    # colnew
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew.o ./colnew.f
    gfortran -c -fPIC -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./colnew_i32.o ./colnew.f
    gfortran -shared -fPIC -o ./colnew.so ./colnew.o
    gfortran -shared -fPIC -o ./colnew_i32.so ./colnew_i32.o

    # bvpm2
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -w -std=legacy -o ./bvp_la-2.o ./bvp_la-2.f
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -std=f2008 -o ./bvp_m-2.o ./bvp_m-2.f90
    gfortran -c -fPIC -fdefault-integer-8 -fdefault-real-8 -fdefault-double-8 -Wall -Wextra -Wimplicit-interface -std=f2008ts -o ./bvp_m_proxy.o ./bvp_m_proxy.f90
    gfortran -shared -fPIC -o ./bvp_m_proxy.so ./bvp_m_proxy.o ./bvp_m-2.o ./bvp_la-2.o
    cp *.so $libdir

    fi
  """
end

sources = [
  GitSource("https://github.com/luchr/ODEInterface.jl", "987d0e2cffd1d7f2fd8e378b7292bed4bf451c06"),
]

script = get_autogenerated_build_script()

platforms = expand_gfortran_versions(supported_platforms())
# platforms = filter(p -> !isa(p, MacOS), platforms)
# platforms = [
#  BinaryBuilder.Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"5.0.0")),
# ]

products = [
  LibraryProduct("dopri5",       :libdopri5),
  LibraryProduct("dopri5_i32",   :libdopri5_i32),
  LibraryProduct("dop853",       :libdop853),
  LibraryProduct("dop853_i32",   :libdop853_i32),
  LibraryProduct("odex",         :libodex),
  LibraryProduct("odex_i32",     :libodex_i32),
  LibraryProduct("radau5",       :libradau5),
  LibraryProduct("radau5_i32",   :libradau5_i32),
  LibraryProduct("radau",        :libradau),
  LibraryProduct("radau_i32",    :libradau_i32),
  LibraryProduct("seulex",       :libseulex),
  LibraryProduct("seulex_i32",   :libseulex_i32),
  LibraryProduct("rodas",        :librodas),
  LibraryProduct("rodas_i32",    :librodas_i32),
  LibraryProduct("ddeabm",       :libddeabm),
  LibraryProduct("ddeabm_i32",   :libddeabm_i32),
  LibraryProduct("ddebdf",       :libddebdf),
  LibraryProduct("ddebdf_i32",   :libddebdf_i32),
  LibraryProduct("bvpsol",       :libbvpsol),
  LibraryProduct("bvpsol_i32",   :libbvpsol_i32),
  LibraryProduct("colnew",       :libcolnew),
  LibraryProduct("colnew_i32",   :libcolnew_i32),
  LibraryProduct("bvp_m_proxy",  :libbvp_m_proxy),
]

dependencies = [
  Dependency("CompilerSupportLibraries_jll")
]

build_tarballs(ARGS, name, version, sources, script,
               platforms, products, dependencies; preferred_gcc_version=v"9.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct QLayout <: MemoryLayout end

MemoryLayout(::Type{<:AbstractQ}) = QLayout()

transposelayout(::QLayout) = QLayout()

copy(M::Lmul{QLayout}) = copyto!(similar(M), M)

function copyto!(dest::AbstractArray{T}, M::Lmul{QLayout}) where T
    A,B = M.A,M.B
    if size(dest,1) == size(B,1)
        copyto!(dest, B)
    else
        copyto!(view(dest,1:size(B,1),:), B)
        zero!(@view(dest[size(B,1)+1:end,:]))
    end
    materialize!(Lmul(A,dest))
end

function copyto!(dest::AbstractArray, M::Ldiv{QLayout})
    A,B = M.A,M.B
    copyto!(dest, B)
    materialize!(Ldiv(A,dest))
end

materialize!(M::Ldiv{QLayout}) = materialize!(Lmul(M.A',M.B))


text:
using Documenter, DataVoyager

makedocs(
	modules=[DataVoyager],
	sitename="DataVoyager.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo="github.com/queryverse/DataVoyager.jl.git"
)


text:
using Random
using LinearAlgebra: rank

function test_gram_schmidt(A::Matrix{T}) where T<:Integer
    B = Cones.gram_schmidt(A)
    @test Cones.check_column_span_same(A, B)
    @test Cones.check_columns_orthogonal(B)
end

Random.seed!(1)

println("\t...gram_schmidt")

N = 10         # Trials
m, n = 5, 3
for counter in 0:(N-1)
    # Randomly generate full column rank matrix
    A = rand(Int, m, n) .% 5
    if rank(A) < n continue end
    test_gram_schmidt(A)
end

# Test normalise_rays
#
function test_normalise_rays(rays::Matrix{T}, basis::Matrix{T}, norm_rays::Matrix{T}) where T<:Integer
    orth_basis = Cones.gram_schmidt(basis)
    norm_rays_out = Cones.normalise_rays(rays, basis)
    @test Cones.check_columns_same(norm_rays, norm_rays_out)
end

println("\t...normalise_rays")

# 0D lineality space
rays = [[2, 3, 1] [-1, 3, 2]]
basis = Array{Int}(undef, 3, 0)
norm_rays = rays
test_normalise_rays(rays, basis, norm_rays)

# 1D lineality space
rays = [[2, 0, 2] [0, 1, 1]]
basis= reshape([0, 0, 1], 3, 1)
norm_rays = [[1, 0, 0] [0, 1, 0]]
test_normalise_rays(rays, basis, norm_rays)

rays = [[0, 2, 1] [3, 1, 1]]
basis = reshape([-3, -3, 0], 3, 1)
norm_rays = [[-1, 1, 1] [1, -1, 1]]
test_normalise_rays(rays, basis, norm_rays)

# 2D lineality space
rays = reshape([1, 1, 1], 3, 1)
basis= [[1, 1, 0] [1, 0, 0]]
norm_rays = reshape([0, 0, 1], 3, 1)
test_normalise_rays(rays, basis, norm_rays)


text:
# General structure for reachability methods

# This function performs layer-by-layer propagation
# It is called by all solvers under reachability
# TODO: also called by ReluVal and FastLin, so move to general utils (or to network.jl)
function forward_network(solver, nnet::Network, input::AbstractPolytope)
    reach = input
    for layer in nnet.layers
        reach = forward_layer(solver, layer, reach)
    end
    return reach
end

# This function checks whether the reachable set belongs to the output constraint
# It is called by all solvers under reachability
# Note vertices_list is not defined for HPolytope: to be defined
function check_inclusion(reach::Vector{<:AbstractPolytope}, output::AbstractPolytope)
    for poly in reach
        issubset(poly, output) || return ReachabilityResult(:UNSAT, reach) # TODO shouldn't this be poly, not the full reach?
    end
    return ReachabilityResult(:SAT, similar(reach, 0))
end

function check_inclusion(reach::P, output::AbstractPolytope) where P<:AbstractPolytope
    if issubset(reach, output)
        return ReachabilityResult(:SAT, P[])
    end
    return ReachabilityResult(:UNSAT, [reach])
end


text:
using Plots
using LaTeXStrings

# backend(:pgfplotsx)
backend(:pyplot)
# backend(:gr)

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(5.00)
#Extracting stored varibles

#Add colors
lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

x0 = 1.0
y0 = 1.0

KL(x,y) = x * log(x / x0) + y * log(y / y0) - x - y + x0 + y0
squared_distance(x,y) = 0.5 * ((x - x0)^2 + (y - y0)^2)

ep = 0.0001

x1 = (0.0 + ep) : 0.05 : (2.0 + ep)
y1 = (0.0 + ep) : 0.05 : (2.0 + ep)

x2 = (0.0 + ep) : 0.10 : (2.0 + ep)
y2 = (0.0 + ep) : 0.10 : (2.0 + ep)

srf = wireframe(x1, y1, squared_distance, camera=(60,30), linecolor=silver, linealpha=0.2, xlabel=L"x", ylabel=L"y", zticks=0.0:0.5:2.0, label=L"\mathrm{true} \ f)", thickness_scaling=1.5, linewidth=0.75)

wireframe!(srf, x2, y2,  KL, linecolor=rust, label=L"\mathbb{D}\left(x\middle|y\right)", legend=:topright)

savefig(srf, "./out/figures/divergences.pdf")
run(`pdfcrop --margins '5 5 5 5' ./out/figures/divergences.pdf ./out/figures/divergences.pdf`)

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
module MCAnalyzer
export mark_start, mark_end, analyze

using LLVM
using LLVM.Interop

const optlevel = Ref{Int}()
const analyzer = Ref{Any}()

function __init__()
    @assert LLVM.InitializeNativeTarget() == false
    optlevel[] = Base.JLOptions().opt_level
    analyzer[] = iaca
end

function iaca(march, objfile, asmfile)
    iaca_path = "iaca"
    if haskey(ENV, "IACA_PATH")
        iaca_path = ENV["IACA_PATH"]
    end
    @assert !isempty(iaca_path)
    Base.run(`$iaca_path -arch $march $objfile`)
end

function llvm_mca(march, objfile, asmfile)
    llvm_mca = "llvm-mca"
    if haskey(ENV, "LLVM_MCA_PATH")
        llvm_mca = ENV["LLVM_MCA_PATH"]
    end
    @assert !isempty(llvm_mca)
    Base.run(`$llvm_mca -mcpu $(llvm_march(march)) $asmfile`)
end

include("irgen.jl")

"""
    analyze(func, tt, march = :SKL)

Analyze a given function `func` with the type signature `tt`.
The specific method needs to be annotated with the `IACA` markers.
Supported `march` are :HSW, :BDW, :SKL, and :SKX.

# Example

```julia
function mysum(A)
    acc = zero(eltype(A))
    for i in eachindex(A)
        mark_start()
        @inbounds acc += A[i]
    end
    mark_end()
    return acc
end

analyze(mysum, Tuple{Vector{Float64}})
```

# Advanced usage
## Switching opt-level

```julia
MCAnalyzer.optlevel[] = 3
analyze(mysum, Tuple{Vector{Float64}}, :SKL)
```

## Changing the optimization pipeline

```julia
myoptimize!(tm, mod) = ...
analyze(mysum, Tuple{Vector{Float64}}, :SKL, myoptimize!)
```

## Changing the analyzer tool

```julia
MCAnalyzer.analyzer[] = MCAnalyzer.llvm_mca
analyze(mysum, Tuple{Vector{Float64}})
```
"""
function analyze(@nospecialize(func), @nospecialize(tt), march=:SKL, optimize!::Core.Function = jloptimize!)
    mktempdir() do dir
        objfile = joinpath(dir, "a.out")
        asmfile = joinpath(dir, "a.S")
        mod, _ = irgen(func, tt)
        target_machine(llvm_march(march)) do tm
            optimize!(tm, mod)
            LLVM.emit(tm, mod, LLVM.API.LLVMAssemblyFile, asmfile)
            LLVM.emit(tm, mod, LLVM.API.LLVMObjectFile, objfile)
        end
        Base.invokelatest(analyzer[], march, objfile, asmfile)
    end
    return nothing
end

function llvm_march(march)
    cpus = Dict(
        :HSW => "haswell",
        :BDW => "broadwell",
        :SKL => "skylake",
        :SKX => "skx",
    )
    @assert haskey(cpus, march) "Arch: $march not supported"
    return cpus[march]
end

nameof(f::Core.Function) = String(typeof(f).name.mt.name)

function target_machine(lambda, cpu, features = "")
    triple = LLVM.triple()
    target = LLVM.Target(triple)
    LLVM.TargetMachine(lambda, target, triple, cpu, features)
end

"""
    jloptimize!(tm, mod)

Runs the Julia optimizer pipeline.
"""
function jloptimize!(tm::LLVM.TargetMachine, mod::LLVM.Module)
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        ccall(:jl_add_optimization_passes, Nothing,
              (LLVM.API.LLVMPassManagerRef, Cint),
               LLVM.ref(pm), optlevel[])
        run!(pm, mod)
    end
end

"""
    mark_start()

Insert `iaca` and `llvm-mca` start markers at this position.
"""
function mark_start()
    @asmcall("""
    movl \$\$111, %ebx
    .byte 0x64, 0x67, 0x90
    # LLVM-MCA-BEGIN
    """, "~{memory},~{ebx}", true)
end

"""
    mark_end()

Insert `iaca` and `llvm-mca` end markers at this position.
"""
function mark_end()
    @asmcall("""
    # LLVM-MCA-END
    movl \$\$222, %ebx
    .byte 0x64, 0x67, 0x90
    """, "~{memory},~{ebx}", true)
end

include("reflection.jl")
end # module


text:
include("RPCond.jl")

struct SimpleFluxCalculator
    rs
    bc
end

function calcFluxes(fc::SimpleFluxCalculator,
                    state::HydroSnapshot,
                    edge_velocities,
                    eos,
                    cached)
    rpcond_list = [RPCond(c.density, c.pressure, c.velocity, ss, e)
                    for (c,ss,e) in zip(state.cells,
                                        cached["sound_speeds"],
                                        cached["energies"])]
    bulk_fluxes = [calcFlux(fc.rs,
                            rpcond_list[n],
                            rpcond_list[n+1],
                            edge_velocities[n+1],
                            eos)
                    for n in 1:length(state.cells)-1]
    left_bc = calcLeftFlux(fc.bc,
                           state,
                           edge_velocities,
                           eos,
                           cached)
    right_bc = calcRightFlux(fc.bc,
                             state,
                             edge_velocities,
                             eos,
                             cached)
    return [[left_bc];
            bulk_fluxes;
            [right_bc]]
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
"""
    blended_pairwise_conditional_gradient(f, grad!, lmo, x0; kwargs...)

Implements the BPCG algorithm from [Tsuji, Tanaka, Pokutta](https://arxiv.org/abs/2110.12650).
The method uses an active set of current vertices.
Unlike away-step, it transfers weight from an away vertex to another vertex of the active set.
"""
function blended_pairwise_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    line_search::LineSearchMethod=Adaptive(),
    L=Inf,
    gamma0=0,
    step_lim=20,
    epsilon=1e-7,
    max_iteration=10000,
    print_iter=1000,
    trajectory=false,
    verbose=false,
    linesearch_tol=1e-7,
    emphasis::Emphasis=memory,
    gradient=nothing,
    callback=nothing,
    timeout=Inf,
    print_callback=print_callback,
    renorm_interval=1000,
    lazy=false,
    lazy_tolerance=2.0,
)
    # add the first vertex to active set from initialization
    active_set = ActiveSet([(1.0, x0)])

    return blended_pairwise_conditional_gradient(
        f,
        grad!,
        lmo,
        active_set,
        line_search=line_search,
        L=L,
        gamma0=gamma0,
        step_lim=step_lim,
        epsilon=epsilon,
        max_iteration=max_iteration,
        print_iter=print_iter,
        trajectory=trajectory,
        verbose=verbose,
        linesearch_tol=linesearch_tol,
        emphasis=emphasis,
        gradient=gradient,
        callback=callback,
        timeout=timeout,
        print_callback=print_callback,
        renorm_interval=renorm_interval,
        lazy=lazy,
        lazy_tolerance=lazy_tolerance,
    )
end

"""
    blended_pairwise_conditional_gradient(f, grad!, lmo, active_set::ActiveSet; kwargs...)

Warm-starts BPCG with a pre-defined `active_set`.
"""
function blended_pairwise_conditional_gradient(
    f,
    grad!,
    lmo,
    active_set::ActiveSet;
    line_search::LineSearchMethod=Adaptive(),
    L=Inf,
    gamma0=0,
    step_lim=20,
    epsilon=1e-7,
    max_iteration=10000,
    print_iter=1000,
    trajectory=false,
    verbose=false,
    linesearch_tol=1e-7,
    emphasis::Emphasis=memory,
    gradient=nothing,
    callback=nothing,
    timeout=Inf,
    print_callback=print_callback,
    renorm_interval=1000,
    lazy=false,
    lazy_tolerance=2.0,
)

    # format string for output of the algorithm
    format_string = "%6s %13s %14e %14e %14e %14e %14e %14i\n"

    t = 0
    primal = Inf
    x = get_active_set_iterate(active_set)
    tt = regular
    traj_data = []
    if trajectory && callback === nothing
        callback = trajectory_callback(traj_data)
    end
    time_start = time_ns()

    d = similar(x)

    if line_search isa Shortstep && L == Inf
        println("WARNING: Lipschitz constant not set. Prepare to blow up spectacularly.")
    end

    if line_search isa FixedStep && gamma0 == 0
        println("WARNING: gamma0 not set. We are not going to move a single bit.")
    end

    if verbose
        println("\nBlended Pairwise Conditional Gradient Algorithm.")
        num_type = eltype(x)
        println(
            "EMPHASIS: $emphasis STEPSIZE: $line_search EPSILON: $epsilon MAXITERATION: $max_iteration TYPE: $num_type",
        )
        grad_type = typeof(gradient)
        println(
            "GRADIENTTYPE: $grad_type LAZY: $lazy lazy_tolerance: $lazy_tolerance",
        )
        if emphasis == memory
            println("WARNING: In memory emphasis mode iterates are written back into x0!")
        end
        headers =
            ("Type", "Iteration", "Primal", "Dual", "Dual Gap", "Time", "It/sec", "#ActiveSet")
        print_callback(headers, format_string, print_header=true)
    end

    # likely not needed anymore as now the iterates are provided directly via the active set
    if gradient === nothing
        gradient = similar(x)
    end

    grad!(gradient, x)
    v = compute_extreme_point(lmo, gradient)
    # if !lazy, phi is maintained as the global dual gap
    phi = max(0, fast_dot(x, gradient) - fast_dot(v, gradient))
    local_gap = zero(phi)
    gamma = 1.0

    while t <= max_iteration && phi >= max(epsilon, eps())

        # managing time limit
        time_at_loop = time_ns()
        if t == 0
            time_start = time_at_loop
        end
        # time is measured at beginning of loop for consistency throughout all algorithms
        tot_time = (time_at_loop - time_start) / 1e9

        if timeout < Inf
            if tot_time ≥ timeout
                if verbose
                    @info "Time limit reached"
                end
                break
            end
        end

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

        # compute current iterate from active set
        x = get_active_set_iterate(active_set)
        grad!(gradient, x)

        _, v_local, v_local_loc, _, a_lambda, a, a_loc, _ ,_ =
            active_set_argminmax(active_set, gradient)

        local_gap = fast_dot(gradient, a) - fast_dot(gradient, v_local)

        if !lazy
            v = compute_extreme_point(lmo, gradient)
            dual_gap = fast_dot(gradient, x) - fast_dot(gradient, v)
            phi = dual_gap
        end
        # minor modification from original paper for improved sparsity
        # (proof follows with minor modification when estimating the step)
        if local_gap ≥ phi / lazy_tolerance
            @. d = a - v_local
            vertex_taken = v_local
            gamma_max = a_lambda
            gamma, L = line_search_wrapper(
                line_search,
                t,
                f,
                grad!,
                x,
                d,
                gradient,
                phi,
                L,
                gamma0,
                linesearch_tol,
                step_lim,
                gamma_max,
            )
            # reached maximum of lambda -> dropping away vertex
            if gamma ≈ gamma_max
                tt = drop
                active_set.weights[v_local_loc] += gamma
                deleteat!(active_set, a_loc)
            else # transfer weight from away to local FW
                tt = pairwise
                active_set.weights[a_loc] -= gamma
                active_set.weights[v_local_loc] += gamma
                @assert active_set_validate(active_set)
            end
            active_set_update_iterate_pairwise!(active_set, gamma, v_local, a)
        else # add to active set
            if lazy # otherwise, v computed above already
                v = compute_extreme_point(lmo, gradient)
            end
            vertex_taken = v
            dual_gap = fast_dot(gradient, x) - fast_dot(gradient, v)
            if (!lazy || dual_gap ≥ phi / lazy_tolerance)
                tt = regular
                @. d = x - v
                gamma_max = one(eltype(x))
                gamma, L = line_search_wrapper(
                    line_search,
                    t,
                    f,
                    grad!,
                    x,
                    d,
                    gradient,
                    dual_gap,
                    L,
                    gamma0,
                    linesearch_tol,
                    step_lim,
                    gamma_max,
                )
                # dropping active set and restarting from singleton
                if gamma ≈ 1.0
                    active_set_initialize!(active_set, v)
                else
                    renorm = mod(t, renorm_interval) == 0
                    active_set_update!(active_set, gamma, v, renorm, nothing)
                end
            else # dual step
                tt = dualstep
                # set to computed dual_gap for consistency between the lazy and non-lazy run.
                # that is ok as we scale with the K = 2.0 default anyways
                phi = dual_gap
            end
        end
        if (
            ((mod(t, print_iter) == 0 || tt == dualstep) == 0 && verbose) ||
            callback !== nothing ||
            !(line_search isa Agnostic || line_search isa Nonconvex || line_search isa FixedStep)
        )
            primal = f(x)
        end
        if callback !== nothing
            state = (
                t=t,
                primal=primal,
                dual=primal - dual_gap,
                dual_gap=phi,
                time=tot_time,
                x=x,
                v=vertex_taken,
                active_set_length=length(active_set),
                gamma=gamma,
            )
            callback(state)
        end

        if verbose && (mod(t, print_iter) == 0 || tt == dualstep)
            if t == 0
                tt = initial
            end
            rep = (
                st[Symbol(tt)],
                string(t),
                Float64(primal),
                Float64(primal - dual_gap),
                Float64(dual_gap),
                tot_time,
                t / tot_time,
                length(active_set),
            )
            print_callback(rep, format_string)
            flush(stdout)
        end
        t += 1
    end

    # recompute everything once more for final verfication / do not record to trajectory though for now!
    # this is important as some variants do not recompute f(x) and the dual_gap regularly but only when reporting
    # hence the final computation.
    # do also cleanup of active_set due to many operations on the same set

    if verbose
        x = get_active_set_iterate(active_set)
        grad!(gradient, x)
        v = compute_extreme_point(lmo, gradient)
        primal = f(x)
        phi = fast_dot(x, gradient) - fast_dot(v, gradient)
        tt = last
        rep = (
            st[Symbol(tt)],
            string(t - 1),
            Float64(primal),
            Float64(primal - dual_gap),
            Float64(dual_gap),
            (time_ns() - time_start) / 1.0e9,
            t / ((time_ns() - time_start) / 1.0e9),
            length(active_set),
        )
        print_callback(rep, format_string)
        flush(stdout)
    end
    active_set_renormalize!(active_set)
    active_set_cleanup!(active_set)
    x = get_active_set_iterate(active_set)
    grad!(gradient, x)
    v = compute_extreme_point(lmo, gradient)
    primal = f(x)
    dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
    if verbose
        tt = pp
        rep = (
            st[Symbol(tt)],
            string(t - 1),
            Float64(primal),
            Float64(primal - dual_gap),
            Float64(dual_gap),
            (time_ns() - time_start) / 1.0e9,
            t / ((time_ns() - time_start) / 1.0e9),
            length(active_set),
        )
        print_callback(rep, format_string)
        print_callback(nothing, format_string, print_footer=true)
        flush(stdout)
    end

    return x, v, primal, dual_gap, traj_data, active_set
end


text:
export wrand, wrandn, wrandnp, wsobol, sobol, soboln, sobolnp, wgrid, xgrid, near

##### random

"""
    wrand(d, k) = Weighted(rand(d, k))
Uniformly distributed vectors in `[0,1]^d`, as columns of a `Weighted{Matrix}`
which knows to `clamp` them to this box. Keyword `weights=true` gives weights ∝ `1 .+ rand(k)` rather than constant.
Default is now `k=1`, making a one-column matrix.

    wrandn(d, k) = π .* Weighted(randn(d, k))
Normally distributed `d`-vectors, of mean zero and std. dev. `scale=π` by default.
Keyword `max=10` clamps absolute values to be less than this.

    wrandnp(d, k)
Absolute value of normally distributed...
"""
wrand(d::Int, k::Int=1; weights=true) = Weighted( rand(d,k), weights ? 1 .+ rand(k) : ones(k) ,0,1)

@doc @doc(wrand)
function wrandn(d::Int, k::Int=1; scale=π, weights=true, max=Inf)
    out = Weighted( scale .* randn(d,k), weights ? 1 .+ rand(k) : ones(k) )
    max==Inf ? out : clamp!(out, -max, max)
end

@doc @doc(wrand)
function wrandnp(d::Int, k::Int=1; scale=π, weights=true, max=Inf)
    out = Weighted( abs.(scale .* randn(d,k)), weights ? 1 .+ rand(k) : ones(k), 0,Inf)
    max==Inf ? out : clamp!(out, 0, max)
end

##### sub-random

import Sobol

"""    sobol(d, k)
First `k` entries from the `d`-dimensional Sobol sequence, as columns of a `Weighted{Matrix}`. """
function sobol(d::Int, k::Int=1)
    s = Sobol.SobolSeq(d)
    cols = hcat([Sobol.next!(s) for i=1:k]...)
    Weighted(cols, ones(k), 0,1)
end

using SpecialFunctions: erfinv

"""    soboln(d, k)
Normally distributed `d`-vectors, of mean zero and std. dev. π, generated from Sobol sequence,
as columns of a `Weighted{Matrix}`. """
function soboln(d::Int, k::Int=1; scale=π)
    boxed = sobol(d,k)
    Weighted( scale .* erfinv.( -1 .+ 2 .*array(boxed)), boxed.weights)
end

function sobolnp(d::Int, k::Int=1; scale=π)
    boxed = sobol(d,k)
    Weighted( scale .* erfinv.(boxed.array), boxed.weights, 0, Inf)
end

##### ordered!

wgrid(d::Int, range::AbstractRange) = Weighted(xgrid(d,range), extrema(range)...)

##### not Weighted

"""    xgrid(d, 0:0.1:5)
Gives the matrix whose colums are `d`-vectors, forming a grid of the given range in all dimensions. """
xgrid(d::Int, n::Int) = xgrid(d, 0:n)

function (xgrid(d::Int, range::AbstractRange{T}; big=false)::Matrix{T}) where T
    d==1 && return Matrix(collect(range)')

    n = length(range)
    !big && n^d > 4_000_000 && @error "xgrid is happiest with n^d = < 4×10^6, you have $n^$d... keyword big=true will disable this check"

    intcols = zeros(Int,d,n^d)
    @inbounds for i=1:n^d
        @views digits!(intcols[:,i], i-1; base=n)
    end

    range[1] .+ step(range) .* intcols
end

function near(x::Matrix, y, dist::Real; verbose=false)
    yesno = minimum( pairwise2(x,array(y)) ,dims=2) .<= (dist^2)
    sum(yesno)==0 && begin @warn("no nearby points, returning complete matrix"); return x end
    verbose && @info string("kept ",sum(yesno)," / ",size(x,2)," columns: ",round(100*sum(yesno)/size(x,2),2)," percent")
    x[:,vec(yesno)]
end

using Distances

"""
    pairwise2(x, y=x) = Distances.pairwise(SqEuclidean(), x, y)
Resulting `mat[i,j]` is distance sqared from `x[:,i]` to `y[:,j]`. Implementation varies.
"""
pairwise2(x::Matrix, y::Matrix) = pairwise(SqEuclidean(), x, y, dims=2)
pairwise2(x::Matrix) = pairwise(SqEuclidean(), x)

pairwise2(x::AbsMat, y::AbsMat) = diag(x'*x) .+ diag(y'*y)' .- 2x'*y
function pairwise2(x::AbsMat) # = diag(x'*x) .+ diag(x'*x)' .- 2x'*x
    mat = x'*x
    vec = diag(mat)
    vec .+ vec' .- 2 .* mat
end

# using Requires
#
# function init_create_req()
#     @require Distances = "???" begin
#
# 	pairwise2(x::Matrix, y::Matrix) = pairwise(SqEuclidean(), x, y)
# 	pairwise2(x::Matrix) = pairwise(SqEuclidean(), x)
#
# 	end # @require
# end