text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@testset "IterationStats" begin

  @testset "max_primal_violation" begin
    # min 0
    # s.t. y == 10
    # z >= 11 (as a linear constraint)
    # -1 <= x <= 1
    lp = FirstOrderLp.linear_programming_problem(
      [-1.0, -Inf, -Inf],  # variable_lower_bound
      [1.0, Inf, Inf],  # variable_upper_bound
      zeros(3),  # objective_vector
      0.0,  # objective_constant
      [0.0 1.0 0.0; 0.0 0.0 1.0],  # constraint_matrix
      [10.0, 11.0],  # right_hand_side
      1,  # num_equalities
    )
    @test FirstOrderLp.max_primal_violation(lp, [0.0, 10.0, 11.0]) == 0.0
    @test FirstOrderLp.max_primal_violation(lp, [-2.0, 10.0, 11.0]) ≈ 1.0
    @test FirstOrderLp.max_primal_violation(lp, [3.0, 10.0, 11.0]) ≈ 2.0
    @test FirstOrderLp.max_primal_violation(lp, [0.0, 11.0, 11.0]) ≈ 1.0
    @test FirstOrderLp.max_primal_violation(lp, [0.0, 9.0, 11.0]) ≈ 1.0
    @test FirstOrderLp.max_primal_violation(lp, [0.0, 11.0, 0.0]) ≈ 11.0
  end

  @testset "primal_obj" begin
    qp = example_qp()
    @test FirstOrderLp.primal_obj(qp, [0.0, 0.0]) == 0.0
    @test FirstOrderLp.primal_obj(qp, [1.0, 1.0]) == 0.5
    @test FirstOrderLp.primal_obj(qp, [1.0, 0.0]) == 1.0
    @test FirstOrderLp.primal_obj(qp, [0.0, 1.0]) == -0.5
    @test FirstOrderLp.primal_obj(qp, [0.0, -1.0]) == 1.5
  end

  @testset "dual_stats" begin
    # Primal:
    # min x + 2y
    # x + y >= 1
    # -1 <= x <= 1
    # Dual:
    # max a - b - c
    # s.t. a + b - c = 1
    #      a = 2
    # a, b, c >= 0
    lp = FirstOrderLp.linear_programming_problem(
      [-1.0, -Inf],  # variable_lower_bound
      [1.0, Inf],  # variable_upper_bound
      [1.0, 2.0],  # objective_vector
      0.0,  # objective_constant
      reshape([1.0, 1.0], 1, 2),  # constraint_matrix
      [1.0],  # right_hand_side
      0,  # num_equalities
    )
    # b and c are reduced costs. The dual vector is [a].
    dual_stats0 = FirstOrderLp.compute_dual_stats(lp, [0.0, 0.0], [0.0])
    @test dual_stats0.dual_objective == -1.0
    @test norm(dual_stats0.dual_residual, Inf) == 2.0
    @test dual_stats0.dual_residual == [0.0; 0.0; 2.0]

    dual_stats1 = FirstOrderLp.compute_dual_stats(lp, [0.0, 0.0], [1.0])
    @test dual_stats1.dual_objective == 1.0
    @test norm(dual_stats1.dual_residual, Inf) == 1.0
    @test dual_stats1.dual_residual == [0.0; 0.0; 1.0]

    dual_stats2 = FirstOrderLp.compute_dual_stats(lp, [0.0, 0.0], [2.0])
    @test dual_stats2.dual_objective == 1.0
    @test norm(dual_stats2.dual_residual, Inf) == 0.0

    dual_stats3 = FirstOrderLp.compute_dual_stats(lp, [0.0, 0.0], [3.0])
    @test dual_stats3.dual_objective == 1.0
    @test norm(dual_stats3.dual_residual, Inf) == 1.0

    dual_stats4 = FirstOrderLp.compute_dual_stats(lp, [0.0, 1.0], [-1.0])
    @test dual_stats4.dual_objective == -3.0
    @test dual_stats4.dual_residual == [1.0, 0.0, 3.0]

    lp = FirstOrderLp.linear_programming_problem(
      [Inf, -Inf],  # variable_lower_bound
      [Inf, Inf],  # variable_upper_bound
      [1.0, 2.0],  # objective_vector
      0.0,  # objective_constant
      reshape([1.0, 1.0], 1, 2),  # constraint_matrix
      [1.0],  # right_hand_side
      0,  # num_equalities
    )
    dual_stats5 = FirstOrderLp.compute_dual_stats(lp, [0.0, 1.0], [-1.0])
    @test dual_stats5.dual_objective == -1.0
    @test dual_stats5.dual_residual == [1.0, 2.0, 3.0]

    qp = example_qp()
    dual_stats6 = FirstOrderLp.compute_dual_stats(qp, [0.0, 0.0], [3.0])
    @test dual_stats6.dual_objective == -3.0
    @test norm(dual_stats6.dual_residual, Inf) == 0.0

    dual_stats7 = FirstOrderLp.compute_dual_stats(qp, [0.0, 0.0], [1.0])
    @test dual_stats7.dual_objective == -1.0
    @test norm(dual_stats7.dual_residual, Inf) == 0.0

    dual_stats8 = FirstOrderLp.compute_dual_stats(qp, [0.5, 0.5], [1.0])
    @test dual_stats8.dual_objective == -1.625
    @test norm(dual_stats8.dual_residual, Inf) == 0.0
  end

  @testset "Primal/dual optimal" begin
    # Primal:
    # min x + 2y
    # x + y >= 1
    # -1 <= x <= 1
    # Dual:
    # max a - b - c
    # s.t. a + b - c = 1
    #      a = 2
    # a, b, c >= 0
    lp = FirstOrderLp.linear_programming_problem(
      [-1.0, -Inf],  # variable_lower_bound
      [1.0, Inf],  # variable_upper_bound
      [1.0, 2.0],  # objective_vector
      0.0,  # objective_constant
      reshape([1.0, 1.0], 1, 2),  # constraint_matrix
      [1.0],  # right_hand_side
      0,  # num_equalities
    )
    # b and c are reduced costs. The dual vector is [a].
    stats = FirstOrderLp.compute_iteration_stats(
      lp,
      FirstOrderLp.cached_quadratic_program_info(lp),
      [1.0, 0.0],  # primal_iterate
      [2.0],  # dual_iterate
      [0.0, 0.0],  # primal_ray_estimate
      [0.0],  # dual_ray_estimate
      5,  # iteration_number
      1.5,  # cumulative_kkt_matrix_passses
      5.0,  # cumulative_time_sec
      1e-6,  # eps_optimal_absolute
      1e-6,  # eps_optimal_relative
      1.0,  # step_size
      1.0,  # primal_weight
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE,
    )
    correct_convergence_info = FirstOrderLp.ConvergenceInformation()
    correct_convergence_info.primal_objective = 1.0
    correct_convergence_info.dual_objective = 1.0
    correct_convergence_info.corrected_dual_objective = 1.0
    correct_convergence_info.l_inf_primal_variable = 1.0
    correct_convergence_info.l2_primal_variable = 1.0
    correct_convergence_info.l_inf_dual_variable = 2.0
    correct_convergence_info.l2_dual_variable = 2.0
    correct_convergence_info.candidate_type =
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_infeas_info = FirstOrderLp.InfeasibilityInformation()
    correct_infeas_info.candidate_type = FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_stats = FirstOrderLp.IterationStats()
    correct_stats.iteration_number = 5
    correct_stats.cumulative_kkt_matrix_passes = 1.5
    correct_stats.cumulative_time_sec = 5.0
    correct_stats.convergence_information = [correct_convergence_info]
    correct_stats.infeasibility_information = [correct_infeas_info]
    correct_stats.step_size = 1.0
    correct_stats.primal_weight = 1.0

    test_fields_equal(stats, correct_stats)
  end

  @testset "Primal infeasible" begin
    # Primal:
    # min x + 2
    # s.t. x = 10
    # 0 <= x <= 1
    # Dual:
    # max 10a - b + 2
    # s.t. a - b <= 1
    # b >= 0.
    lp = FirstOrderLp.linear_programming_problem(
      [0.0],  # variable_lower_bound
      [1.0],  # variable_upper_bound
      [1.0],  # objective_vector
      2.0,  # objective_constant
      reshape([1.0], 1, 1),  # constraint_matrix
      [10.0],  # right_hand_side
      1,   # num_equalities
    )
    # b is a reduced cost. The dual vector is [a].
    stats = FirstOrderLp.compute_iteration_stats(
      lp,
      FirstOrderLp.cached_quadratic_program_info(lp),
      [2.0],  # primal_iterate
      [1.0],  # dual_iterate
      [0.0],  # primal_ray_estimate
      [1.0],  # dual_ray_estimate
      5,  # iteration_number
      1.5,  # cumulative_kkt_matrix_passses
      5.0,  # cumulative_time_sec
      1e-6,  # eps_optimal_absolute
      1e-6,  # eps_optimal_relative
      1.0,  # step_size
      1.0,  # primal_weight
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE,
    )
    correct_convergence_info = FirstOrderLp.ConvergenceInformation()
    correct_convergence_info.primal_objective = 4.0
    correct_convergence_info.dual_objective = 10.0 + 2.0
    correct_convergence_info.corrected_dual_objective = 12.0
    correct_convergence_info.l_inf_primal_residual = 8.0
    correct_convergence_info.l2_primal_residual = norm([8.0, 1.0], 2)
    correct_convergence_info.relative_l_inf_primal_residual = 8.0 / (1.0 + 10.0)
    correct_convergence_info.relative_l2_primal_residual =
      norm([8.0, 1.0], 2) / (1.0 + 10.0)
    correct_convergence_info.relative_optimality_gap = 8.0 / (1.0 + 16.0)
    correct_convergence_info.l_inf_primal_variable = 2.0
    correct_convergence_info.l2_primal_variable = 2.0
    correct_convergence_info.l_inf_dual_variable = 1.0
    correct_convergence_info.l2_dual_variable = 1.0
    correct_convergence_info.candidate_type =
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_infeas_info = FirstOrderLp.InfeasibilityInformation()
    correct_infeas_info.dual_ray_objective = 9.0
    correct_infeas_info.candidate_type = FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_stats = FirstOrderLp.IterationStats()
    correct_stats.iteration_number = 5
    correct_stats.cumulative_kkt_matrix_passes = 1.5
    correct_stats.cumulative_time_sec = 5.0
    correct_stats.convergence_information = [correct_convergence_info]
    correct_stats.infeasibility_information = [correct_infeas_info]
    correct_stats.step_size = 1.0
    correct_stats.primal_weight = 1.0

    test_fields_equal(stats, correct_stats)
  end

  @testset "Dual infeasible" begin
    # Primal:
    # min -x
    # s.t. x >= 10
    # Dual:
    # max 10a
    # s.t. a <= -1
    # a >= 0.
    lp = FirstOrderLp.linear_programming_problem(
      [-Inf],  # variable_lower_bound
      [Inf],  # variable_upper_bound
      [-1.0],  # objective_vector
      0.0,  # objective_constant
      reshape([1.0], 1, 1),  # constraint_matrix
      [10.0],  # right_hand_side
      0,  # num_equalities
    )
    stats = FirstOrderLp.compute_iteration_stats(
      lp,
      FirstOrderLp.cached_quadratic_program_info(lp),
      [10.0],  # primal_iterate
      [0.0],  # dual_iterate
      [1.0],  # primal_ray_estimate
      [0.0],  # dual_ray_estimate
      5,  # iteration_number
      1.5,  # cumulative_kkt_matrix_passses
      5.0,  # cumulative_time_sec
      1e-6,  # eps_optimal_absolute
      1e-6,  # eps_optimal_relative
      1.0,  # step_size
      1.0,  # primal_weight
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE,
    )
    correct_convergence_info = FirstOrderLp.ConvergenceInformation()
    correct_convergence_info.primal_objective = -10.0
    correct_convergence_info.corrected_dual_objective = -Inf
    correct_convergence_info.l_inf_dual_residual = 1.0
    correct_convergence_info.l2_dual_residual = 1.0
    correct_convergence_info.relative_l_inf_dual_residual = 1.0 / (1.0 + 1.0)
    correct_convergence_info.relative_l2_dual_residual = 1.0 / (1.0 + 1.0)
    correct_convergence_info.relative_optimality_gap = 10.0 / (1.0 + 10.0)
    correct_convergence_info.l_inf_primal_variable = 10.0
    correct_convergence_info.l2_primal_variable = 10.0
    correct_convergence_info.candidate_type =
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_infeas_info = FirstOrderLp.InfeasibilityInformation()
    correct_infeas_info.primal_ray_linear_objective = -1.0
    correct_infeas_info.candidate_type = FirstOrderLp.POINT_TYPE_CURRENT_ITERATE

    correct_stats = FirstOrderLp.IterationStats()
    correct_stats.iteration_number = 5
    correct_stats.cumulative_kkt_matrix_passes = 1.5
    correct_stats.cumulative_time_sec = 5.0
    correct_stats.convergence_information = [correct_convergence_info]
    correct_stats.infeasibility_information = [correct_infeas_info]
    correct_stats.step_size = 1.0
    correct_stats.primal_weight = 1.0

    test_fields_equal(stats, correct_stats)
  end

  @testset "print_to_screen_this_iteration" begin
    termination_evaluation_frequency = convert(Int32, 10)
    @test FirstOrderLp.print_to_screen_this_iteration(
      false, # termination_reason
      1, # iteration,
      2, # verbosity
      termination_evaluation_frequency,
    )
    @test FirstOrderLp.print_to_screen_this_iteration(
      false, # termination_reason
      101, # iteration,
      5, # verbosity
      termination_evaluation_frequency,
    )
    @test !FirstOrderLp.print_to_screen_this_iteration(
      false, # termination_reason
      31, # iteration,
      5, # verbosity
      termination_evaluation_frequency,
    )
    @test !FirstOrderLp.print_to_screen_this_iteration(
      false, # termination_reason
      531, # iteration,
      5, # verbosity
      termination_evaluation_frequency,
    )
    @test FirstOrderLp.print_to_screen_this_iteration(
      true, # termination_reason
      124, # iteration,
      5, # verbosity
      termination_evaluation_frequency,
    )
  end
