text:
maxNbWorkers = 1
using Distributed
while nworkers() < maxNbWorkers
  addprocs(maxNbWorkers - nworkers())
end
println(nworkers())

@everywhere using MSM
@everywhere using DataStructures
@everywhere using OrderedCollections
@everywhere using Random
@everywhere using Distributions
@everywhere using Statistics
@everywhere using LinearAlgebra
using Test
using Pkg
using Optim
using BlackBoxOptim
using DataFrames
using StatsBase

# Uncomment to make sure Travis CI works as expected
# @test 1 == 2

# OPTIONS
do_plots = false #To create "visual tests". Set to false when using Travis CI.

@everywhere function eye(n::Int)
    Diagonal(ones(n))
end

@testset "MSM.jl" begin

    @testset "testing Types" begin

        @testset "testing MSMOptions" begin

            t = MSMOptions()

            # Testing default values
            #-----------------------------------------------------------------------
            @test t.globalOptimizer == :dxnes
            @test t.localOptimizer == :LBFGS
            @test t.maxFuncEvals == 1000
            @test t.showDistance == false
            @test t.minBox == false
            @test t.populationSize == 50
            @test t.penaltyValue == 999999.0
            @test t.gridType == :LHC
            @test t.saveStartingValues == false
            @test t.maxTrialsStartingValues == 20
            @test t.thresholdStartingValue == t.penaltyValue/10.0

        end

        @testset "testing MSMOptions" begin

            t = MSMProblem()

            # When initialized, iter is equal to 0
            @test typeof(t.priors) == OrderedDict{String,Array{Float64,1}}
            @test typeof(t.empiricalMoments) == OrderedDict{String,Array{Float64,1}}
            @test typeof(t.simulatedMoments) == OrderedDict{String, Float64}
            @test typeof(t.distanceEmpSimMoments) == Float64
            # the functions t.simulate_empirical_moments are initialized with x->x
            @test t.simulate_empirical_moments(1.0) == 1.0
            @test t.objective_function(1.0) == 1.0
            @test typeof(t.options) == MSMOptions

        end

    end #end "testing Types"

    @testset "Latin hypercube sampling" begin

        # Test the home-made function
        a = zeros(3)
        b = ones(3)
        nums = 100
        nums = 10
        points = latin_hypercube_sampling(a, b, nums)

        @test size(points, 1) == nums
        @test size(points, 2) == length(a)

        for i=1:size(points, 1)
            for j=1:size(points, 2)
                @test points[i,j] >= a[j]
                @test points[i,j] <= b[j]
            end
        end

        # Test the function using the package LatinHypercubeSampling.jl
        myProblem = MSMProblem(options = MSMOptions());
        dictPriors = OrderedDict{String,Array{Float64,1}}()
        # Of the form: [initial_guess, lower_bound, upper_bound]
        dictPriors["alpha"] = [0.5, 0.0, 1.0]
        dictPriors["beta1"] = [0.5, 0.0, 1.0]
        dictPriors["beta2"] = [0.5, 0.0, 1.0]
        set_priors!(myProblem, dictPriors)
        points_2 = latin_hypercube_sampling(generate_bbSearchRange(myProblem), nums)

        @test size(points_2, 1) == nums
        @test size(points_2, 2) == length(a)

        for i=1:size(points_2, 1)
            for j=1:size(points_2, 2)
                @test points_2[i,j] >= a[j]
                @test points_2[i,j] <= b[j]
            end
        end

        if do_plots == true
            using Plots
            gr()
            p1 = scatter(points[:,1], points[:,2], points[:,3], title="Home made")
            p2 = scatter(points_2[:,1], points_2[:,2], points_2[:,3], title="LatinHypercubeSampling.jl")
            plot(p1,p2 )
        end

    end

    @testset "Sobol sampling" begin

        # Test the home-made function
        a = zeros(3)
        b = ones(3)
        nums = 100
        points = sobol_sampling(a, b, nums)

        @test size(points, 1) == nums
        @test size(points, 2) == length(a)

        for i=1:size(points, 1)
            for j=1:size(points, 2)
                @test points[i,j] >= a[j]
                @test points[i,j] <= b[j]
            end
        end

        # For fun, let's compare to random sampling
        points_2 = rand(nums, length(a))

        if do_plots == true
            using Plots
            gr()
            p1 = scatter(points_2[:,1], points_2[:,2], points_2[:,3], title="Random")
            p2 = scatter(points[:,1], points[:,2], points[:,3], title="Sobol")
            plot(p1,p2 )
        end

    end

    @testset "testing checks on global algo" begin

        listValidGlobalOptimizers = [:dxnes, :adaptive_de_rand_1_bin_radiuslimited, :xnes,
                         :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin,
                         :generating_set_search, :de_rand_1_bin,
                         :separable_nes, :resampling_inheritance_memetic_search,
                         :probabilistic_descent, :resampling_memetic_search,
                         :de_rand_2_bin_radiuslimited, :de_rand_2_bin,
                         :random_search, :simultaneous_perturbation_stochastic_approximation]

        for globalOptim in listValidGlobalOptimizers

            @test is_global_optimizer(globalOptim) == true

        end

    end

    @testset "testing checks on local algo" begin

        listValidLocalOptimizers = [:NelderMead, :SimulatedAnnealing, :ParticleSwarm,
                            :BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
                            :MomentumGradientDescent, :AcceleratedGradientDescent]

        for localOptim in listValidLocalOptimizers

            @test is_local_optimizer(localOptim) == true

        end

    end

    @testset "testing Optim" begin

        atolOptim = 0.5
        # Algorthims working with FminBox()
        # * GradientDescent
        # * BFGS
        # * LBFGS
        # * ConjugateGradient
        listValidLocalOptimizers = [:GradientDescent, :NelderMead, :SimulatedAnnealing,
                            :BFGS, :LBFGS, :ConjugateGradient,
                            :MomentumGradientDescent, :AcceleratedGradientDescent]

        f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
        x0 = [0.0, 0.0]
        lower = [-2.0; -2.0]
        upper = [2.0; 2.0]

        for localOptim in listValidLocalOptimizers

            results = optimize(f, x0, convert_to_optim_algo(localOptim), Optim.Options(iterations = 2000))

            @test Optim.minimizer(results)[1] ≈ 1.0 atol = atolOptim
            @test Optim.minimizer(results)[2] ≈ 1.0 atol = atolOptim

        end

    end

    @testset "testing FminBox, non binding" begin

        atolOptim = 0.5
        # Algorthims working with FminBox()
        # * GradientDescent
        # * BFGS
        # * LBFGS
        # * ConjugateGradient
        listValidLocalOptimizers = [:GradientDescent, :BFGS, :LBFGS, :ConjugateGradient]

        f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
        x0 = [0.0, 0.0]
        lower = [-2.0; -2.0]
        upper = [2.0; 2.0]

        for localOptim in listValidLocalOptimizers

            results = optimize(f, lower, upper, x0, convert_to_fminbox(localOptim), Optim.Options(iterations = 2000))

            @test Optim.minimizer(results)[1] ≈ 1.0 atol = atolOptim
            @test Optim.minimizer(results)[2] ≈ 1.0 atol = atolOptim

        end

    end

    @testset "testing FminBox binding" begin

        atolOptim = 1e-1

        listValidLocalOptimizers = [:GradientDescent, :BFGS, :LBFGS, :ConjugateGradient]

        f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
        x0 = [0.0, 0.0]
        lower = [-2.0; -2.0]
        upper = [0.5; 0.5]

        for localOptim in listValidLocalOptimizers

            results = optimize(f, lower, upper, x0, convert_to_fminbox(localOptim), Optim.Options(iterations = 2000))

            @test Optim.minimizer(results)[1] < 0.5
            @test Optim.minimizer(results)[1] > -2.0
            @test Optim.minimizer(results)[2] < 0.5
            @test Optim.minimizer(results)[2] > -2.0

        end

    end

    @testset "testing loading priors and empirical moments" begin

       @testset "testing read_priors" begin

            dictPriors = read_priors(joinpath(Pkg.dir("MSM"), "test/priorsTest.csv"))

            @test typeof(dictPriors) == OrderedDict{String,Array{Float64,1}}
            # First component stores the value
            @test dictPriors["alpha"][1] == 0.5
            # Second component stores the lower bound:
            @test dictPriors["alpha"][2] == 0.01
            # Third component stores the upper bound:
            @test dictPriors["alpha"][3] == 0.9

       end

       @testset "testing read_empirical_moments" begin

            dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))

            @test typeof(dictEmpiricalMoments) == OrderedDict{String,Array{Float64,1}}

            # First component stores the value
            @test dictEmpiricalMoments["meanU"][1] == 0.05
            # Second component stores the weight associated to
            @test dictEmpiricalMoments["meanU"][2] == 0.05

       end

    end

    @testset "set_priors!, set_empirical_moments!, set_weight_matrix!" begin

        t = MSMProblem();

        @testset "testing set_priors!" begin

            dictPriors = read_priors(joinpath(Pkg.dir("MSM"), "test/priorsTest.csv"))

            set_priors!(t, dictPriors)

            @test t.priors == dictPriors

        end

        @testset "set_empirical_moments!" begin

            dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))

            set_empirical_moments!(t, dictEmpiricalMoments)

            @test t.empiricalMoments == dictEmpiricalMoments

        end

        @testset "set_weight_matrix!" begin

            N = 5
            W = collect(1:N).*ones(N,N)

            set_weight_matrix!(t, W)

            @test t.W == W

        end
    end

    @testset "testing the construction of the objective function" begin

       @testset "set_simulate_empirical_moments!" begin

            function functionTest(x::Vector)

                output = OrderedDict{String,Float64}()
                output["mom1"] = x[1]
                output["mom2"] = x[2]

                return output
            end

            t = MSMProblem();

            set_simulate_empirical_moments!(t, functionTest)

            x1Value = 1.0
            x2Value = 2.0
            simulatedMoments = t.simulate_empirical_moments([x1Value; x2Value])
            @test simulatedMoments["mom1"]  == x1Value
            @test simulatedMoments["mom2"]  == x2Value

       end

       @testset "Testing construct_objective_function!" begin

            Random.seed!(1234)
            tol1dMean = 0.01

            function functionTest(x::Vector)

                output = OrderedDict{String,Float64}()
                d = Normal(x[1])
                output["meanU"] = mean(rand(d, 100000))

                return output
            end

            t = MSMProblem();

            set_simulate_empirical_moments!(t, functionTest)

            # For the test to make sense, we need to set the field
            # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
            #------------------------------------------------------
            dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))
            set_empirical_moments!(t, dictEmpiricalMoments)

            # A. Set the function: parameter -> simulated moments
            set_simulate_empirical_moments!(t, functionTest)

            # A'. Attach weight marix
            W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
            #Special case: diagonal matrix
            #(you may choose something else)
            for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
            end

            set_weight_matrix!(t, W)

            # B. Construct the objective function, using the function: parameter -> simulated moments
            # and moments' weights:
            construct_objective_function!(t)

            # The objective function should be very close to 0 when
            # evaluated at the true value (modulo randomness)
            @test t.objective_function([dictEmpiricalMoments["meanU"][1]]) ≈ 0. atol = tol1dMean

        end

        @testset "Testing generate_bbSearchRange" begin

            # A.
            #----
            t = MSMProblem();

            dictPriors = read_priors(joinpath(Pkg.dir("MSM"), "test/priorsTest.csv"))

            set_priors!(t, dictPriors)

            testSearchRange = generate_bbSearchRange(t)

            @test testSearchRange[1][1] == 0.01
            @test testSearchRange[1][2] == 0.9
            @test testSearchRange[2][1] == 0.0
            @test testSearchRange[2][2] == 1.0

            # B.
            #---
            t = MSMProblem()

            dictPriors = OrderedDict{String,Array{Float64,1}}()
            dictPriors["mu1"] = [0., -5.0, 5.0]
            dictPriors["mu2"] = [0., -15.0, -10.0]
            dictPriors["mu3"] = [0., -20.0, -15.0]

            set_priors!(t, dictPriors)

            testSearchRange = generate_bbSearchRange(t)

            @test testSearchRange[1][1] == -5.0
            @test testSearchRange[1][2] == 5.0
            @test testSearchRange[2][1] == -15.0
            @test testSearchRange[2][2] == -10.0
            @test testSearchRange[3][1] == -20.0
            @test testSearchRange[3][2] == -15.0

        end

    end

    @testset "Testing msm_optimize!" begin

        # 1d problem
        #-----------
        @testset "Testing msm_optimize! with 1d" begin

            # Rermark:
            # It is important NOT to use Random.seed!()
            # within the function simulate_empirical_moments!
            # Otherwise, BlackBoxOptim does not find the solution
            #----------------------------------------------------
            tol1dMean = 0.2

            @everywhere function functionTest1d(x)

                d = Normal(x[1])
                output = OrderedDict{String,Float64}()

                output["meanU"] = mean(rand(d, 1000000))

                return output
            end

            t = MSMProblem(options = MSMOptions(maxFuncEvals=1000))

            # For the test to make sense, we need to set the field
            # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
            #------------------------------------------------------
            dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))
            set_empirical_moments!(t, dictEmpiricalMoments)

            dictPriors = OrderedDict{String,Array{Float64,1}}()
            dictPriors["mu1"] = [0., -2.0, 2.0]
            set_priors!(t, dictPriors)

            # A. Set the function: parameter -> simulated moments
            #----------------------------------------------------
            set_simulate_empirical_moments!(t, functionTest1d)

            # A'. Attach weight marix
            W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
            #Special case: diagonal matrix
            #(you may choose something else)
            for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
            end

            set_weight_matrix!(t, W)

            # B. Construct the objective function, using the function: parameter -> simulated moments
            # and moments' weights:
            #----------------------------------------------------
            construct_objective_function!(t)

            msm_optimize!(t, verbose = true)

            @test best_candidate(t.bbResults)[1] ≈ 0.05 atol = tol1dMean

            # C. Testing refinement of the global max using a local routine
            #--------------------------------------------------------------
            @test msm_refine_globalmin!(t, verbose = true)[1] ≈ 0.05 atol = tol1dMean

            @test msm_local_minimizer(t)[1] ≈ 0.05 atol = tol1dMean

        end

        @testset "Testing mutlistart 1d" begin

          tol1dMean = 0.1

          @everywhere function functionTest1d(x)

              # When using one of the deterministic methods of Optim,
              # we can safely "control" for randomness
              #-----------------------------------------------------
              Random.seed!(1234)
              d = Normal(x[1])
              output = OrderedDict{String,Float64}()

              output["meanU"] = mean(rand(d, 1000000))

              return output
          end

          t = MSMProblem(options = MSMOptions(maxFuncEvals=1000))

          # For the test to make sense, we need to set the field
          # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
          #------------------------------------------------------
          dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))
          set_empirical_moments!(t, dictEmpiricalMoments)

          dictPriors = OrderedDict{String,Array{Float64,1}}()
          dictPriors["mu1"] = [0., -2.0, 2.0]
          set_priors!(t, dictPriors)

          # A. Set the function: parameter -> simulated moments
          #----------------------------------------------------
          set_simulate_empirical_moments!(t, functionTest1d)

          # A'. Attach weight marix
          W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
          #Special case: diagonal matrix
          #(you may choose something else)
          for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
              W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
          end

          set_weight_matrix!(t, W)

          # B. Construct the objective function, using the function: parameter -> simulated moments
          # and moments' weights:
          #----------------------------------------------------
          construct_objective_function!(t)

          msm_multistart!(t, nums = maxNbWorkers, verbose = true)

          @test msm_local_minimum(t) ≈ 0.0 atol = tol1dMean

          @test msm_local_minimizer(t)[1] ≈ 0.05 atol = tol1dMean

        end

        @testset "Testing local to global 1d with FminBox" begin

          tol1dMean = 0.1

          @everywhere function functionTest1d(x)

              # When using one of the deterministic methods of Optim,
              # we can safely "control" for randomness
              #-----------------------------------------------------
              Random.seed!(1234)
              d = Normal(x[1])
              output = OrderedDict{String,Float64}()

              output["meanU"] = mean(rand(d, 1000000))

              return output
          end

          # Let's try with the optim minBox on
          #-----------------------------------
          t = MSMProblem(options = MSMOptions(maxFuncEvals=1000, localOptimizer = :LBFGS, minBox = true))

          # For the test to make sense, we need to set the field
          # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
          #------------------------------------------------------
          dictEmpiricalMoments = read_empirical_moments(joinpath(Pkg.dir("MSM"), "test/empiricalMomentsTest.csv"))
          set_empirical_moments!(t, dictEmpiricalMoments)

          dictPriors = OrderedDict{String,Array{Float64,1}}()
          dictPriors["mu1"] = [0., -2.0, 2.0]
          set_priors!(t, dictPriors)

          # A. Set the function: parameter -> simulated moments
          #----------------------------------------------------
          set_simulate_empirical_moments!(t, functionTest1d)

          # A'. Attach weight marix
          W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
          #Special case: diagonal matrix
          #(you may choose something else)
          for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
              W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
          end

          set_weight_matrix!(t, W)

          # B. Construct the objective function, using the function: parameter -> simulated moments
          # and moments' weights:
          #----------------------------------------------------
          construct_objective_function!(t)

          msm_multistart!(t, nums = maxNbWorkers, verbose = true)

          @test msm_local_minimum(t) ≈ 0.0 atol = tol1dMean

          @test msm_local_minimizer(t)[1] ≈ 0.05 atol = tol1dMean

        end

        # 2d problem
        #-----------
        @testset "Testing smmoptimize with 2d and same magnitude" begin

            # Rermark:
            # It is important NOT to use Random.seed!()
            # within the function simulate_empirical_moments!
            # Otherwise, BlackBoxOptim does not find the solution
            #----------------------------------------------------
            tol2dMean = 0.2

            @everywhere function functionTest2d(x)

                d = MvNormal( [x[1]; x[2]], eye(2))
                output = OrderedDict{String,Float64}()

                draws = rand(d, 1000000)
                output["mean1"] = mean(draws[1,:])
                output["mean2"] = mean(draws[2,:])

                return output
            end

            t = MSMProblem(options = MSMOptions(maxFuncEvals=1000))

            # For the test to make sense, we need to set the field
            # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
            #------------------------------------------------------
            dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
            dictEmpiricalMoments["mean1"] = [1.0; 1.0]
            dictEmpiricalMoments["mean2"] = [-1.0; -1.0]
            set_empirical_moments!(t, dictEmpiricalMoments)

            dictPriors = OrderedDict{String,Array{Float64,1}}()
            dictPriors["mu1"] = [0., -5.0, 5.0]
            dictPriors["mu2"] = [0., -5.0, 5.0]
            set_priors!(t, dictPriors)

            # A. Set the function: parameter -> simulated moments
            set_simulate_empirical_moments!(t, functionTest2d)

            # A'. Attach weight marix
            W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
            #Special case: diagonal matrix
            #(you may choose something else)
            for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
            end

            set_weight_matrix!(t, W)

            # B. Construct the objective function, using the function: parameter -> simulated moments
            # and moments' weights:
            construct_objective_function!(t)

            # C. Run the optimization
            # This function first modifies t.bbSetup
            # and then modifies t.bbResults
            msm_optimize!(t, verbose = true)

            @test best_candidate(t.bbResults)[1] ≈ 1.0 atol = tol2dMean
            @test best_candidate(t.bbResults)[2] ≈ -1.0 atol = tol2dMean

        end

        # 2d problem
        #-----------
        @testset "Testing msm_multistart! with 2d and same magnitude" begin

            # Rermark:
            # It is important NOT to use Random.seed!()
            # within the function simulate_empirical_moments!
            # Otherwise, BlackBoxOptim does not find the solution
            #----------------------------------------------------
            tol2dMean = 0.2

            @everywhere function functionTest2d(x)

                d = MvNormal( [x[1]; x[2]], eye(2))
                output = OrderedDict{String,Float64}()

                # When using one of the deterministic methods of Optim,
                # we can safely "control" for randomness
                #-----------------------------------------------------
                Random.seed!(1234)
                draws = rand(d, 1000000)
                output["mean1"] = mean(draws[1,:])
                output["mean2"] = mean(draws[2,:])

                return output
            end

            t = MSMProblem(options = MSMOptions(maxFuncEvals=1000))

            # For the test to make sense, we need to set the field
            # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
            #------------------------------------------------------
            dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
            dictEmpiricalMoments["mean1"] = [1.0; 1.0]
            dictEmpiricalMoments["mean2"] = [-1.0; -1.0]
            set_empirical_moments!(t, dictEmpiricalMoments)

            dictPriors = OrderedDict{String,Array{Float64,1}}()
            dictPriors["mu1"] = [0., -5.0, 5.0]
            dictPriors["mu2"] = [0., -5.0, 5.0]
            set_priors!(t, dictPriors)

            # A. Set the function: parameter -> simulated moments
            set_simulate_empirical_moments!(t, functionTest2d)

            # A'. Attach weight marix
            W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
            #Special case: diagonal matrix
            #(you may choose something else)
            for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
            end

            set_weight_matrix!(t, W)

            # B. Construct the objective function, using the function: parameter -> simulated moments
            # and moments' weights:
            construct_objective_function!(t)

            msm_multistart!(t, nums = nworkers(), verbose = true)

            @test msm_local_minimum(t) ≈ 0.0 atol = tol2dMean

            @test msm_local_minimizer(t)[1] ≈ 1.0 atol = tol2dMean
            @test msm_local_minimizer(t)[2] ≈ - 1.0 atol = tol2dMean

          end

          @testset "Testing msm_multistart! with 2d, same magnitude and minBox = true" begin

              # Rermark:
              # It is important NOT to use Random.seed!()
              # within the function simulate_empirical_moments!
              # Otherwise, BlackBoxOptim does not find the solution
              #----------------------------------------------------
              tol2dMean = 0.2

              @everywhere function functionTest2d(x)

                  d = MvNormal( [x[1]; x[2]], eye(2))
                  output = OrderedDict{String,Float64}()

                  # When using one of the deterministic methods of Optim,
                  # we can safely "control" for randomness
                  #-----------------------------------------------------
                  Random.seed!(1234)
                  draws = rand(d, 1000000)
                  output["mean1"] = mean(draws[1,:])
                  output["mean2"] = mean(draws[2,:])

                  return output
              end

              t = MSMProblem(options = MSMOptions(maxFuncEvals=1000, localOptimizer = :GradientDescent, minBox = true))

              # For the test to make sense, we need to set the field
              # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
              #------------------------------------------------------
              dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
              dictEmpiricalMoments["mean1"] = [1.0; 1.0]
              dictEmpiricalMoments["mean2"] = [-1.0; -1.0]
              set_empirical_moments!(t, dictEmpiricalMoments)

              dictPriors = OrderedDict{String,Array{Float64,1}}()
              dictPriors["mu1"] = [0., -5.0, 5.0]
              dictPriors["mu2"] = [0., -5.0, 5.0]
              set_priors!(t, dictPriors)

              # A. Set the function: parameter -> simulated moments
              set_simulate_empirical_moments!(t, functionTest2d)

              # A'. Attach weight marix
              W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
              #Special case: diagonal matrix
              #(you may choose something else)
              for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                  W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
              end

              set_weight_matrix!(t, W)

              # B. Construct the objective function, using the function: parameter -> simulated moments
              # and moments' weights:
              construct_objective_function!(t)

              msm_multistart!(t, nums = nworkers(), verbose = true)

              @test msm_local_minimum(t) ≈ 0.0 atol = tol2dMean

              @test msm_local_minimizer(t)[1] ≈ 1.0 atol = tol2dMean
              @test msm_local_minimizer(t)[2] ≈ - 1.0 atol = tol2dMean

            end

        # 2d problem
        #-----------
        @testset "Testing smmoptimize with 2d with a 1-order magnitude difference" begin

            # Rermark:
            # It is important NOT to use Random.seed!()
            # within the function simulate_empirical_moments!
            # Otherwise, BlackBoxOptim does not find the solution
            #----------------------------------------------------
            # The difference of magniture make it more difficult to find the minimum
            tol2dMean = 2.0

            function functionTest2d(x)

                d = MvNormal( [x[1]; x[2]], eye(2))
                output = OrderedDict{String,Float64}()

                draws = rand(d, 1000000)
                output["mean1"] = mean(draws[1,:])
                output["mean2"] = mean(draws[2,:])

                return output
            end

            t = MSMProblem(options = MSMOptions(maxFuncEvals=1000))

            # For the test to make sense, we need to set the field
            # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
            #------------------------------------------------------
            dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
            dictEmpiricalMoments["mean1"] = [ 1.0; 1.0]
            dictEmpiricalMoments["mean2"] = [-12.0; 12.0]
            set_empirical_moments!(t, dictEmpiricalMoments)

            dictPriors = OrderedDict{String,Array{Float64,1}}()
            dictPriors["mu1"] = [0., -5.0, 5.0]
            dictPriors["mu2"] = [0., -15.0, -10.0]
            set_priors!(t, dictPriors)

            # A. Set the function: parameter -> simulated moments
            set_simulate_empirical_moments!(t, functionTest2d)

            # A'. Attach weight marix
            W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
            #Special case: diagonal matrix
            #(you may choose something else)
            for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
            end

            set_weight_matrix!(t, W)

            # B. Construct the objective function, using the function: parameter -> simulated moments
            # and moments' weights:
            construct_objective_function!(t)

            # C. Run the optimization
            # This function first modifies t.bbSetup
            # and then modifies t.bbResults
            msm_optimize!(t, verbose = true)

            @test best_candidate(t.bbResults)[1] ≈  1.0 atol = tol2dMean
            @test best_candidate(t.bbResults)[2] ≈ -12.0 atol = tol2dMean

        end

    end

    @testset "Testing minimizing a function that may fail" begin

              #---------------------------------------------------
              tol2dMean = 0.5
              @everywhere d_Uni = Uniform(0,1)
              @everywhere uniform_draws = rand(d_Uni, 10000)

              @everywhere function functionTest2d(x, uniform_draws::Array{Float64,1})

                  # function that fails when the first input
                  # is smaller than minus 1:
                  #------------------------------------------
                  if x[1] < -1.0
                    error("I failed")
                  end

                  d = Normal(x[1], 1.0)
                  # Inverse cdf (i.e. quantile)
                  draws = quantile.(d, uniform_draws)
                  output = OrderedDict{String,Float64}()
                  output["mean1"] = mean(draws)

                  return output
              end

              t = MSMProblem(options = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes,
                            localOptimizer = :NelderMead, penaltyValue = 100.0, showDistance=true))

              #---------------------------------------------------------------------
              # Using multistart algo
              #---------------------------------------------------------------------
              # For the test to make sense, we need to set the field
              # t.empiricalMoments::OrderedDict{String,Array{Float64,1}}
              #------------------------------------------------------
              dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
              dictEmpiricalMoments["mean1"] = [1.0; 1.0]
              set_empirical_moments!(t, dictEmpiricalMoments)

              dictPriors = OrderedDict{String,Array{Float64,1}}()
              dictPriors["mu1"] = [0., -5.0, 5.0]
              set_priors!(t, dictPriors)

              # A. Set the function: parameter -> simulated moments
              set_simulate_empirical_moments!(t, x -> functionTest2d(x,uniform_draws))

              # A'. Attach weight marix
              W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
              #Special case: diagonal matrix
              #(you may choose something else)
              for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                  W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
              end

              set_weight_matrix!(t, W)

              # B. Construct the objective function, using the function: parameter -> simulated moments
              # and moments' weights:
              construct_objective_function!(t)

              msm_multistart!(t, nums = nworkers(), verbose = true)
              @test msm_local_minimum(t) ≈ 0.0 atol = tol2dMean
              @test msm_local_minimizer(t)[1] ≈ 1.0 atol = tol2dMean

              #---------------------------------------------------------------------
              # Using BlackBoxOptim
              #---------------------------------------------------------------------
              t = MSMProblem(options = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes))

              dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
              dictEmpiricalMoments["mean1"] = [1.0; 1.0]
              set_empirical_moments!(t, dictEmpiricalMoments)

              dictPriors = OrderedDict{String,Array{Float64,1}}()
              dictPriors["mu1"] = [0., -5.0, 5.0]
              set_priors!(t, dictPriors)

              # A. Set the function: parameter -> simulated moments
              set_simulate_empirical_moments!(t, x -> functionTest2d(x, uniform_draws))

              # A'. Attach weight marix
              W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
              #Special case: diagonal matrix
              #(you may choose something else)
              for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
                  W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
              end

              set_weight_matrix!(t, W)

              # B. Construct the objective function, using the function: parameter -> simulated moments
              # and moments' weights:
              construct_objective_function!(t)

              msm_optimize!(t, verbose = true)

              @test best_candidate(t.bbResults)[1] ≈  1.0 atol = tol2dMean

    end

    @testset "Testing Var-Covariance estimation" begin

        d = Normal(0, 0.2)
        T=10000
        y1 = rand(d, T);
        y2 = rand(d, T);
        y3 = rand(d, T);
        data = hcat(y1, y2, y3);

        #When l=0, cov_NW(data) should be cov(data)
        @test maximum(abs.(cov_NW(data, l=0) .- cov(data))) < 1.0e-10

        # Even when including lags, with no serial correlation the difference should
        # be small
        @test maximum(abs.(cov_NW(data, l=1) .- cov(data))) < 1.0e-2

    end

    @testset "Testing Inference" begin

      tolLinear = 0.05

      # Inference in the linear model
      #------------------------------
      Random.seed!(1234)         #for replicability reasons
      T = 100000          #number of periods
      P = 2               #number of dependent variables
      beta0 = [2.0; 3.0]     #choose true coefficients by drawing from a uniform distribution on [0,1]
      alpha0 = 1.0  #intercept
      theta0 = 0.0        #coefficient to create serial correlation in the error terms
      println("True intercept = $(alpha0)")
      println("True coefficient beta0 = $(beta0)")
      println("Serial correlation coefficient theta0 = $(theta0)")

      # Simulation of a sample:
      # Generation of error terms
      #--------------------------
      # row = individual dimension
      # column = time dimension
      U = zeros(T)
      d = Normal()
      U[1] = rand(d, 1)[] #first error term
      # loop over time periods
      for t = 2:T
          U[t] = rand(d, 1)[] + theta0*U[t-1]
      end
      # Let's simulate x_t
      #-------------------
      x = zeros(T, P)

      d = Uniform(0, 5)
      for p = 1:P
              x[:,p] = rand(d, T)
      end

      # Let's calculate the resulting y_t
      y = zeros(T)

      for t=1:T
          y[t] = alpha0 + x[t,1]*beta0[1] + x[t,2]*beta0[2] + U[t]
      end

      myProblem = MSMProblem(options = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes, localOptimizer = :LBFGS));

      # Empirical moments
      #------------------
      dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
      dictEmpiricalMoments["mean"] = [mean(y); mean(y)] #informative on the intercept
      dictEmpiricalMoments["mean_x1y"] = [mean(x[:,1] .* y); mean(x[:,1] .* y)] #informative on betas
      dictEmpiricalMoments["mean_x2y"] = [mean(x[:,2] .* y); mean(x[:,2] .* y)] #informative on betas
      dictEmpiricalMoments["mean_x1y^2"] = [mean((x[:,1] .* y).^2); mean((x[:,1] .* y).^2)] #informative on betas
      dictEmpiricalMoments["mean_x2y^2"] = [mean((x[:,2] .* y).^2); mean((x[:,2] .* y).^2)] #informative on betas

      set_empirical_moments!(myProblem, dictEmpiricalMoments)

      dictPriors = OrderedDict{String,Array{Float64,1}}()
      dictPriors["alpha"] = [0.5, 0.001, 1.0]
      dictPriors["beta1"] = [0.5, 0.001, 1.0]
      dictPriors["beta2"] = [0.5, 0.001, 1.0]

      set_priors!(myProblem, dictPriors)

      # A'. Attach weight marix
      W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
      #Special case: diagonal matrix
      #(you may choose something else)
      for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
          W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
      end

      set_weight_matrix!(myProblem, W)

      # x[1] corresponds to the intercept
      # x[1] corresponds to beta1
      # x[3] corresponds to beta2
      @everywhere function functionLinearModel(x; uniform_draws::Array{Float64,1}, simX::Array{Float64,2}, nbDraws::Int64 = length(uniform_draws), burnInPerc::Int64 = 10)

          T = nbDraws
          P = 2       #number of dependent variables

          alpha = x[1]
          beta = x[2:end]
          theta = 0.0     #coefficient to create serial correlation in the error terms

          # Creation of error terms
          # row = individual dimension
          # column = time dimension
          U = zeros(T)
          d = Normal()
          # Inverse cdf (i.e. quantile)
          gaussian_draws = quantile.(d, uniform_draws)
          U[1] = gaussian_draws[1] #first error term

          # loop over time periods
          for t = 2:T
              U[t] = gaussian_draws[t] + theta*U[t-1]
          end

          # Let's calculate the resulting y_t
          y = zeros(T)

          for t=1:T
              y[t] = alpha + simX[t,1]*beta[1] + simX[t,2]*beta[2] + U[t]
          end

          # Get rid of the burn-in phase:
          #------------------------------
          startT = div(nbDraws, burnInPerc)

          # Moments:
          #---------
          output = OrderedDict{String,Float64}()
          output["mean"] = mean(y[startT:nbDraws])
          output["mean_x1y"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws])
          output["mean_x2y"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws])
          output["mean_x1y^2"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2)
          output["mean_x2y^2"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)

          return output
      end

      # Let's freeze the randomness during the minimization
      @everywhere d_Uni = Uniform(0,1)
      @everywhere nbDraws = 100000 #number of draws in the simulated data
      @everywhere uniform_draws = rand(d_Uni, nbDraws)
      simX = zeros(length(uniform_draws), 2)
      d = Uniform(0, 5)
      for p = 1:2
          simX[:,p] = rand(d, length(uniform_draws))
      end

      set_simulate_empirical_moments!(myProblem, x -> functionLinearModel(x, uniform_draws = uniform_draws, simX = simX))

      # Construct the objective function using:
      #* the function: parameter -> simulated moments
      #* emprical moments values
      #* emprical moments weights
      construct_objective_function!(myProblem)

      # Run the optimization in parallel using n different starting values
      # where n is equal to the number of available workers
      #--------------------------------------------------------------------
      listOptimResults = msm_multistart!(myProblem, verbose = true)

      # Remark: it would not be appropriate to use BlackBoxOptim because the
      # No big deal here, because we use Optim
      minimizer = msm_local_minimizer(myProblem)

      # The minimizer should not be too far from the true values
      #---------------------------------------------------------
      @test minimizer[1] ≈  alpha0[1] atol = tolLinear
      @test minimizer[2] ≈  beta0[1] atol = tolLinear
      @test minimizer[3] ≈  beta0[2] atol = tolLinear

      # Empirical Distance matrix
      #--------------------------
      X = zeros(T, 5)

      X[:,1] = y
      X[:,2] = (x[:,1] .* y)
      X[:,3] = (x[:,2] .* y)
      X[:,4] = (x[:,1] .* y).^2
      X[:,5] = (x[:,2] .* y).^2

      Sigma0 = cov(X)

      set_Sigma0!(myProblem, Sigma0)

      @test myProblem.Sigma0 == Sigma0

      calculate_Avar!(myProblem, minimizer, tau = T/nbDraws)

      # The asymptotic variance should be
      # * symmetric
      # * positive semi-definite
      #-----------------------------------
      @test issymmetric(myProblem.Avar) == true

      # test for positive semi-definiteness
      # all eigenvalues should be non-negative
      eigv = eigvals(myProblem.Avar)

      counterNegativeEigVals = 0
      for i=1:length(eigv)
        if eigv[i] < 0
          counterNegativeEigVals += 1
        end
      end

      @test counterNegativeEigVals == 0

      # summary table:
      #---------------
      df = summary_table(myProblem, minimizer, T, 0.05)
      @test typeof(df) == CoefTable

      # Convert coeftable to a DataFrame
      df = DataFrame(df)

      # first column : point estimates
      @test df[:, "Coef."] == minimizer

      # 2nd column : std error
      for i =1:size(df,1)
        @test df[i, "Std. Error"] > 0.
      end

      # confidence interval
      for i =1:size(df,1)
        @test df[i, "CI Lower"] <= df[i, "CI Upper"]
      end

      D = calculate_D(myProblem, minimizer)
      @test rank(D) == size(df,1)

      # Slice the objective functions
      nb_p = 7
      vXGrid, vYGrid = msm_slices(myProblem, minimizer, nbPoints = nb_p)

      @test size(vXGrid,1) == nb_p
      @test size(vXGrid,2) == size(df,1)

      @test size(vYGrid,1) == nb_p
      @test size(vYGrid,2) == size(df,1)

      # Plot the results
      if do_plots == true

          gr()
          list_plots = []
          for (keyIndex, keyValue) in enumerate(keys(myProblem.priors))

              p = plot(vXGrid[:,keyIndex], vYGrid[:,keyIndex], title = "$(keyValue)", label = "")
              plot!([minimizer[keyIndex]], seriestype = :vline, label = "")
              push!(list_plots, p)

          end

          #Let's combine all the plots in a single plot
          s0 = ""
          for i = 1:length(keys(myProblem.priors))
              if i==1
                  s0 = string("list_plots[$(i)]" )
              else
                  s0 = string(s0, ", ", "list_plots[$(i)]" )
              end
          end

          plot_combined = eval(Meta.parse(string("plot(", s0, ")")))
          display(plot_combined)

      end

    end

