text:
# Let's increase the number of workers
#-------------------------------------
maxNbWorkers = 3
while nworkers() < maxNbWorkers
  addprocs(1)
end

@everywhere using SMM
@everywhere using Distributions
using Optim

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@static if VERSION < v"0.7.0"
    @everywhere using DataStructures
else
    @everywhere using OrderedCollections
    # Because of this issue (https://github.com/JuliaIO/JLD2.jl/issues/107)
    # we also need to import BlackBoxOptim and Optim to load and save
    #----------------------------------------------------------------------
    # (no need when using Julia 0.6)
    using BlackBoxOptim

end

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

@testset "testing Types" begin

    @testset "testing SMMOptions" begin

        saveName = get_now()

        t = SMMOptions()

        # Testing default values
        #-----------------------------------------------------------------------
        @test t.globalOptimizer == :dxnes
	      @test t.localOptimizer == :LBFGS
        @test t.maxFuncEvals == 1000
        @test t.saveSteps == 100
        @test t.saveName == saveName
        @test t.showDistance == false
        @test t.minBox == false
        @test t.populationSize == 50
        @test t.penaltyValue == 999999.0
        @test t.gridType == :latin
        @test t.saveStartingValues == true
        @test t.maxTrialsStartingValues == 1000
        @test t.thresholdStartingValue == 99999.0

    end

    @testset "testing SMMOptions" begin

        t = SMMProblem()

        # When initialized, iter is equal to 0
        @test t.iter == 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) == SMMOptions

    end

end

@testset "testing creation of Cartesian grids" begin

  a = zeros(3)
  b = ones(3)
  nums = 2

  points = create_grid(a, b, nums)

  @test points[1,:] == [0.0; 0.0; 0.0]
  @test points[2,:] == [1.0, 0.0, 0.0]
  @test points[3,:] == [0.0, 1.0, 0.0]
  @test points[4,:] == [1.0, 1.0, 0.0]
  @test points[5,:] == [0.0, 0.0, 1.0]
  @test points[6,:] == [1.0, 0.0, 1.0]
  @test points[7,:] == [0.0, 1.0, 1.0]
  @test points[8,:] == [1.0, 1.0, 1.0]

end

@testset "testing creation of random grids" begin

  a = zeros(3)
  b = [1.0; 2.0; 3.0]
  numPoints = 100

  @testset "drawing from a multivariate normal" begin

    points = create_grid_stochastic(a, b, numPoints, gridType = :normal, alpha = 0.05)

    @test size(points, 1) == numPoints
    @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

  end

  @testset "drawing from a multivariate uniform" begin

    points = create_grid_stochastic(a, b, numPoints, gridType = :uniform)

    @test size(points, 1) == numPoints
    @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

  end

  @testset "Latin hypercube sampling" begin

    points = latin_hypercube_sampling(a, b, numPoints)

    @test size(points, 1) == numPoints
    @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

  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 FminBox, non binding" begin

  atolOptim = 0.5
  # ParticleSwarm is not working here
  # I don't know why
  listValidLocalOptimizers = [:NelderMead, :SimulatedAnnealing,
                            :BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
                            :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, lower, upper, iterations = 2000, convert_to_fminbox(localOptim))

      @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
  # ParticleSwarm is not working here
  # I don't know why
  listValidLocalOptimizers = [:NelderMead, :SimulatedAnnealing,
                            :BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
                            :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 = [0.5; 0.5]

    for localOptim in listValidLocalOptimizers

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

      @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("SMM"), "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("SMM"), "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!" begin

    t = SMMProblem()

    @testset "testing set_priors!" begin

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

        set_priors!(t, dictPriors)

        @test t.priors == dictPriors

    end

    @testset "set_empirical_moments!" begin

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

        set_empirical_moments!(t, dictEmpiricalMoments)

        @test t.empiricalMoments == dictEmpiricalMoments

    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 = SMMProblem()

        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

        srand(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 = SMMProblem()

        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("SMM"), "test/empiricalMomentsTest.csv"))
        set_empirical_moments!(t, dictEmpiricalMoments)

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

        # 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 = SMMProblem()

        dictPriors = read_priors(joinpath(Pkg.dir("SMM"), "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 = SMMProblem()

        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 loading and saving an optimization" begin

    srand(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 = SMMProblem()

    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("SMM"), "test/empiricalMomentsTest.csv"))
    set_empirical_moments!(t, dictEmpiricalMoments)

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

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

    saveSMMOptim(t, saveName = "iamatest")
    t2 = loadSMMOptim("iamatest")

    # Test the objective function is correctly loaded
    #------------------------------------------------
    # [TODO] test other fields
    @test t.objective_function([dictEmpiricalMoments["meanU"][1]]) ≈ t2.objective_function([dictEmpiricalMoments["meanU"][1]]) atol = tol1dMean

end

@testset "Testing smmoptimize" begin

    # 1d problem
    #-----------
    @testset "Testing smmoptimize with 1d" begin

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

        @everywhere function functionTest1d(x)

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

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

            return output
        end

        t = SMMProblem(options = SMMOptions(maxFuncEvals=200,saveSteps = 100))

        # 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("SMM"), "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)

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

        smmoptimize!(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 smm_refine_globalmin!(t, verbose = true)[1] ≈ 0.05 atol = tol1dMean

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

    end

    @testset "Testing local to global 1d" begin

      tol1dMean = 0.1

      @everywhere function functionTest1d(x)

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

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

          return output
      end

      t = SMMProblem(options = SMMOptions(maxFuncEvals=200,saveSteps = 100))

      # 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("SMM"), "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)

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

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

      @test smm_local_minimum(t) ≈ 0.0 atol = tol1dMean

      @test smm_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
          #-----------------------------------------------------
          srand(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 = SMMProblem(options = SMMOptions(maxFuncEvals=200, saveSteps = 100, 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("SMM"), "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)

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

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

      @test smm_local_minimum(t) ≈ 0.0 atol = tol1dMean

      @test smm_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 srand()
        # 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 = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 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)

        # 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
        smmoptimize!(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 local_to_global! with 2d and same magnitude" begin

        # Rermark:
        # It is important NOT to use srand()
        # 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
            #-----------------------------------------------------
            srand(1234)
            draws = rand(d, 1000000)
            output["mean1"] = mean(draws[1,:])
            output["mean2"] = mean(draws[2,:])

            return output
        end

        t = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 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)

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

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

        @test smm_local_minimum(t) ≈ 0.0 atol = tol2dMean

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

      end

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

          # Rermark:
          # It is important NOT to use srand()
          # 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
              #-----------------------------------------------------
              srand(1234)
              draws = rand(d, 1000000)
              output["mean1"] = mean(draws[1,:])
              output["mean2"] = mean(draws[2,:])

              return output
          end

          t = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 1000, localOptimizer = :NelderMead, 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)

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

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

          @test smm_local_minimum(t) ≈ 0.0 atol = tol2dMean

          @test smm_local_minimizer(t)[1] ≈ 1.0 atol = tol2dMean
          @test smm_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 srand()
        # 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 = 0.5

        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 = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 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)

        # 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
        smmoptimize!(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 function functionTest2d(x)

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

              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 = SMMProblem(options = SMMOptions(maxFuncEvals=2000,saveSteps = 2000, globalOptimizer = :dxnes, localOptimizer = :NelderMead, minBox = false))

          #---------------------------------------------------------------------
          # 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]
          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)

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

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

          @test smm_local_minimum(t) ≈ 0.0 atol = tol2dMean

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

          #---------------------------------------------------------------------
          # Using BlackBoxOptim
          #---------------------------------------------------------------------
          t = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 1000, globalOptimizer = :dxnes, localOptimizer = :NelderMead, minBox = false))

          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)

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

          smmoptimize!(t, verbose = true)

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

end

@testset "Testing Inference" begin

  tolLinear = 0.05

  # Inference in the linear model
  #------------------------------
  srand(1234)         #for replicability reasons
  T = 100000          #number of periods
  P = 2               #number of dependent variables
  beta0 = rand(P)     #choose true coefficients by drawing from a uniform distribution on [0,1]
  alpha0 = rand(1)[]  #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 = SMMProblem(options = SMMOptions(maxFuncEvals=500, saveSteps = 500, globalOptimizer = :dxnes, localOptimizer = :LBFGS, minBox = false, showDistance = false));

  # 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)

  # x[1] corresponds to the intercept
  # x[1] corresponds to beta1
  # x[3] corresponds to beta2
  @everywhere function functionLinearModel(x; nbDraws::Int64 = 1000000, burnInPerc::Int64 = 10)

      # Structural Model
      #-----------------
      srand(1234) #for replicability reasons
      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()
      U[1] = rand(d, 1)[] #first error term
      # loop over time periods
      for t = 2:T
          U[t] = rand(d, 1)[] + theta*U[t-1]
      end

      simX = zeros(T, P)
      d = Uniform(0, 5)
      for p = 1:P
              simX[:,p] = rand(d, T)
      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

  set_simulate_empirical_moments!(myProblem, functionLinearModel)

  # 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
  #--------------------------------------------------------------------
  @time listOptimResults = local_to_global!(myProblem, verbose = true)

  minimizer = smm_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

  nbDraws = 1000000 #number of draws in the simulated data
  calculate_Avar!(myProblem, minimizer, 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)

  # first column : point estimates
  @test df[:Estimate] == minimizer

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

  # confidence interval
  for i =1:size(df,1)
    @test df[:ConfIntervalLower][i] <= df[:ConfIntervalUpper][i]
  end

end


text:
f = x->0;
u_exact = x->x;

a = 0;
b = 1;
n = 9;

problem = FiniteDifferenceBVPProblem(0, 1, n, f, DirichletBC(0), DirichletBC(1));
assemble_system!(problem);
u = solve_bvp(problem);
err = @. problem.x - u;
# check error
norm(err, Inf)< 1e-14

ext: jl


text:
using Statistics, Dates