end


text:
function get_mesh_1d(i)

    coordinates_min = (-1,)
    coordinates_max = ( 1,)

    mesh1_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=4,
                    n_cells_max=10_000
                    )

    coordinates_min = (-1,)
    coordinates_max = ( 1,)
    #refinement_patches = (
    #  (type="box", coordinates_min=(-1.0), coordinates_max=(1.0)),
    #)
    mesh2_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=5,
                    #refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )

    coordinates_min = (0,)
    coordinates_max = ( 2,)
    refinement_patches = (
    (type="box", coordinates_min=(0.5), coordinates_max=(1.5)),
    )
    mesh3_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=6,
                    refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )

    coordinates_min = (-1,)
    coordinates_max = ( 1,)
    refinement_patches = (
    (type="box", coordinates_min=(-0.5), coordinates_max=(0.5)),
    (type="box", coordinates_min=(-0.2), coordinates_max=(0.2)),
    )
    mesh4_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=5,
                    refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )
    if i==1
        mesh = mesh1_1d
    elseif i==2
        mesh = mesh2_1d
    elseif i==3
        mesh = mesh3_1d
    elseif i==4
        mesh = mesh4_1d
    end
    return mesh
end

ext: jl


text:
begin
	using Functors:functor
	using Flux

	using Transformers
	using Transformers.Basic #for loading the positional embedding
	using Flux: @functor, Params
	using Flux: update!

	#flux function for update parameters
	using Flux: gradient
	using Statistics: norm, mean, std

	labels = collect(1:10)
	startsym = 11
	endsym = 12
	unksym = 0
	labels = [unksym, startsym, endsym, labels...]
	vocab = Vocabulary(labels, unksym)
	#function for generate training datas
	sample_data() = rand(1:10, 10)
	#function for adding start & end symbol
	preprocess(x) = [startsym, x..., endsym]

	sample = preprocess(sample_data())

	batch=hcat(preprocess.((sample_data() for i = 1:32))...)
end

begin
	mutable struct  simpleLayer
		embd::Embed
		pe::PositionEmbedding
		linear::Positionwise
	end

	@functor simpleLayer

	function simpleLayer()
		embd=Embed(512, length(vocab))
		pe = PositionEmbedding(512, trainable=true)
		linear = Positionwise(Dense(512, length(vocab)))
		simpleLayer(embd, pe, linear)
	end

	function (t::simpleLayer)(x)
		we = t.embd(x, inv(sqrt(512)))
		e = we .+ t.pe(we)
		p = t.linear(e)
	end

	model=simpleLayer()

	a, b = functor(model)

	gs = gradient(m->mean(m(batch)), model)[1]
  #diff way to get the gradient, define the fuction inside do, don't pass the
	#use model parameters
	gs2=gradient(params(model)) do
		mean(model(batch))
	end

   #PROBLEM:  We still get nothing for the gradient of the embedding layer
	for p in params(model)
		println(typeof(gs2[p]))
	end
	model.pe.embedding
end

begin
	struct GPTPositionEmbedding{F, W <: AbstractArray{F}}
	    embedding::W
	end

	@functor GPTPositionEmbedding

	function GPTPositionEmbedding(size::Int, max_len::Int = 1024)
	    GPTPositionEmbedding( zeros(Float32,size, max_len))
	end

	function (pe::GPTPositionEmbedding{F})(x::AbstractArray{F}) where F
	    pe.embedding[:, 1:size(x,2)]  # just return the needed position embedding
	end

end

begin

   struct  GPTTestLayer
		embd::Embed
		pe::GPTPositionEmbedding
		linear::Positionwise
	end
		@functor GPTTestLayer

		function (t::GPTTestLayer)(x)
			we = t.embd(x, inv(sqrt(512)))
			e = we .+ t.pe(we)
			p = t.linear(e)
		end

	function GPTTestLayer()
		embd=Embed(512, length(vocab))
		pe = GPTPositionEmbedding(512, length(vocab))
		linear = Positionwise(Dense(512, length(vocab), initW=ones,initb=ones))
		GPTTestLayer(embd, pe, linear)
	end
	model=GPTTestLayer()

	for i in 1:10
		ps = params(model)

		gs2=gradient(ps) do
			mean(model(batch))
		end

		for p in ps
			println(size(p), mean(p), mean(gs2[p]), typeof(gs2[p]), size(gs2[p]))
		end
		opt = ADAM(1e-4)
		update!(opt, ps, gs2)
	end
end

#Try simple model with straight array as possition embedding.
begin

	mutable struct  simpleLayerArrayPosEmd
		embd::Embed
		pe
		#linear::Positionwise
	end
		@functor simpleLayerArrayPosEmd

		function (t::simpleLayerArrayPosEmd)(x)
			we = t.embd(x, inv(sqrt(512)))
			@show size(we)

			e = we .+ t.pe[:, 1:size(we,2)]
			#p = t.linear(e)
		end

	function simpleLayerArrayPosEmd()
		embd=Embed(512, length(vocab))
		pe = zeros(512,100)   ## pos embedding is larger that training, so some will remain zero
		#linear = Positionwise(Dense(512, length(vocab)))
		simpleLayerArrayPosEmd(embd, pe) #, linear)
	end
	model=simpleLayerArrayPosEmd()

    model(batch)
	ps = params(model)
	gs3=gradient(ps) do
		mean(model(batch))
	end

	for p in ps
		println(typeof(gs3[p]))
	end

	for i in 1:10

		gs4=gradient(ps) do
			mean(model(batch))
		end

		for p in ps
			println(size(p), mean(p), mean(gs4[p]), typeof(gs4[p]), size(gs4[p]))
		end
		opt = ADAM(1e-4)
		update!(opt, ps, gs4)
	end

end

max_stars_repo_name: daryoush/CudaAndFluxExperiments

max_issues_repo_name: daryoush/CudaAndFluxExperiments

max_forks_repo_name: daryoush/CudaAndFluxExperiments


text:
using NKLandscapes
using Base.Test

const BOUNDARIES = ["Nearest-neighbor", "Random"]
const LC = 100 # Landscape count
const PAIRS = [(8, 0), (8, 2), (8, 4), (8, 8),
(16, 0), (16, 2), (16, 4), (16, 8), (16, 16),
(24, 0), (24, 2), (24, 4), (24, 8), (24, 16), (24, 24),
(48, 0), (48, 2), (48, 4), (48, 8), (48, 16), (48, 24), (48, 48),
(96, 0), (96, 2), (96, 4), (96, 8), (96, 16), (96, 24), (96, 48), (96, 96)]

# Mean fitness of local optima (nearest-neighbor interactions) & Mean walk
# lengths to local optima (nearest-neighbor interactions) Kauffman, p. 55-56

# Mean fitness of local optima (random interactions) & Mean walk lengths to
# local optima (random interactions) Kauffman, p. 57

println("Length of Paths to Optima and Fitness of Optima")
for b = [1, 2]
  println(BOUNDARIES[b], " Interactions")
  println("N\tK\tμ len\tσ len\tμ fit\tσ fit")
  for (n, k) = PAIRS
    if n == k
      k = k - 1
      # Periodic boundary requires even K
      if b == 1
        k = k - 1
      end
    end
    lengths = zeros(LC)
    fitnesses = zeros(LC)
    for i = 1:LC
      l = NKLandscape(n, k, near=b == 1)
      g = rand(Genotype, l)
      w = random_adaptive_walk(g)
      lengths[i] = w.length
      fitnesses[i] = fitness(w.history_list[end])
    end
    @printf("%d\t%d\t%.2f\t%.2f\t%.3f\t%.3f\n", n, k, mean(lengths), std(lengths), mean(fitnesses), std(fitnesses))
  end
end

# Number of optima. Kauffman, p. 60.

MAX_FOUND = 10_000
MAX_FAILS = 100

println("Number of Optima")
for b = [1, 2]
  println(BOUNDARIES[b], " Interactions")
  println("N\tK\tμ\tσ")
  for (n, k) = PAIRS[1:9]  # Only runs tests for N = 8 and N = 16.  Modify to run more tests (time consuming).
    if n == k
      k = k - 1
      if b == 1
        k = k - 1
      end
    end
    counts = zeros(LC)
    for i = 1:LC
      since_last = 0
      found_opts = Set{Genotype}()
      l = NKLandscape(n, k, near=b == 1)
      while length(found_opts) <= MAX_FOUND && since_last <= MAX_FAILS
        g = rand(Genotype, l)
        w = random_adaptive_walk(g)
        opt = w.history_list[end]
        if !(opt in found_opts)
          push!(found_opts, opt)
          since_last = 0
        else
          since_last += 1
        end
      end
      counts[i] = length(found_opts)
    end
    @printf("%d\t%d\t%.2f\t%.2f\n", n, k, mean(counts), std(counts))
  end
end


text:
"""
    MMALagOptions

A struct that stores all the options of the MMA algorithms. Th following are the fields of `MMALagOptions`:
 - `maxiter`: the maximum number of inner iterations. For `MMA87`, there is 1 inner iteration per outer iteration.
 - `outer_maxiter`: the maximum number of outer iterations.
 - `maxinner`: the maximum number of inner iterations per outer iteration of [`MMA02`](@ref). Not applicable for [`MMA87`](@ref).
 - `tol`: a tolerance struct of type [`Tolerance`](@ref).
 - `s_init`: defined in the original [`MMA02`](@ref) paper.
 - `s_incr`: defined in the original [`MMA02`](@ref) paper.
 - `s_decr`: defined in the original [`MMA02`](@ref) paper.
 - `store_trace`: if true, a trace will be stored.
 - `auto_scale`: if set to `true`, the objective will be auto-scaled to make the norm of the objective and its gradient similar to the constriants'.
 - `dual_options`: the options passed to the dual optimizer from [`Optim.jl`](https://github.com/JuliaNLSolvers/Optim.jl).
 - `tonormalize`: if set to `true`, the dual optimizer will normalize the dual variables and perform gradient descent on the hyper-sphere.
 - `norm`: the norm type used to normalize the dual variables if `tonormalize` is set to `true`.
 - `callback_steps`: the [`MMALagCallback`](@ref) is called every `callback_steps` outer iterations. The callback increments the quadratic penalty factor and the linear weights are updated.
 - `quad_incr`: the factor by which the quadratic penalty is incremented every time the [`MMALagCallback`](@ref) is called.
 - `dual_step`: the dual solution will be incremented according to a gradient descent scheme which is `Descent(dual_step * nconstr)` by default where `nconstr` is the number of constraints in the problem.
 - `dual_opt`: the dual gradient descent scheme that the dual solution is incremented by. This defaults to `Descent(dual_step * nconstr)` where `dual_step` is the keyword argument passed in and `nconstr` is the number of constraints. Other examples include [`Nesterov`](@ref), [`Momentum`](@ref) and [`ExpDecay`](@ref).
 - `verbose`: when set to `true`, more information is printed to the screen.
"""
@with_kw mutable struct MMALagOptions{T, Ttol <: Tolerance, TSubOptions <: Optim.Options}
    maxiter::Int = 1000
    outer_maxiter::Int = 10^8
    maxinner::Int = 10
    tol::Ttol = Tolerance()
    s_init::T = 0.1
    s_incr::T = 1.2
    s_decr::T = 0.7
    store_trace::Bool = false
    show_trace::Bool = false
    auto_scale::Bool = false
    keep_best::Bool = false
    dual_options::TSubOptions = Optim.Options(allow_f_increases = false, iterations = 1000, outer_iterations=1000)
    tonormalize::Bool = true
    norm = 2
    callback_steps::Int = 3
    quad_incr::T = 1.2
    dual_step::T = 0.01
    dual_opt = nothing
    verbose::Bool = false
end

@params struct MMALag <: AbstractOptimizer
    parent::AbstractOptimizer
end
MMALag() = MMALag(MMA02())

default_options(model::MMALagModel, alg::MMA87) = MMALagOptions()
default_options(model::MMALagModel, alg::MMA02) = MMALagOptions()
default_options(model::VecModel, alg::MMALag) = MMALagOptions()

@params mutable struct MMALagCallback <: Function
    previter::Int
    prevf::Real
    stochopt
end
function MMALagCallback(options, nconstr::Int=1)
    @unpack dual_opt, dual_step = options
    #opt = ExpDecay(0.1, 0.9, 1, 1e-6)
    #opt = Descent(nconstr * 0.01)
    nconstr = 1
    if dual_opt === nothing
        dual_opt = Descent(nconstr * dual_step)
    end
    return MMALagCallback(0, Inf, dual_opt)
end

function (callback::MMALagCallback)(workspace::MMAWorkspace)
    @unpack dualmodel, optimizer, options, solution = workspace
    @unpack convstate = solution
    @unpack tol = options
    infeastol = tol.infeas
    model = dualmodel |> getparent |> getparent
    callbackrun = workspace.outer_iter == 2 ||
        workspace.iter >= options.maxiter ||
        workspace.outer_iter >= options.outer_maxiter ||
        convstate.converged && convstate.infeas > infeastol ||
        workspace.outer_iter - callback.previter > options.callback_steps

    if callbackrun
        mmalagupdate!(
            model,
            workspace,
            solution.prevf,
            callback.stochopt,
            first = workspace.outer_iter == 2,
        )
        callback.previter = workspace.outer_iter
    end
    return callbackrun
end