end


text:
import Base: read, write, close
import MQTT: read_len, Message
using Base.Test, MQTT

#include("smoke.jl")
include("mocksocket.jl")
include("packet.jl")
include("unittests.jl")

ext: jl

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
#=

lucas_stokey.jl

@author: Shunsuke Hori

=#

module LS

using QuantEcon

using NLsolve
using NLopt

mutable struct Para{TF <: AbstractFloat,
                    TM <: AbstractMatrix{TF},
                    TV <: AbstractVector{TF}}
    beta::TF
    Pi::TM
    G::TV
    Theta::TV
    transfers::Bool
    U::Function
    Uc::Function
    Ucc::Function
    Un::Function
    Unn::Function
    n_less_than_one::Bool
end

"""
Class returns planner's allocation as a function of the multiplier
on the implementability constraint mu
"""
struct Planners_Allocation_Sequential{TP <: Para,
                                      TI <: Integer,
                                      TV <: AbstractVector}
    para::TP
    mc::MarkovChain
    S::TI
    cFB::TV
    nFB::TV
    XiFB::TV
    zFB::TV
end

"""
Initializes the class from the calibration Para
"""
function Planners_Allocation_Sequential(para::Para)
    beta, Pi, G, Theta =
        para.beta, para.Pi, para.G, para.Theta
    mc = MarkovChain(Pi)
    S = size(Pi, 1) # number of states
    #now find the first best allocation
    cFB, nFB, XiFB, zFB = find_first_best(para, S, 1)

    return Planners_Allocation_Sequential(para, mc, S, cFB, nFB, XiFB, zFB)