function test_ach()

	@eval using Plots
	gr()

	nn = 1
	ts = 10
	ACh = fill(0.0, nn)
	synt = 3 # or 8
	Φ = 0.5
	spikesequence = repeat([1, 0, 0, 0, 0], outer=2)
	@show spikesequence
	spt = sim_spikes(nn, ts, spikesequence)
	p = scatter()

	for t = 1:ts

		update_ACh!(ACh, synt, Φ, t, spt[t])
		scatter!(p, [t], [ACh[1]])

	end

	display(p)
end

function sim_spikes(nn, ts, spikesequence)

	lastspike = 1
	out = []

	for i = 1 : ts

		if spikesequence[i] == 1
			push!(out, i)
			lastspike = i
		else
			push!(out, lastspike)
		end

	end
	return out
end

function test_da()

	@eval using Plots
	gr()

	ts = 40
	τ = 20
	da = 0
	BA = repeat([0.5, 0.0, 0.0, 0.0],inner=8, outer=5)
	p = Plots.scatter()

	for t = 1:ts
		da = update_da(da, BA[t], τ; δt=1)
		@show da, BA[t]
		Plots.scatter!(p, [t], [da], show=true, legend=false)
	end
end

function test_transmission()

	# get2 layers, with one synapse between them, and see how activity in the 1st translates to activity in the second.
	nn = [1, 1]
	ts = 60
	ACh = [0.0]
	w = 4.0
	spt1 = sim_spikes(1, 60, repeat([1, 1, 0],inner=20))
	spt2 = 0
	synt = 3 # or 8
	Φ = 0.5
	v = [-60.0]
	rec = [0.0]
	sp = [false]
	C = 100
	a = 0.3
	b = -0.3
	c = -65
	d = 8
	k = 2

	p = Plots.scatter()

	for t = 1:ts

		update_ACh!(ACh, synt, Φ, t, spt1[t])
		# calc output.
		input = (-v) .* w .* ACh
		update_activation!(1, v, -60, sp, spt2, t, -40, rec, input, C, a, b, c, d, k;)
		Plots.scatter!(p,[t], v, legend=false, color=RGB(0.5, 0.5, 1))

	end

	display(p)
end

function test_weight()

	@eval using Plots
	gr()

	# get2 layers, with one synapse between them, and see how activity in the 1st translates to activity in the second.
	nn = [1, 1]
	ts = 60
	ACh = hcat(0.0)
	w = hcat(4.0)
	γ = hcat(0.0)
	connected = [1]
	spt1 = sim_spikes(1, 60, repeat([1, 1, 0],inner=20))
	spt2 = hcat(0)
	synt = 3 # or 8
	Φ = 0.5
	v = hcat(-60.0)
	rec = hcat(0.0)
	sp = hcat(false)
	C = 100
	a = 0.3
	b = -0.3
	c = -65
	d = 8
	k = 2

	p = scatter()

	for t = 1:ts
		update_ACh!(ACh, synt, Φ, t, spt1[t])
		input = (-v) .* w .* ACh
		update_activation!(1, v, -60, sp, spt2, t, -40, rec, input, C, a, b, c, d, k;)
		update_weights!(w, γ, connected, t, spt1[t], spt2, 1, 40; δt=1, A₋=-1, t₋=15)

		scatter!(p, [t], w, legend=false, color=RGB(0.5, 0.5, 1))
	end
	display(p)
end

const stimtype = (ColorInput,)
const sstages = [20, 5, 50]
const inputstages = Bool[1,1,1]
const dastages = Bool[0,1,0]
const reportvar = "spiked"
const init_da=0

struct normdata
	mean::Float64
	std::Float64
end
normdata(intuple::Tuple{Float64,Float64}) = normdata(intuple[1],intuple[2])
function normdata(norm::Vector{normdata})

	mean = 0
	std = 0

	for n in norm
		mean += n.mean
		std += n.std
	end

	return normdata(mean/length(norm),std/length(norm))
end
normdata(a::Array) = normdata(mean(a), std(a))

function setup()

	p = get_parameters()
	m = MatrixTypes(initialise_matrices(p.nn, p)...)
	return (p, m)
end

function setup(p::ParameterTypes)

	m = MatrixTypes(initialise_matrices(p.nn, p)...)
	return (p, m)
end

function initialise_reporter(sensory_input, m)
	[initialise_return_variable(duration(s), m, reportvar) for (i,s) in enumerate(sensory_input)]
end

function train_input!(p, m, sensory_input::AbstractInput)
	numsteps = duration(sensory_input)
	da = Dopamine(init_da)
	run_all_steps(p.nn, numsteps, m, p, sensory_input, da, savevars=nothing, update=true, normweights=true)
	reset!(m,p)
	return
end

InputArrays = Union{InputSequence, Vector{<:AbstractInput}}

function train_input!(p, m, sensory_input::InputArrays)

	for s in sensory_input
		train_input!(p, m, s)
	end
end

function test_input(p, m, sensory_input::AbstractInput)
	reset!(m,p)
	numsteps = duration(sensory_input)
	da = Dopamine(init_da)
	return run_all_steps(p.nn, numsteps, m, p, sensory_input, da, savevars=nothing, update=false, normweights=false, reportvar="spiked")
end

function test_input(p, m, sensory_input::InputArrays)
	report = initialise_reporter(sensory_input, m)
	for (i, s) in enumerate(sensory_input)
		report[i] = test_input(p, m, s)
	end
	return report
end

function spiked_statistics(reporter::Array{NeuronLayers{Bool,1}}; layer=3)

	meanspikes = zeros(length(reporter))
	for (i,timestep) in enumerate(reporter)
		meanspikes[i] = mean(timestep.layers[layer])
	end

	return normdata(mean(meanspikes), std(meanspikes))
end

function spiked_statistics(reporter::Array{Array{NeuronLayers{Bool,1},1},1}; layer=3)
	stats = Vector{normdata}(undef, length(reporter))
	for (i, epoch) in enumerate(reporter)
		stats[i] = normdata(spiked_statistics(epoch, layer=layer))
	end
	return stats
end

function kenyon_cell_representation(p, m, sensory_input::AbstractInput)
	spiked = test_input(p,m,sensory_input)
	representation = spike_representation(spiked, p.nn[2], layer=2)
end

function spike_representation(spiked, neurons_in_layer; layer=2)

	representation = zeros(neurons_in_layer)
	for timestep in spiked
		representation += timestep.layers[layer]
	end
	representation ./ length(spiked)
end

function spike_learning_measure(p,m, train, test)

    for sensory_train in train
        train_input!(p,m,sensory_train)
    end

    test_spikes = Vector{normdata}(undef, length(test))
    for (i,sensory_test) in enumerate(test)
        test_spikes[i] = spiked_statistics(test_input(p,m,sensory_test))
    end

    train_spikes = similar(test_spikes)
    for (i,sensory_train) in enumerate(train)
        train_spikes[i] = spiked_statistics(test_input(p,m,sensory_train))
    end

    return (train_spikes, test_spikes)
end

function manhatten_distance(array1, array2)

	sum( abs.( array1 .- array2 ) )
end

function nearest_neighbour(array1, all_other_arrays)

	distances = zeros(length(all_other_arrays))

	for (i,a) in enumerate(all_other_arrays)

		distances[i] = manhatten_distance(array1, a)

	end

	return findmin(distances)
end


text:
__precompile__()

module SoilTracers

include("consts.jl")  # Constants
include("physchem/physchem.jl")  # PhysChem
include("soilgrid/soilgrid.jl")  # SoilGrid
include("tracers/tracers.jl")  # Tracers

# expose submodule public members to the main module namespace
using SoilTracers.PhysChem
using SoilTracers.SoilGrid

end # module


text:
(pwd() != @__DIR__) && cd(@__DIR__) # allow starting app from bin/ dir

using HelloWorldJl
push!(Base.modules_warned_for, Base.PkgId(HelloWorldJl))
HelloWorldJl.main()

max_stars_repo_name: QuestofIranon/Hello-World-Julia

max_issues_repo_name: QuestofIranon/Hello-World-Julia

max_forks_repo_name: QuestofIranon/Hello-World-Julia


text:
#----------------------------------------------------------------------
# Humanize.jl    https://github.com/IainNZ/Humanize.jl
# Based on jmoiron's humanize Python library (MIT licensed):
#  https://github.com/jmoiron/humanize/
# All original code is (c) Iain Dunning and MIT licensed.
#----------------------------------------------------------------------

isdefined(Base, :__precompile__) && __precompile__()

module Humanize

export datasize, timedelta, digitsep

#---------------------------------------------------------------------
const dec_suf = [" B", " kB", " MB", " GB", " TB", " PB", " EB", " ZB", " YB"]
const bin_suf = [" B", " KiB", " MiB", " GiB", " TiB", " PiB", " EiB", " ZiB", " YiB"]
const gnu_suf = ["", "K", "M", "G", "T", "P", "E", "Z", "Y"]
"""
datasize(value::Number; style=:dec, format="%.1f")

Format a number of bytes in a human-friendly format (eg. 10 kB).

style=:dec - default, decimal suffixes (kB, MB), base 10^3
style=:bin - binary suffixes (KiB, MiB), base 2^10
style=:gnu - GNU-style (ls -sh style, K, M), base 2^10
"""
function datasize(value::Number; style=:dec, format="%.1f")
    suffix  = style == :gnu ? gnu_suf : (style == :bin ? bin_suf : dec_suf)
    base    = style == :dec ? 1000.0 : 1024.0
    bytes   = float(value)
    format  = "$format%s"
    fmt_str = @eval (v,s)->@sprintf($format,v,s)
    unit    = base
    s       = suffix[1]
    for (i,s) in enumerate(suffix)
        unit = base ^ (i)
        bytes < unit && break
    end
    return fmt_str(base * bytes / unit, s)
end