mmalagupdate!(model::VecModel, workspace, prevf, stochopt; first = false) = model
function mmalagupdate!(model::MMALagModel, workspace, prevf, stochopt; first = false)
    @unpack dualmodel, solution, optimizer, options, convcriteria = workspace
    @unpack convstate = solution
    @unpack tol, tonormalize = options
    normn = options.norm
    ftol, infeastol = tol.f, tol.infeas
    model = dualmodel |> getparent |> getparent

    f = solution.f
    g = getorigconstrval(model)

    N = length(getlinweights(model)) ÷ 2
    if tonormalize
        dims = getdim.(getineqconstraints(getparent(model)).fs)
        nc = length(dims)

        λ = getlinweights(model)[1:N]
        start = 1
        for i in 1:nc
            inds = start:start+dims[i]-1
            λ[inds] .= getlinweights(model)[inds] .* solution.λ[i] .+ getlinweights(model)[inds .+ N] .* solution.λ[i .+ nc]
            start += dims[i]
        end
        normalize!(λ, normn)

        direction = copy(g)
        start = 1
        for i in 1:nc
            inds = start:start+dims[i]-1
            normval = norm(direction[inds], normn)
            if normval != 0
                direction[inds] .= direction[inds] ./ normval
            end
            start += dims[i]
        end
        apply!(stochopt, λ, direction)

        newλ = similar(λ)
        start = 1
        for i in 1:nc
            inds = start:start+dims[i]-1
            if length(inds) == 1
                newλ[start] = 1
            else
                temp = max.(0, (λ[inds] .+ direction[inds]))
                newλ[inds] = temp
                normval = norm(temp, normn)
                if normval != 0
                    newλ[inds] .= newλ[inds] ./ normval
                end
            end
            start += dims[i]
        end
        setlinweights!(model, [λ; newλ])
    else
        dims = getdim.(getineqconstraints(getparent(model)).fs)
        nc = length(dims)

        λ = getlinweights(model)[1:N]
        start = 1
        for i in 1:nc
            inds = start:start+dims[i]-1
            λ[inds] .= getlinweights(model)[inds] .* solution.λ[i] .+ getlinweights(model)[inds .+ N] .* solution.λ[i + nc]
            start += dims[i]
        end
        direction = copy(g)#.* (1 .+ rand.())
        apply!(stochopt, λ, direction)
        newλ = max.(0, (λ .+ direction))
        setlinweights!(model, [λ; newλ])
    end

    c = getquadweight(model)
    if convstate.infeas <= infeastol
        newc = c
    else
        newc = min(options.quad_incr * c, 1e6)
    end
    setquadweight!(model, newc)
    if debugging[]
        @show getquadweight(model)

        println()
        #@show maximum(g)
        #@show convstate.infeas
        #@show getoptimalx(dualmodel), solution.f
        #@show newλ, newc
    end

    approxmodel = getparent(dualmodel)
    optimalx = getoptimalx(dualmodel)
    aggfg, agg∇fg = value_jacobian(
        approxmodel.objective_ineq_constraints,
        optimalx,
    )
    workspace.fcalls += 1

    # Update the convex approximation at the current x
    updateapprox!(dualmodel, optimalx, aggfg, agg∇fg)

    # Update the objective and constraint values and gradients in `solution`
    updatefg!(solution, aggfg, agg∇fg)

    correctsolution!(solution, model, options)

    return model
end

function init!(model::MMALagModel)
    set_objective_multiple!(model, 1)
    setquadweight!(model, 1e-5)
    setlinweights!(model, 1)
    return model
end

function scalequadweight!(model::MMALagModel, s::Real)
    setquadweight!(model, getquadweight(model) * s)
    return model
end

function correctsolution!(solution::Solution, model::MMALagModel, options::MMALagOptions)
    @unpack convstate = solution
    @unpack tol = options
    infeastol = tol.infeas
    f = getorigobjval(model)
    if debugging[]
        @show f
    end
    g = getorigconstrval(model)
    infeas = max(maximum(g), 0)
    if debugging[]
        @show infeas
    end
    convstate.infeas = infeas
    convstate.converged = convstate.converged && infeas <= infeastol
    solution.f = f
    return solution
end

function Workspace(model::VecModel, alg::MMALag, x0::AbstractVector, args...; kwargs...)
    return MMAWorkspace(MMALagModel(model; kwargs...), alg.parent, x0, args...; kwargs...)
end


text:
using Base.Test, CodeTools

@test CodeTools.getthing("Base.fft") == fft
@test CodeTools.getthing(Base, [:fft]) == fft

@test CodeTools.filemodule(Pkg.dir("CodeTools", "src", "module.jl")) == "CodeTools"


text:
# Precompute flags for the B matrix
@enum PrecomputeFlags begin
  FULL = 1
  TENSOR = 2
  LINEAR = 3
  POLYNOMIAL = 4
end

# Timing functions that allow for timing parts of an NFFT

mutable struct TimingStats
  pre::Float64
  conv::Float64
  fft::Float64
  apod::Float64
  conv_adjoint::Float64
  fft_adjoint::Float64
  apod_adjoint::Float64
end

TimingStats() = TimingStats(0.0,0.0,0.0,0.0,0.0,0.0,0.0)

function Base.println(t::TimingStats)
  print("Timing: ")
  @printf "pre = %.4f s apod = %.4f / %.4f s fft = %.4f / %.4f s conv = %.4f / %.4f s\n" t.pre t.apod t.apod_adjoint t.fft t.fft_adjoint t.conv t.conv_adjoint

  total = t.apod + t.fft + t.conv
  totalAdj = t.apod_adjoint + t.fft_adjoint + t.conv_adjoint
  @printf "                       apod = %.4f / %.4f %% fft = %.4f / %.4f %% conv = %.4f / %.4f %%\n" 100*t.apod/total 100*t.apod_adjoint/totalAdj 100*t.fft/total 100*t.fft_adjoint/totalAdj 100*t.conv/total 100*t.conv_adjoint/totalAdj

end

function reltolToParams(reltol)
  w = ceil(Int, log(10,1/reltol)) + 1
  m = (w)÷2
  return m, 2.0
end

function paramsToReltol(m::Int, σ)
  w = 2*m
  return 10.0^(-(w-1))
end

"""
  accuracyParams(; [m, σ, reltol]) -> m, σ, reltol

Calculate accuracy parameters m, σ, reltol based on either
* reltol
or
* m, σ

TODO: Right now, the oversampling parameter is not taken into account, i.e. σ=2.0 is assumed
"""
function accuracyParams(; kargs...)

  if haskey(kargs, :reltol)
    reltol = kargs[:reltol]
    m, σ = reltolToParams(reltol)
  elseif haskey(kargs, :m) && haskey(kargs, :σ)
    m = kargs[:m]
    σ = kargs[:σ]
    reltol = paramsToReltol(m, σ)
  else
    reltol = 1e-9
    m, σ = reltolToParams(reltol)
  end

  return m, σ, reltol
end


text:
# Test rank correlation

using StatsBase
using Base.Test

X = Float64[1 0; 2 1; 3 0; 4 1; 5 10]

x1 = X[:,1]
x2 = X[:,2]
y = [5, 3, 4, 2, 5]

# corspearman

@test_approx_eq corspearman(x1, y) -0.102597835208515
@test_approx_eq corspearman(x2, y) -0.081110710565381

@test_approx_eq corspearman(X, y) [-0.102597835208515, -0.081110710565381]
@test_approx_eq corspearman(y, X) [-0.102597835208515 -0.081110710565381]

c11 = corspearman(x1, x1)
c12 = corspearman(x1, x2)
c22 = corspearman(x2, x2)
@test_approx_eq c11 1.0
@test_approx_eq c22 1.0
@test_approx_eq corspearman(X, X) [c11 c12; c12 c22]
@test_approx_eq corspearman(X)    [c11 c12; c12 c22]

# corkendall

@test_approx_eq corkendall(x1, y) -0.105409255338946
@test_approx_eq corkendall(x2, y) -0.117851130197758

@test_approx_eq corkendall(X, y) [-0.105409255338946, -0.117851130197758]
@test_approx_eq corkendall(y, X) [-0.105409255338946 -0.117851130197758]

c11 = corkendall(x1, x1)
c12 = corkendall(x1, x2)
c22 = corkendall(x2, x2)

@test_approx_eq c11 1.0
@test_approx_eq c22 1.0
@test_approx_eq corkendall(X, X) [c11 c12; c12 c22]
@test_approx_eq corkendall(X)    [c11 c12; c12 c22]


text:
using RandomMatrices
using LinearAlgebra: I, tr
using Test

@testset "Haar" begin

N=5
A=randn(N,N)
B=randn(N,N)
Q=rand(Haar(1), N)