end

"""
Find the first best allocation
"""
function find_first_best(para::Para, S::Integer, version::Integer)
    if version != 1 && version != 2
        throw(ArgumentError("version must be 1 or 2"))
    end
    beta, Theta, Uc, Un, G, Pi =
        para.beta, para.Theta, para.Uc, para.Un, para.G, para.Pi
    function res!(z, out)
        c = z[1:S]
        n = z[S+1:end]
        out[1:S] = Theta.*Uc(c, n)+Un(c, n)
        out[S+1:end] = Theta.*n - c - G
    end
    res = nlsolve(res!, 0.5*ones(2*S))

    if converged(res) == false
        error("Could not find first best")
    end

    if version == 1
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        XiFB = Uc(cFB, nFB) #multiplier on the resource constraint.
        zFB = vcat(cFB, nFB, XiFB)
        return cFB, nFB, XiFB, zFB
    elseif version == 2
        cFB = res.zero[1:S]
        nFB = res.zero[S+1:end]
        IFB = Uc(cFB, nFB).*cFB + Un(cFB, nFB).*nFB
        xFB = \(eye(S) - beta*Pi, IFB)
        zFB = [vcat(cFB[s], xFB[s], xFB) for s in 1:S]
        return cFB, nFB, IFB, xFB, zFB
    end