#---------------------------------------------------------------------
"""
timedelta(secs::Integer)
timedelta{T<:Integer}(years::T,months::T,days::T,hours::T,mins::T,secs::T)
timedelta(dt_diff::Dates.Millisecond)
timedelta(d_diff::Dates.Day)

Format a time length in a human friendly format.

timedelta(secs::Integer)
    e.g. 3699   -> 'An hour'
timedelta{T<:Integer}(years::T,months::T,days::T,hours::T,mins::T,secs::T)
    e.g. days=1,hours=4,...     -> 'A day'
         hours=4,mins=2,...     -> '4 hours'
         years=1,months=2,...   -> '1 year, 2 months'
timedelta(dt_diff::Dates.Millisecond)
    e.g. DateTime(2014,2,3) - DateTime(2013,3,7) -> '11 months'
timedelta(d_diff::Dates.Day)
    e.g. Date(2014,3,7) - Date(2013,2,4) -> '1 year, 1 month'
"""
function timedelta(secs::Integer)
    mins   = div(  secs, 60); secs   -= 60*mins
    hours  = div(  mins, 60); mins   -= 60*hours
    days   = div( hours, 24); hours  -= 24*days
    months = div(  days, 30); days   -= 30*months
    years  = div(months, 12); months -= 12*years

    if years == 0
        if days == 0 && months == 0
            hours >= 1 && return hours == 1 ? "an hour"  : "$hours hours"
            mins  >= 1 && return mins  == 1 ? "a minute" : "$mins minutes"
            secs  >= 1 && return secs  == 1 ? "a second" : "$secs seconds"
            return "a moment"
        end
        months == 0 && return days == 1 ? "a day" : "$days days"
        return months == 1 ? "a month" : "$months months"
    elseif years == 1
        months == 0 && return "a year"
        months == 1 && return "1 year, 1 month"
        return "1 year, $months months"
    else
        return "$years years"
    end
end
# Assume nothing about magnitudes of inputs, so cast to seconds first
timedelta{T<:Integer}(years::T,months::T,days::T,hours::T,mins::T,secs::T) =
    timedelta(((((years*12+months)*30+days)*24+hours)*60+mins)*60+secs)
timedelta(dt_diff::Dates.Millisecond) = timedelta(div(Int(dt_diff),1000))
timedelta(d_diff::Dates.Day) = timedelta(Int(d_diff)*24*3600)

#---------------------------------------------------------------------
"""
digitsep(value::Integer, sep = ",", k = 3)

Convert an integer to a string, separating each 'k' digits by 'sep'.  'k'
defaults to 3, separating by thousands.  The default "," for 'sep' matches the
commonly used digit separator in the US.

digitsep(value::Integer)
    e.g. 12345678 -> "12,345,678"
digitsep(value::Integer, sep = "'")
    e.g. 12345678 -> "12'345'678"
digitsep(value::Integer, sep = "'", k = 4)
    e.g. 12345678 -> "1234'5678"
"""
function digitsep(value::Integer, sep = ",", k = 3)
    value = string(value)
    n = length(value)
    starts = reverse(collect(n:-k:1))
    groups = [value[max(x-k+1, 1):x] for x in starts]
    return join(groups, sep)
end

end


text:
using KernelFunctions
using Test

# minimal definition of a custom kernel
struct MyKernel <: Kernel end

KernelFunctions.kappa(::MyKernel, d2::Real) = exp(-d2)
KernelFunctions.metric(::MyKernel) = SqEuclidean()

@test kappa(MyKernel(), 3) == kappa(SqExponentialKernel(), 3)
@test kappa(MyKernel(), 1, 3) == kappa(SqExponentialKernel(), 1, 3)
@test kappa(MyKernel(), [1, 2], [3, 4]) == kappa(SqExponentialKernel(), [1, 2], [3, 4])
@test kernelmatrix(MyKernel(), [1 2; 3 4], [5 6; 7 8]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4], [5 6; 7 8])
@test kernelmatrix(MyKernel(), [1 2; 3 4]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4])

# some syntactic sugar
(κ::MyKernel)(d::Real) = kappa(κ, d)
(κ::MyKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
(κ::MyKernel)(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}; obsdim = 2) = kernelmatrix(κ, X, Y; obsdim = obsdim)
(κ::MyKernel)(X::AbstractMatrix{<:Real}; obsdim = 2) = kernelmatrix(κ, X; obsdim = obsdim)

@test MyKernel()(3) == SqExponentialKernel()(3)
@test MyKernel()([1, 2], [3, 4]) == SqExponentialKernel()([1, 2], [3, 4])
@test MyKernel()([1 2; 3 4], [5 6; 7 8]) == SqExponentialKernel()([1 2; 3 4], [5 6; 7 8])
@test MyKernel()([1 2; 3 4]) == SqExponentialKernel()([1 2; 3 4])


text:
__precompile__(true)

module KernelDensityEstimate

using Gadfly, Colors, Cairo, Fontconfig

import Base: promote_rule, *, rand

export
    kde!,
    getPoints,
    getBW,
    root,
    Npts,
    Ndim,
    getWeights,
    marginal,
    sample,
    rand,
    resample,
    evaluateDualTree,
    BallTree,
    BallTreeDensity,
    getKDERange,
    getKDEMax,
    getKDEMean,
    # approximate intersection volume
    intersIntgAppxIS,

    # product operation with Multiscale Gibbs sampling
    prodAppxMSGibbsS,

    # Gadfly plotting functions
    plotKDE,
    stackMarginals,
    vstackedPlots,
    drawHorDens,
    toggleYTicks,

    # add * operator for kde product approximate
    *

include("BallTree01.jl")
include("BallTreeDensity01.jl")
include("DualTree01.jl")
include("KDE01.jl")
include("KDEPlotting01.jl")
include("MSGibbs01.jl")

"""
A Julia package for Kernel Density Estimation and approximations of their products

"""

end


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

if length(ARGS) == 0
  results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-7-4"
else
  results_dir = ARGS[1]  # path to results directory
end

"""
Returns a list of monitors, where each element is the result from
one temperature.
"""
function getsamples(path_to_output)
  output = BSON.load(path_to_output)
  samples = output[:samples]
  num_tempers = length(output[:tempers])
  num_monitors = length(samples)

  return [let
     monitor = [let
                  num_samples = length(output[:samples][m])
                  [output[:samples][m][i][t] for i in 1:num_samples]
                end for m in 1:num_monitors]
   end for t in 1:num_tempers]
end

using Distributed
rmprocs(filter(w -> w > 1, workers()))
addprocs(20)

@everywhere include("../../../PlotUtils/PlotUtils.jl")
@everywhere include("../../../PlotUtils/imports.jl")

@everywhere function plot_params(samples, simdat, imgdir)
  # Create a directory for images / txt if needed.
  mkpath("$(imgdir)/txt/")

  # Extract samples
  extract(s) = map(o -> o[s], samples[1])
  Zs = extract(:theta__Z)
  Ws = extract(:theta__W)
  lams = extract(:theta__lam)
  alphas = extract(:theta__alpha)
  vs = extract(:theta__v)
  deltas = extract(:theta__delta)
  sig2s = extract(:theta__sig2)
  etas = extract(:theta__eta)

  # Number of samples
  I = length(sig2s[1])

  # Print R
  computeR(w) = vec(sum(w .> 0, dims=2))
  open("$(imgdir)/txt/Rcounts.txt", "w") do io
    I, K = size(Ws[1])
    Rs = map(computeR, Ws)
    Rcounts = zeros(Int, I, K)
    for R in Rs
      for i in 1:I
        Rcounts[i, R[i]] += 1
      end
    end
    writedlm(io, Rcounts)
  end

  # Plot parameters
  PlotUtils.plot_W(Ws, imgdir=imgdir, W_true=simdat[:W])
  PlotUtils.plot_v(vs, imgdir)
  PlotUtils.plot_alpha(alphas, imgdir)
  PlotUtils.plot_mus(deltas, imgdir)
  PlotUtils.plot_sig2(sig2s, imgdir, sig2_true=simdat[:sig2])

  # Plot y / Z
  PlotUtils.make_yz(simdat[:y], Zs, Ws, lams, imgdir, vlim=(-4,4),
                    Z_true=simdat[:Z], w_thresh=0.0)
end

function makeplots(path_to_output, imgdir)
  mkpath("$(imgdir)/txt")

  output_dir = splitdir(path_to_output)[1]
  path_to_simdat = joinpath(output_dir, "simdat.bson")
  simdat = BSON.load(path_to_simdat)[:simdat]

  println(path_to_output)
  output = BSON.load(path_to_output)
  samples = output[:samples][1]

  # Swap props
  println("Making swap proprs ...")
  if :swapcounts in keys(output) && :paircounts in keys(output)
    swapcounts = output[:swapcounts]
    paircounts = output[:paircounts]
    tempers = output[:tempers]
    ntempers = length(tempers)
    plt.imshow(swapcounts ./ (paircounts .+ 1e-6))
    plt.colorbar()
    plt.xticks(0:(ntempers-1), round.(tempers, digits=2), rotation=45)
    plt.yticks(0:(ntempers-1), round.(tempers, digits=2))
    plt.savefig("$(imgdir)/swapprops.pdf", bbox_inches="tight")
    plt.close()
  end

  # Number of burn in iterations
  nburn = output[:nburn]

  # Plot loglikelihood
  println("Loglikes ...")
  llkey = if :loglike in keys(output)
    PlotUtils.plot_loglike(output[:loglike], imgdir, fname="loglike.pdf")
    PlotUtils.plot_loglike(output[:loglike][(nburn + 1):end],
                           imgdir, fname="loglike_postburn.pdf")
  elseif :lls in keys(output)
    lls = output[:lls]
    for t in 1:ntempers
      temper = tempers[t]
      title = "temperature: $(temper)"
      ll_path = joinpath(imgdir, "loglike")
      mkpath(ll_path)
      PlotUtils.plot_loglike(lls[t], ll_path, fname="loglike_$(t).pdf",
                             title=title)
      PlotUtils.plot_loglike(lls[t][(nburn + 1):end],
                             ll_path, fname="loglike_postburn_$(t).pdf",
                             title=title)
    end
  end

  # Extract samples
  println("Extracting samples ...")
  extract(s) = map(o -> o[s], samples[1])
  Zs = extract(:theta__Z)
  Ws = extract(:theta__W)
  deltas = extract(:theta__delta)
  sig2s = extract(:theta__sig2)
  etas = extract(:theta__eta)

  # Number of samples
  I = length(sig2s[1])

  # Plot missing mechanism
  println("Plotting missmechs ...")
  for i in 1:I
    PlotUtils.plot_missmech(output[:c].constants.beta, i, xlim=[-5, 0])
    plt.savefig("$(imgdir)/missmech_$(i).pdf", bbox_inches="tight")
    plt.close()
  end

  # Print missing mechanism
  open(joinpath(imgdir, "txt/beta.txt"), "w") do io
    writedlm(io, output[:c].constants.beta, ',')
  end

  # Plot data density
  println("Plotting dden ...")
  PlotUtils.plot_dden(ddens=output[:dden],
                      etas=etas, Ws=Ws, Zs=Zs, sig2s=sig2s, deltas=deltas,
                      ygrid=output[:c].constants.y_grid,
                      imgdir=imgdir, simdat=simdat)