#Test case where there are no symbolic Haar matrices
@test_broken eval(expectation(:(A*B))) ≈ A*B
#Test case where there is one pair of symbolic Haar matrices
@test_broken tr(eval(expectation(:(A*Q*B*Q')))) ≈ tr(A)*tr(B)/N

println("Case 3")
@test_broken println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))

for elty in (Float64, ComplexF64)
	A = Stewart(elty, N)
        @test A'A ≈ Matrix{elty}(I, N, N)
end

end # testset


text:
using Koala
using Base.Test
using DataFrames

# loss functions:
y = [1, 2, 3, 4]
yhat = y + 1
@test isapprox(rms(y, yhat), 1.0)
@test isapprox(rmsl(y, yhat),
               sqrt((log(1/2)^2 + log(2/3)^2 + log(3/4)^2 + log(4/5)^2)/4))
@test isapprox(rmslp1(y, yhat),
               sqrt((log(2/3)^2 + log(3/4)^2 + log(4/5)^2 + log(5/6)^2)/4))
@test isapprox(rmsp(y, yhat), sqrt((1 + 1/4 + 1/9 + 1/16)/4))

# helpers:
dict = Dict{String,Int}()
dict["c"] = 3
dict["b"] = 2
dict["a"] = 1
@test keys_ordered_by_values(dict) == ["a", "b", "c"]
bootstrap_resample_of_mean(randn(100))
v = [1,2,5,4,missing,7,missing]
@test hasmissing(v)
@test ismissingtype(eltype(v))
@test Koala.leadingtype(eltype(v)) == Int
@test countmissing(v) == 2

v[5] = 42; v[7] = 42
@test eltype(purify(v)) == Int

X, y = load_boston();
train, test = split(eachindex(y), 0.8); # 80:20 split
describe(X)
get_meta(X)

transformer = Koala.FeatureSelector(features=[:Indus, :Chas])
transformerM = Machine(transformer, X)
@test transform(transformerM, X) == X[[:Indus, :Chas]]

rgs = ConstantRegressor()
mach = Machine(rgs, X, y, train, transformer_X=transformer)
showall(mach)
fit!(mach, train)
score = err(mach, test)
println("score = $score")
fit!(mach) # fit again without recomputing cache
@test 9.5 < score && 9.6 > score

@test score == err(mach, test, raw=true)
@test predict(mach, X, test) ==
    predict(mach.model, mach.predictor, mach.Xt[test,:], false, false)

learning_curve(mach, train, test, [2, 4, 8, 1000], raw=false)
learning_curve(mach, train, test, [2000, 3000], restart=false)
cv(mach, vcat(test, train))

u,v = @curve r linspace(0,10,50) (r^2 + 1)
u,v = @pcurve r linspace(0,10,50) (r^2 + 1)
u,v,w =@curve r linspace(0,10,5) s linspace(0,5,4) r*s^2

# test split_seen_unseen
v = ['a', 'b', 'b', 'c',
     'a', 'd', 'a', 'b', 'e']
trainrows = 1:4
testrows = 5:9
@test split_seen_unseen(v, trainrows, testrows) == ([5, 7, 8], [6, 9])

w = ["log", "house", "house", "house",
     "brick", "house", "log", "log", "log"]
df = DataFrame(v=v, w=w)
@test split_seen_unseen(df, trainrows, testrows) == ([7, 8], [5, 6, 9])

# time tests:
v=rand(UInt32, 10^7)
trainrows, testrows = split(eachindex(v), 0.9)
@time map(length, split_seen_unseen(v, trainrows, testrows))

# test drop_unseen capability
X[4] = map(Char, X[4])
X[9] = map(Char, X[9])
X[10] = map(Char, X[10])
train = 1:length(y) - 6
test = length(y) - 5 : length(y)
model = ConstantRegressor()
mach = Machine(model, X, y, train, drop_unseen=true)
fit!(mach, train)
err(mach, test)
learning_curve(mach, train, test, [2, 4, 8, 1000])
learning_curve(mach, train, test, [2000, 3000], restart=false)

t = Koala.RowsTransformer(10)
tM  = Machine(t, [2, 4, 5, 9])
@test transform(tM, 3:7) == [2, 3, 4]
df = DataFrame(letters=['a','b','c','d','e','f','g','h','i','j'])
df_shrunk = df[[1, 3, 6, 7, 8], :]
@test df_shrunk[transform(tM, 3:7), :letters] == ['c', 'f', 'g']

rgs = ConstantRegressor()
train, test = split(eachindex(y), 0.8); # 80:20 split
mach = Machine(rgs, X, y, train, drop_unseen=true)
showall(mach)
fit!(mach, train)
score = err(mach, test)
unseen = mach.rows_with_unseen
seen = filter(test) do i
    !(i in unseen)
end
@test mean(y[train]) ≈ mach.predictor
yhat = mean(y[train])*ones(length(y))
@test rms(yhat, y, seen) ≈ score
println("score = $score")
fit!(mach) # fit again without recomputing cache
@test 9 < score && 10 > score
@test score == err(mach, test, raw=true)
learning_curve(mach, train, test, [2, 4, 8, 1000], raw=false)
learning_curve(mach, train, test, [2000, 3000], restart=false)
cv(mach, vcat(test, train))

# test compete function:
@test compete(randn(1000) + 5, randn(1000))[1] == '1'

load_iris()
load_ames()


text:
using KitBase.PyCall

try
    using Conda
    Conda.add_channel("conda-forge")
    Conda.add("tensorflow")
    tf = pyimport("tensorflow")
catch
    cmd = `pip3 install tensorflow --user`
    run(cmd)
    tf = pyimport("tensorflow")
end


text:
uci_datasets = ["abalone", "blood-transfusion", "breast-cancer-wisconsin", "breast-tissue",
	"cardiotocography", "ecoli", "glass", "haberman", "ionosphere", "iris", "isolet",
	"letter-recognition", "libras", "magic-telescope", "miniboone", "multiple-features",
	"page-blocks", "parkinsons", "pendigits", "pima-indians", "sonar", "spect-heart", "statlog-satimage",
	"statlog-segment", "statlog-shuttle", "statlog-vehicle", "synthetic-control-chart",
	"wall-following-robot", "waveform-1", "waveform-2", "wine", "yeast"]
	# "gisette", "madelon"] - no easy anomalies + very large in size
	# "vertebral-column"] - no easy and medium anomalies

"""
	load_uci_data(dataset::String)

Loads basic UCI data. For a list of available datasets, check `GenerativeAD.Datasets.uci_datasets`.
"""
function load_uci_data(dataset::String)
	# I have opted for the original Loda datasets, use of multiclass problems in all vs one case
	# does not necessarily represent a good anomaly detection scenario
	data, _, _ = UCI.get_loda_data(dataset)
	# return only easy and medium anomalies
	UCI.normalize(data.normal, hcat(data.easy, data.medium)) # data (standardized)
end

other_datasets = ["annthyroid", "arrhythmia", "htru2", "kdd99", "kdd99_small", "spambase", "mammography", "har", "seismic"]
function __init__()
	register(
		DataDep(
			"annthyroid",
			"""
			Dataset: ANNThyroid
			Authors: Ross Quinlan
			Website: https://archive.ics.uci.edu/ml/datasets/Thyroid+Disease

			[Quinlan, 1987]
				Quinlan, John Ross, et al.
				"Inductive knowledge acquisition: a case study."
				Proceedings of the Second Australian Conference on Applications of expert systems. 1987.

			The original thyroid disease (ann-thyroid) dataset from UCI machine learning repository is
			a classification dataset, which is suited for training ANNs. It has 3772 training instances
			and 3428 testing instances. It has 15 categorical and 6 real attributes. The problem is to
			determine whether a patient referred to the clinic is hypothyroid. Therefore three classes
			are built: normal (not hypothyroid), hyperfunction and subnormal functioning.
			""",
			[
				"https://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/ann-train.data",
				"https://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/ann-test.data"
			],
			"fbecd00c2ab44e0c95fb135e9c2f39d597fb8074dfaa8581ac2f889c927d40ad"
		))
	# this ensures eager loading (upon first usage of the package)
	# datadep"annthyroid"
	# we dont do that no more

	register(
		DataDep(
			"arrhythmia",
			"""
			Dataset: Arrhythmia
			Authors: H. Altay Guvenir
			Website: https://archive.ics.uci.edu/ml/datasets/Arrhythmia

			[Guvenir, 1997]
				Guvenir, H. Altay, et al.
				"A supervised machine learning algorithm for arrhythmia analysis."
				Computers in Cardiology 1997. IEEE, 1997.

			The aim is to distinguish between the presence and absence of cardiac arrhythmia and to classify
			it in one of the 16 groups. Class 01 refers to 'normal' ECG classes 02 to 15 refers to different
			classes of arrhythmia and class 16 refers to the rest of unclassified ones. For the time being,
			there exists a computer program that makes such a classification. However there are differences
			between the cardiolog's and the programs classification. Taking the cardiolog's as a gold standard
			we aim to minimise this difference by means of machine learning tools.
			""",
			"https://archive.ics.uci.edu/ml/machine-learning-databases/arrhythmia/arrhythmia.data",
			"a7f0f4ca289a4c58b5ed85a9adb793189acd38425828ce3dfbb70adb45f30169"
		))

	register(
		DataDep(
			"htru2",
			"""
			Dataset: HTRU2
			Authors: Dr. Robert Lyon
			Website: https://archive.ics.uci.edu/ml/datasets/HTRU2

			[Keith, 2010]
				M. J. Keith et al.
				"The High Time Resolution Universe Pulsar Survey - I. System Configuration and Initial Discoveries."
				onthly Notices of the Royal Astronomical Society, 2010.

			HTRU2 is a data set which describes a sample of pulsar candidates collected during
			the High Time Resolution Universe Survey (South).
			""",
			"https://archive.ics.uci.edu/ml/machine-learning-databases/00372/HTRU2.zip",
			"ba442c076dd22a8952700f26e38499fc1806037dcf7bea0e125e6bfba393f379",
			post_fetch_method = unpack
		))

	register(
		DataDep(
			"kdd99",
			"""
			Dataset: KDD Cup 1999
			Authors: Stolfo et al.
			Website: http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html

			[Stolfo, 2000]
				 S. J. Stolfo, W. Fan, W. Lee, A. Prodromidis, and P. K. Chan.
				"Costbased modeling for fraud and intrusion detection: Results from the jam project."
				discex, 2000.

			This is the data set used for The Third International Knowledge Discovery and Data Mining Tools
			Competition, which was held in conjunction with KDD-99 The Fifth International Conference on
			Knowledge Discovery and Data Mining. The competition task was to build a network intrusion
			detector, a predictive model capable of distinguishing between ``bad'' connections, called
			intrusions or attacks, and ``good'' normal connections. This database contains a standard set
			of data to be audited, which includes a wide variety of intrusions simulated in a military
			network environment.

			!!!!!!!!!!!!!!!!!!!!

			WARNING --- for some users, the automatic unpacking may fail. Run this instead:

			mkdir ~/.julia/datadeps/kdd99
			cd ~/.julia/datadeps/kdd99
			wget http://kdd.ics.uci.edu/databases/kddcup99/kddcup.data.gz
			wget http://kdd.ics.uci.edu/databases/kddcup99/kddcup.data_10_percent.gz
			gzip -d kddcup.data.gz
			gzip -d kddcup.data_10_percent.gz

			!!!!!!!!!!!!!!!!!!!!
			""",
			["http://kdd.ics.uci.edu/databases/kddcup99/kddcup.data.gz",
			"http://kdd.ics.uci.edu/databases/kddcup99/kddcup.data_10_percent.gz"],
			"bb29388a787b1cea818a6fd427d14ad45b556176b5fbf13257ca33c1d2dad7f3",
			post_fetch_method = unpack
		))

	register(
		DataDep(
			"spambase",
			"""
			Dataset: Spambase
			Authors: Mark Hopkins, Erik Reeber, George Forman, Jaap Suermondt
			Website: https://archive.ics.uci.edu/ml/datasets/spambase

			Our collection of spam e-mails came from our postmaster and individuals who had filed spam.
			Our collection of non-spam e-mails came from filed work and personal e-mails, and hence
			the word 'george' and the area code '650' are indicators of non-spam. These are useful when
			constructing a personalized spam filter. One would either have to blind such non-spam indicators
			or get a very wide collection of non-spam to generate a general purpose spam filter.
			""",
			"https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.data",
			"b1ef93de71f97714d3d7d4f58fc9f718da7bbc8ac8a150eff2778616a8097b12"
		))

	register(
		DataDep(
			"mammography",
			"""
			Dataset: Mammography
			Authors: Tobias Kuehn
			Website: https://www.openml.org/d/310

			Mammography dataset.
			""",
			"https://www.openml.org/data/get_csv/52214/phpn1jVwe",
			"d39362668aa89b2b48aa0e6c2d12277850ace4c390e15b89862b7717f1525f0c"
			))

	register(
		DataDep(
			"har",
			"""
			Dataset: Human Activity Recognition
			Authors: Jorge L. Reyes
			Website: https://archive.ics.uci.edu/ml/datasets/human+activity+recognition+using+smartphones

			Human Activity Recognition dataset. Normal class = "WALKING", anomalous class = other.
			""",
			"https://archive.ics.uci.edu/ml/machine-learning-databases/00240/UCI%20HAR%20Dataset.zip",
			"2045e435c955214b38145fb5fa00776c72814f01b203fec405152dac7d5bfeb0",
			post_fetch_method = unpack
			))

	register(
		DataDep(
			"seismic",
			"""
			Dataset: Seismic activity
			Authors: Marek Sikora
			Website: https://archive.ics.uci.edu/ml/datasets/seismic-bumps

			Seismic bumps from mining activity.
			""",
			"https://archive.ics.uci.edu/ml/machine-learning-databases/00266/seismic-bumps.arff",
			"aabe512fab65b36d1dfb462650b75cfd8d99d8cc2723e8ecb4e6f5e1caccd5a7",
			))
end

"""
	init_tabular()

If not present, downloads tabular datasets.
"""
init_tabular() = map(d->(@datadep_str d), filter(x->x != "kdd99_small", other_datasets))

"""
	load_other_dataset(dataset[; standardize=false])

Load a dataset that is not a part of the UCI.jl package. For a list of available datasets, check
`GenerativeAD.Datasets.other_datasets`.
"""
function load_other_data(dataset; standardize = false)
	load_f = eval(Symbol("load_"*dataset))
	data_normal, data_anomalous = load_f()
	if standardize
		return UCI.normalize(data_normal, data_anomalous)
	else
		return data_normal, data_anomalous
	end
end

"""
	load_arrhythmia()
"""
function load_arrhythmia()
	data_path = datadep"arrhythmia"
	f = joinpath(data_path, readdir(data_path)[1])
	raw_data = readdlm(f, ',')

	# trow away features containing missing values
	proc_data = raw_data[:,filter(x->!(x in [11,12,13,14]), 1:size(raw_data,2))]
	# now throw away one observation with a missing value
	proc_data = proc_data[filter(x->!(x in [5]), 1:size(proc_data,1)), :]
	data = Array{Float32,2}(transpose(proc_data[:,1:end-1]))
	labels = proc_data[:,end]

	data[:, labels.==1], data[:, labels.!=1]
end

"""
	load_annthyroid()
"""
function load_annthyroid()
	data_path = datadep"annthyroid"
	fs = joinpath.(data_path, readdir(data_path))
	raw_data = readdlm.(fs, Float32)
	raw_data = Array(transpose(vcat(raw_data...)))
	labels = raw_data[end,:]
	data = raw_data[1:end-1,:]
	data[:, labels.==3], data[:, labels.!=3]
end

"""
	load_htru2()
"""
function load_htru2()
	data_path = datadep"htru2"
	f = joinpath(data_path, "HTRU_2.csv")
	raw_data = CSV.read(f, DataFrame, header = false)
	data = Array{Float32,2}(transpose(Array(raw_data[:,1:end-1])))
	labels = raw_data[:,end]
	data[:, labels.==0], data[:, labels.==1]
end

function _load_kdd99(f::String)
	raw_data = readdlm(f, ',')
	M,N = size(raw_data)
	labels = raw_data[:,end]
	# the 2nd, 3rd and 4th columns contain nonnumericalvalues and must be one hot encoded
	# 2nd col = 3 unique vals
	# 3 = 66
	# 4 = 11
	unss = [unique(raw_data[:,2]), unique(raw_data[:,3]), unique(raw_data[:,4])]
	ls = length.(unss)
	data = zeros(Float32,N+sum(ls)-4, M) # this is the final output array
	# copy first col
	data[1, :] = raw_data[:,1]
	# onehot encode the 3 following
	start_ind = 1
	for (uns,l,j) in zip(unss,ls,2:4)
		for i in 1:M
			data[(start_ind+1):(start_ind+l),i] = Flux.onehot(raw_data[i,j], uns)
		end
		start_ind = start_ind + l
	end
	# now copy the rest
	data[(sum(ls)+2):end,:] = transpose(raw_data[:,5:end-1])

	data[:, labels.=="normal."], data[:, labels.!="normal."]
end

"""
	load_kdd99_small()
"""
function load_kdd99_small()
	data_path = datadep"kdd99"
	f = joinpath(data_path, "kddcup.data_10_percent")
	_load_kdd99(f)
end

"""
	load_kdd99()
"""
function load_kdd99()
	data_path = datadep"kdd99"
	f = joinpath(data_path, "kddcup.data")
	_load_kdd99(f)
end

"""
	load_spambase()
"""
function load_spambase()
	data_path = datadep"spambase"
	f = joinpath(data_path, "spambase.data")
	raw_data = readdlm(f, ',', Float32)
	data = Array{Float32,2}(transpose(Array(raw_data[:,1:end-1])))
	labels = raw_data[:,end]
	data[:, labels.==0], data[:, labels.==1]
end

"""
	load_mammography()
"""
function load_mammography()
	data_path = datadep"mammography"
	f = joinpath(data_path, readdir(data_path)[1])
	raw_data = CSV.read(f, DataFrame)
	labels = parse.(Int,replace.(raw_data[!,:class], "'" => ""))
	data = Array{Float32,2}(transpose(Array(raw_data[:,1:end-1])))
	data[:, labels.==-1], data[:, labels.==1]
end

"""
	load_har()
"""
function load_har()
	data_path = datadep"har"
	X_test = readdlm(joinpath(data_path, "UCI HAR Dataset/test/X_test.txt"), Float32)
	X_train = readdlm(joinpath(data_path, "UCI HAR Dataset/train/X_train.txt"), Float32)
	y_test = readdlm(joinpath(data_path, "UCI HAR Dataset/test/y_test.txt"), Float32)
	y_train = readdlm(joinpath(data_path, "UCI HAR Dataset/train/y_train.txt"), Float32)

	data = Array(transpose(vcat(X_test, X_train)))
	labels = vec(vcat(y_test, y_train))
	data[:, labels.!=6], data[:, labels.==6] # other x laying
end

"""
	load_seismic()
"""
function load_seismic()
	data_path = datadep"seismic"
	f = joinpath(data_path, "seismic-bumps.arff")
	M,N = (2584,19)
	raw_data = Array{Any,2}(zeros(M,N))
	open(f) do file
		for (i,ln) in enumerate(eachline(file))
			(i > 154) ? (raw_data[i-154,:] = split(ln, ",")) : nothing
		end
	end
	labels = parse.(Int,raw_data[:,end])

	# the 1st, 2nd, 3rd and 8th columns contain nonnumericalvalues and must be one hot encoded
	unss = [unique(raw_data[:,1]), unique(raw_data[:,2]), unique(raw_data[:,3]), unique(raw_data[:,8])]
	ls = length.(unss)
	data = zeros(Float32,N+sum(ls)-5, M) # this is the final output array
	start_ind = 0
	for (uns,l,j) in zip(unss[1:3],ls[1:3],1:3)
		for i in 1:M
			data[(start_ind+1):(start_ind+l),i] = Flux.onehot(raw_data[i,j], uns)
		end
		start_ind = start_ind + l
	end

	# copy columns 4:7
	data[start_ind+1:start_ind+4,:] = transpose(parse.(Float32,raw_data[:,4:7]))

	# encode and fill col 8
	start_ind = start_ind+5
	uns = unss[4]
	l = ls[4]
	for i in 1:M
		data[start_ind:start_ind+l-1,i] = Flux.onehot(raw_data[i,8], uns)
	end

	# now copy the rest
	data[(sum(ls)+5):end,:] = transpose(parse.(Float32, raw_data[:,9:end-1]))

	data[:,labels.==0], data[:,labels.==1]
end


text:
using TicTacToe
using Test

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


text:
using EmpiricalBayes
using Distributions

m=LikelihoodModel(xs=collect(linspace(0,1,100)), phi=(x->x), measerr=Normal(0,.1), ndata=5, prior=Uniform(0,1), zmult = 10)

npmle(m)
refprior(m)
@time mple(m, 0.5)
@time mple(m, 0.5, config=OptConfig(METHOD=:auglag))
@time EmpiricalBayes.mpb(m, .5)
#@time EmpiricalBayes.jump(m, 0.5)

#em(m, 1000)
# plots


text:
#!/usr/bin/env julia

using OhMyJulia
using JSON

read_topo = true
node_list = String[]
node_dict = Dict{String, Int}()
topo = Set{Tuple{Int, Int}}()
demands = Tuple{Int, Int}[]

# requires switches and hosts numbers to be a 1 to 1 map
for line in readlines()
    if line == "***"
        read_topo = false
        continue
    end

    src, dst = split(line, ' ')
    src, dst = cdr(src), cdr(dst)
    src == dst && continue

    if read_topo
        for n in (src, dst)
            if n ∉ keys(node_dict)
                push!(node_list, n)
                node_dict[n] = length(node_list)
            end
        end
        push!(topo, (node_dict[src], node_dict[dst]))
    else
        push!(demands, (node_dict[src], node_dict[dst]))
    end
end

nnode = length(node_list)

if length(demands) == 0
    for i in 1:nnode, j in 1:nnode @when i != j
        push!(demands, (i, j))
    end
end

const lib = rel"a.out"

groups = []

tic()

for (src, dst) in demands
    mf = ccall((:init, lib), *(Void), (Cint,), nnode)

    for (s, d) in topo @when s < d
        ccall((:addEdge, lib), Void, (*(Void), Cint, Cint), mf, s, d)
    end

    max_flow = ccall((:dinitz, lib), Cint, (*(Void), Cint, Cint), mf, src, dst)
    npath = ccall((:npath, lib), Cint, (*(Void),), mf)

    group = []

    for i in 1:npath
        len = ccall((:path_len, lib), Cint, (*(Void), Cint), mf, i)
        path = Vector{Cint}(len)
        ccall((:get_path, lib), Void, (*(Void), Cint, *(Cint)), mf, i, path)
        push!(group, [src] ++ path)
    end

    push!(groups, group)
end

println(STDERR, "TIME: ", toq())

# println(STDERR, groups)

to_pair(x) = "h$(node_list[x[][]]) h$(node_list[x[][end]])"

to_path(x) = begin
    path = ["s$(node_list[x[i]]) s$(node_list[x[i+1]])" for i in 1:length(x)-1]
    hsrc = "h$(node_list[x[]]) s$(node_list[x[]])"
    hdst = "s$(node_list[x[end]]) h$(node_list[x[end]])"
    [hsrc, path..., hdst]
end

println(JSON.json(Dict( to_pair(group) => map(to_path, group) for group in groups )))

max_stars_repo_name: ylxdzsw/TEexp

max_issues_repo_name: ylxdzsw/TEexp

max_forks_repo_name: ylxdzsw/TEexp


text:
### A Pluto.jl notebook ###
# v0.12.11

using Markdown
using InteractiveUtils

# ╔═╡ 2ebffa2a-feac-11ea-27d0-3b7cf030672a
using Pkg, DrWatson

# ╔═╡ 2ec03170-feac-11ea-3cd8-e9c1d16cd34c
begin
  #@quickactivate "StatisticalRethinkingStan"
  using StanSample
  using StatisticalRethinking
end

# ╔═╡ 37e0328a-feab-11ea-031f-df82e9c87154
md"## Clip-06-13-15s.jl"

# ╔═╡ 2ec0c25c-feac-11ea-2439-07e617727eb5
begin
  N = 100
  df = DataFrame(
    :h0 => rand(Normal(10,2 ), N),
    :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
  );
  df.fungus = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
  df.h1 = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]
  scale!(df, [:h0, :treatment, :fungus, :h1])
end;

# ╔═╡ 2ed0d892-feac-11ea-1ad3-5dfeacd40f38
Text(precis(df, io=String))

# ╔═╡ 2ed17b7e-feac-11ea-2d49-3ba0f6783ba4
sim_p = DataFrame(:sim_p => rand(LogNormal(0, 0.25), 10000));

# ╔═╡ 2edbc4b2-feac-11ea-0977-37cb5da52497
Text(precis(sim_p; io=String))

# ╔═╡ 2edc81fe-feac-11ea-3f3c-2ff7b2c13458
stan6_6 = "
data {
  int <lower=1> N;
  vector[N] h0;
  vector[N] h1;
}
parameters{
  real<lower=0> p;
  real<lower=0> sigma;
}
model {
  vector[N] mu;
  p ~ lognormal(0, 0.25);
  sigma ~ exponential(1);
  mu = h0 * p;
  h1 ~ normal(mu, sigma);
}
";

# ╔═╡ 2ee7bd6a-feac-11ea-1120-776954c32bd8
begin
	m6_6s = SampleModel("m6.6s", stan6_6)
	m6_6_data = Dict(:N => nrow(df), :h0 => df.h0, :h1 => df.h1)
	rc6_6s = stan_sample(m6_6s; data=m6_6_data)
	success(rc6_6s) && (post6_6s_df = read_samples(m6_6s, :dataframe))
end;

# ╔═╡ 2ef5aa9c-feac-11ea-1da4-67356630f549
if success(rc6_6s)
	part6_6s = Particles(post6_6s_df)
end

# ╔═╡ 2f01856c-feac-11ea-0b6b-537d35689cfc
success(rc6_6s) && (Text(precis(post6_6s_df; io=String)))

# ╔═╡ 2f093c3a-feac-11ea-2325-37124b7e6bdf
md"## End of clip-06-13-15s.jl"

# ╔═╡ Cell order:
# ╠═37e0328a-feab-11ea-031f-df82e9c87154
# ╠═2ebffa2a-feac-11ea-27d0-3b7cf030672a
# ╠═2ec03170-feac-11ea-3cd8-e9c1d16cd34c
# ╠═2ec0c25c-feac-11ea-2439-07e617727eb5
# ╠═2ed0d892-feac-11ea-1ad3-5dfeacd40f38
# ╠═2ed17b7e-feac-11ea-2d49-3ba0f6783ba4
# ╠═2edbc4b2-feac-11ea-0977-37cb5da52497
# ╠═2edc81fe-feac-11ea-3f3c-2ff7b2c13458
# ╠═2ee7bd6a-feac-11ea-1120-776954c32bd8
# ╠═2ef5aa9c-feac-11ea-1da4-67356630f549
# ╠═2f01856c-feac-11ea-0b6b-537d35689cfc
# ╟─2f093c3a-feac-11ea-2325-37124b7e6bdf


text:
# Addition operator

Base.:(+)(::AbstractSet, d::SetDifference) = union(d.comvals, d.addvals)
Base.:(+)(d::SetDifference, a::AbstractSet) = +(a, d)

function Base.:(+)(a::AbstractDict, d::DictDifference)
    result = d.addvals
    modkeys = keys(d.modvals)
    for k in modkeys
        typeof(a[k]) != typeof(d.modvals[k]) &&
            throw(ArgumentError("type of values of common keys in `a` and `d` must match"))
        result = merge(result, Dict(k=>a[k] + d.modvals[k]))
    end
    return result
end
Base.:(+)(d::DictDifference, a::AbstractDict) = +(a, d)


text:
# descriptor

mutable struct FilterDesc
  ptr::cudnnFilterDescriptor_t
end

unsafe_free!(fd::FilterDesc) = cudnnDestroyFilterDescriptor(fd.ptr)

Base.unsafe_convert(::Type{cudnnFilterDescriptor_t}, fd::FilterDesc) = fd.ptr

function createFilterDesc()
  d = Ref{cudnnFilterDescriptor_t}()
  cudnnCreateFilterDescriptor(d)
  return d[]
end

function FilterDesc(T::Type, size::Tuple; format = CUDNN_TENSOR_NCHW)
    # The only difference of a FilterDescriptor is no strides.
    sz = Cint.(size) |> reverse |> collect
    d = createFilterDesc()
    version() >= v"5" ?
        cudnnSetFilterNdDescriptor(d, cudnnDataType(T), format, length(sz), sz) :
    version() >= v"4" ?
        cudnnSetFilterNdDescriptor_v4(d, cudnnDataType(T), format, length(sz), sz) :
        cudnnSetFilterNdDescriptor(d, cudnnDataType(T), length(sz), sz)
    this = FilterDesc(d)
    finalizer(unsafe_free!, this)
    return this
end

FilterDesc(a::CuArray; format = CUDNN_TENSOR_NCHW) = FilterDesc(eltype(a), size(a), format = format)

function Base.size(f::FilterDesc)
  typ = Ref{Cuint}()
  format = Ref{Cuint}()
  ndims = Ref{Cint}()
  dims = Vector{Cint}(undef, 8)
  cudnnGetFilterNdDescriptor(f, 8, typ, format, ndims, dims)
  @assert ndims[] ≤ 8
  return (dims[1:ndims[]]...,) |> reverse
end


text:
mutable struct EstimatedLandmark <: AbstractLandmark
    pos_::Vector{Float64}
    id::Int64
    cov_::Union{Matrix{Float64},Nothing}
    function EstimatedLandmark(id::Int64)
        # this cov_ initial value is for debugging
        new([0.0, 0.0], id, nothing)
    end
end

#=
function Base.copy(lm::EstimatedLandmark)
    if lm.cov_ == nothing
        return EstimatedLandmark(copy(lm.pos_), lm.id, nothing)
    else
        return EstimatedLandmark(copy(lm.pos_), lm.id, copy(lm.cov_))
    end
end
=#

function draw(mark::EstimatedLandmark, p::Plot{T}) where {T}
    if mark.cov_ == nothing
        return
    else
        p = scatter!(
            [mark.pos_[1]],
            [mark.pos_[2]],
            markershape = :star,
            markersize = 10,
            color = "blue",
        )
        p = annotate!(
            mark.pos_[1] + 0.5,
            mark.pos_[2] + 0.5,
            text("id: $(mark.id)", 10, :blue),
        )
        p = myellipse!(mark.pos_, mark.cov_[1:2, 1:2], n_std = 3, aspect_ratio = 1)
    end
end

mutable struct MapParticle <: AbstractParticle
    pose_::Vector{Float64}
    weight_::Float64
    map_::Map
    function MapParticle(init_pose::Vector{Float64}, weight::Float64, landmark_num::Int64)
        envmap = Map()
        for i = 1:landmark_num
            push!(envmap, EstimatedLandmark(i - 1))
        end
        new(copy(init_pose), weight, envmap)
    end
end

function Base.copy(p::MapParticle)
    p_ = MapParticle(copy(p.pose_), p.weight_, 0)
    p_.map_ = copy(p.map_)
    return p_
end

function motion_update2(
    particle::MapParticle,
    v::Float64,
    ω::Float64,
    dt::Float64,
    motion_noise_stds::Dict{String,Float64},
    observation::Vector{Vector{Float64}},
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    M = matM(v, ω, dt, motion_noise_stds)
    A = matA(v, ω, dt, particle.pose_[3])
    Rₜ = A * M * transpose(A)
    x̂ = state_transition(particle.pose_, v, ω, dt)

    N = size(observation)[1]
    for i = 1:N
        z = observation[i][1:2]
        idx = convert(Int64, observation[i][3])
        zₜ̂, Qzₜ, Hxₜ = drawing_params(
            x̂,
            particle.map_.landmarks_[idx+1],
            distance_dev_rate,
            direction_dev,
        )
        Σzₜ = Hxₜ * Rₜ * transpose(Hxₜ) + Qzₜ
        Σzₜ = (Σzₜ + transpose(Σzₜ)) / 2.0
        particle.weight_ = pdf(MvNormal(zₜ̂, Σzₜ), z)
    end

    for i = 1:N
        obsv = observation[i]
        z = obsv[1:2]
        idx = convert(Int64, obsv[3])
        x̂, Rₜ = gauss_for_drawing(
            x̂,
            Rₜ,
            z,
            particle.map_.landmarks_[idx+1],
            distance_dev_rate,
            direction_dev,
        )
    end

    Rₜ = (Rₜ + transpose(Rₜ)) / 2.0
    particle.pose_ = rand(MvNormal(x̂, Rₜ + Matrix(1.0e-10 * I, 3, 3)))
end

function drawing_params(
    x̂::Vector{Float64},
    lm::EstimatedLandmark,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    d = hypot((x̂[1:2] - lm.pos_)...)
    Qzₜ̂ = matQ(distance_dev_rate * d, direction_dev)
    zₜ̂ = observation_function(x̂, lm.pos_)
    Hm = -matH(x̂, lm.pos_)[1:2, 1:2]
    Hxₜ = matH(x̂, lm.pos_)

    Qzₜ = Hm * lm.cov_ * transpose(Hm) + Qzₜ̂

    return zₜ̂, Qzₜ, Hxₜ
end

function gauss_for_drawing(
    x̂::Vector{Float64},
    Rₜ::Matrix{Float64},
    z::Vector{Float64},
    lm::EstimatedLandmark,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    zₜ̂, Qzₜ, Hxₜ = drawing_params(x̂, lm, distance_dev_rate, direction_dev)
    K = Rₜ * transpose(Hxₜ) * inv(Qzₜ + Hxₜ * Rₜ * transpose(Hxₜ))

    return K * (z - zₜ̂) + x̂, (Matrix(1.0I, 3, 3) - K * Hxₜ) * Rₜ
end

function observation_update(
    particle::MapParticle,
    observation::Vector{Vector{Float64}},
    envmap::Map,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    N = size(observation)[1]
    for i = 1:N
        obsv = observation[i]
        d, ϕ = obsv[1], obsv[2]
        lm_id = convert(Int64, obsv[3])
        landmark = particle.map_.landmarks_[lm_id+1]
        if landmark.cov_ == nothing
            init_landmark_estimation(
                particle,
                landmark,
                d,
                ϕ,
                distance_dev_rate,
                direction_dev,
            )
        else
            observation_update_landmark(
                particle,
                landmark,
                d,
                ϕ,
                distance_dev_rate,
                direction_dev,
            )
        end
    end
end

function init_landmark_estimation(
    particle::MapParticle,
    lm::EstimatedLandmark,
    d::Float64,
    ϕ::Float64,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    θ = particle.pose_[3]
    lm.pos_ = particle.pose_[1:2] + d * [cos(ϕ + θ), sin(ϕ + θ)]
    H = matH(particle.pose_, lm.pos_)[1:2, 1:2]
    Q = matQ(distance_dev_rate * d, direction_dev)
    lm.cov_ = inv(transpose(H) * inv(Q) * H)
end

function observation_update_landmark(
    particle::MapParticle,
    lm::EstimatedLandmark,
    d::Float64,
    ϕ::Float64,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    # estimated position [d, ϕ] from estimated position of landmark
    estm_z = observation_function(particle.pose_, lm.pos_)
    if estm_z[1] < 0.01
        # too close
        return
    end

    H = -matH(particle.pose_, lm.pos_)[1:2, 1:2]
    Q = matQ(distance_dev_rate * estm_z[1], direction_dev)
    K = lm.cov_ * transpose(H) * inv(Q + H * lm.cov_ * transpose(H))

    # update weight
    # Q_z = H * lm.cov_ * transpose(H) + Q
    # Q_z = (Q_z + transpose(Q_z)) / 2.0
    # particle.weight_ *= pdf(MvNormal(estm_z, Q_z), [d, ϕ])

    # update landmark position esimate
    lm.pos_ = K * ([d, ϕ] - estm_z) + lm.pos_
    lm.cov_ = (Matrix(1.0I, 2, 2) - K * H) * lm.cov_
end

mutable struct FastSlam1 <: AbstractMcl
    particles_::Vector{MapParticle}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::MapParticle
    pose_::Vector{Float64}
    function FastSlam1(
        init_pose::Vector{Float64},
        particle_num::Int64,
        landmark_num::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05,
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])
        new(
            [
                MapParticle(init_pose, 1.0 / particle_num, landmark_num) for
                i = 1:particle_num
            ],
            MvNormal([0.0, 0.0, 0.0, 0.0], cov),
            distance_dev_rate,
            direction_dev,
            MapParticle(init_pose, 1.0, 0),
            init_pose,
        )
    end
end

mutable struct FastSlam2 <: AbstractMcl
    particles_::Vector{MapParticle}
    motion_noise_stds::Dict{String,Float64}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::MapParticle
    pose_::Vector{Float64}
    function FastSlam2(
        init_pose::Vector{Float64},
        particle_num::Int64,
        landmark_num::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05,
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])
        new(
            [
                MapParticle(init_pose, 1.0 / particle_num, landmark_num) for
                i = 1:particle_num
            ],
            copy(v),
            MvNormal([0.0, 0.0, 0.0, 0.0], cov),
            distance_dev_rate,
            direction_dev,
            MapParticle(init_pose, 1.0, 0),
            init_pose,
        )
    end
end

function motion_update(slam::FastSlam1, v::Float64, ω::Float64, dt::Float64; kwargs...)
    N = length(slam.particles_)
    for i = 1:N
        motion_update(slam.particles_[i], v, ω, dt, slam.motion_noise_rate_pdf)
    end
end

function motion_update(slam::FastSlam2, v::Float64, ω::Float64, dt::Float64; kwargs...)
    observation = kwargs[:observation] # Vector{Vector{Float64}}
    not_first_observation = Vector{Vector{Float64}}(undef, 0)
    M = size(observation)[1]
    for i = 1:M
        obsv = observation[i]
        idx = convert(Int64, obsv[3])
        if slam.particles_[1].map_.landmarks_[idx+1].cov_ != nothing
            push!(not_first_observation, obsv)
        end
    end

    N = length(slam.particles_)
    if size(not_first_observation)[1] > 0
        for i = 1:N
            motion_update2(
                slam.particles_[i],
                v,
                ω,
                dt,
                slam.motion_noise_stds,
                not_first_observation,
                slam.distance_dev_rate,
                slam.direction_dev,
            )
            #motion_update(slam.particles_[i], v, ω, dt, slam.motion_noise_rate_pdf)
        end
    else
        for i = 1:N
            motion_update(slam.particles_[i], v, ω, dt, slam.motion_noise_rate_pdf)
        end
    end
end

function observation_update(
    slam::Union{FastSlam1,FastSlam2},
    observation::Vector{Vector{Float64}},
    envmap::Map;
    kwargs...,
)
    resample = true
    if haskey(kwargs, :resample)
        resample = convert(Bool, kwargs[:resample])
    end
    # currently same as the one for mcl::Mcl
    N = length(slam.particles_)
    for i = 1:N
        observation_update(
            slam.particles_[i],
            observation,
            envmap,
            slam.distance_dev_rate,
            slam.direction_dev,
        )
    end
    set_ml(slam)
    if resample
        resampling(slam)
    end
end

function draw(slam::Union{FastSlam1,FastSlam2}, p::Plot{T}) where {T}
    xs = [p.pose_[1] for p in slam.particles_]
    ys = [p.pose_[2] for p in slam.particles_]
    vxs = [
        cos(p.pose_[3]) * 0.5 * p.weight_ * length(slam.particles_) for
        p in slam.particles_
    ]
    vys = [
        sin(p.pose_[3]) * 0.5 * p.weight_ * length(slam.particles_) for
        p in slam.particles_
    ]
    p = quiver!(xs, ys, quiver = (vxs, vys), color = "blue", alpha = 0.5)
    pose = slam.pose_
    annota = "($(round(pose[1], sigdigits=3)), $(round(pose[2], sigdigits=3)), $(round(pose[3], sigdigits=3)))"
    p = annotate!(pose[1] + 1.0, pose[2] + 1.0, text(annota, 10))
    # draw the belief map of maximum likelihood particle
    draw(slam.ml_.map_, p)
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
if abspath(PROGRAM_FILE) == @__FILE__
    values = Vector{UInt8}()
    @assert(isa(values, bytearray) == true)
end


text:
using DataDeps

register(DataDep(
  "breast-cancer-wisconsin-prognostic",
  "http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+%28Prognostic%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wpbc.data",
  "f2bd5ccdb6086fa832370dc043319907430cf6bfcfd3b96a3961a357be345c8e",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=2,
        feature_indices=4:35,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# This file contains methods that are dispatched to when using neural networks
# with a preallocated cache (called CachedNetworks). Those are used to greatly
# improve performance, especially in

export cached, weights, grad_cache

"""
    NNCache{N}

The base abstract type holding a cache for the network type `N<:NeuralNetwork`.
"""
abstract type NNCache{N}
    #valid::Bool
end

"""
    invalidate_cache!(net)

Invalidates the cache of a cached Neural Network. Usually called when weights
are changed to signal the fact that lookup tables must be cleared.
"""
@inline invalidate_cache!(net::NNCache) = net.valid = false
@inline invalidate_cache!(net::Nothing) = nothing # if not implemented

"""
    CachedNet{N, NC}

A Cached Version of Neural Network `N` using Cache `NC`. It behaves the same as
the network `N`, but uses a preallocated cache to speed up computation of the
gradient and the function.
"""
struct CachedNet{N, NC} <: NeuralNetwork
    net::N
    cache::NC
end

trainable(cnet::CachedNet) = trainable(cnet.net)

"""
    cached(net)

Constructs a cached version of the neural network `net`, preallocating several
arrays for intermediate results for improved performance. This object will
behave identically to a standard network `net`.
"""
cached(net::NeuralNetwork) = CachedNet(net, cache(net))
cached(net::CachedNet) = CachedNet(net.net, cache(net))

# When copying shallow copy the net but deepcopy the der_vec
"""
    copy(cnet::CachedNet)

Copy a cached network, building a shallow copy of the network and a deep-copy
of the cache.
"""
Base.copy(cnet::CachedNet) = CachedNet(cnet.net, deepcopy(cnet.cache))

"""
    cache(net)

Constructs the `NNCache{typeof(net)}` object that holds the cache for this network.
If it has not been implemented returns nothing.
"""
cache(net) = nothing
cache(net::CachedNet) = cache(net.net)

"""
    weights(net)

Access the underlying structure holding the weights of the network. This is usefull
when dealing with Symmetrized networks or cached networks, as it gives access to
the wrapped data structure.
"""
weights(net) = trainable(net)

@specialize_vararg 5 @inline (cnet::CachedNet)(σ...) = logψ(cnet, σ...)
# When you call logψ on a cached net use the cache to compute the net
@inline logψ(cnet::CachedNet, σ::NTuple{N,<:AbstractVector}) where N =
    cnet.net(cnet.cache, σ...)
@inline logψ(cnet::CachedNet, σ::ADoubleState) where {N,V} =
    cnet.net(cnet.cache, σ...)
@inline logψ(cnet::CachedNet, σ::AState) where N =
    cnet.net(cnet.cache, σ)
@inline logψ(cnet::CachedNet, σr::T, σc::T) where {T<:Union{AState}} =
    cnet.net(cnet.cache, σr, σc)

function logψ_and_∇logψ(n::CachedNet, σ::Vararg{N,V}) where {N,V}
    #@warn "Inefficient calling logψ_and_∇logψ for cachedNet"
    ∇lnψ   = grad_cache(n)
    lψ, _g = logψ_and_∇logψ!(∇lnψ, n, σ...);
    return (lψ, ∇lnψ)
end

# Declare the two functions, even if config(blabla)=blabla, because of a shitty
# Julia's performance bug #32761
# see https://github.com/JuliaLang/julia/issues/32761
@inline function logψ_and_∇logψ!(der, n::CachedNet, σ::NTuple{N,AbstractArray}) where N
    lψ = logψ_and_∇logψ!(der, n.net, n.cache, σ...)
    return (lψ, der)
end
@inline function logψ_and_∇logψ!(der, n::CachedNet, σ::Vararg{<:AbstractArray,N}) where N
    lψ = logψ_and_∇logψ!(der, n.net, n.cache, σ...);
    return (lψ, der)
end
# this is required for networks with only one state (KetNet)
@inline function logψ_and_∇logψ!(der, n::CachedNet, σ::AbstractArray)
    lψ = logψ_and_∇logψ!(der, n.net, n.cache, σ);
    return (lψ, der)
end

## Inplace stuff
function logψ!(out, net::NeuralNetwork, cache::NNCache, vv::Union{AStateBatch,ADoubleStateBatch})
    for i=1:num_batches(vv)
        v = state_uview(vv, i)
        out[i] = logψ(net, cache, v)
    end
    return out
end

function logψ!(out, net::NeuralNetwork, cache::NNCache, vvl::AStateBatch, vvr::AStateBatch)
    for i=1:num_batches(vvr)
        vl = state_uview(vvl, i)
        vr = state_uview(vvr, i)
        out[i] = logψ(net, cache, vl, vr)
    end
    return out
end

## Optimisation of cachednet
update!(opt, cnet::CachedNet, Δ, state=nothing) = begin
    update!(opt, weights(cnet), weights(Δ), state)
    invalidate_cache!(cnet.cache)
    return nothing
end

apply!(opt, val1::Union{NeuralNetwork, CachedNet}, val2::Union{NeuralNetwork, CachedNet}, args...) =
    apply!(weights(val1), weights(val2), args...)

# The common operations are forwarded to the underlying network.
out_type(cnet::CachedNet) = out_type(cnet.net)
is_analytic(cnet::CachedNet) = is_analytic(cnet.net)
num_params(cnet::CachedNet) = num_params(cnet.net)

# Visualization
function Base.show(io::IO, m::CachedNet)
    print(io, "CachedNet{$(m.net)}")
    try
        bsz = batch_size(m)
        print(io, ": batch_size = $bsz")
    catch
    end
end

function Base.show(io::IO, ::MIME"text/plain", m::CachedNet)
    print(io, "CachedNet{$(m.net)}")
    try
        bsz = batch_size(m)
        print(io, ": batch_size = $bsz")
    catch
    end
end

## Particoular kinds of networks
# Thing for Matrix Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
const MatrixNet   = Union{MatrixNeuralNetwork, CachedNet{<:MatrixNeuralNetwork}}
# Diagonal
@inline (net::MatrixNeuralNetwork)(cache::NNCache, σ::AState) = net(cache, σ, σ)
@inline log_prob_ψ(net::CachedNet{<:MatrixNeuralNetwork},
                    σ::ADoubleStateOrBatchOrVec) = 2.0*real(net(σ))
# For Matrix network evaluated on diagonal
@inline log_prob_ψ(net::CachedNet{<:MatrixNeuralNetwork},
                        σ::AStateOrBatchOrVec) = real(net(σ))

# Thing for Pure Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
const PureNet   = Union{PureNeuralNetwork, CachedNet{<:PureNeuralNetwork}}

# Thing for Pure-state (Closed system) Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
const KetNet   = Union{KetNeuralNetwork, CachedNet{<:KetNeuralNetwork}}

# This overrides the standard behaviour of net(σ...) because Vector unpacking
# should not happen
#@inline logψ(cnet::CachedNet{<:KetNeuralNetwork}, σ) = cnet.net(cnet.cache, config(σ))

function logψ_and_∇logψ(n::CachedNet{<:KetNeuralNetwork}, σ::AbstractVector)
    if σ isa AbstractVector
        ∇lnψ = grad_cache(n)
    elseif σ isa AbstractMatrix
        ∇lnψ = grad_cache(n, batch_size(σ))
    end
    lnψ = logψ_and_∇logψ!(∇lnψ, n.net, n.cache, σ);
    return (lnψ, ∇lnψ)
end


text:
"""
# Row-merging (vcat/rbind)

```julia
vcat(ts1::TS, ts2::TS; colmerge::Symbol=:union)
```

Concatenate rows of two TS objects, append `ts2` to `ts1`.

The `colmerge` keyword argument specifies the column merge
strategy. The value of `colmerge` is directly passed to `cols`
argument of `DataFrames.vcat`.

Currently, `DataFrames.vcat` supports four types of column-merge strategies:

1. `:setequal`: only merge if both objects have same column names, use the order of columns in `ts1`.
2. `:orderequal`: only merge if both objects have same column names and columns are in the same order.
3. `:intersect`: only merge the columns which are common to both objects, ignore the rest.
4. `:union`: merge even if columns differ, the resulting object has all the columns filled with `missing`, if necessary.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;

julia> random(x) = rand(MersenneTwister(123), x);

julia> dates1 = collect(Date(2017,1,1):Day(1):Date(2017,1,10));

julia> dates2 = collect(Date(2017,1,11):Day(1):Date(2017,1,30));

julia> ts1 = TS([randn(length(dates1)) randn(length(dates1))], dates1)
julia> show(ts1)
(10 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
────────────────────────
 2017-01-01  -0.420348
 2017-01-02   0.109363
 2017-01-03  -0.0702014
 2017-01-04   0.165618
 2017-01-05  -0.0556799
 2017-01-06  -0.147801
 2017-01-07  -2.50723
 2017-01-08  -0.099783
 2017-01-09   0.177526
 2017-01-10  -1.08461

julia> df = DataFrame(x1 = randn(length(dates2)), y1 = randn(length(dates2)))
julia> ts2 = TS(df, dates2)
julia> show(ts2)
(20 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
────────────────────────
 2017-01-11   2.15087
 2017-01-12   0.9203
 2017-01-13  -0.0879142
 2017-01-14  -0.930109
 2017-01-15   0.061117
 2017-01-16   0.0434627
 2017-01-17   0.0834733
 2017-01-18  -1.52281
     ⋮           ⋮
 2017-01-23  -0.756143
 2017-01-24   0.491623
 2017-01-25   0.549672
 2017-01-26   0.570689
 2017-01-27  -0.380011
 2017-01-28  -2.09965
 2017-01-29   1.37289
 2017-01-30  -0.462384
          4 rows omitted

julia> vcat(ts1, ts2)
(30 x 3) TS with Date Index

 Index       x1          x2              y1
 Date        Float64     Float64?        Float64?
─────────────────────────────────────────────────────────
 2017-01-01  -0.524798        -1.4949    missing
 2017-01-02  -0.719611        -1.1278    missing
 2017-01-03   0.0926092        1.19778   missing
 2017-01-04   0.236237         1.39115   missing
 2017-01-05   0.369588         1.21792   missing
 2017-01-06   1.65287         -0.930058  missing
 2017-01-07   0.761301         0.23794   missing
 2017-01-08  -0.571046        -0.480486  missing
 2017-01-09  -2.01905         -0.46391   missing
 2017-01-10   0.193942        -1.01471   missing
 2017-01-11   0.239041   missing              -0.473429
 2017-01-12   0.286036   missing              -0.90377
 2017-01-13   0.683429   missing              -0.128489
 2017-01-14  -1.51442    missing              -2.39843
 2017-01-15  -0.581341   missing              -0.12265
 2017-01-16   1.07059    missing              -0.916064
 2017-01-17   0.859396   missing               0.0162969
 2017-01-18  -1.93127    missing               2.11127
 2017-01-19   0.529477   missing               0.636964
 2017-01-20   0.817429   missing              -0.34038
 2017-01-21  -0.682296   missing              -0.971262
 2017-01-22   1.36232    missing              -0.236323
 2017-01-23   0.143188   missing              -0.501722
 2017-01-24   0.621845   missing              -1.20016
 2017-01-25   0.076199   missing              -1.36616
 2017-01-26   0.379672   missing              -0.555395
 2017-01-27   0.494473   missing               1.05389
 2017-01-28   0.278259   missing              -0.358983
 2017-01-29   0.0231765  missing               0.712526
 2017-01-30   0.516704   missing               0.216855

julia> vcat(ts1, ts2; colmerge=:intersect)
(30 x 1) TS with Date Index

 Index       x1
 Date        Float64
────────────────────────
 2017-01-01  -0.524798
 2017-01-02  -0.719611
 2017-01-03   0.0926092
 2017-01-04   0.236237
 2017-01-05   0.369588
 2017-01-06   1.65287
 2017-01-07   0.761301
 2017-01-08  -0.571046
 2017-01-09  -2.01905
 2017-01-10   0.193942
 2017-01-11   0.239041
 2017-01-12   0.286036
 2017-01-13   0.683429
 2017-01-14  -1.51442
 2017-01-15  -0.581341
 2017-01-16   1.07059
 2017-01-17   0.859396
 2017-01-18  -1.93127
 2017-01-19   0.529477
 2017-01-20   0.817429
 2017-01-21  -0.682296
 2017-01-22   1.36232
 2017-01-23   0.143188
 2017-01-24   0.621845
 2017-01-25   0.076199
 2017-01-26   0.379672
 2017-01-27   0.494473
 2017-01-28   0.278259
 2017-01-29   0.0231765
 2017-01-30   0.516704

```
"""
function Base.vcat(ts1::TS, ts2::TS; colmerge=:union)
    result_df = DataFrames.vcat(ts1.coredata, ts2.coredata; cols=colmerge)
    return TS(result_df)
end
# alias
rbind = vcat


text:
const QNIndexSet = IndexSet{QNIndex}

const QNIndices = Union{QNIndexSet,Tuple{Vararg{QNIndex}}}

# Get a list of the non-zero blocks given a desired flux
# TODO: make a fillqns(inds::Indices) function that makes all indices
# in inds have the same qns. Then, use a faster comparison:
#   ==(flux(inds,block; assume_filled=true), qn; assume_filled=true)
function nzblocks(qn::QN, inds::Indices)
  N = length(inds)
  blocks = Block{N}[]
  for block in eachblock(inds)
    if flux(inds, block) == qn
      push!(blocks, block)
    end
  end
  return blocks
end

function nzdiagblocks(qn::QN, inds::Indices)
  N = length(inds)
  blocks = NTuple{N,Int}[]
  for block in eachdiagblock(inds)
    if flux(inds, block) == qn
      push!(blocks, Tuple(block))
    end
  end
  return blocks
end

removeqns(is::QNIndices) = map(i -> removeqns(i), is)

anyfermionic(is::Indices) = any(isfermionic, is)

allfermionic(is::Indices) = all(isfermionic, is)

# TODO: specialize for Vector vs. Tuple
function Base.show(io::IO, is::QNIndices)
  print(io, "$(length(is)) Indices\n")
  for i in is
    print(io, i)
    println(io)
  end
end


text:
"""
    mutable struct Context
A contextual data object.
"""
Base.@kwdef mutable struct Context
    id::Int=new_id()
    Window::Ptr{GLFWwindow} = Ptr{GLFWwindow}(C_NULL)
    ClientApi::GlfwClientApi = GlfwClientApi_Unknown
    Time::Cfloat = Cfloat(0.0f0)
    MouseJustPressed::Vector{Bool} = fill(false, Int(ImGuiMouseButton_COUNT))
    MouseCursors::Vector{Ptr{GLFWcursor}} = fill(Ptr{GLFWcursor}(C_NULL), Int(ImGuiMouseCursor_COUNT))
    KeyOwnerWindows::Vector{Ptr{GLFWwindow}} = fill(Ptr{GLFWwindow}(C_NULL), 512)
    InstalledCallbacks::Bool = true
    WantUpdateMonitors::Bool = true
    PrevUserCallbackMousebutton::Ptr{Cvoid} = C_NULL
    PrevUserCallbackScroll::Ptr{Cvoid} = C_NULL
    PrevUserCallbackKey::Ptr{Cvoid} = C_NULL
    PrevUserCallbackChar::Ptr{Cvoid} = C_NULL
    PrevUserCallbackMonitor::Ptr{Cvoid} = C_NULL
end

Base.show(io::IO, x::Context) = print(io, "Context(id=$(x.id))")

const __GLFW_CONTEXTS = Dict{Int,Context}()
const __GLFW_CONTEXT_COUNTER = Threads.Atomic{Int}(0)

reset_counter() = Threads.atomic_sub!(__GLFW_CONTEXT_COUNTER,  __GLFW_CONTEXT_COUNTER[])
add_counter() = Threads.atomic_add!(__GLFW_CONTEXT_COUNTER, 1)
get_counter() = __GLFW_CONTEXT_COUNTER[]
new_id() = (add_counter(); get_counter();)

"""
    create_context(window=C_NULL; install_callbacks=true, client_api=GlfwClientApi_OpenGL)
Return a GLFW backend contextual data object.
"""
function create_context(window=C_NULL; install_callbacks::Bool=true, client_api::GlfwClientApi=GlfwClientApi_OpenGL)
    client_api == GlfwClientApi_Unknown && throw(ArgumentError("Expected input argument `client_api`: `GlfwClientApi_OpenGL` or `GlfwClientApi_Vulkan`, got `GlfwClientApi_Unknown`."))
    client_api == GlfwClientApi_Vulkan && throw(ArgumentError("Backend support for Vulkan has not been implemented yet, please use `GlfwClientApi_OpenGL` instead."))

    if window == C_NULL
        ctx = create_default_context()
    else
        ctx = Context(; Window=window)
    end

    ctx.InstalledCallbacks = install_callbacks
    ctx.ClientApi = client_api

    # store the ctx in a global variable to prevent it from being GC-ed
    __GLFW_CONTEXTS[ctx.id] = ctx

    # set window userdata
    glfwSetWindowUserPointer(ctx.Window, pointer_from_objref(ctx))

    return ctx
end

"""
    release_context(ctx::Context)
Relese the ctx so it can be GC-ed.
"""
release_context(ctx::Context) = delete!(__GLFW_CONTEXTS, ctx.id)

function create_default_context()
    glfwDefaultWindowHints()
    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3)
    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2)
    @static if Sys.isapple()
        glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE) # 3.2+ only
        glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, true) # required on Mac
    end
    glfwSetErrorCallback(@cfunction(ImGui_ImplGlfw_ErrorCallback, Cvoid, (Cint, Ptr{Cchar})))
    w, h = GLFW_WINDOW_DEFAULT_SIZE[]
    window = glfwCreateWindow(w, h, GLFW_WINDOW_DEFAULT_TITLE[], C_NULL, C_NULL)
    @assert window != C_NULL "failed to create GLFW window."
    glfwMakeContextCurrent(window)
    glfwSwapInterval(1)  # enable vsync
    return Context(; Window=window)
end

set_default_glfw_window_size(width::Integer, height::Integer) = (GLFW_WINDOW_DEFAULT_SIZE[] = (width, height);)
set_default_glfw_window_title(title::AbstractString) = (GLFW_WINDOW_DEFAULT_TITLE[] = title;)

get_window(ctx::Context) = ctx.Window


text:
using GaussQuadrature

T = Float64

print("\nFloating point data type is ", T, '\n')
@printf("\teps = %0.2e\n\n", eps(T))

n = 10
d = rand(T, n)
e = rand(T, n)
z = zero(d)

A = SymTridiagonal(d, e[1:end-1])
D, V = eig(A)

steig!(d, e, z, 30)
idx = sortperm(d)
d = d[idx]
z = z[idx]

@printf("\nEigenvalues:\n\n")
@printf("%20s  %20s  %12s\n\n", "eig", "steig!", "difference")
for i=1:length(d)
    @printf("%20.15f  %20.15f  %12.2e\n", d[i], D[i], d[i]-D[i])
end

@printf("\nFirst component of eigenvectors:\n\n")
@printf("%20s  %20s  %12s\n\n", "eig", "steig!", "difference")
for i=1:length(d)
    @printf("%20.15f  %20.15f  %12.2e\n", z[i], V[1,i],
            abs(z[i])-abs(V[1,i]))
end


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Sequentially moves each vertex to its best point on the tour.
Repeats until no more moves can be found
"""
function moveopt!(tour::Array{Int64, 1}, dist::Array{Int64, 2}, sets::Array{Any, 1},
				  member::Array{Int64,1}, setdist::Distsv)
    improvement_found = true
    number_of_moves = 0
    start_position = 1

    @inbounds while improvement_found && number_of_moves < 10
        improvement_found = false
        for i = start_position:length(tour)
            select_vertex = tour[i]
            delete_cost = removal_cost(tour, dist, i)
			set_ind = member[select_vertex]
            splice!(tour, i)  # remove vertex from tour

            # find the best place to insert the vertex
            v, pos, cost = insert_cost_lb(tour, dist, sets[set_ind], set_ind, setdist,
										  select_vertex, i, delete_cost)
            insert!(tour, pos, v)
            # check if we found a better position for vertex i
            if cost < delete_cost
                improvement_found = true
                number_of_moves += 1
                start_position = min(pos, i) # start looking for swaps where tour change began
                break
            end
        end
    end
end

function moveopt_rand!(tour::Array{Int64, 1}, dist::Array{Int64, 2}, sets::Array{Any, 1},
				  member::Array{Int64,1}, iters::Int64, setdist::Distsv)
	tour_inds = collect(1:length(tour))
	@inbounds for i = 1:iters # i = rand(1:length(tour), iters)
		i = incremental_shuffle!(tour_inds, i)
		select_vertex = tour[i]

		# first check if this vertex should be moved
		delete_cost = removal_cost(tour, dist, i)
	    set_ind = member[select_vertex]
		splice!(tour, i)  # remove vertex from tour
	    v, pos, cost = insert_cost_lb(tour, dist, sets[set_ind], set_ind, setdist,
								      select_vertex, i, delete_cost)
		insert!(tour, pos, v)
    end
end

"""
compute the cost of inserting vertex v into position i of tour
"""
@inline function insert_cost_lb(tour::Array{Int64,1}, dist::Array{Int64,2}, set::Array{Int64, 1}, setind::Int64,
						   setdist::Distsv, bestv::Int64, bestpos::Int, best_cost::Int64)
    @inbounds for i = 1:length(tour)
		v1 = prev_tour(tour, i) # first check lower bound
		lb = setdist.vert_set[v1, setind] + setdist.set_vert[setind, tour[i]] - dist[v1, tour[i]]
		lb > best_cost && continue

		for v in set
	        insert_cost = dist[v1, v] + dist[v, tour[i]] - dist[v1, tour[i]]
	        if insert_cost < best_cost
				best_cost = insert_cost
				bestv = v
				bestpos = i
	        end
		end
    end
    return bestv, bestpos, best_cost
end

"""
determine the cost of removing the vertex at position i in the tour
"""
@inline function removal_cost(tour::Array{Int64, 1}, dist::Array{Int64, 2}, i::Int32)
    if i == 1
        return dist[tour[end], tour[i]] + dist[tour[i], tour[i+1]] - dist[tour[end], tour[i+1]]
    elseif i == length(tour)
        return dist[tour[i-1], tour[i]] + dist[tour[i], tour[1]] - dist[tour[i-1], tour[1]]
    else
        return dist[tour[i-1], tour[i]] + dist[tour[i], tour[i+1]] - dist[tour[i-1], tour[i+1]]
	end
end

""" repeatedly perform moveopt and reopt_tour until there is no improvement """
function opt_cycle!(current::Tour, dist::Array{Int64,2}, sets::Array{Any, 1},
					member::Array{Int64,1}, param::Dict{Symbol, Any}, setdist::Distsv, use)
	current.cost = tour_cost(current.tour, dist)
	prev_cost = current.cost
	for i=1:5
		if i % 2 == 1
			current.tour = reopt_tour(current.tour, dist, sets, member, param)
		elseif param[:mode] == "fast" || use == "partial"
			moveopt_rand!(current.tour, dist, sets, member, param[:max_removals], setdist)
		else
			moveopt!(current.tour, dist, sets, member, setdist)
		end
		current.cost = tour_cost(current.tour, dist)
		if i > 1 && (current.cost >= prev_cost || use == "partial")
			return
		end
		prev_cost = current.cost
	end
end

"""
Given an ordering of the sets, this alg performs BFS to find the
optimal vertex in each set
"""
function reopt_tour(tour::Array{Int64,1}, dist::Array{Int64,2}, sets::Array{Any, 1},
					member::Array{Int64,1}, param::Dict{Symbol, Any})
    best_tour_cost = tour_cost(tour, dist)
	new_tour = copy(tour)
	min_index = min_setv(tour, sets, member, param)
    tour = [tour[min_index:end]; tour[1:min_index-1]]

    prev = zeros(Int64, param[:num_vertices])   # initialize cost_to_come
    cost_to_come = zeros(Int64, param[:num_vertices])
    @inbounds for start_vertex in sets[member[tour[1]]]
 		relax_in!(cost_to_come, dist, prev, Int64[start_vertex], sets[member[tour[2]]])
        for i = 3:length(tour)  # cost to get to ith set on path through (i-1)th set
            relax_in!(cost_to_come, dist, prev, sets[member[tour[i-1]]], sets[member[tour[i]]])
        end
        # find the cost back to the start vertex.
        tour_cost, start_prev = relax(cost_to_come, dist, sets[member[tour[end]]], start_vertex)
        if tour_cost < best_tour_cost   # reconstruct the path
			best_tour_cost = tour_cost
            new_tour = extract_tour(prev, start_vertex, start_prev)
        end
    end
	return new_tour
end

""" Find the set with the smallest number of vertices """
function min_setv(tour::Array{Int64, 1}, sets::Array{Any, 1}, member::Array{Int64, 1},
				  param::Dict{Symbol, Any})
	min_set = param[:min_set]
	@inbounds for i = 1:length(tour)
		member[tour[i]] == min_set && return i
	end
	return 1
end

"""
extracting a tour from the prev pointers.
"""
function extract_tour(prev::Array{Int64,1}, start_vertex::Int64, start_prev::Int64)
    tour = [start_vertex]
    vertex_step = start_prev
    while prev[vertex_step] != 0
        push!(tour, vertex_step)
        vertex_step = prev[vertex_step]
    end
    return reverse(tour)
end

"""
outputs the new cost and prev for vertex v2 after relaxing
does not actually update the cost
"""
@inline function relax(cost::Array{Int64, 1}, dist::Array{Int64, 2}, set1::Array{Int64, 1}, v2::Int64)
	v1 = set1[1]
	min_cost = cost[v1] + dist[v1, v2]
	min_prev = v1
    @inbounds for i = 2:length(set1)
		v1 = set1[i]
		newcost = cost[v1] + dist[v1, v2]
        if min_cost > newcost
            min_cost, min_prev = newcost, v1
        end
    end
    return min_cost, min_prev
end

"""
relaxes the cost of each vertex in the set set2 in-place.
"""
@inline function relax_in!(cost::Array{Int64, 1}, dist::Array{Int64, 2}, prev::Array{Int64, 1},
				 set1::Array{Int64, 1}, set2::Array{Int64, 1})
	@inbounds for v2 in set2
		v1 = set1[1]
		cost[v2] = cost[v1] + dist[v1, v2]
		prev[v2] = v1
        for i = 2:length(set1)
			v1 = set1[i]
			newcost = cost[v1] + dist[v1, v2]
            if cost[v2] > newcost
				cost[v2], prev[v2] = newcost, v1
            end
        end
    end
end


text:
function theAnswer()
    return 42
end

max_stars_repo_name: takbal/linter-julia

max_issues_repo_name: takbal/linter-julia

max_forks_repo_name: takbal/linter-julia


text:
"""otherIncome
other unmentioned means of Temporary Income

"""
mutable struct otherIncome:>
    nonfixedIncome

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
## Helper functions

"""
    remove!(a, item)

Removes `item` from the collection `a`.
"""
function remove!(a, item)
    deleteat!(a, findfirst(==(item), a))
end

"""
    prepare_sums(results, depth)

Returns the DataFrame for the calculation of the majority of votes, as well as the sum of
evaluations up to a given rank.
"""
function prepare_sums(results, depth)
    calculation = DataFrame(Number = results[!, :Number])
    sum_of_eval = DataFrame(Number = results[!, :Number])
    for i = 1:depth
        insertcols!(calculation, Symbol("1-$(i)") => vec(count(<=(i), results[!, Not(:Number)] |> Array, dims = 2)))
        if i == 1
            insertcols!(sum_of_eval, Symbol("1-$(i)") => calculation[:, i+1])
        else
            insertcols!(sum_of_eval, Symbol("1-$(i)") => i .* (calculation[:, i+1] - calculation[:, i]) + sum_of_eval[:, i])
        end
    end
    return calculation, sum_of_eval
end


text:
include("../benchmark.jl")

using JuMP #, Ipopt

m = Model(solver=OnePhaseSolver(init!mu_scale=0.01))
@variable(m, x >= 0)
@variable(m, y >= 0)
@objective(m, Min, x-y)
@NLconstraint(m, x^2 + y^2 <= 1)

solve(m)

getobjectivevalue(m)


text:
"""
	dynamics(y, env, x, u, w)

	evaluates an environment's dynamics (in-place)

	y: state at next time step
	env: Environment
	x: state at current time step
	u: input control (impulse in N.s)
	w: system parameters
"""
function dynamics(y, env::Environment, x, u, w;
	gradients=false,
	attitude_decompress=false)

	timestep = env.mechanism.timestep
	step(env, x, u/timestep,
		gradients=gradients,
		attitude_decompress=attitude_decompress)[1]
    y .= env.state
end

"""
	dynamics_jacobian_state(dx, env, x, u, w)

	evaluates an environment's dynamics Jacobian wrt state (in-place)

	dx: state Jacobian at next time step
	env: Environment
	x: state at current time step
	u: input
	w: system parameters
"""
function dynamics_jacobian_state(dx, env::Environment, x, u, w;
	attitude_decompress=false)

	timestep = env.mechanism.timestep
	step(env, x, u/timestep,
		gradients=true,
		attitude_decompress=attitude_decompress)
    dx .= env.dynamics_jacobian_state
end

"""
	dynamics_jacobian_input(du, env, x, u, w)

	evaluates an environment's dynamics Jacobian wrt input (in-place)

	du: input Jacobian at next time step
	env: Environment
	x: state at current time step
	u: input
	w: system parameters
"""
function dynamics_jacobian_input(du, env::Environment, x, u, w;
	attitude_decompress=false)
	# step(env, x, u, diff=true) # this is run in dynamics_jacobian_state
	du .= env.dynamics_jacobian_input
end


text:
## Could tidy up this HTML to make it look nicer
html_templates = Dict()

grading_partial = """
  if(correct) {
    msgBox.innerHTML = "
👍  {{#:CORRECT}}{{{:CORRECT}}}{{/:CORRECT}}{{^:CORRECT}}Correct{{/:CORRECT}}
";
  } else {
    msgBox.innerHTML = "
👎  {{#:INCORRECT}}{{{:INCORRECT}}}{{/:INCORRECT}}{{^:INCORRECT}}Incorrect{{/:INCORRECT}}
";
  }
"""

## Basic question
## has label and hint option.
## Hint is put with label when present; otherwise, it appears at bottom of form.
## this is overridden with input widget in how show method is called
html_templates["question_tpl"] = mt"""
<form class="mx-2 my-3" name='WeaveQuestion' data-id='{{:ID}}' data-controltype='{{:TYPE}}'>

    {{#:LABEL}}
        {{{:LABEL}}}{{#:HINT}} 🎁{{/:HINT}}

    {{/:LABEL}}

    {{{:FORM}}}
    {{^:LABEL}}{{#:HINT}} 🎁{{/:HINT}}{{/:LABEL}}

</form>

"""

html_templates["input_grading_script"] = jmt"""
document.getElementById("{{:ID}}").addEventListener("change", function() {
  var correct = {{{:CORRECT_ANSWER}}};
  var msgBox = document.getElementById('{{:ID}}_message');
  $(grading_partial)
});
"""

##
html_templates["inputq_form"] = mt"""

    <input id="{{:ID}}" type="{{:TYPE}}" class="form-control" placeholder="{{:PLACEHOLDER}}">

  {{#:UNITS}}{{{:UNITS}}}{{/:UNITS}}{{#:HINT}} 🎁{{/:HINT}}

"""

## Multiple choice (one of many)
## XXX add {{INLINE}}
html_templates["Radioq"] = mt"""
{{#:ITEMS}}

    <input class="form-check-input" type="radio" name="radio_{{:ID}}"
              id="radio_{{:ID}}_{{:VALUE}}" value="{{:VALUE}}">

        {{{:LABEL}}}

    </input>

{{/:ITEMS}}
"""

html_templates["radio_grading_script"] = """
document.querySelectorAll('input[name="radio_{{:ID}}"]').forEach(function(rb) {
rb.addEventListener("change", function() {
    var correct = rb.value == {{:CORRECT_ANSWER}};
    var msgBox = document.getElementById('{{:ID}}_message');
    $(grading_partial)
})});
"""

html_templates["Multiq"] = mt"""
{{#:ITEMS}}

    <input class="form-check-input" type="checkbox" name="check_{{:ID}}"
              id="check_{{:ID}}_{{:VALUE}}" value="{{:VALUE}}">

        {{{:LABEL}}}

    </input>

{{/:ITEMS}}
"""

html_templates["multi_grading_script"] = """
document.querySelectorAll('input[name="check_{{:ID}}"]').forEach(function(rb) {
rb.addEventListener("change", function() {
    var choice_buttons = document.getElementsByName("check_{{:ID}}");
    selected = [];
    for (var i=0; i < choice_buttons.length; i++) {
        if (choice_buttons[i].checked) {
           selected.push(i+1)
        }
    }
    var a = selected;
    b = {{{:CORRECT_ANSWER}}};
    // https://stackoverflow.com/questions/7837456/how-to-compare-arrays-in-javascript
    correct =  (a.length === b.length && a.find((v,i) => v !== b[i]) === undefined)
    var msgBox = document.getElementById('{{:ID}}_message');
    $(grading_partial)
})});
"""