end

"""
Computes optimal allocation for time ``t\geq 1`` for a given ``\mu``
"""
function time1_allocation(pas::Planners_Allocation_Sequential, mu::Real)
    para, S = pas.para, pas.S
    Theta, beta, Pi, G, Uc, Ucc, Un, Unn =
        para.Theta, para.beta, para.Pi, para.G,
        para.Uc, para.Ucc, para.Un, para.Unn
    function FOC!(z::Vector, out)
        c = z[1:S]
        n = z[S+1:2S]
        Xi = z[2S+1:end]
        out[1:S] = Uc(c,n) - mu*(Ucc(c,n).*c+Uc(c,n)) -Xi #foc c
        out[S+1:2S] = Un(c,n) - mu*(Unn(c,n).*n+Un(c,n)) + Theta.*Xi #foc n
        out[2S+1:end] = Theta.*n - c - G #resource constraint
        return out
    end
    #find the root of the FOC
    res = nlsolve(FOC!, pas.zFB)
    if res.f_converged == false
        error("Could not find LS allocation.")
    end
    z = res.zero
    c, n, Xi = z[1:S], z[S+1:2S], z[2S+1:end]
    #now compute x
    I  = Uc(c,n).*c +  Un(c,n).*n
    x = \(eye(S) - beta*para.Pi, I)
    return c, n, x, Xi
end

"""
Finds the optimal allocation given initial government debt `B_` and state `s_0`
"""
function time0_allocation(pas::Planners_Allocation_Sequential,
                          B_::AbstractFloat, s_0::Integer)
    para = pas.para
    Pi, Theta, G, beta =
        para.Pi, para.Theta, para.G, para.beta
    Uc, Ucc, Un, Unn =
        para.Uc, para.Ucc, para.Un, para.Unn
    #first order conditions of planner's problem
    function FOC!(z, out)
        mu, c, n, Xi = z[1], z[2], z[3], z[4]
        xprime = time1_allocation(pas, mu)[3]
        out .= vcat(
        Uc(c, n).*(c-B_) + Un(c, n).*n + beta*dot(Pi[s_0, :], xprime),
        Uc(c, n) - mu*(Ucc(c, n).*(c-B_) + Uc(c, n)) - Xi,
        Un(c, n) - mu*(Unn(c, n).*n+Un(c, n)) + Theta[s_0].*Xi,
        (Theta.*n - c - G)[s_0]
        )
    end
    #find root
    res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.XiFB[s_0]])
    if res.f_converged == false
        error("Could not find time 0 LS allocation.")
    end
    return (res.zero...)
end

"""
Find the value associated with multiplier `mu`
"""
function time1_value(pas::Planners_Allocation_Sequential, mu::Real)
    para = pas.para
    c, n, x, Xi = time1_allocation(pas, mu)
    U_val = para.U.(c, n)
    V = \(eye(pas.S) - para.beta*para.Pi, U_val)
    return c, n, x, V
end

"""
Computes Tau given `c`, `n`
"""
function Tau(para::Para, c::Union{Real,Vector}, n::Union{Real,Vector})
    Uc, Un = para.Uc.(c, n), para.Un.(c, n)
    return 1+Un./(para.Theta .* Uc)
end

"""
Simulates planners policies for `T` periods
"""
function simulate(pas::Planners_Allocation_Sequential,
                  B_::AbstractFloat, s_0::Integer, T::Integer,
                  sHist::Union{Vector, Void}=nothing)
    para = pas.para
    Pi, beta, Uc = para.Pi, para.beta, para.Uc
    if sHist == nothing
        sHist = QuantEcon.simulate(pas.mc, T, init=s_0)
    end
    cHist = zeros(T)
    nHist = zeros(T)
    Bhist = zeros(T)
    TauHist = zeros(T)
    muHist = zeros(T)
    RHist = zeros(T-1)
    #time0
    mu, cHist[1], nHist[1], _  = time0_allocation(pas, B_, s_0)
    TauHist[1] = Tau(pas.para, cHist[1], nHist[1])[s_0]
    Bhist[1] = B_
    muHist[1] = mu
    #time 1 onward
    for t in 2:T
        c, n, x, Xi = time1_allocation(pas,mu)
        u_c = Uc(c,n)
        s = sHist[t]
        TauHist[t] = Tau(pas.para, c, n)[s]
        Eu_c = dot(Pi[sHist[t-1],:], u_c)
        cHist[t], nHist[t], Bhist[t] = c[s], n[s], x[s]/u_c[s]
        RHist[t-1] = Uc(cHist[t-1], nHist[t-1])/(beta*Eu_c)
        muHist[t] = mu
    end
    return cHist, nHist, Bhist, TauHist, sHist, muHist, RHist
end

"""
Bellman equation for the continuation of the Lucas-Stokey Problem
"""
mutable struct BellmanEquation{TP <: Para,
                               TI <: Integer,
                               TV <: AbstractVector,
                               TM <: AbstractMatrix{TV},
                               TVV <: AbstractVector{TV}}
    para::TP
    S::TI
    xbar::TV
    time_0::Bool
    z0::TM
    cFB::TV
    nFB::TV
    xFB::TV
    zFB::TVV
end

"""
Initializes the class from the calibration `Para`
"""
function BellmanEquation(para::Para, xgrid::AbstractVector, policies0::Vector)
    S = size(para.Pi, 1) # number of states
    xbar = [minimum(xgrid), maximum(xgrid)]
    time_0 = false
    cf, nf, xprimef = policies0
    z0 = [vcat(cf[s](x), nf[s](x), [xprimef[s, sprime](x) for sprime in 1:S])
                        for x in xgrid, s in 1:S]
    cFB, nFB, IFB, xFB, zFB = find_first_best(para, S, 2)
    return BellmanEquation(para, S, xbar, time_0, z0, cFB, nFB, xFB, zFB)
end