end

### MAIN ###

# Name of output file
OUTPUT_FILE = "output.bson"

# PATH TO ALL OUTPUT FILES
output_paths = [joinpath(root, OUTPUT_FILE)
                for (root, _, files) in walkdir(results_dir)
                if OUTPUT_FILE in files]

# Plot the posterior distributions of parameters
simname = "maxtemp2-ntempts20-degree2-N500"
path_to_output = joinpath(results_dir, simname, "output.bson")
path_to_simdat = joinpath(results_dir, simname, "simdat.bson")
samples = getsamples(path_to_output);
simdat = BSON.load(path_to_simdat)[:simdat];
@everywhere simdat = $simdat
genimgdir(t) = joinpath(results_dir, simname, "img/temper_$(t)")
@everywhere plot_params(s, imgdir) = plot_params(s, simdat, imgdir)
_ = pmap(plot_params, samples, genimgdir.(1:length(samples)))

# Make general plots
imgdir = joinpath(results_dir, simname, "img")
makeplots(path_to_output, imgdir)

# Remove extra processors
# rmprocs(filter(w -> w > 1, workers()))

println("DONE!")

# Detailed Analysis
output = BSON.load(path_to_output);
last_states = output[:all_last_states];
tempers = output[:tempers]
num_tempers = length(tempers)

ll_fn(s, t) = CytofRepFAM.Model.compute_marg_loglike(s,
                                                     output[:c].constants,
                                                     output[:d].data, t)
i, j = 6, 7
ll_fn(last_states[i].theta, tempers[i]*1.001) -
ll_fn(last_states[i].theta, tempers[i])
ll_fn(last_states[j].theta, tempers[i]) -
ll_fn(last_states[j].theta, tempers[i]*1.001)

log_accept_mat = [if i >= j
                MCMC.WSPT.compute_log_accept_ratio(
                  ll_fn, (last_states[i].theta, last_states[j].theta),
                  (tempers[i], tempers[j]), verbose=0)
              else
                0.0
              end for i in 1:num_tempers, j in 1:num_tempers]
for i in 1:num_tempers, j in 1:i
  log_accept_mat[j, i] = log_accept_mat[i, j]
end

plt.imshow(exp.(log_accept_mat), vmin=.01, vmax=1.0)
plt.colorbar()
plt.savefig(joinpath(results_dir, simname, "img/ll_compare.pdf"),
            bbox_inches="tight")
plt.close()

plt.plot(tempers, log_accept_mat[:, 20], marker="o")
plt.xlabel("temperatures")
plt.ylabel("log accept ratio (t1 vs others)")
plt.savefig(joinpath(results_dir, simname, "img/log_accept_compare.pdf"),
            bbox_inches="tight")
plt.close()

ll1_vary_tempers = [ll_fn(last_states[1].theta, t) for t in tempers]
plt.plot(tempers, ll1_vary_tempers, marker="o")
plt.xlabel("temperatures")
plt.ylabel("log-likelihood")
plt.savefig(joinpath(results_dir, simname, "img/ll1_vary_tempers.pdf"),
            bbox_inches="tight")
plt.close()


text:
# Autogenerated wrapper script for EarCut_jll for armv7l-linux-gnueabihf
export libearcut

JLLWrappers.@generate_wrapper_header("EarCut")
JLLWrappers.@declare_library_product(libearcut, "libearcut.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libearcut,
        "lib/libearcut.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StochasticBlockModel
using Test

@testset "OptResults struct" begin

    @test_throws ArgumentError StochasticBlockModel.OptResults(2.0, 1.0, :Optimal, 2.0, 0, 0, 0) # Check LB, UB
    @test_throws ArgumentError StochasticBlockModel.OptResults(2.0, 4.0, :Optimal, 2.0, 0, 0, 0) # Check optimal status
    #@test_throws DomainError StochasticBlockModel.OptResults(2.0, 2.0, :Optimal, 0.0, 0, 0, 0) # Check solvetime domain
    @test_throws DomainError StochasticBlockModel.OptResults(2.0, 2.0, :Optimal, -2.0, 0, 0, 0) # Check solvetime domain
    @test_throws DomainError StochasticBlockModel.OptResults(2.0, 2.0, :Optimal, 2.0, -1, 0, 0) # Check iterations domain
    @test_throws DomainError StochasticBlockModel.OptResults(2.0, 2.0, :Optimal, 2.0, 0, -1, 0) # Check nodecount domain
    @test_throws DomainError StochasticBlockModel.OptResults(2.0, 2.0, :Optimal, 2.0, 0, 0, -1) # Check nodecount domain
end


text:
function flat1(arr)
    rst = Any[]
    grep(v) = for x in v
        if isa(x, Array) grep(x) else push!(rst, x) end
    end
    grep(arr)
    rst
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
pfile.jl
read GE MRI k-space data from .p file
based on https://gitlab.com/fMRI/toppe/+toppe/+utils/loadpfile.m

(c) 2016 The Regents of the University of Michigan
Jon-Fredrik Nielsen, jfnielse@umich.edu

Comments from Matlab version at
https://gitlab.com/fMRI/toppe/+toppe/+utils/loadpfile.m

This file was derived from part of the TOPPE development environment
for platform-independent MR pulse programming.

TOPPE is free software: you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as published by
the Free Software Foundation version 2.0 of the License.

TOPPE is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Library General Public License for more details.

You should have received a copy of the GNU Library General Public License
along with TOPPE.
If not, see <http://www.gnu.org/licenses/old-licenses/lgpl-2.0.html>.

2019-05-22 Julia version by Jeff Fessler
=#

export loadpfile

# using MIRTio: read_rdb_hdr, header_init, header_size, header_string, header_write

"""
    (dat, rdb_hdr) = loadpfile(fid::IOStream ; ...)

Load data for one or more echoes from GE MRI scan Pfile.

in
- `fid::IOStream`

option
- `coils::AbstractVector{Int}`
 only get data for these coils; default: all coils
- `echoes::AbstractVector{Int}`
 only get data for these echoes; default: all echoes
- `slices::AbstractVector{Int}`
 only get data for these slices; default: `2:nslices` (NB!)
 because first slice (dabslice=0 slot) may contain corrupt data.
- `views::AbstractVector{Int}`
 only get data for these views; default: all views
- `quiet::Bool` non-verbosity, default `false`

out
- `dat::Array{Complex{Int16}}` `[ndat, ncoil, nslice, necho, nview]`
- `rdb_hdr::NamedTuple` header information

To save memory the output type is complex-valued Int16.
"""
function loadpfile(
    fid::IOStream;
    filesize::Int = 0,
    coils::AbstractVector{Int} = Int[],
    echoes::AbstractVector{Int} = Int[],
    slices::AbstractVector{Int} = Int[],
    views::AbstractVector{Int} = Int[],
    quiet::Bool = false,
)

    rdb_hdr = read_rdb_hdr(fid) # read header

    # Header parameters
    ndat = Int(rdb_hdr.frame_size)
    nslices = rdb_hdr.nslices
    ptsize = rdb_hdr.point_size # 2: data stored in short int format, int16)
    ptsize != 2 && # 4 (extended precision) unsupported here
        throw("ptsize = $ptsize unsupported")
    nechoes = rdb_hdr.nechoes
    nviews = rdb_hdr.nframes
    ncoils = rdb_hdr.dab[2] - rdb_hdr.dab[1] + 1

    # Calculate size of data chunks
    # See pfilestruct.jpg, and rhrawsize calculation in .e file.
    # number of data points per 'echo' loaddab slot. Includes baseline (0) view:
    # Apparent file data order is: [ndat, nview, necho, nslice, ncoil]
    echores = ndat * (nviews + 1)
    sliceres = nechoes * echores # number of data points per 'slice'
    coilres = nslices * sliceres # number of data points per receive coil

    # This size should match the Pfile size exactly:
    pfilesize =
        rdb_hdr.off_data + 2 * ptsize * ncoils * nslices * nechoes * (nviews + 1) * ndat

    # File size check
    (filesize != 0) &&
        (pfilesize != filesize) &&
        throw("Expected $(pfilesize/1e6) MB but see $(filesize/1e6) MB file")

    # Check if second view is empty
    # This happens when there's only one view, but the scanner sets nviews = 2
    if nviews == 2
        # Seek to view slice 2, echo 1, view 2, coil 1:
        seek(fid, rdb_hdr.off_data + (2 * ptsize * (sliceres + 2 * ndat)))
        view2tmp = Array{Int16}(undef, 2 * ndat)
        view2tmp = read!(fid, view2tmp)
        if all(view2tmp .== Int16(0))
            nviews = 1 # set nviews to 1 so we don't read in all the empty data
            !quiet && @warn("View 2 appears to be empty, only loading view 1.")
        end
    end

    # Determine which portions to load
    # NB: default is to skip first slice (sometimes contains corrupted data)
    isempty(coils) && (coils = 1:ncoils)
    isempty(echoes) && (echoes = 1:nechoes)
    isempty(slices) && (slices = 2:nslices)
    isempty(views) && (views = 1:nviews)

    # output dimensions:
    dims = (ndat, length(coils), length(slices), length(echoes), length(views))

    # Warn if output size will exceed system memory.
    # Exceeding free RAM may cause a system freeze due to attempted mem swap
    memneeded = 4 * prod(dims) # 4 bytes per complex value because Int16
    memfree = Sys.total_memory() # system total memory
    mempercentuse = 100 * memneeded / memfree
    (mempercentuse > 100) &&
        @warn("Loading data ($(memneeded/1e9) GB) may exceed available RAM and freeze system")
    # fprintf('Press enter to continue anyway...');
    # input('');
    (mempercentuse > 90) && # Warn if we will we use 90% of memory
        @warn("Loading data ($(memneeded/1e9) GB) will use $mempercentuse % of your available RAM. Proceed with caution!")

    !quiet && (@info("ndat = $ndat, memory = $(memneeded/1e9) GB"))
    !quiet && (@info("dims = $dims"))
    !quiet && (@info("slices $(slices[1]):$(slices[end]) / $nslices"))
    !quiet && (@info("echoes $(echoes[1]):$(echoes[end]) / $nechoes"))
    !quiet && (@info("views $(views[1]):$(views[end]) / $nviews"))
    !quiet && (@info("coils $(coils[1]):$(coils[end]) / $ncoils"))

    # Read data from file
    # Julia stores complex data as real1,imag1, real2,imag2, ... unlike matlab!
    # This simplifies the IO here!
    data = Array{Complex{Int16}}(undef, dims)
    dtmp = Array{Complex{Int16}}(undef, ndat) # one readout
    # !quiet && textprogressbar('Loading data: ')
    for icoil in coils
        !quiet && (@show icoil)
        # !quiet && textprogressbar(icoil/ncoils*100)
        for islice in slices
            for iecho in echoes
                for iview in views
                    # !quiet && (@info "icoil=$icoil islice=$islice iecho=$iecho iview=$iview")
                    offsetres =
                        (icoil - 1) * coilres +
                        (islice - 1) * sliceres +
                        (iecho - 1) * echores +
                        iview * ndat
                    offsetbytes = 2 * ptsize * offsetres
                    seek(fid, rdb_hdr.off_data + offsetbytes)
                    read!(fid, dtmp)
                    data[
                        :,
                        icoil-coils[1]+1,
                        islice-slices[1]+1,
                        iecho-echoes[1]+1,
                        iview-views[1]+1,
                    ] .= dtmp
                end
            end
        end
    end

    # !quiet && textprogressbar(' done.')

    return (data, rdb_hdr)
end

"""
    (dat, rdb_hdr) = loadpfile(file::String ; ...)

Load from `file`.
"""
function loadpfile(file::String; kwargs...)
    open(file, "r") do fid
        return loadpfile(fid; filesize = filesize(file), kwargs...)
    end
end

"""
    (dat, rdb_hdr) = loadpfile(file, echo::Integer ; ...)

Load a single echo from `file`.
"""
function loadpfile(pfile::String, echo::Integer; kwarg...)
    return loadpfile(pfile, echoes = [echo]; kwarg...)
end


text:
# TODO: move this next line to MLJBase:
MLJBase.iterator(r::NominalRange, ::Nothing) = iterator(r)

"""
    MLJTuning.grid([rng, ] prototype, ranges, resolutions)

Given an iterable `ranges` of `ParamRange` objects, and an iterable
`resolutions` of the same length, return a vector of models generated
by cloning and mutating the hyperparameters (fields) of `prototype`,
according to the Cartesian grid defined by the specifed
one-dimensional `ranges` (`ParamRange` objects) and specified
`resolutions`. A resolution of `nothing` for a `NominalRange`
indicates that all values should be used.

Specification of an `AbstractRNG` object `rng` implies shuffling of
the results. Otherwise models are ordered, with the first
hyperparameter referenced cycling fastest.

"""
grid(rng::AbstractRNG, prototype::Model, ranges, resolutions) =
    shuffle(rng, grid(prototype, ranges, resolutions))

function grid(prototype::Model, ranges, resolutions)

    iterators = broadcast(iterator, ranges, resolutions)

    A = MLJBase.unwind(iterators...)

    N = size(A, 1)
    map(1:N) do i
        clone = deepcopy(prototype)
        for k in eachindex(ranges)
            field = ranges[k].field
            recursive_setproperty!(clone, field, A[i,k])
        end
        clone
    end
end

"""
    process_user_range(user_specified_range, resolution, verbosity)

Utility to convert user-specified range (see [`Grid`](@ref)) into a
pair of tuples `(ranges, resolutions)`.

For example, if `r1`, `r2` are `NumericRange`s and `s` is a
NominalRange` with 5 values, then we have:

    julia> MLJTuning.process_user_range([(r1, 3), r2, s], 42, 1) ==
                            ((r1, r2, s), (3, 42, 5))
    true

If `verbosity` > 0, then a warning is issued if a `Nominal` range is
paired with a resolution.

"""
process_user_range(user_specified_range, resolution, verbosity) =
    process_user_range([user_specified_range, ], resolution, verbosity)
function process_user_range(user_specified_range::AbstractVector,
                    resolution, verbosity)
    stand(r) = throw(ArgumentError("Unsupported range. "))
    stand(r::NumericRange) = (r, resolution)
    stand(r::NominalRange) = (r, length(r.values))
    stand(t::Tuple{NumericRange,Integer}) = t
    function stand(t::Tuple{NominalRange,Integer})
        verbosity < 0 ||
            @warn  "Ignoring a resolution specified for a `NominalRange`. "
        return (first(t), length(first(t).values))
    end

    ret = zip(stand.(user_specified_range)...) |> collect
    return first(ret), last(ret)
end


text:
using ConvexBodyProximityQueries: proj
@testset "projection operation" begin
    a⃗ = SVector{2}([1.0, 2.0])
    b⃗ = SVector{2}([-1.0, 1.0])
    @test isapprox(proj(a⃗,b⃗), [0.2, 0.4], rtol=1e-3)
    @test isapprox(proj(b⃗,a⃗), [-0.5, 0.5], rtol=1e-3)
    a⃗ = SVector{3}([1.0, 2.0, 2.0])
    b⃗ = SVector{3}([-1.0, 1.0, -2.0])
    @test isapprox(proj(a⃗,b⃗), [-1.0, -2.0, -2.0]/3, rtol=1e-3)
    @test isapprox(proj(b⃗,a⃗), [0.5, -0.5, 1.0], rtol=1e-3)
end # projection operation

#  TODO: fix tests
#  using GJK: findsimplex
#  @testset "find simplex" begin
#      @testset "find line simplex" begin
#          spx = SMatrix{2,2}([1.0 -1.0; 2.0 1.0])
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [1, 2])
#          @test all(ret[2] .≈ [0.6, -1.2])
#          @test ret[3] == false
#          spx = SMatrix{2,2}([1.0 0.5; 2.0 1.5])
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2])
#          @test all(ret[2] .≈ [-0.5, -1.5])
#          @test ret[3] == false
#
#          spx = SMatrix{3,2}([1.0 -1.0; 2.0 1.0; 2.0 -2.0])
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [1, 2])
#          @test all(ret[2] .≈ [0.1428571428571429, -1.4285714285714286, 0.2857142857142858])
#          @test ret[3] == false
#          spx = SMatrix{3,2}([1.0 0.0; 2.0 1.5; 2.0 1.0])
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2])
#          @test all(ret[2] .≈ [0.0, -1.5, -1.0])
#          @test ret[3] == false
#
#          @testset "degeneracies" begin
#              rvec = rand(2, 1) .+ 1.0
#              spx = [rvec -rvec]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#              rvec = rand(2, 1) .+ 1.0
#              spx = [rvec zero(rvec)]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#              rvec = rand(2, 1) .+ 1.0
#              spx = [zero(rvec) rvec]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#
#              rvec = rand(3, 1) .+ 1.0
#              spx = [rvec -rvec]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#              rvec = rand(3, 1) .+ 1.0
#              spx = [rvec zero(rvec)]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#              rvec = rand(3, 1) .+ 1.0
#              spx = [zero(rvec) rvec]
#              ret = findsimplex(spx)
#              @test ret[3] .== true
#          end # degneracies
#      end # find line simplex
#
#      @testset "find triangle simplex" begin
#          spx = [1.0 -1.0 0.0; 2.0 1.0 -1.0]
#          ret = findsimplex(spx)
#          @test ret[3] == true
#          spx = [0.0 1.0 -1.0; 3.0 2.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2, 3])
#          @test all(ret[2] .≈ [0.6, -1.2])
#          @test ret[3] == false
#          spx = [1.0 0.0 -1.0; 2.0 3.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [1, 3])
#          @test all(ret[2] .≈ [0.6, -1.2])
#          @test ret[3] == false
#          spx = [1.0 -1.0 0.0; 2.0 1.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [3])
#          @test all(ret[2] .≈ [0.0, -1.0])
#          @test ret[3] == false
#
#          spx = [1.0 -1.0 0.0; 2.0 1.0 -1.0; 2.0 1.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2, 1, 3])
#          @test all(ret[2] .≈ [0.4, 0.2, -1.0])
#          @test ret[3] == false
#          spx = [0.0 1.0 -1.0; 3.0 2.0 1.0; 2.0 1.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2, 3])
#          @test all(ret[2] .≈ [0.6, -1.2, -1.0])
#          @test ret[3] == false
#          spx = [1.0 -1.0 0.0; 2.0 1.0 1.0; 2.0 1.0 0.5]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [3])
#          @test all(ret[2] .≈ [0.0, -1.0, -0.5])
#          @test ret[3] == false
#
#          @testset "degeneracies" begin
#              x = rand() + 1.0
#              y = rand() + 1.0
#              z = rand() + 1.0
#              spx = [x 0 -x; 0 y -y; z z -2z]
#              ret = findsimplex(spx)
#              @test ret[3] == true
#              spx = [x 0 0; 0 y 0; z -z 0]
#              ret = findsimplex(spx)
#              @test ret[3] == true
#
#              spx = [1.0 -1.0 0.0; -0.5 0.5 0.0]
#              ret = findsimplex(spx)
#              @test ret[3] == true
#              spx = [1.0 1.0 1.0; -3.0 2.0 -0.5]
#              ret = findsimplex(spx)
#              @test all(ret[1] .== [2, 3])
#              @test all(ret[2] .≈ [-1.0, 0.0])
#              @test ret[3] == false
#              spx = [1.0 1.0 1.0; -3.0 -2.0 -0.5]
#              ret = findsimplex(spx)
#              @test all(ret[1] .== [3])
#              @test all(ret[2] .≈ [-1.0, 0.5])
#              @test ret[3] == false
#
#              spx = [1.0 -1.0 0.0; -0.5 0.5 0.0; 1.0 -1.0 0.0]
#              ret = findsimplex(spx)
#              @test ret[3] == true
#              spx = [1.0 0.0 0.4; 0.0 1.0 0.6; 1.0 -1.0 -0.2]
#              ret = findsimplex(spx)
#              @test all(ret[1] .== [3])
#              @test all(ret[2] .≈ [-0.4, -0.6, 0.2])
#              @test ret[3] == false
#              spx = [1.0 0.0 0.5; 0.0 1.0 0.5; 1.0 0.0 0.5]
#              ret = findsimplex(spx)
#              @test all(ret[1] .== [2, 3])
#              @test all(ret[2] .≈ [-1.0, -2.0, -1.0]/3.0)
#              @test ret[3] == false
#          end # degneracies
#      end # find triangle simplex
#
#      @testset "find tetrahedron simplex" begin
#          spx = [1.0 0.0 -1.0 0.0; 2.0 -1.0 1.0 0.0; 2.0 1.0 1.0 -2.0]
#          ret = findsimplex(spx)
#          @test ret[3] == true
#          spx = [1.0 0.0 -1.0 0.0; 2.0 -1.0 1.0 0.0; 2.0 1.0 1.0 1.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [2, 3, 4])
#          @test all(ret[2] .≈ [0.0, 0.0, -1.0])
#          @test ret[3] == false
#          spx = [-2.0 -2.0 2.0 1.0; -1.0 1.0 0.0 0.0; 1.0 1.0 1.0 0.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [1, 2, 4])
#          @test all(ret[2] .≈ [-0.1, 0.0, -0.3])
#          @test ret[3] == false
#          spx = [-2.0 2.0 -2.0 1.0; -1.0 0.0 1.0 0.0; -1.0 -1.0 -1.0 0.0]
#          ret = findsimplex(spx)
#          @test all(ret[1] .== [3, 1, 4])
#          @test all(ret[2] .≈ [-0.1, 0.0, 0.3])
#          @test ret[3] == false
#
#          @testset "degeneracies" begin
#              x = rand() + 1.0
#              y = rand() + 1.0
#              z = rand() + 1.0
#              spx = [x 0 -x 0; 0 y -y 0; z z -z 0]
#              ret = findsimplex(spx)
#              @test ret[3] == true
#          end # degneracies
#      end # find tetrahedron simplex
#  end # find simplex