"""
Finds the optimal policies
"""
function get_policies_time1(T::BellmanEquation,
                        i_x::Integer, x::AbstractFloat,
                        s::Integer, Vf::AbstractArray)
    para, S = T.para, T.S
    beta, Theta, G, Pi = para.beta, para.Theta, para.G, para.Pi
    U, Uc, Un = para.U, para.Uc, para.Un

    function objf(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + beta * dot(Pi[s, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n=c+G[s]
        return x - Uc(c, n)*c-Un(c, n)*n - beta*dot(Pi[s, :], xprime)
    end
    lb = vcat(0, T.xbar[1]*ones(S))
    ub = vcat(1-G[s], T.xbar[2]*ones(S))
    opt = Opt(:LN_COBYLA, length(T.z0[i_x, s])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.z0[i_x, s][1], T.z0[i_x, s][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    T.z0[i_x, s] = vcat(minx[1], minx[1]+G[s], minx[2:end])
    return vcat(-minf, T.z0[i_x, s])
end
"""
Finds the optimal policies
"""
function get_policies_time0(T::BellmanEquation,
                        B_::AbstractFloat, s0::Integer, Vf::Array)
    para, S = T.para, T.S
    beta, Theta, G, Pi = para.beta, para.Theta, para.G, para.Pi
    U, Uc, Un = para.U, para.Uc, para.Un
    function objf(z, grad)
        c, xprime = z[1], z[2:end]
        n = c+G[s0]
        Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S]
        return -(U(c, n) + beta*dot(Pi[s0, :], Vprime))
    end
    function cons(z::Vector, grad)
        c, xprime = z[1], z[2:end]
        n = c+G[s0]
        return -Uc(c, n)*(c-B_)-Un(c, n)*n - beta*dot(Pi[s0, :], xprime)
    end
    lb = vcat(0, T.xbar[1]*ones(S))
    ub = vcat(1-G[s0], T.xbar[2]*ones(S))
    opt = Opt(:LN_COBYLA, length(T.zFB[s0])-1)
    min_objective!(opt, objf)
    equality_constraint!(opt, cons)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 300)
    maxtime!(opt, 10)
    init = vcat(T.zFB[s0][1], T.zFB[s0][3:end])
    for (i, val) in enumerate(init)
        if val > ub[i]
            init[i] = ub[i]
        elseif val < lb[i]
            init[i] = lb[i]
        end
    end
    (minf, minx, ret) = optimize(opt, init)
    return vcat(-minf, vcat(minx[1], minx[1]+G[s0], minx[2:end]))
end

"""
Compute the planner's allocation by solving Bellman
equation.
"""
struct Planners_Allocation_Bellman{TP <: Para, TI <: Integer,
                                   TVg <: AbstractVector, TVv <: AbstractVector,
                                   TVp <: AbstractArray}
    para::TP
    mc::MarkovChain
    S::TI
    T::BellmanEquation
    mugrid::TVg
    xgrid::TVg
    Vf::TVv
    policies::TVp
end

"""
Initializes the class from the calibration `Para`
"""
function Planners_Allocation_Bellman(para::Para, mugrid::AbstractArray)
    mc = MarkovChain(para.Pi)
    G = para.G
    S = size(para.Pi, 1) # number of states
    #now find the first best allocation
    Vf, policies, T, xgrid = solve_time1_bellman(para, mugrid)
    T.time_0 = true #Bellman equation now solves time 0 problem
    return Planners_Allocation_Bellman(para, mc, S, T, mugrid, xgrid, Vf, policies)
end

"""
Solve the time 1 Bellman equation for calibration `Para` and initial grid `mugrid0`
"""
function solve_time1_bellman{TF <: AbstractFloat}(para::Para{TF}, mugrid::AbstractArray)
    mugrid0 = mugrid
    S = size(para.Pi, 1)
    #First get initial fit
    PP = Planners_Allocation_Sequential(para)
    c = Matrix{TF}(length(mugrid), 2)
    n = Matrix{TF}(length(mugrid), 2)
    x = Matrix{TF}(length(mugrid), 2)
    V = Matrix{TF}(length(mugrid), 2)
    for (i, mu) in enumerate(mugrid0)
        c[i, :], n[i, :], x[i, :], V[i, :] = time1_value(PP, mu)
    end
    Vf = Vector{LinInterp}(2)
    cf = Vector{LinInterp}(2)
    nf = Vector{LinInterp}(2)
    xprimef = Array{LinInterp}(2, S)
    for s in 1:2
        cf[s] = LinInterp(x[:, s][end:-1:1], c[:, s][end:-1:1])
        nf[s] = LinInterp(x[:, s][end:-1:1], n[:, s][end:-1:1])
        Vf[s] = LinInterp(x[:, s][end:-1:1], V[:, s][end:-1:1])
        for sprime in 1:S
            xprimef[s, sprime] = LinInterp(x[:, s][end:-1:1], x[:, s][end:-1:1])
        end
    end
    policies = [cf, nf, xprimef]
    #create xgrid
    xbar = [maximum(minimum(x, 1)), minimum(maximum(x, 1))]
    xgrid = linspace(xbar[1], xbar[2], length(mugrid0))
    #Now iterate on bellman equation
    T = BellmanEquation(para, xgrid, policies)
    diff = 1.0
    while diff > 1e-6
        if T.time_0 == false
            Vfnew, policies =
                fit_policy_function(PP,
                (i_x, x, s) -> get_policies_time1(T, i_x, x, s, Vf), xgrid)
        elseif T.time_0 == true
            Vfnew, policies =
                fit_policy_function(PP,
                (i_x, B_, s0) -> get_policies_time0(T, i_x, B_, s0, Vf), xgrid)
        else
            error("T.time_0 is $(T.time_0), which is invalid")
        end
        diff = 0.0
        for s in 1:S
            diff = max(diff, maximum(abs, (Vf[s].(xgrid)-Vfnew[s].(xgrid))./Vf[s].(xgrid)))
        end
        print("diff = $diff \n")
        Vf = Vfnew
    end
    # store value function policies and Bellman Equations
    return Vf, policies, T, xgrid
end

"""
Fits the policy functions PF using the points `xgrid` using interpolation
"""
function fit_policy_function(PP::Planners_Allocation_Sequential,
                            PF::Function, xgrid::AbstractArray)
    S = PP.S
    Vf = Vector{LinInterp}(S)
    cf = Vector{LinInterp}(S)
    nf = Vector{LinInterp}(S)
    xprimef = Array{LinInterp}(S, S)
    for s in 1:S
        PFvec = Array{typeof(PP.para).parameters[1]}(length(xgrid), 3+S)
        for (i_x, x) in enumerate(xgrid)
            PFvec[i_x, :] = PF(i_x, x, s)
        end
        Vf[s] = LinInterp(xgrid, PFvec[:, 1])
        cf[s] = LinInterp(xgrid, PFvec[:, 2])
        nf[s] = LinInterp(xgrid, PFvec[:, 3])
        for sprime in 1:S
            xprimef[s, sprime] = LinInterp(xgrid, PFvec[:, 3+sprime])
        end
    end
    return Vf, [cf, nf, xprimef]
end

"""
Finds the optimal allocation given initial government debt `B_` and state `s_0`
"""
function time0_allocation(pab::Planners_Allocation_Bellman,
                          B_::AbstractFloat, s0::Integer)
    xgrid = pab.xgrid
    if pab.T.time_0 == false
        z0 = get_policies_time1(pab.T, i_x, x, s, pab.Vf)
    elseif pab.T.time_0 == true
        z0 = get_policies_time0(pab.T, B_, s0, pab.Vf)
    else
        error("T.time_0 is $(T.time_0), which is invalid")
    end
    c0, n0, xprime0 = z0[2], z0[3], z0[4:end]
    return c0, n0, xprime0
end

"""
Simulates Ramsey plan for `T` periods
"""
function simulate(pab::Planners_Allocation_Bellman,
                B_::AbstractFloat, s_0::Integer, T::Integer,
                sHist::Vector=QuantEcon.simulate(mc, s_0, T))
    para, S, policies = pab.para, pab.S, pab.policies
    beta, Pi, Uc = para.beta, para.Pi, para.Uc
    cf, nf, xprimef = policies[1], policies[2], policies[3]
    TF = typeof(para).parameters[1]
    cHist = Vector{TF}(T)
    nHist = Vector{TF}(T)
    Bhist = Vector{TF}(T)
    TauHist = Vector{TF}(T)
    muHist = Vector{TF}(T)
    RHist = Vector{TF}(T-1)
    #time0
    cHist[1], nHist[1], xprime = time0_allocation(pab, B_, s_0)
    TauHist[1] = Tau(pab.para, cHist[1], nHist[1])[s_0]
    Bhist[1] = B_
    muHist[1] = 0.0
    #time 1 onward
    for t in 2:T
        s, x = sHist[t], xprime[sHist[t]]
        n = nf[s](x)
        c = [cf[shat](x) for shat in 1:S]
        xprime = [xprimef[s, sprime](x) for sprime in 1:S]
        TauHist[t] = Tau(pab.para, c, n)[s]
        u_c = Uc(c, n)
        Eu_c = dot(Pi[sHist[t-1], :], u_c)
        muHist[t] = pab.Vf[s](x)
        RHist[t-1] = Uc(cHist[t-1], nHist[t-1])/(beta*Eu_c)
        cHist[t], nHist[t], Bhist[t] = c[s], n, x/u_c[s]
    end
    return cHist, nHist, Bhist, TauHist, sHist, muHist, RHist
end
end


text:
@testset "Backend Unit Tests:" begin

    @testset "Constructor" begin
        b = BlockOptBackend(test_model, test_driver)

        @testset "inhereted model methods" begin

            @test model(b) ≡ test_model

            @test b.xₖ ≈ initial_iterate(test_model)

            @test b.n ≡ dimension(test_model)

            @test obj(b, b.xₖ) ≈ obj(test_model, initial_iterate(test_model))

            @test grad!(b, similar(b.xₖ), b.xₖ) ≈ grad!(test_model, similar(b.xₖ), b.xₖ)

            @test grad(b, b.xₖ) ≈ grad(test_model, b.xₖ)

            @test dimension(b) ≡ dimension(test_model)
        end

        @testset "inhereted driver methods" begin

            @test driver(b) ≡ test_driver

            @test b.w ≡ Int(samples(test_driver) / 2)

            @test QN_update(b) ≡ QN_update(test_driver)

            @test S_update(b) ≡ S_update(test_driver)

            @test pflag(b) ≡ pflag(test_driver)

            @test samples(b) ≡ samples(test_driver)

            @test Δ_max(b) ≡ Δ_max(test_driver)

            @test δ_tol(b) ≡ δ_tol(test_driver)

            @test ϵ_tol(b) ≡ ϵ_tol(test_driver)

            @test max_iterations(b) ≡ max_iterations(test_driver)
        end

        @testset "backend accessors" begin

            @test fₖ(b) ≡ ∇fₖ_norm(b) ≡ pₖ_norm(b) ≡ Δₖ(b) ≡ ρ(b) ≡ 0.0

            b.fₖ = 1.0

            @test fₖ(b) ≡ 1.0
            b.fₖ = 0.0

            b.∇fₖ_norm = 1.0

            @test ∇fₖ_norm(b) ≡ 1.0
            b.∇fₖ_norm = 0.0

            b.pₖ_norm = 1.0

            @test pₖ_norm(b) ≡ 1.0
            b.pₖ_norm = 0.0

            b.Δₖ = 1.0

            @test Δₖ(b) ≡ 1.0
            b.Δₖ = 0.0

            b.ρ = 1.0

            @test ρ(b) ≡ 1.0
            b.ρ = 0.0
        end

        @testset "initialization of fields" begin

            @test size(b.C) ≡ size(b.P) ≡ (samples(b) + 1, samples(b) + 1)

            @test size(b.b) ≡ size(b.aₖ) ≡ (samples(b) + 1,)

            @test size(b.Qₖ) ≡ (dimension(b), samples(b) + 1)

            @test size(b.xₖ) ≡ size(b.∇fₖ) ≡ size(b.hₖ) ≡ (dimension(b),)

            @test size(b.pₖ) ≡ size(b.qₖ) ≡ size(b.xₜ) ≡ (dimension(b),)

            @test size(b.Sₖ) ≡ size(b.Yₖ) ≡ (dimension(b), samples(b) - 1)

            @test size(b.Uₖ) ≡ size(b.Vₖ) ≡ (dimension(b), samples(b))

            @test size(b.Hₖ) ≡ (dimension(b), dimension(b))

            @test map(issymmetric, (b.C, b.P, b.Hₖ)) ≡ (true, true, true)
        end
    end

    # @testset "Memory manipulation" begin
    #     b = BlockOptBackend(test_model, test_driver)

    #     # manipulates fₖ, Δₖ, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, ∇fₖ, hₖ, ∇fₖ_norm
    #     @testset "initialize subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         initialize(b)

    #         @test !(fₖ, Δₖ, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, ∇fₖ, hₖ, ∇fₖ_norm ≈
    #         b.fₖ, b.Δₖ, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ, b.Hₖ, b.∇fₖ, b.hₖ, b.∇fₖ_norm)

    #         @test xₖ, pₖ, xₜ, fₜ, ρ, pₖ_norm, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.pₖ_norm, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # manipulates P, b, C
    #     @testset "build_trs subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         build_trs(b)

    #         @test !(P, b, C ≈ P, b, C)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.aₖ, b.qₖ
    #     end

    #     # manipulates aₖ, qₖ, pₖ, pₖ_norm
    #     @testset "solve_trs subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         solve_trs(b)

    #         @test !(aₖ, qₖ, pₖ, pₖ_norm ≈ b.aₖ, b.qₖ, b.pₖ, b.pₖ_norm)

    #         @test xₖ, ∇fₖ, hₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C
    #     end

    #     # manipulates xₜ, fₜ, ρ
    #     @testset "build_trial subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         build_trial(b)

    #         @test !(xₜ, fₜ, ρ ≈ b.xₜ, b.fₜ, b.ρ)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, fₖ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.fₖ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # manipulates Sₖ
    #     @testset "update_Sₖ subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         update_Sₖ(b)

    #         @test !(Sₖ ≈ b.Sₖ)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # manipulates ∇fₖ_norm, ∇fₖ, hₖ, Yₖ Uₖ, Vₖ
    #     @testset "gHS subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         gHS(b)

    #         @test !(∇fₖ_norm, ∇fₖ, hₖ, Yₖ, Uₖ, Vₖ ≈ b.∇fₖ_norm, b.∇fₖ, b.hₖ, b.Yₖ, b.Uₖ, b.Vₖ)

    #         @test xₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, pₖ_norm, Sₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.pₖ_norm, b.Sₖ, b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # manipulates Hₖ
    #     @testset "blockQN subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         blockQN(b)

    #         @test !(Hₖ ≈ b.Hₖ)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # manipulates Hₖ
    #     @testset "secantQN subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         secantQN(b)

    #         @test !(Hₖ ≈ b.Hₖ)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # may manipulate xₜ, fₜ
    #     @testset "accept_trial subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         accept_trial(b)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, fₖ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.fₖ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # may manipulate Δₖ
    #     @testset "update_Δₖ subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         update_Δₖ(b)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # no memory manipulation
    #     @testset "terminal subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         terminal(b, 0)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end

    #     # no memory manipulation
    #     @testset "pflag subroutine" begin
    #         @unpack xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ = b

    #         update_Δₖ(b)

    #         @test xₖ, ∇fₖ, hₖ, pₖ, xₜ, fₖ, fₜ, ρ, Δₖ, ∇fₖ_norm, pₖ_norm, Sₖ, Yₖ, Uₖ, Vₖ, Hₖ, Qₖ, P, b, C, aₖ, qₖ ≈
    #         b.xₖ, b.∇fₖ, b.hₖ, b.pₖ, b.xₜ, b.fₖ, b.fₜ, b.ρ, b.Δₖ, b.∇fₖ_norm, b.pₖ_norm, b.Sₖ, b.Yₖ, b.Uₖ, b.Vₖ,
    #         b.Hₖ, b.Qₖ, b.P, b.b, b.C, b.aₖ, b.qₖ
    #     end
    # end

    @testset "Edge cases" begin

        @testset "terminal subroutine" begin

            b = BlockOptBackend(test_model, test_driver)

            b.∇fₖ_norm = ϵ_tol(b) + 1.0

            @test terminal(b, max_iterations(b) - 1) ≡ false # terminal iteration

            @test terminal(b, max_iterations(b)) ≡ true      # exit condition

            b.∇fₖ_norm = ϵ_tol(b)

            @test terminal(b, max_iterations(b) - 1) ≡ false  # not good enough

            b.∇fₖ_norm = ϵ_tol(b) - eps(Float64)

            @test terminal(b, max_iterations(b) - 1) ≡ true  # convergence reached
        end

        @testset "terminal subroutine" begin

            b = BlockOptBackend(test_model, test_driver)

            # TODO: Return b and pipeline above operations.
            initialize(b)
            build_trs(b)
            solve_trs(b)
            build_trial(b)

            b.ρ = 0.0

            @test !accept_trial(b) && !(b.xₖ ≈ b.xₜ) && !(b.fₖ ≈ b.fₜ)

            b.ρ = 1.0

            @test accept_trial(b) && (b.xₖ ≈ b.xₜ) && (b.fₖ ≈ b.fₜ)
        end

        @testset "update_Δₖ subroutine" begin

            b = BlockOptBackend(test_model, test_driver)

            b.ρ = 0.0
            b.pₖ_norm = b.Δₖ = 1.0  # test for shrinking Δ

            update_Δₖ(b)
            @test b.Δₖ ≈ 0.25

            b.ρ = 0.25
            b.pₖ_norm = b.Δₖ = 1.0 # inclusive lower bound of ρ for keeping Δ

            update_Δₖ(b)
            @test b.Δₖ ≈ 1.0

            b.ρ = 0.75
            b.pₖ_norm = b.Δₖ = 1.0 # non-inclusive lower bound of ρ for increasing Δ

            update_Δₖ(b)
            @test b.Δₖ ≈ 1.0

            b.ρ = b.pₖ_norm = 0.9
            b.Δₖ = 1.0  # increasing Δ boundary requirement

            update_Δₖ(b)
            @test b.Δₖ ≈ 1.0
        end
    end
end

function backend_behavior_test(driver)

    b = BlockOptBackend(test_model, driver)

    initialize(b)

    @testset "initialize subroutine" begin

        @test b.fₖ ≈ obj(b, b.xₖ)

        @test b.∇fₖ ≈ grad(b, b.xₖ)

        @test b.∇fₖ_norm ≈ norm(grad(b, b.xₖ))

        @test b.hₖ ≈ hess_sample(test_model, b.xₖ, b.∇fₖ)

        @test b.Yₖ ≈ hess_sample(test_model, b.xₖ, b.Sₖ)

        @test rank(b.Hₖ) ≡ dimension(b)

        @test b.Hₖ ≈ b.Hₖ'

        @test b.Yₖ' * b.Sₖ ≈ b.Sₖ' * b.Yₖ

        @test b.Uₖ' * b.Vₖ ≈ b.Uₖ' * b.Vₖ

        @test b.Sₖ' * ∇²f * b.Sₖ ≈ b.Sₖ' * b.Yₖ

        # @test b.Yₖ ≈ b.Hₖ * b.Sₖ # line 231 Y₀ = H₀ S₀ (typo - test condition below)

        @test b.Yₖ ≈ inv(b.Hₖ) * b.Sₖ

        @test b.Hₖ * b.Yₖ ≈ b.Sₖ

        @test b.Uₖ ≈ [b.Sₖ b.∇fₖ / b.∇fₖ_norm]

        @test b.Vₖ ≈ [b.Yₖ b.hₖ / b.∇fₖ_norm]
    end

    iterations = 0

    while !terminal(b, iterations)

        iterations += 1

        build_trs(b)

        Mₖ = [b.∇fₖ / b.∇fₖ_norm (b.Hₖ * b.∇fₖ) / b.∇fₖ_norm b.Sₖ]

        Qₖ = [b.hₖ / b.∇fₖ_norm b.∇fₖ / b.∇fₖ_norm b.Yₖ]

        Hₖ⁻¹ = inv(b.Hₖ)

        # The tests below fail occosionally, due to error accumulation,
        # or possibly ill-conditioning. It only occurs when pflag = true,
        # so stability of including secant update should be investigated.
        # However, when it fails, convergence still occurs and the tests
        # break down in the terminal phase of the iteration.
        #     There is a similar problem in block blockQN subroutine tests

        # @test Hₖ⁻¹ * b.∇fₖ ≈ b.hₖ # lhs col of 5.3

        # @test Hₖ⁻¹ * b.∇fₖ ≈ b.hₖ # middle col of 5.3

        # @test Hₖ⁻¹ * b.Sₖ ≈ b.Yₖ  # rhs col of 5.3

        # @test Hₖ⁻¹'* Mₖ  ≈  Qₖ    # equivalence of 5.3

        """
        Equation 5.1: ``pₖ = arg min model_p(p) : ||p|| < Δ,``
        for p ∈ ℜⁿ; the standard n-dimensional trust-region subproblem.
        """
        mp(p) = 0.5 * p' * Hₖ⁻¹ * p + b.∇fₖ'p

        """
        Equation 5.2: ``qₖ = arg min model_q(q) : ||q|| < Δ,``
        for q ∈ ℜⁿ; when Hₖ is full rank, pₖ = Hₖqₖ.
        """
        mq(q) = 0.5 * q' * b.Hₖ * q + (b.Hₖ * b.∇fₖ)'q

        """
        Equation 5.3: ``aₖ = arg min model_a(a) : ||Hₖ Qₖ aₖ|| < Δ,``
        for a ∈ ℜ²ʷ⁺¹; solve to obtain step pₖ = HₖQₖaₖ and
        """
        ma(a) = 0.5 * a' * b.P * a + b.b'a

        @testset "build_trs subroutine" begin

            @test b.P ≈ b.Qₖ' * b.Hₖ * b.Qₖ

            @test b.b ≈ b.Qₖ' * b.Hₖ * b.∇fₖ

            @test b.C ≈ b.Qₖ' * b.Hₖ^2 * b.Qₖ

            @test b.Qₖ ≈ orth(Qₖ)

            @test b.Qₖ'b.Qₖ ≈ zeros(samples(b) + 1, samples(b) + 1) + I

            @test cond(b.Qₖ) ≈ 1.0 # this is redundant

            @test mp(zeros(dimension(b))) ≈ 0.0

            @test mq(zeros(dimension(b))) ≈ 0.0

            @test ma(zeros(samples(b) + 1)) ≈ 0.0

            a = randn(2w + 1)

            @test ma(a) ≈ mq(b.Qₖ * a)

            @test ma(a) ≈ mp(b.Hₖ * b.Qₖ * a) # line 186? p ≠ Mₖa

            @test sqrt(a' * b.C * a) ≈ norm(b.Hₖ * b.Qₖ * a)
        end

        solve_trs(b)

        @testset "solve_trs subroutine" begin

            @test b.qₖ ≈ b.Qₖ * b.aₖ

            @test b.pₖ ≈ b.Hₖ * b.qₖ

            @test ma(b.aₖ) ≈ mq(b.qₖ)

            @test ma(b.aₖ) ≈ mp(b.pₖ)

            @test ma(b.aₖ) ≤ 0.0

            @test b.∇fₖ' * b.pₖ ≤ 0.0  # descent direction

            @test abs(sqrt(b.aₖ' * b.C * b.aₖ) - b.pₖ_norm) ≤ sqrt(eps(Float64))

            @test abs(sqrt(b.aₖ' * b.C * b.aₖ) - norm(b.Hₖ * b.qₖ)) ≤ sqrt(eps(Float64))
        end

        build_trial(b)

        @testset "build_trial subroutine" begin

            @test b.xₜ ≈ b.xₖ + b.pₖ

            @test b.fₜ ≈ obj(b, b.xₜ)

            model_reduction = mq(b.qₖ)

            actual_reduction = b.fₜ - b.fₖ

            @test model_reduction ≤ 0.0

            @test b.ρ ≈ actual_reduction / model_reduction
        end

        successful_trial = accept_trial(b)

        @test successful_trial ≡ (b.ρ > 0.0)

        @test pflag(b) ≡ pflag(driver)

        successful_trial && pflag(b) && @testset "secantQN subroutine" begin

            sₖ = b.pₖ

            yₖ = grad(b, b.xₜ) - b.∇fₖ

            @test sₖ' * yₖ > 0.0    # curvature condition

            secantQN(b)

            @test rank(b.Hₖ) ≡ n

            @test b.Hₖ * yₖ ≈ sₖ
        end

        successful_trial && @testset "update_Sₖ subroutine" begin

            Sₖ = b.Sₖ

            update_Sₖ(b)

            @test b.Sₖ' * b.Sₖ ≈ I

            S_update(b) ≡ S_update_b && @test norm(b.Sₖ' * Sₖ) ≤ sqrt(eps(Float64))

            S_update(b) ≡ S_update_c && @test norm(b.Sₖ' * Sₖ) ≤ sqrt(eps(Float64))
        end

        successful_trial && @testset "gHS subroutine" begin

            gHS(b)

            @test b.∇fₖ ≈ grad(b, b.xₜ)

            @test b.∇fₖ_norm ≈ norm(grad(b, b.xₜ))

            @test b.hₖ ≈ hess_sample(test_model, b.xₖ, b.∇fₖ)

            @test b.Yₖ ≈ hess_sample(test_model, b.xₖ, b.Sₖ)

            @test b.Yₖ' * b.Sₖ ≈ b.Sₖ' * b.Yₖ

            @test b.Uₖ ≈ [b.Sₖ b.∇fₖ / b.∇fₖ_norm]

            @test b.Vₖ ≈ [b.Yₖ b.hₖ / b.∇fₖ_norm]

            @test b.Uₖ' * b.Vₖ ≈ b.Uₖ' * b.Vₖ
        end

        successful_trial && @testset "blockQN subroutine" begin

            blockQN(b)

            @test rank(b.Hₖ) ≡ n

            @test b.Hₖ ≈ b.Hₖ'

            # The tests below fail occosionally, but it is not expected that
            # it will hold in terminal phase.

            # @test b.Hₖ * b.hₖ ≈ b.∇fₖ

            # @test b.Hₖ * b.Vₖ ≈ b.Uₖ

            # @test b.Hₖ * b.Yₖ ≈ b.Sₖ
        end

        update_Δₖ(b)

    end

    @test b.∇fₖ_norm ≤ ϵ_tol(driver)

    nothing
end

@testset "Backend Behavior Tests:" begin

    backend_behavior_test(Driver(S_update = S_update_a))

    backend_behavior_test(Driver(S_update = S_update_b))

    backend_behavior_test(Driver(S_update = S_update_c))

    backend_behavior_test(Driver(S_update = S_update_d))

    backend_behavior_test(Driver(S_update = S_update_e))

    backend_behavior_test(Driver(S_update = S_update_f))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_a))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_b))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_c))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_d))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_e))

    backend_behavior_test(Driver(pflag = true, S_update = S_update_f))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_a))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_b))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_c))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_d))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_e))

    backend_behavior_test(Driver(QN_update = PSB, S_update = S_update_f))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_a))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_b))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_c))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_d))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_e))

    backend_behavior_test(Driver(QN_update = PSB, pflag = true, S_update = S_update_f))