#  using GJK: findcombination
#  @testset "find combination" begin
#      @testset "find line combination" begin
#          spx = [1.0 1.0; -3.0 2.0]
#          vec = [1.0, 0.5]
#          @test all(findcombination(spx, vec) .≈ [0.3, 0.7])
#          vec = [1.0, 2.5]
#          @test all(findcombination(spx, vec) .≈ [0.0, 1.0])
#      end # find line combination
#
#      @testset "find triangle combination" begin
#          spx = [1.0 1.0 -1.0; 1.0 2.0 1.0; 2.0 2.0 1.0]
#          vec = [1.0, 4.0, 5.0]/3
#          @test all(findcombination(spx, vec) .≈ ones(3)/3.0)
#          vec = [-3.0, 0.0, 0.0]
#          @test all(findcombination(spx, vec) .≈ [0.0, 0.0, 1.0])
#          vec = [0.5, 0.0, 0.0]
#          @test all(findcombination(spx, vec) .≈ [0.4, 0.0, 0.6])
#          vec = [0.0, 2.0, 1.0]
#          @test all(findcombination(spx, vec) .≈ [0.0, 0.5, 0.5])
#      end # find triangle combination
#  end # find combination


text:
mutable struct Frequency
    type::DataType
    freq::Int
end

Base.isless(a::Frequency, b::Frequency) = isless(a.freq, b.freq)

function increment!(c::Frequency)
    c.freq += 1
end

struct DataTypeFrequencies
    d::Dict{Int, Frequency}
    DataTypeFrequencies() = new(Dict())
end

function increment!(tf::DataTypeFrequencies, t::DataType)
    f = get!(() -> Frequency(t, 0), tf.d, t.hash)
    increment!(f)
end

Base.values(tf::DataTypeFrequencies) = values(tf.d)
Base.empty!(tf::DataTypeFrequencies) = empty!(tf.d)


text:
function simulate()#(DTC, SR)
    global net, pᵢ, CFᵢ, CFₒ, states

    # initialization
    pᵢ = zeros(nlinks, T, nsinks, nclasses); # inflow proportions
    CFᵢ = zeros(nlinks, T); # cumulative inflows
    CFₒ = zeros(nlinks, T); # cumulative outflows
    for (rno,r) in enumerate(sources(net))
        tmp = @view DTC[:,rno,..]
        pᵢ[rno,..] .= tmp./max.(Ref(1e-16), sum(tmp, dims=(2,3)))
        CFᵢ[r,:] .= round.(cumsum(squeezesum(tmp .* expand(trips[rno,..], dims=1), dims=(2,3))), digits=10)
    end

    rlinkids = [outlinkids(net, r)[1] for r in srcs]
    states = Matrix{Union{Nothing, LinkState}}(fill(nothing, nlinks, T+1))
    for i in 1:nlinks
        l = length(link(net, i))
        for j in 1:l
            states[i,j] = LinkState(j-1-l, 0., (i in rlinkids) ? Inf : Q)
        end
    end

    for t in 1:T
        tracs = zeros(nlinks)

        for k in snks
            li = inlinkids(net,k)[1]
            l = length(link(net, li))

            CFₒ[li,t] = (t > l) ? CFᵢ[li,t-l] : 0.
            tracs[li] = tracker(states[li,t]) + 1
        end

        for m in mrgs
            #@show t, m
            # sending and receiving capacities
            ili = inlinkids(net,m)
            oli = outlinkids(net,m)[1]
            sval = svalue.(states[ili,t])
            rval = rvalue(states[oli,t])

            ls = length.(link.(Ref(net), ili))

            # cuminflows, props, sval, rval, tracs, maxt
            newcf, p, newtracs = mflows([@view(CFᵢ[ili[1],:]), @view(CFᵢ[ili[2],:])], [@view(pᵢ[ili[1],..]), @view(pᵢ[ili[2],..])], sval, rval, tracker.(states[ili,t]), t.-ls)

            CFᵢ[oli,t] = round(sum(newcf), digits=ROUND_DIGITS)
            CFₒ[ili,t] .= round.(newcf, digits=ROUND_DIGITS)
            pᵢ[oli,t,..] .= p
            tracs[ili] .= newtracs
        end

        for d in divs
            #@show t, d
            # sending and receiving capacities
            ili = inlinkids(net,d)[1]
            oli = outlinkids(net,d)
            sval = svalue(states[ili,t])
            rval = rvalue.(states[oli,t])

            l = length(link(net, ili))

            # cuminflows, props, sr, sval, rval, trac, maxt
            if t > l
                sr = [SR[d][ti][t,..] for ti in 1:2]
                newcf, p, newtrac, newtf = dflows((@view(CFᵢ[ili,:])), @view(pᵢ[ili,..]), sr, sval, rval, tracker(states[ili,t]), t-l)

                CFᵢ[oli,t] .= round.(newtf .+ CFᵢ[oli,t-1], digits=ROUND_DIGITS)
                CFₒ[ili,t] = round(newcf, digits=ROUND_DIGITS)
                for (i,li) in enumerate(oli)
                    pᵢ[li,t,..] .= p[i]
                end
                tracs[ili] = newtrac
            end
        end

        for (li,lnk) in enumerate(net.links)
            l = length(lnk)
            if t >= l
                sval = round(min(Q, CFᵢ[li, t+1-l] - CFₒ[li,t]), digits=ROUND_DIGITS)
                rval = if li in rlinkids
                    Inf
                else
                    round(min(Q, ((t >= Int(round(l/δ))) ? CFₒ[li,t+1-Int(round(l/δ))] : 0.) - CFᵢ[li,t] + N*l), digits=ROUND_DIGITS)
                end

                sval = max.(Ref(0.), sval)
                rval = max(0., rval)

                states[li,t+1] = LinkState(tracs[li], sval, rval)
            end
        end
    end

    return CFᵢ, CFₒ, states
end


text:
module Bijective
import ..BijectiveHilbert: encode_hilbert_zero, decode_hilbert_zero!, HilbertAlgorithm, axis_type

struct Simple2D{T} <: HilbertAlgorithm{T}
end

Simple2D(::Type{T}) where {T} = Simple2D{T}()

axis_type(::Simple2D{T}) where {T} = Int