end


text:
export pinvfact, PInv, PInvFact, rank, adjoint

abstract type PInv{T} <: Factorization{T} end

"""
Structure to hold pseudo-inverse factorization.
"""
mutable struct PInvFact{T} <: PInv{T}
    m::Int                  # number of rows of input matrix
    n::Int                  # number of columns
    k::Int                  # rank
    F::QRWrapper   # first QR
    G::QRWrapper   # second QR (in case k < n)
end

Base.propertynames(PI::PInvFact) = (:R, :Q, :prow, :pcol)

import LinearAlgebra: rank, adjoint
rank(psi::PInv) = rank(psi.F)
rank(psi::Adjoint{<:Number,<:PInv}) = rank(psi.parent)

adjoint(psi::PInv) = Adjoint(psi)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, apsi::Adjoint{T,<:PInv}) where T
    print(io, typeof(apsi))
    print(io, apsi.parent)
end

"""
    pinvfact(A::AbstractMatrix; tolrel, tolabs)

Create a Moore-Penrose pseudo-inverse factorization of input matrix A.
A Moore-Penrose pseudo-inverse is represented as one or two pivoted
QR-factorizations. The first factorization determines also the rank of the
matrix. The second one is applied to the transposed of the results of the
first factorization.

The usual way to calculate a pseudo inverse (pinv) makes use of a "SVD"
(singular value decomposition) of `A`, which in turn requires an eigenvalue
decomposition of `A'A`. As the pseudo-inverse itself is typically dense, also
if the input matrix is sparse, it is expensive to calculate it as a matrix.
This implementation stores only the results of two QR-decompositions, which
remain sparse, which makes the method suitable for big sparse matrices as well.

Usage:
psi = pinvfact(A)
x = psi * b
"""
function pinvfact(A::AbstractMatrix; tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0)
    m, n = size(A)
    # TODO decide if that makes sense:
    # 1. probable use case m >= n
    # 2. better error bounds on rank-determining pivots if m < n ???
    if m >= n
        tol = tolerance(A, true, tolrel, tolabs)
        F = qrfactors(A, tolabs=tol)
        k = rank(F, tolabs=tol)
        G = qrfactors(copy(adjoint(adjustsize(F.R, k, n))), tolabs=tol/2)
        @assert k == rank(G) "rank defect $(rank(G)) < $k during second QR factorization"
        PInvFact{eltype(A)}(m, n, k, F, G)
    else
        adjoint(pinvfact(copy(adjoint(A)), tolrel=tolrel, tolabs=tolabs))
    end