function encode_hilbert_zero(::Simple2D{T}, X::Vector{A})::T where {A, T}
    x = X[1]
    y = X[2]
    z = zero(T)
    if x == zero(A) && y == zero(A)
        return z
    end
    rmin = convert(Int, floor(log2(max(x, y))) + 1)
    w = one(A) << (rmin - 1)
    while rmin > 0
        z <<= 2
        if rmin&1 == 1  # odd
            if x < w
                if y >= w
                    # 1
                    x, y = (y - w, x)
                    z += one(T)
                end  # else x, y remain the same.
            else
                if y < w
                    x, y = ((w<<1) - x - one(w), w - y - one(w))
                    # 3
                    z += T(3)
                else
                    x, y = (y - w, x - w)
                    z += T(2)
                    # 2
                end
            end
        else  # even
            if x < w
                if y >= w
                    # Quadrant 3
                    x, y = (w - x - one(w), (w << 1) - y - one(w))
                    z += T(3)
                end  # else do nothing for quadrant 0.
            else
                if y < w
                    # 1
                    x, y = (y, x-w)
                    z += one(T)
                else
                    # 2
                    x, y = (y-w, x-w)
                    z += T(2)
                end
            end
        end
        rmin -= 1
        w >>= 1
    end
    z
end

function decode_hilbert_zero!(::Simple2D{T}, X::Vector{A}, z::T) where {A,T}
    r = z & T(3)
    x, y = [(zero(A), zero(A)), (zero(A), one(A)), (one(A), one(A)), (one(A), zero(A))][r + 1]
    z >>= 2
    rmin = 2
    w = one(A) << 1
    while z > zero(T)
        r = z & T(3)
        parity = 2 - rmin&1
        if rmin & 1 != 0
            # Nothing to do for quadrant 0.
            if r == 1
                x, y = (y, x + w)
            elseif r == 2
                x, y = (y + w, x + w)
            elseif r == 3
                x, y = ((w << 1) - x - one(A), w - y - one(A))
            end
        else
            if r == 1
                x, y = (y + w, x)
            elseif r == 2
                x, y = (y + w, x + w)
            elseif r == 3
                x, y = (w - x - one(A), (w << 1) - y - one(A))
            end
        end
        z >>= 2
        rmin += 1
        w <<= 1
    end
    X[1] = x
    X[2] = y
end

end


text:
module AdvancedProbabilisticCircuits

using MLStyle
using StatsFuns
using Crayons
using Crayons.Box
using SpecialFunctions
using LogExpFunctions

import Base.show

export Prod, Sum, Normal, VISum, TruncatedNormal, Indicator
export Node, SumNode, VISumNode, ProductNode, NodeType
export RealInterval
export Decomposability, Smoothness, Determinism, Invertability, NodeProperties
export support, scope, loglikelihood, params, descendants
export partitionfunction, logpdf, children, leaves

abstract type NodeType end
struct ProductNode <: NodeType end
struct SumNode <: NodeType end
struct VISumNode <: NodeType end

printnodetype(t::Type{<:NodeType}) = string(t)
printnodetype(::Type{SumNode}) = "(+)"
printnodetype(::Type{VISumNode}) = "(vi+)"
printnodetype(::Type{ProductNode}) = "(×)"

abstract type NodeProperties end
struct Decomposability <: NodeProperties end
struct Smoothness <: NodeProperties end
struct Determinism <: NodeProperties end
struct Invertability <: NodeProperties end

abstract type AbstractNode end
abstract type AbstractLeaf <: AbstractNode end

struct RealInterval{T<:Real}
    lb::T
    ub::T
end
RealInterval(lb::T1, ub::T2) where {T1<:Real,T2<:Real} = RealInterval(lb, T1(ub))

struct Node{V<:Union{Vector{Int},Int},S<:Union{Vector,RealInterval},T<:AbstractVector{<:Real},P<:Tuple,N<:NodeType} <: AbstractNode
    scope::V
    support::S
    children::Vector{AbstractNode} # this might lead to inefficiencies
    params::T
    properties::P
end

children(n::Node) = n.children
scope(n::T) where {T<:AbstractNode} = n.scope
params(n::T) where {T<:AbstractNode} = n.params

partitionfunction(n::T) where {T<:AbstractLeaf} = 1.0
partitionfunction(n::N) where {N<:Node} = logpdf(n, missing)

support(n::Node{Int,S,T,P,N}) where {T,S,P,N} = (n.scope => n.support)
support(n::Node{Vector{Int},S,T,P,N}) where {T,S,P,N} = [s => sup for (s, sup) in zip(n.scope, n.support)]

logpdf(n::T, x::Missing) where {T<:AbstractLeaf} = 1.0
logpdf(n::T, x::AbstractVector{<:Real}) where {T<:AbstractLeaf} = logpdf(n, x[n.scope])
logpdf(n::T, x::AbstractMatrix{<:Real}) where {T<:AbstractLeaf} = @inbounds @views logpdf.(Ref(n), x[:,n.scope])

loglikelihood(n::T, x::X) where {T<:AbstractNode, X<:Real} = logpdf(n, x) - partitionfunction(n)
loglikelihood(n::T, x::X) where {T<:AbstractNode, X<:AbstractArray} = logpdf(n, x) .- partitionfunction(n)

# includes
include("macros.jl")
include("nodes.jl")
include("properties.jl")
include("optimize.jl")
include("util.jl")
include("print.jl")
include("math.jl")

end


text:
@doc raw"""
    asymptote_export_S2_signals(filename; points, curves, tVectors, colors, options...)
Export given `points`, `curves`, and `tVectors` on the sphere $\mathbb S^2$
to Asymptote.

# Input
* `filename` – a file to store the Asymptote code in.

# Optional Arguments (Data)
* `colors` - dictionary of color arrays (indexed by symbols `:points`, `:curves`
  and `:tvector`) where each entry has to provide as least as many colors as
  the length of the corresponding sets.
* `curves` – an `Array` of `Arrays` of points on the sphere, where each inner array is
  interpreted as a curve and is accompanied by an entry within `colors`
* `points` – an `Array` of `Arrays` of points on the sphere where each inner array is
  itnerpreted as a set of points and is accompanied by an entry within `colors`
* `tVectors` – an `Array` of `Arrays` of tuples, where the first is a points, the second a
  tangent vector and each set of vectors is accompanied by an entry
  from within `colors`

# Optional Arguments (Asymptote)
* `arrowHeadSize` - (`6.0`) size of the arrowheads of the tangent vectors
* `cameraPosition` - (`(1., 1., 0.)`) position of the camera in the Asymptote
  szene
* `lineWidth` – (`1.0`) size of the lines used to draw the curves.
* `lineWidths` – overrides the previous value to specify a value per curve and tVector set.
* `dotSize` – (`1.0`) size of the dots used to draw the points.
* `dotSizes` – overrides the previous value to specify a value per point set.
* `sphereColor` – (`RGBA{Float64}(0.85, 0.85, 0.85, 0.6)`) color of the sphere
  the data is drawn on
* `sphereLineColor` –  (`RGBA{Float64}(0.75, 0.75, 0.75, 0.6)`) color of the lines on the sphere
* `sphereLineWidth` – (`0.5`) line width of the lines on the sphere
* `target` – (`(0.,0.,0.)`) position the camera points at
"""
function asymptote_export_S2_signals(filename::String;
    points::Array{Array{T,1},1} where T = Array{Array{Float64,1},1}(undef,0),
    curves::Array{Array{T,1},1} where T = Array{Array{Float64,1},1}(undef,0),
    tVectors::Array{Array{Tuple{T,T},1},1} where T = Array{Array{Tuple{Float64,Float64},1},1}(undef,0),
    colors::Dict{Symbol, Array{RGBA{Float64},1} } = Dict{Symbol,Array{RGBA{Float64},1}}(),
    arrowHeadSize::Float64 = 6.,
    cameraPosition::Tuple{Float64,Float64,Float64} = (1., 1., 0.),
    lineWidth::Float64 = 1.0,
    lineWidths::Array{Float64,1} = fill(lineWidth,length(curves)+length(tVectors)),
    dotSize::Float64 = 1.0,
    dotSizes::Array{Float64,1} = fill(dotSize,length(points)),
    sphereColor::RGBA{Float64} = RGBA{Float64}(0.85, 0.85, 0.85, 0.6),
    sphereLineColor::RGBA{Float64} = RGBA{Float64}(0.75, 0.75, 0.75, 0.6),
    sphereLineWidth::Float64 = 0.5,
    target::Tuple{Float64,Float64,Float64} = (0.,0.,0.),
    )
    io = open(filename,"w")
    try
        #
        # Header
        # ---
        write(io,string("import settings;\nimport three;\nimport solids;",
            "unitsize(4cm);\n\n",
            "currentprojection=perspective( ",
            "camera = $(cameraPosition), ",
            "target = $(target) );\n",
            "currentlight=nolight;\n\n",
            "revolution S=sphere(O,1);\n",
            "pen SpherePen = rgb($(red(sphereColor)),",
            "$(green(sphereColor)),$(blue(sphereColor)))",
            "+opacity($(alpha(sphereColor)));\n",
            "pen SphereLinePen = rgb($(red(sphereLineColor)),",
            "$(green(sphereLineColor)),$(blue(sphereLineColor)))",
            "+opacity($(alpha(sphereLineColor)))+linewidth($(sphereLineWidth)pt);\n",
            "draw(surface(S), surfacepen=SpherePen, meshpen=SphereLinePen);\n")
        );
        write(io,"\n/*\n  Colors\n*/\n")
        j=0
        for (key,value) in colors # colors for all keys
            penPrefix = "$(j)"
            sets = 0
            if key==:points
                penPrefix="point"
                sets = length(points)
            elseif key==:curves
                penPrefix="curve"
                sets = length(curves)
            elseif key==:tvectors
                penPrefix="tVector"
                sets = length(tVectors)
            end
            if length(value) < sets
                throw( ErrorException(
                    "Not enough colors ($(length(value))) provided for $(sets) sets in $(key)."
                ))
            end
            i=0
            # export all colors
            for c in value
                i=i+1;
                if i>sets
                    # avoid access errors in lineWidth or dotSizes if more colors then sets are given
                    break
                end
                write(io,string("pen $(penPrefix)Style$(i) = ",
                    "rgb($(red(c)),$(green(c)),$(blue(c)))",
                    (key==:curves || key==:tvectors) ? "+linewidth($(lineWidths[i])pt)" : "",
                    (key==:tvectors) ? "+linewidth($(lineWidths[length(curves)+i])pt)" : "",
                    (key==:points) ? "+linewidth($(dotSizes[i])pt)" : "",
                    "+opacity($(alpha(c)));\n"));
            end
        end
        if length(points)>0
            write(io,"\n/*\n  Exported Points\n*/\n")
        end
        i=0
        for pSet in points
            i=i+1
            for point in pSet
                write(io,string("dot( (",string([string(v,",") for v in point]...)[1:end-1],"), pointStyle$(i));\n"));
            end
        end
        i=0
        if length(curves)>0
            write(io,"\n/*\n  Exported Curves\n*/\n")
        end
        for curve in curves
            i=i+1
            write(io,"path3 p$(i) = ");
            j=0
            for point in curve
                j=j+1
                pString = "("*string([ "$v," for v in point]...)[1:end-1]*")"
                write(io, j>1 ? " .. $(pString)" : pString)
            end
            write( io,string(";\n draw(p$(i), curveStyle$(i));\n") );
        end
        i=0
        if length(tVectors)>0
            write(io,"\n/*\n  Exported tangent vectors\n*/\n")
        end
        for tVecs in tVectors
            i=i+1
            j=0
            for vector in tVecs
                j=j+1
                base = vector[1]
                endPoints = base + vector[2]
                write(io,string("draw( (",
                    string( [string(v,",") for v in base]...)[1:end-1],")--(",
                    string( [string(v,",") for v in endPoints]...)[1:end-1],
                    "), tVectorStyle$(i),Arrow3($(arrowHeadSize)));\n"));
            end
        end
    finally
        close(io)
    end