end

function pinvfact(A::Adjoint{<:AbstractMatrix};
                  tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0)

    adjoint(pinvfact(adjoint(A), tolrel=tolrel, tolabs=tolabs))
end

Base.size(PI::PInvFact, args...) = size(PI.F, args...)
Base.size(PIA::Adjoint{<:Number,<:PInv}) = reverse(size(PIA.parent))

import Base.\

"""
    Apply the pseudo-inverse factorization to the rhs `A`.
"""
function \(PI::PInv{T}, A::StridedVecOrMat{T}) where {T}
    m, n = size(PI)
    if m != size(A, 1)
        throw(ArgumentError("Dimension mismatch"))
    end
    k = PI.k
    F = PI.F
    G = PI.G
    nr = size(A, 2)
    tmpm = similar(A, ntuple(i-> i == 1 ? m : nr, ndims(A)))
    tmpm[:,:] = view(A, F.prow, :)
    lmul!(F.Q', tmpm)
    tmpx = view(tmpm, ntuple(i-> 1:(i == 1 ? k : nr), ndims(A))...)
    tmpx = view(tmpx, G.pcol, :)
    y1 = adjustsize((G.R' \ tmpx), n, nr)
    y1 = (G.Q * y1)[invperm(G.prow), :]
    y1[invperm(F.pcol), :]
end

"""
    Apply the adjoint of pseudoinverse factorization to rhs.
"""
function \(PIT::Adjoint{T,PS}, A::StridedVecOrMat{T}) where {T,PS<:PInv{T}}
    PI = PIT.parent
    m, n = size(PI)
    if n != size(A, 1)
        throw(ArgumentError("Dimension mismatch"))
    end
    k = PI.k
    F = PI.F
    G = PI.G
    nr = size(A, 2)
    x1 = similar(A)
    x1[invperm(G.prow),:] = A[F.pcol, :]
    lmul!(G.Q', x1)
    x2 = x1[1:k,:]
    x2[G.pcol,:] = G.R \ x2
    (F.Q * adjustsize(x2, m, nr))[invperm(F.prow),:]
end

function \(PI::PInv{T}, A::StridedVecOrMat{S}) where {T,S}
    P = promote_type(T, S)
    PI \ P.(A)
end

function \(PI::Adjoint{T,PS}, A::StridedVecOrMat{S}) where {T,S,PS<:PInv{T}}
    P = promote_type(T, S)
    PI \ P.(A)
end


text:
using JLLWrappers
using Pkg
using Test

module TestJLL end

@testset "JLLWrappers.jl" begin
    mktempdir() do dir
        Pkg.activate(dir)
        # Package with a FileProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "Vulkan_Headers_jll")))
        @test_nowarn @eval TestJLL using Vulkan_Headers_jll
        @test @eval TestJLL Vulkan_Headers_jll.is_available()
        @test isfile(@eval TestJLL vk_xml)
        @test isfile(@eval TestJLL Vulkan_Headers_jll.vk_xml_path)
        @test isfile(@eval TestJLL Vulkan_Headers_jll.get_vk_xml_path())
        @test isdir(@eval TestJLL Vulkan_Headers_jll.artifact_dir)
        @test isempty(@eval TestJLL Vulkan_Headers_jll.PATH[])
        @test occursin(Sys.BINDIR, @eval TestJLL Vulkan_Headers_jll.LIBPATH[])
        # Package with an ExecutableProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "HelloWorldC_jll")))
        @test_nowarn @eval TestJLL using HelloWorldC_jll
        if Sys.isfreebsd()
            @test @eval TestJLL !HelloWorldC_jll.is_available()
        else
            @test @eval TestJLL HelloWorldC_jll.is_available()
            @test "Hello, World!" == @eval TestJLL hello_world(h->readchomp(`$h`))
            @test isfile(@eval TestJLL HelloWorldC_jll.hello_world_path)
            @test isfile(@eval TestJLL HelloWorldC_jll.get_hello_world_path())
            @test isdir(@eval TestJLL HelloWorldC_jll.artifact_dir)
            @test !isempty(@eval TestJLL HelloWorldC_jll.PATH[])
            @test occursin(Sys.BINDIR, @eval TestJLL HelloWorldC_jll.LIBPATH[])
        end
        # Package with a LibraryProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "OpenLibm_jll")))
        @test_nowarn @eval TestJLL using OpenLibm_jll
        @test @eval TestJLL OpenLibm_jll.is_available()
        @test exp(3.14) ≈ @eval TestJLL ccall((:exp, libopenlibm), Cdouble, (Cdouble,), 3.14)
        @test isfile(@eval TestJLL OpenLibm_jll.libopenlibm_path)
        @test isfile(@eval TestJLL OpenLibm_jll.get_libopenlibm_path())
        @test isdir(@eval TestJLL OpenLibm_jll.artifact_dir)
        @test isempty(@eval TestJLL OpenLibm_jll.PATH[])
        @test occursin(Sys.BINDIR, @eval TestJLL OpenLibm_jll.LIBPATH[])

        # Issue #20
        if Sys.iswindows()
            @test Sys.BINDIR ∈ JLLWrappers.get_julia_libpaths()
        end
    end
end


text:
using FCSFiles
using FileIO
using Test, HTTP

@testset "FCSFiles test suite" begin

    # test loading an FCS 2.0 file
    @testset "Loading an FCS 2.0 file" begin
        # download the FCS 2.0 file
        cwd = pwd()
        cd(cwd*"/testdata")
        @info "Downloading FCS 2.0 file ..."
        io = open("testFCS2.fcs", "w")
        r = HTTP.request("GET", "https://flowrepository.org/experiments/4/fcs_files/326/download", response_stream=io)
        close(io)
        cd(cwd)
        @info "Done."

        # load the FCS 2.0 file
        @test_throws ErrorException @test_warn "FSC2.0 files are not guaranteed to work" flowrun = load("testdata/testFCS2.fcs")

        # cleanup
        rm("testdata/testFCS2.fcs", force=true)
        @info "FCS 2.0 file removed"
    end

    # test the size of the file
    @testset "SSC-A size" begin
        flowrun = load("testdata/BD-FACS-Aria-II.fcs")

        @test length(flowrun["SSC-A"]) == 100000
    end

    # test the loading of a large FCS file
    @testset "Loading of large FCS file" begin
        # download the large FCS file
        cwd = pwd()
        cd(cwd*"/testdata")
        @info "Downloading large FCS file ..."
        io = open("testLargeFile.fcs", "w")
        r = HTTP.request("GET", "https://flowrepository.org/experiments/1177/fcs_files/113151/download", response_stream=io)
        close(io)
        cd(cwd)
        @info "Done."

        # load the large file
        flowrun = load("testdata/testLargeFile.fcs")
        @test length(flowrun.data) == 50
        @test length(flowrun.params) == 268

        # cleanup
        rm("testdata/testLargeFile.fcs", force=true)
        @info "Large file removed"
    end
end


text:
# mpirun -n 16 julia MPI_forward.jl
using Revise
using ADSeismic
using ADCME
using JLD2
ADCME.options.customop.verbose = false
mpi_init()
M = N = Int(sqrt(mpi_size()))
n = 100

param = MPIElasticPropagatorParams(NX= n * M, NY= n * N, n=n, NSTEP=2000, DELTAT=1e-4/2,
                        DELTAX=1.0, DELTAY=1.0)
compute_PML_Params!(param)
vp = 3300.
vs = 3300. / 1.732
rho = 2800.
λ, ρ, μ = compute_default_properties(param, vp, vs, rho)
source = Ricker(param, 50.0, 200.0)

srci = [div(param.NX,5)]
srcj = [div(param.NY,2)]
srctype = [0]
srcv = repeat(source, 1, 1)
src = MPIElasticSource(param, srci, srcj, srctype, srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
rcvtype = zeros(Int64, length(rcvj))
receiver =  MPIElasticReceiver(param, rcvi, rcvj, rcvtype)

propagator = MPIElasticPropagatorSolver(param, src, ρ, λ, μ)
MPISimulatedObservation!(propagator, receiver)

###### second source ######

srci = [div(param.NX,5)]
srcj = [div(param.NY,2)]
srctype = [0]
srcv = repeat(source, 1, 2)
src = MPIElasticSource(param, srci, srcj, srctype, srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
rcvtype = zeros(Int64, length(rcvj))
receiver2 =  MPIElasticReceiver(param, rcvi, rcvj, rcvtype)

propagator2 = MPIElasticPropagatorSolver(param, src, ρ, λ, μ)
MPISimulatedObservation!(propagator2, receiver2)

function is_not_noop(c)
    try
        length(c)
    catch
        return false
    end
    return true
end

loss = constant(0.0)
if is_not_noop(receiver)
    loss = sum(receiver.rcvv^2) + sum(receiver2.rcvv^2)
end
loss = mpi_sum(loss)

sess = Session(); init(sess)
@info "starting..."
L = run(sess, loss)

if mpi_size()>1
    mpi_finalize()
end
# if mpi_rank()==0 && mpi_size()==1
#     close("all")
#     pcolormesh(rcvv, cmap="gray")
#     colorbar()
#     xlabel("Location")
#     ylabel("Time")
#     savefig("receiver.png")

#     close("all")
#     visualize_wavefield(Vx, param)

# end


text:
module Player

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Player (Spawn Point)" => Ahorn.EntityPlacement(
        Maple.Player
    )
)

sprite = "characters/player/sitDown00.png"

function Ahorn.selection(entity::Maple.Player)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Player) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)

end


text:
module AnalyticalEOR

using Base:Real
using ForwardDiff
using Roots
using Plots

include("types.jl")
include("flowFunctions.jl")
include("buckleyLeverett.jl")
include("polymerFlooding.jl")
include("tracers.jl")
include("plottingFunctions.jl")

export RelPerms, WaterFlooding
export water_rel_perm, oil_rel_perm, krw_derivative, kro_derivative
export fractional_flow, fw_derivative
export solve_waterflooding, saturation_profile
export solve_tracer, tracer_profile
export plot_fw, plot_sw_profile, animate_sw_profile
export plot_tracer_fw
export solve_polymerflood

end


text:
module RecurrentTests

using Photon, Test

resetContext()
ctx = getContext()

gethistoric(steps=100, features=10) = KorA(randn(ctx.dtype,features,steps,4))

function basic_layer(T)

    d = gethistoric(100)

    l = T(20)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20,2)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20;dropout=0.5, bidirectional=true)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20;dropout=0.5, last_only=false)
    pred = l(d)
    @test size(pred) == (20,100,4)

    l = T(20; initw=zeros, initb=zeros, initf=zeros)
    pred = l(d)
    @test size(pred) == (20,4)
    @test sum(pred) == 0.0

end

function lstm_model()

    model = Sequential(
        LSTM(20),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)

    pred = model(gethistoric(200))
    @test size(pred) == (10,4)

end

function lstm_model2()

    model = Sequential(
        LSTM(20, last_only=false),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)

    @test_throws DimensionMismatch model(gethistoric(200))
end

function gru_model()

    model = Sequential(
        GRU(20,2;dropout=0.5),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)
end

function rnn_model()

    model = Sequential(
        RNN(20,2;activation=:relu, dropout=0.5),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)
end

@testset "Recurrent" begin
    basic_layer(LSTM)
    basic_layer(GRU)
    basic_layer(RNN)
    lstm_model()
    lstm_model2()
    gru_model()
    rnn_model()
end

end


text:
#!/usr/local/bin/julia

#=
Julia has a built-in package manager for installing add-on
functionality written in Julia.

The package manager is
  - declarative rather than imperative
=#

# Print what packages are installed, their version and state
Pkg.status()

# Declare a new requirement
Pkg.add("Distributions")

#=
Pkg.add first adds "Distributions" to your
  ~/.julia/v0.6/REQUIRE file
and then runs Pkg.resolve() to determine the
best set of packages
=#

# Verify Distributions was installed
Pkg.status()

# Package directory
Pkg.dir("Distributions")

# Remove a requirement
Pkg.rm("Distributions")

# Update your packages
Pkg.update()

# Checkout, Pin, and Free
Pkg.checkout("Distributions")
Pkg.pin("Distributions")
Pkg.free("Distributions")

#=
PkgDev for creating your own packages:
  - PkgDev.generate("PackageName", "License")
  - PkgDev.register("PackageName")
  - PkgDev.publish("PackageName")
=#

max_stars_repo_name: ellisvalentiner/julia-talk

max_issues_repo_name: ellisvalentiner/julia-talk

max_forks_repo_name: ellisvalentiner/julia-talk


text:
module RoadefRenaultBrazil

using Random

include("parser.jl")
include("solution.jl")
include("functions.jl")
include("constants.jl")
include("greedy.jl")
include("ils_hprc.jl")
include("vns_lprc.jl")
include("repair.jl")
include("vns_pcc.jl")

end # module


text:
"""
    FileBuffer <: IO

A generic buffer type to provide an IO interface for none IO based path types.

NOTES:
- All `read` operations will read the entire file into the internal buffer at once.
  Subsequent calls to `read` will only operate on the internal buffer and will not access
  the filepath.
- All `write` operations will only write to the internal buffer and `flush`/`close` are
  required to update the filepath contents.
"""
struct FileBuffer <: IO
    path::AbstractPath
    io::IOBuffer
    read::Bool
    write::Bool
    create::Bool
end

function FileBuffer(
    fp::AbstractPath; read=true, write=false, create=false, truncate=false, append=false
)
    buffer = FileBuffer(fp, IOBuffer(), read, write, create)

    # If we're wanting to append data then we we need to prepopulate the internal buffer
    if write && append
        _read(buffer)
        seekend(buffer)
    end

    return buffer
end

Base.isreadable(buffer::FileBuffer) = buffer.read
Base.iswritable(buffer::FileBuffer) = buffer.write
Base.seek(buffer::FileBuffer, n::Integer) = seek(buffer.io, n)
Base.seekstart(buffer::FileBuffer) = seekstart(buffer.io)
Base.seekend(buffer::FileBuffer) = seekend(buffer.io)
Base.position(buffer::FileBuffer) = position(buffer.io)
function Base.eof(buffer::FileBuffer)
    if position(buffer) == 0
        _read(buffer)
        seekstart(buffer)
    end
    return eof(buffer.io)
end

function _read(buffer::FileBuffer)
    # If our IOBuffer is empty then populate it with the
    # filepath contents
    if buffer.io.size == 0
        write(buffer.io, read(buffer.path))
    end
end

function Base.read(buffer::FileBuffer)
    isreadable(buffer) || throw(ArgumentError("read failed, FileBuffer is not readable"))
    _read(buffer)
    seekstart(buffer)
    read(buffer.io)
end

function Base.read(buffer::FileBuffer, ::Type{String})
    isreadable(buffer) || throw(ArgumentError("read failed, FileBuffer is not readable"))
    _read(buffer)
    seekstart(buffer)
    read(buffer.io, String)
end

function Base.read(buffer::FileBuffer, ::Type{UInt8})
    if buffer.io.size == 0
        write(buffer.io, read(buffer.path))
        seekstart(buffer)
    end
    read(buffer.io, UInt8)
end

#=
NOTE: We need to define multiple methods because of ambiguity error with base IO methods.
=#
function Base.write(buffer::FileBuffer, x::Vector{UInt8})
    iswritable(buffer) || throw(ArgumentError("write failed, FileBuffer is not writeable"))
    write(buffer.io, x)
end

function Base.write(buffer::FileBuffer, x::String)
    iswritable(buffer) || throw(ArgumentError("write failed, FileBuffer is not writeable"))
    write(buffer.io, x)
end

function Base.write(buffer::FileBuffer, x::UInt8)
    iswritable(buffer) || throw(ArgumentError("write failed, FileBuffer is not writeable"))
    write(buffer.io, x)
end

function Base.flush(buffer::FileBuffer)
    if iswritable(buffer)
        seekstart(buffer)
        write(buffer.path, read(buffer.io))
    end
end

function Base.close(buffer::FileBuffer)
    flush(buffer)
    close(buffer.io)
end


text:
module ScatteringLayers

using Optim
using NNLS
using Distributions
using DiffEqOperators
using Images
using ProgressMeter

export AbstractSpaceTimePoint,
    SpacePoint,
    SpaceTimePoint,
    Echogram,
    Ping,
    depths,
    spacetime,
    getping,
    LayerMixture,
    getpeaks,
    guesslayers,
    mixpdf,
    getbasis!,
    getbasis,
    refine_heights!,
    refine_σ!,
    refine_μσ!,
    fit!,
    fit,
    fitlayers

include("echograms.jl")
include("layerfitting.jl")

end # module

max_stars_repo_name: ElOceanografo/ScatteringLayers

max_issues_repo_name: ElOceanografo/ScatteringLayers

max_forks_repo_name: ElOceanografo/ScatteringLayers