end
@doc raw"""
    asymptote_export_S2_data(filename)
Export given `data` as an array of points on the sphere, i.e. one-, two-
or three-dimensional data with points on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$.

# Input
* `filename` – a file to store the Asymptote code in.

# Optional Arguments (Data)
* `data` – a point representing the 1-,2-, or 3-D array of points
* `elevationColorScheme` - A `ColorScheme` for elevation
* `scaleAxes` - (`(1/3,1/3,1/3)`) move spheres closer to each other by a factor
  per direction

# Optional Arguments (Asymptote)
* `arrowHeadSize` - (`1.8`) size of the arrowheads of the vectors (in mm)
* `cameraPosition` - position of the camrea (default: centered above xy-plane)
  szene
* `target` - position the camera points at (default: center of xy-plane within
  data).
"""
function asymptote_export_S2_data(filename::String;
    data = fill([0.,0.,1.],0,0),
    arrowHeadSize::Float64 = 1.8,
    scaleAxes = (1/3.,1/3.,1/3.),
    cameraPosition::Tuple{Float64,Float64,Float64} = scaleAxes.*( (size(data,1)-1)/2 ,(size(data,2)-1)/2, max(size(data,3),0)+10),
    target::Tuple{Float64,Float64,Float64} = scaleAxes.*( (size(data,1)-1)/2 ,(size(data,2)-1)/2, 0.),
    elevationColorScheme = ColorSchemes.viridis,
    )
    io = open(filename,"w")
    try
      write(io,string("import settings;\nimport three;\n",
        "size(7cm);\n",
        "DefaultHead.size=new real(pen p=currentpen) {return $(arrowHeadSize)mm;};\n",
        "currentprojection=perspective( ",
        "camera = $(cameraPosition), up=Y,",
        "target = $(target) );\n\n"));
      dims = [size(data,i) for i=[1,2,3]]
      for x=1:dims[1]
        for y=1:dims[2]
          for z=1:dims[3]
            v = Tuple(data[x,y,z]) #extract value
            el = asin( min(1,max(-1,v[3])) ); # since 3 is between -1 and 1 this yields a value between 0 and pi
            # map elevation to colormap
            c = get(elevationColorScheme,el+π/2, (0.,Float64(π)) );
            # write arrow in this colormap
            # transpose image to comply with image adresses (first index column downwards, second rows)
            write(io,string("draw( $(scaleAxes.*(x-1,y-1,z-1))",
              "--$(scaleAxes.*(x-1,y-1,z-1).+v),",
              " rgb($(red(c)),$(green(c)),$(blue(c))), Arrow3);\n"));
          end
        end
      end
    finally
      close(io)
    end
end
@doc raw"""
    asymptote_export_SPD(filename)

export given `data` as a point on a `Power{SPDPoint}` manifold, i.e. one-, two-
or three-dimensional data with points on the manifold of symmetric positive
definite matrices.

# Input
* `filename` – a file to store the Asymptote code in.

# Optional Arguments (Data)
* `data` – a point representing the 1-,2-, or 3-D array of `SPDPoints`
* `colorScheme` - A `ColorScheme` for Geometric Anisotropy Index
* `scaleAxes` - (`(1/3,1/3,1/3)`) move symmetric positive definite matrices
  closer to each other by a factor per direction compared to the distance
  esimated by the maximal eigenvalue of all involved SPD points

# Optional Arguments (Asymptote)
* `cameraPosition` - position of the camrea (default: centered above xy-plane)
  szene.
* `target` - position the camera points at (default: center of xy-plane within data).

Both values `cameraPosition` and `target` are scaled by `scaledAxes*EW`, where
`EW` is the maximal eigenvalue in the `data`.
"""
function asymptote_export_SPD(filename::String;
    data = fill(Matrix{Float64}(I,3,3),0,0),
    scaleAxes = (1/3.,1/3.,1/3.) #multiplied with the maximal eigenvalue if data is present
      .* ( length(data) > 0 ? maximum(maximum(eigvals.( data ))) : 1 ),
    cameraPosition::Tuple{Float64,Float64,Float64} = ( (size(data,1)-1)/2 ,(size(data,2)-1)/2, max(size(data,3),0.)+10.),
    target::Tuple{Float64,Float64,Float64} = ( (size(data,1)-1)/2 ,(size(data,2)-1)/2, 0.),
    colorScheme = ColorSchemes.viridis,
    )
    io = open(filename,"w")
    try
      write(io,string("import settings;\nimport three;\n",
      "surface ellipsoid(triple v1,triple v2,triple v3,real l1,real l2, real l3, triple pos=O) {\n",
      "  transform3 T = identity(4);\n",
      "  T[0][0] = l1*v1.x;\n  T[1][0] = l1*v1.y;\n  T[2][0] = l1*v1.z;\n",
      "  T[0][1] = l2*v2.x;\n  T[1][1] = l2*v2.y;\n  T[2][1] = l2*v2.z;\n",
      "  T[0][2] = l3*v3.x;\n  T[1][2] = l3*v3.y;\n  T[2][2] = l3*v3.z;\n",
      "  T[0][3] = pos.x;\n  T[1][3] = pos.y;\n  T[2][3] = pos.z;\n",
      "  return T*unitsphere;\n}\n\n",
      "size(200);\n\n",
      "real gDx=$(scaleAxes[1]);\n",
      "real gDy=$(scaleAxes[2]);\n",
      "real gDz=$(scaleAxes[3]);\n\n",
      "currentprojection=perspective(up=Y, ",
      "camera = (gDx*$(cameraPosition[1]),gDy*$(cameraPosition[2]),gDz*$(cameraPosition[3])), ",
      "target = (gDx*$(target[1]),gDy*$(target[2]),gDz*$(target[3])) );\n",
      "currentlight=Viewport;\n\n",
      ));
      dims = [size(data,1) size(data,2) size(data,3) ];
      for x=1:dims[1]
        for y=1:dims[2]
          for z=1:dims[3]
            A = data[x,y,z] #extract matrix
            F = eigen(A)
            if maximum(abs.(A)) > 0. # a nonzero matrix (exclude several pixel
              # Following Moakher & Batchelor: Geometric Anisotropic Index:
              λ = F.values
              V = F.vectors
              Lλ = log.(λ)
              GAI = sqrt( 2/3*sum( Lλ.^2 )  - 2/3 * sum(sum( tril( Lλ * Lλ',-1), dims=1), dims=2)[1])
              c = get(colorScheme, GAI/(1+GAI), (0,1) )
              write(io,string("  draw(  ellipsoid( ($(V[1,1]),$(V[2,1]),$(V[3,1])),",
              " ($(V[1,2]),$(V[2,2]),$(V[3,2])), ($(V[1,3]),$(V[2,3]),$(V[3,3])),",
              " $(λ[1]), $(λ[2]), $(λ[3]), ",
              " (gDx*$(x-1), gDy*$(y-1), gDz*$(z-1))),",
              " rgb($(red(c)),$(green(c)),$(blue(c)))  );\n"));
            end
          end
        end
      end
    finally
      close(io)
    end
end

export renderAsymptote, asyExportS2Signals, asyExportS2Data, asyExportSPDData
"""
    render_asymptote(filename; render=4, format="png", ...)
render an exported asymptote file specified in the `filename`, which can also
be given as a relative or full path

# Input

* `filename` – filename of the exported `asy` and rendered image

# Keyword Arguments

the default values are given in brackets

* `render` – (`4`) render level of asymptote, i.e. its `-render` option
* `format` – (`"png"`) final rendered format, i.e. asymptote's `-f` option
* `export_file` - (the filename with format as ending) specify the export filename
"""
function render_asymptote(
    filename;
    render::Int=4,
    format="png",
    exportFolder = string( filename[1:( [findlast(".",filename)...][1])], format),
)
    renderCmd = `asy -render $(render) -f $(format) -globalwrite  -o "$(relpath(exportFolder))" $(filename)`
    run(renderCmd)
end