text:
include("utilities.jl")

using SEAL
using Printf
using Random

"""
    bfv_performance_test(context)

Perform multiple performance tests with the BFV scheme for the SEAL encryption context given in
`context`.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref),
          [`example_ckks_performance_default`](@ref), [`example_ckks_performance_custom`](@ref),
          [`bfv_performance_test`](@ref)
"""
function bfv_performance_test(context)
  print_parameters(context)
  println()

  enc_params = parms(first_context_data(context))
  plain_modulus_ = plain_modulus(enc_params)
  poly_modulus_degree_ = poly_modulus_degree(enc_params)

  print("Generating secret/public keys: ")
  keygen = KeyGenerator(context)
  println("Done")

  secret_key_ = secret_key(keygen)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)

  if using_keyswitching(context)
    print("Generating relinearization keys: ")
    time_diff = @elapsedus begin
      relin_keys_ = RelinKeys()
      create_relin_keys!(relin_keys_, keygen)
    end
    println("Done [", time_diff, " microseconds]")

    if !using_batching(qualifiers(key_context_data(context)))
      println("Given encryption parameters do not support batching.")
      return
    end

    print("Generating Galois keys: ")
    time_diff = @elapsedus begin
      galois_keys_ = GaloisKeys()
      create_galois_keys!(galois_keys_, keygen)
    end
    println("Done [", time_diff, " microseconds]")
  end

  encryptor = Encryptor(context, public_key_)
  decryptor = Decryptor(context, secret_key_)
  evaluator = Evaluator(context)
  batch_encoder = BatchEncoder(context)

  time_batch_sum = 0
  time_unbatch_sum = 0
  time_encrypt_sum = 0
  time_decrypt_sum = 0
  time_add_sum = 0
  time_multiply_sum = 0
  time_multiply_plain_sum = 0
  time_square_sum = 0
  time_relinearize_sum = 0
  time_rotate_rows_one_step_sum = 0
  time_rotate_rows_random_sum = 0
  time_rotate_columns_sum = 0
  time_serialize_sum = 0
  time_serialize_zlib_sum = 0
  time_serialize_zstd_sum = 0

  count = 10

  slot_count_ = slot_count(batch_encoder)
  pod_vector = rand(RandomDevice(), UInt64, slot_count_) .% value(plain_modulus_)

  print("Running tests ")
  for i in 1:count
    plain = Plaintext(poly_modulus_degree_, 0)
    plain1 = Plaintext(poly_modulus_degree_, 0)
    plain2 = Plaintext(poly_modulus_degree_, 0)
    time_batch_sum += @elapsedus encode!(plain, pod_vector, batch_encoder)

    pod_vector2 = Vector{UInt64}(undef, slot_count_)
    time_unbatch_sum += @elapsedus decode!(pod_vector2, plain, batch_encoder)

    if pod_vector2 != pod_vector
      error("Batch/unbatch failed. Something is wrong.")
    end

    encrypted = Ciphertext(context)
    time_encrypt_sum += @elapsedus encrypt!(encrypted, plain, encryptor)

    plain2 = Plaintext(poly_modulus_degree(enc_params), 0)
    time_decrypt_sum += @elapsedus decrypt!(plain2, encrypted, decryptor)

    if plain2 != plain
      error("Encrypt/decrypt failed. Something is wrong.")
    end

    encrypted1 = Ciphertext(context)
    encode!(plain1, fill(UInt64(i), slot_count_), batch_encoder)
    encrypt!(encrypted1, plain1, encryptor)
    encrypted2 = Ciphertext(context)
    encode!(plain2, fill(UInt64(i+1), slot_count_), batch_encoder)
    encrypt!(encrypted2, plain2, encryptor)
    time_add_sum += @elapsedus begin
      add_inplace!(encrypted1, encrypted1, evaluator)
      add_inplace!(encrypted2, encrypted2, evaluator)
      add_inplace!(encrypted1, encrypted2, evaluator)
    end

    reserve!(encrypted1, 3)
    time_multiply_sum += @elapsedus multiply_inplace!(encrypted1, encrypted2, evaluator)

    time_multiply_plain_sum += @elapsedus multiply_plain_inplace!(encrypted2, plain, evaluator)

    time_square_sum += @elapsedus square_inplace!(encrypted2, evaluator)

    if using_keyswitching(context)
      time_relinearize_sum += @elapsedus relinearize_inplace!(encrypted1, relin_keys_, evaluator)

      time_rotate_rows_one_step_sum += @elapsedus begin
        rotate_rows_inplace!(encrypted,  1, galois_keys_, evaluator)
        rotate_rows_inplace!(encrypted, -1, galois_keys_, evaluator)
      end

      row_size = div(slot_count(batch_encoder), 2)
      random_rotation = rand(RandomDevice(), Int64) & (row_size - 1)
      time_rotate_rows_random_sum += @elapsedus begin
        rotate_rows_inplace!(encrypted, random_rotation, galois_keys_, evaluator)
      end

      time_rotate_columns_sum += @elapsedus begin
        rotate_columns_inplace!(encrypted, galois_keys_, evaluator)
      end
    end

    buf_size = save_size(ComprModeType.none, encrypted)
    buf = Vector{UInt8}(undef, buf_size)
    time_serialize_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.none, encrypted)
    end

    buf_size = save_size(ComprModeType.zlib, encrypted)
    resize!(buf, buf_size)
    time_serialize_zlib_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.zlib, encrypted)
    end

    buf_size = save_size(ComprModeType.zstd, encrypted)
    resize!(buf, buf_size)
    time_serialize_zstd_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.zstd, encrypted)
    end

    print(".")
    flush(stdout)
  end

  println(" Done")
  println()

  avg_batch = div(time_batch_sum, count)
  avg_unbatch = div(time_unbatch_sum, count)
  avg_encrypt = div(time_encrypt_sum, count)
  avg_decrypt = div(time_decrypt_sum, count)
  avg_add = div(time_add_sum, (3 * count))
  avg_multiply = div(time_multiply_sum, count)
  avg_multiply_plain = div(time_multiply_plain_sum, count)
  avg_square = div(time_square_sum, count)
  avg_relinearize = div(time_relinearize_sum, count)
  avg_rotate_rows_one_step = div(time_rotate_rows_one_step_sum, (2 * count))
  avg_rotate_rows_random = div(time_rotate_rows_random_sum, count)
  avg_rotate_columns = div(time_rotate_columns_sum, count)
  avg_serialize = div(time_serialize_sum, count)
  avg_serialize_zlib = div(time_serialize_zlib_sum, count)
  avg_serialize_zstd = div(time_serialize_zstd_sum, count)

  println("Average batch: ", avg_batch, " microseconds")
  println("Average unbatch: ", avg_unbatch, " microseconds")
  println("Average encrypt: ", avg_encrypt, " microseconds")
  println("Average decrypt: ", avg_decrypt, " microseconds")
  println("Average add: ", avg_add, " microseconds")
  println("Average multiply: ", avg_multiply, " microseconds")
  println("Average multiply plain: ", avg_multiply_plain, " microseconds")
  println("Average square: ", avg_square, " microseconds")
  if using_keyswitching(context)
    println("Average relinearize: ", avg_relinearize, " microseconds")
    println("Average rotate rows one step: ", avg_rotate_rows_one_step, " microseconds")
    println("Average rotate rows random: ", avg_rotate_rows_random, " microseconds")
    println("Average rotate columns: ", avg_rotate_columns, " microseconds")
  end
  println("Average serialize ciphertext: ", avg_serialize, " microseconds")
  println("Average compressed (ZLIB) serialize ciphertext: ", avg_serialize_zlib, " microseconds")
  println("Average compressed (Zstandard) serialize ciphertext: ", avg_serialize_zstd, " microseconds")
end

"""
    ckks_performance_test(context)

Perform multiple performance tests with the CKKS scheme for the SEAL encryption context given in
`context`.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref),
          [`example_ckks_performance_default`](@ref), [`example_ckks_performance_custom`](@ref),
          [`bfv_performance_test`](@ref)
"""
function ckks_performance_test(context)
  print_parameters(context)
  println()

  enc_params = parms(first_context_data(context))
  poly_modulus_degree_ = poly_modulus_degree(enc_params)

  print("Generating secret/public keys: ")
  keygen = KeyGenerator(context)
  println("Done")

  secret_key_ = secret_key(keygen)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)

  if using_keyswitching(context)
    print("Generating relinearization keys: ")
    time_diff = @elapsedus begin
      relin_keys_ = RelinKeys()
      create_relin_keys!(relin_keys_, keygen)
    end
    println("Done [", time_diff, " microseconds]")

    if !using_batching(qualifiers(key_context_data(context)))
      println("Given encryption parameters do not support batching.")
      return
    end

    print("Generating Galois keys: ")
    time_diff = @elapsedus begin
      galois_keys_ = GaloisKeys()
      create_galois_keys!(galois_keys_, keygen)
    end
    println("Done [", time_diff, " microseconds]")
  end

  encryptor = Encryptor(context, public_key_)
  decryptor = Decryptor(context, secret_key_)
  evaluator = Evaluator(context)
  ckks_encoder = CKKSEncoder(context)

  time_encode_sum = 0
  time_decode_sum = 0
  time_encrypt_sum = 0
  time_decrypt_sum = 0
  time_add_sum = 0
  time_multiply_sum = 0
  time_multiply_plain_sum = 0
  time_square_sum = 0
  time_relinearize_sum = 0
  time_rescale_sum = 0
  time_rotate_one_step_sum = 0
  time_rotate_random_sum = 0
  time_conjugate_sum = 0
  time_serialize_sum = 0
  time_serialize_zlib_sum = 0
  time_serialize_zstd_sum = 0

  count = 10

  pod_vector = 1.001 * collect(1:slot_count(ckks_encoder))

  print("Running tests ")
  for i in 1:count
    plain = Plaintext(poly_modulus_degree(enc_params) * length(coeff_modulus(enc_params)), 0)

    scale_ = sqrt(value(last(coeff_modulus(enc_params))))
    time_encode_sum += @elapsedus encode!(plain, pod_vector, scale_, ckks_encoder)

    pod_vector2 = Vector{Float64}(undef, slot_count(ckks_encoder))
    time_decode_sum += @elapsedus decode!(pod_vector2, plain, ckks_encoder)

    encrypted = Ciphertext(context)
    time_encrypt_sum += @elapsedus encrypt!(encrypted, plain, encryptor)

    plain2 = Plaintext(poly_modulus_degree_, 0)
    time_decrypt_sum += @elapsedus decrypt!(plain2, encrypted, decryptor)

    encrypted1 = Ciphertext()
    encode!(plain, i + 1, ckks_encoder)
    encrypt!(encrypted1, plain, encryptor)
    encrypted2 = Ciphertext()
    encode!(plain2, i + 1, ckks_encoder)
    encrypt!(encrypted2, plain2, encryptor)
    time_add_sum += @elapsedus begin
      add_inplace!(encrypted1, encrypted1, evaluator)
      add_inplace!(encrypted2, encrypted2, evaluator)
      add_inplace!(encrypted1, encrypted2, evaluator)
    end

    reserve!(encrypted1, 3)
    time_multiply_sum += @elapsedus multiply_inplace!(encrypted1, encrypted2, evaluator)

    time_multiply_plain_sum += @elapsedus multiply_plain_inplace!(encrypted2, plain, evaluator)

    time_square_sum += @elapsedus square_inplace!(encrypted2, evaluator)

    if using_keyswitching(context)
      time_relinearize_sum += @elapsedus relinearize_inplace!(encrypted1, relin_keys_, evaluator)

      time_rescale_sum += @elapsedus rescale_to_next_inplace!(encrypted1, evaluator)

      time_rotate_one_step_sum += @elapsedus begin
        rotate_vector_inplace!(encrypted,  1, galois_keys_, evaluator)
        rotate_vector_inplace!(encrypted, -1, galois_keys_, evaluator)
      end

      random_rotation = rand(RandomDevice(), Int64) & (slot_count(ckks_encoder) - 1)
      time_rotate_random_sum += @elapsedus begin
        rotate_vector_inplace!(encrypted, random_rotation, galois_keys_, evaluator)
      end

      time_conjugate_sum += @elapsedus complex_conjugate_inplace!(encrypted, galois_keys_, evaluator)
    end

    buf_size = save_size(ComprModeType.none, encrypted)
    buf = Vector{UInt8}(undef, buf_size)
    time_serialize_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.none, encrypted)
    end

    buf_size = save_size(ComprModeType.zlib, encrypted)
    resize!(buf, buf_size)
    time_serialize_zlib_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.zlib, encrypted)
    end

    buf_size = save_size(ComprModeType.zstd, encrypted)
    resize!(buf, buf_size)
    time_serialize_zstd_sum += @elapsedus begin
      save!(buf, buf_size, ComprModeType.zstd, encrypted)
    end

    print(".")
    flush(stdout)
  end

  println(" Done")
  println()

  avg_encode = div(time_encode_sum, count)
  avg_decode = div(time_decode_sum, count)
  avg_encrypt = div(time_encrypt_sum, count)
  avg_decrypt = div(time_decrypt_sum, count)
  avg_add = div(time_add_sum, (3 * count))
  avg_multiply = div(time_multiply_sum, count)
  avg_multiply_plain = div(time_multiply_plain_sum, count)
  avg_square = div(time_square_sum, count)
  avg_relinearize = div(time_relinearize_sum, count)
  avg_rescale = div(time_rescale_sum, count)
  avg_rotate_one_step = div(time_rotate_one_step_sum, (2 * count))
  avg_rotate_random = div(time_rotate_random_sum, count)
  avg_conjugate = div(time_conjugate_sum, count)
  avg_serialize = div(time_serialize_sum, count)
  avg_serialize_zlib = div(time_serialize_zlib_sum, count)
  avg_serialize_zstd = div(time_serialize_zstd_sum, count)

  println("Average encode: ", avg_encode, " microseconds")
  println("Average decode: ", avg_decode, " microseconds")
  println("Average encrypt: ", avg_encrypt, " microseconds")
  println("Average decrypt: ", avg_decrypt, " microseconds")
  println("Average add: ", avg_add, " microseconds")
  println("Average multiply: ", avg_multiply, " microseconds")
  println("Average multiply plain: ", avg_multiply_plain, " microseconds")
  println("Average square: ", avg_square, " microseconds")
  if using_keyswitching(context)
    println("Average relinearize: ", avg_relinearize, " microseconds")
    println("Average rescale: ", avg_rescale, " microseconds")
    println("Average rotate vector one step: ", avg_rotate_one_step, " microseconds")
    println("Average rotate vector random: ", avg_rotate_random, " microseconds")
    println("Average complex conjugate: ", avg_conjugate, " microseconds")
  end
  println("Average serialize ciphertext: ", avg_serialize, " microseconds")
  println("Average compressed (ZLIB) serialize ciphertext: ", avg_serialize_zlib, " microseconds")
  println("Average compressed (Zstandard) serialize ciphertext: ", avg_serialize_zstd, " microseconds")
  flush(stdout)
end

"""
    example_bfv_performance_default()

Perform multiple performance tests of the BFV scheme for several polynomial modulus degrees.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref), [`example_bfv_performance_custom`](@ref),
          [`bfv_performance_test`](@ref),
          [`example_ckks_performance_default`](@ref)
"""
function example_bfv_performance_default()
  print_example_banner("BFV Performance Test with Degrees: 4096, 8192, and 16384")

  enc_parms = EncryptionParameters(SchemeType.bfv)
  poly_modulus_degree = 4096
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  set_plain_modulus!(enc_parms, 786433)
  bfv_performance_test(SEALContext(enc_parms))

  println()
  poly_modulus_degree = 8192
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  set_plain_modulus!(enc_parms, 786433)
  bfv_performance_test(SEALContext(enc_parms))

  println()
  poly_modulus_degree = 16384
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  set_plain_modulus!(enc_parms, 786433)
  bfv_performance_test(SEALContext(enc_parms))

  # Comment out the following to run the biggest example
  #
  # println()
  # poly_modulus_degree = 32768
  # set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  # set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  # set_plain_modulus!(enc_parms, 786433)
  # bfv_performance_test(SEALContext(enc_parms))
end

"""
    example_bfv_performance_custom()

Show a terminal menu for selecting a custom polynomial modulus degree and use it for performance
tests with the BFV scheme.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref), [`example_bfv_performance_default`](@ref),
          [`bfv_performance_test`](@ref),
          [`example_ckks_performance_custom`](@ref)
"""
function example_bfv_performance_custom()
  poly_modulus_degree = 0
  println()
  print("Set poly_modulus_degree (1024, 2048, 4096, 8192, 16384, or 32768): ")
  input = readline()
  if !isopen(stdin)
    return
  end
  try
    poly_modulus_degree = parse(Int, input)
  catch
    println("Invalid option")
    return
  end
  if (poly_modulus_degree < 1024 || poly_modulus_degree > 32768 || (poly_modulus_degree &
                                                                    (poly_modulus_degree - 1)) != 0)
    println("Invalid option")
    return
  end

  print_example_banner("BFV Performance Test with Degree: $poly_modulus_degree")

  enc_parms = EncryptionParameters(SchemeType.bfv)
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  if poly_modulus_degree == 1024
    set_plain_modulus!(enc_parms, 12289)
  else
    set_plain_modulus!(enc_parms, 786433)
  end
  bfv_performance_test(SEALContext(enc_parms))
end

"""
    example_ckks_performance_default()

Perform multiple performance tests of the CKKS scheme for several polynomial modulus degrees.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref), [`example_ckks_performance_custom`](@ref),
          [`ckks_performance_test`](@ref),
          [`example_bfv_performance_default`](@ref)
"""
function example_ckks_performance_default()
  print_example_banner("CKKS Performance Test with Degrees: 4096, 8192, and 16384")

  # BFV primes are not recommended for CKKS, but are good enough for performance testing
  enc_parms = EncryptionParameters(SchemeType.ckks)
  poly_modulus_degree = 4096
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  ckks_performance_test(SEALContext(enc_parms))

  println()
  poly_modulus_degree = 8192
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  ckks_performance_test(SEALContext(enc_parms))

  println()
  poly_modulus_degree = 16384
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  ckks_performance_test(SEALContext(enc_parms))

  # Comment out the following to run the biggest example
  #
  # println()
  # poly_modulus_degree = 32768
  # set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  # set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  # ckks_performance_test(SEALContext(enc_parms))
end

"""
    example_ckks_performance_custom()

Show a terminal menu for selecting a custom polynomial modulus degree and use it for performance
tests with the CKKS scheme.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_performance_test`](@ref), [`example_ckks_performance_default`](@ref),
          [`ckks_performance_test`](@ref),
          [`example_bfv_performance_custom`](@ref)
"""
function example_ckks_performance_custom()
  poly_modulus_degree = 0
  println()
  print("Set poly_modulus_degree (1024, 2048, 4096, 8192, 16384, or 32768): ")
  input = readline()
  if !isopen(stdin)
    return
  end
  try
    poly_modulus_degree = parse(Int, input)
  catch
    println("Invalid option")
    return
  end
  if (poly_modulus_degree < 1024 || poly_modulus_degree > 32768 || (poly_modulus_degree &
                                                                    (poly_modulus_degree - 1)) != 0)
    println("Invalid option")
    return
  end

  print_example_banner("CKKS Performance Test with Degree: $poly_modulus_degree")

  enc_parms = EncryptionParameters(SchemeType.ckks)
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)
  set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(poly_modulus_degree))
  ckks_performance_test(SEALContext(enc_parms))
end

"""
    example_performance_test()

Show a terminal menu for selecting different performance tests.

This function is based on the file `native/examples/7_performance.cpp` of the original
SEAL library and should yield the exact same output, except (of course) for the actual performance
numbers.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/7_performance.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/7_performance.cpp)

See also: [`example_bfv_performance_default`](@ref), [`example_bfv_performance_custom`](@ref),
          [`bfv_performance_test`](@ref),
          [`example_ckks_performance_default`](@ref), [`example_ckks_performance_custom`](@ref),
          [`ckks_performance_test`](@ref)
"""
function example_performance_test()
  print_example_banner("Example: Performance Test")

  while true
    println()
    println("Select a scheme (and optionally poly_modulus_degree):")
    println("  1. BFV with default degrees")
    println("  2. BFV with a custom degree")
    println("  3. CKKS with default degrees")
    println("  4. CKKS with a custom degree")
    println("  0. Back to main menu")

    selection = 0
    println()
    print("> Run performance test (1 ~ 4) or go back (0): ")
    input = readline()
    if !isopen(stdin)
      println()
      return
    end
    try
      selection = parse(Int, input)
    catch
      println("Invalid option")
      continue
    end

    if selection == 1
      example_bfv_performance_default()
    elseif selection == 2
      example_bfv_performance_custom()
    elseif selection == 3
      example_ckks_performance_default()
    elseif selection == 4
      example_ckks_performance_custom()
    elseif selection == 0
      println()
      return
    else
      println("Invalid option")
    end
  end
end


text:
# Copyright (c) 2021, Rolfe Power

using Test
using Pat

#  These exist mainly to warn if there's been an unintentional change
@test Pat.GRAVITATIONAL_CONSTANT == 6.6743e-20
@test Pat.SPEED_OF_LIGHT == 299792458e-3
@test Pat.EARTH_STANDARD_GRAVITY == 9.80665e-03
@test Pat.BOLTZMANN_CONSTANT == 1.380649e-23
@test Pat.ASTRONOMICAL_UNIT == 149597870700e-3
@test Pat.ELECTRON_MASS == 9.109_383_7015e-31
@test Pat.NEUTRON_MASS == 1.674_927_498_04e-27
@test Pat.PROTON_MASS == 1.672_621_923_69e-27

ext: jl


text:
module TestReadWriteLock

using ConcurrentUtils
using Test

using ..Utils: poll_until, unfair_sleep

function test_no_blocks()
    lck = ReadWriteLock()

    @sync begin
        lock_read(lck)
        lock_read(lck)
        Threads.@spawn begin
            lock_read(lck)
            unlock_read(lck)
        end
        unlock_read(lck)
        unlock_read(lck)
    end

    lock(lck)
    unlock(lck)
end

function check_minimal_lock_interface(lck)
    phase = Threads.Atomic{Int}(0)
    lock(lck)
    @sync begin
        Threads.@spawn begin
            phase[] = 1
            lock(lck)
            unlock(lck)
            phase[] = 2
        end

        @test poll_until(() -> phase[] != 0)
        @test phase[] == 1
        sleep(0.01)
        @test phase[] == 1

        unlock(lck)
    end
    @test phase[] == 2
end

test_wlock() = check_minimal_lock_interface(ReadWriteLock())

function test_a_writer_blocks_a_reader()
    lck = ReadWriteLock()
    locked = Threads.Atomic{Bool}(false)
    @sync lock(lck) do
        Threads.@spawn begin
            lock_read(lck)
            locked[] = true
            unlock_read(lck)
        end

        sleep(0.01)
        @test !locked[]
    end
    @test locked[]
end

function test_a_writer_blocks_a_writer()
    lck = ReadWriteLock()

    locked = Threads.Atomic{Bool}(false)
    @sync lock(lck) do
        Threads.@spawn begin
            lock(lck)
            locked[] = true
            unlock(lck)
        end

        sleep(0.01)
        @test !locked[]
    end
    @test locked[]
end

function test_a_reader_blocks_a_writer()
    lck = ReadWriteLock()

    locked = Threads.Atomic{Bool}(false)
    @sync lock_read(lck) do
        Threads.@spawn begin
            lock(lck)
            locked[] = true
            unlock(lck)
        end

        sleep(0.01)
        @test !locked[]
    end
    @test locked[]
end

function check_concurrent_mutex(nreaders, nwriters, ntries)
    lck = ReadWriteLock()

    limit = nwriters * ntries
    # nreads = Threads.Atomic{Int}(0)
    ref = Ref(0)
    @sync begin
        for _ in 1:nreaders
            Threads.@spawn while true
                lock_read(lck) do
                    # Threads.atomic_add!(nreads, 1)
                    ref[] < limit
                end || break
                yield()
            end
        end
        sleep(0.01)

        for _ in 1:nwriters
            Threads.@spawn for _ in 1:ntries
                lock(lck) do
                    local x = ref[]

                    # sleep about 3 μs
                    unfair_sleep(100)

                    ref[] = x + 1
                end
            end
        end
    end
    # @show nreads[]

    return ref[]
end

function test_concurrent_mutex()
    @testset for nwriters in [Threads.nthreads(), 64 * Threads.nthreads()]

        nwriters = 2
        ntries = 1000
        actual = check_concurrent_mutex(nwriters, nwriters, ntries)
        @test actual == nwriters * ntries
    end
end

end  # module


text:
# ---
# title: 848. Shifting Letters
# id: problem848
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/shifting-letters/description/>
# hidden: true
# ---
#
# We have a string `S` of lowercase letters, and an integer array `shifts`.
#
# Call the _shift_ of a letter, the next letter in the alphabet, (wrapping
# around so that `'z'` becomes `'a'`).
#
# For example, `shift('a') = 'b'`, `shift('t') = 'u'`, and `shift('z') = 'a'`.
#
# Now for each `shifts[i] = x`, we want to shift the first `i+1` letters of `S`,
# `x` times.
#
# Return the final string after all such shifts to `S` are applied.
#
# **Example 1:**
#
#
#
#     Input: S = "abc", shifts = [3,5,9]
#     Output: "rpl"
#     Explanation:
#     We start with "abc".
#     After shifting the first 1 letters of S by 3, we have "dbc".
#     After shifting the first 2 letters of S by 5, we have "igc".
#     After shifting the first 3 letters of S by 9, we have "rpl", the answer.
#
#
# **Note:**
#
#   1. `1 <= S.length = shifts.length <= 20000`
#   2. `0 <= shifts[i] <= 10 ^ 9`
#
#
## @lc code=start
using LeetCode

function shifting_letters(s::String, shift::Vector{Int})::String
    for i in length(shift):-1:2
        shift[i - 1] += shift[i]
    end
    res = codeunits(s)[:]
    for i in 1:length(shift)
        res[i] += shift[i]
    end
    return String(res)
end
## @lc code=end


text:
# Autogenerated wrapper script for MbedTLS_jll for x86_64-w64-mingw32
export libmbedcrypto, libmbedtls, libmbedx509

JLLWrappers.@generate_wrapper_header("MbedTLS")
JLLWrappers.@declare_library_product(libmbedcrypto, "libmbedcrypto.dll")
JLLWrappers.@declare_library_product(libmbedtls, "libmbedtls.dll")
JLLWrappers.@declare_library_product(libmbedx509, "libmbedx509.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmbedcrypto,
        "bin\\libmbedcrypto.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmbedtls,
        "bin\\libmbedtls.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmbedx509,
        "bin\\libmbedx509.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, Unitful

makedocs()


text:
@polly function kernel_covariance(float_n, data, cov, mean)
    n,m = size(data)

    for j = 1:m
        mean[j] = zero(eltype(mean))
        for i = 1:n
            mean[j] += data[i,j]
        end
        mean[j] /= float_n
    end

    for i = 1:n, j = 1:m
        data[i,j] -= mean[j]
    end

    for i = 1:m, j = i:m
        cov[i,j] = zero(eltype(cov))
        for k = 1:n
            cov[i,j] += data[k,i] * data[k,j]
        end
        cov[i,j] /= (float_n - one(eltype(cov)))
        cov[j,i] = cov[i,j]
    end
end

let
    m = 1200
    n = 1400

    data = zeros(Float32,n,m)
    cov = zeros(Float32,m,m)
    mean = zeros(Float32,m)

    float_n = Float32(n)

    for i = 1:n, j = 1:m
        data[i,j] = (i*j) / m
    end

    SUITE["covariance"] = @benchmarkable kernel_covariance($float_n, data, cov, mean) setup = (data = copy($data); cov = copy($cov); mean = copy($mean))
end


text:
module SpikingNeuralNetworks

const SNN = SpikingNeuralNetworks

using SparseArrays
try
    using Reexport
    using Parameters
    using Requires

catch
    import Pkg; Pkg.add("Reexport")
    import Pkg; Pkg.add("Parameters")
    import Pkg; Pkg.add("Requires")
    using Reexport
    using Parameters
    using Requires
end
#using Unitful
#using Unitful.DefaultSymbols
#@reexport using Utils

const SNNInt = Int32
const SNNFloat = Float32
# srand(1000)

include("units.jl")
include("main.jl")
include("utils.jl")

include("neuron/if.jl")
include("neuron/if2.jl")
include("neuron/noisy_if.jl")
include("neuron/poisson.jl")
include("neuron/iz.jl")
include("neuron/iz_more_param.jl")
include("neuron/adexp.jl")

include("neuron/hh.jl")
include("neuron/rate.jl")

include("synapse/spiking_synapse.jl")
include("synapse/rate_synapse.jl")
include("synapse/pinning_full_synapse.jl")
# include("synapse/pinning_synapse.jl")
include("synapse/fl_full_synapse.jl")
# include("synapse/fl_synapse.jl")

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
end

end

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
using LockandKeyLookups
using Test, DataFrames

#DEVELOPERS: Useful snippit for debugging
# for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a, X -> X.b )
#      println( i )
#      println(b[i[2][2],:b])
# end

@testset "LockandKeyLookups - iterations" begin
    # examples...
    endrng = 3
    a = DataFrame( Dict( :a => 1:1:endrng, :b => 1:1:endrng ) )
    b = DataFrame( Dict( :b => 1:1:endrng, :c => 1:1:endrng ) )
    #Test single stream complete match
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a, X -> X.b ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ) , 2 => ( 1, 2 ), 3 => ( 1, 3 ) ] )
    #Test two streams complete match
    c = DataFrame( Dict( :b => 1:1:endrng, :c => 1:1:endrng ) )
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b, c ] ), X -> X.a, X -> X.b ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ), 1 => ( 2, 1 ), 2 => ( 1, 2 ), 2 => ( 2, 2 ), 3 => ( 1, 3 ), 3 => ( 2, 3 ) ] )
    #Test 3 streams with heterogeniety of matches
    endrng = 10
    a = DataFrame( Dict( :a => 1:1:endrng, :b => 1:1:endrng ) )
    b = DataFrame( Dict( :b => 1:5:endrng, :c => 1:5:endrng ) )
    c = DataFrame( Dict( :b => 6:7:endrng, :c => 6:7:endrng ) )
    d = DataFrame( Dict( :b => 6:7:endrng, :c => 6:7:endrng ) )
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b, c, d ] ), X -> X.a, X -> X.b ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ) , 6 => ( 1, 2 ), 6 => ( 2, 1 ), 6 => ( 3, 1 ) ] )

    endrng = 5
    a = DataFrame( Dict( :a => 1:2:endrng, :b => 1:2:endrng ) )
    b = DataFrame( Dict( :b => 1:1:endrng, :c => 1:1:endrng ) )
    #Test sparse key stream!
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a, X -> X.b ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ) , 2 => ( 1, 3 ), 3 => ( 1, 5 ) ] )

    c = DataFrame( Dict( :b => 1:2:endrng, :c => 1:2:endrng ) )
    #Test sparse key stream with 2 pins!
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b, c ] ), X -> X.a, X -> X.b ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ), 1 => ( 2, 1 ), 2 => ( 1, 3 ), 2 => ( 2, 2 ), 3 => ( 1, 5 ), 3 => ( 2, 3 ) ] )

    b = DataFrame( Dict( :b => 2:2:endrng, :c => 2:2:endrng ) )
    #Test sparse key stream with NO matches!
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a, X -> X.b ) ]
    @test isempty( LKL_result )

    c = DataFrame( Dict( :b => 2:2:20, :c => 2:2:20 ) )
    #Test sparse key stream with 2 streams that have NO matches!
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b, c ] ), X -> X.a, X -> X.b ) ]
    @test isempty( LKL_result )

    a = DataFrame( Dict( :a => 101:2:111, :b => 101:2:111 ) )
    b = DataFrame( Dict( :b => 20:1:25, :c => 20:1:25 ) )
    c = DataFrame( Dict( :b => 107:1:115, :c => 107:1:115 ) )
    #Sparse stream some matching some not
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b,c ] ), X -> X.a, X -> X.c ) ]
    @test all( LKL_result .== [ 4 => ( 2, 1 ) , 5 => ( 2, 3 ), 6 => ( 2, 5 ) ] )
end

@testset "LockandKeyLookups - Arbitrary functions" begin
    a = DataFrame( Dict( :a => 7:1:11, :b => 7:1:11 ) )
    b = DataFrame( Dict( :b => 1:1:22, :c => 1:1:22 ) )

    #test changing a value
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a,
                        X -> X.b/2.0 ) ]
    @test all( LKL_result .== [ 1 => ( 1, 14 ) , 2 => ( 1, 16 ), 3 => ( 1, 18 ), 4 => ( 1, 20 ), 5 => ( 1, 22 ) ] )

    #test changing a comparator but not really
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a,
                        X -> X.b / 2.0, (k, t) -> (1.0 .* k) == (t ./ 1.0) ) ]
    println(LKL_result)
    @test all( LKL_result .== [ 1 => ( 1, 14 ) , 2 => ( 1, 16 ), 3 => ( 1, 18 ), 4 => ( 1, 20 ), 5 => ( 1, 22 ) ] )

end

@testset "LockandKeyLookups - Internal code tests..." begin
    a = DataFrame( Dict( :a => 1:1:11, :b => 1:1:11 ) )
    b = DataFrame( Dict( :b => 1:1:11, :c => 1:1:11 ) )
    #test changing a value
    LKL_iterator = LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X.a, X -> X.b )
    @test length( LKL_iterator ) == Base.SizeUnknown()
    @test Base.IteratorSize( LKL_iterator ) == Base.SizeUnknown()
    #
    @test all( LockandKeyLookups.unpack_2tuples( [(1,2), (3,4)] ) .== [[1,3], [2,4] ])
    @test all( LockandKeyLookups.unpack_2tuples( [(1,2), nothing, (3,4)] ) .== [[1,nothing,3], [2,nothing,4] ])
end

@testset "LockandKeyLookups - run simple test with arrays" begin
    # examples...
    a = [ [ 1 2 ];
            [ 3 4 ];
            [ 5 6] ]
    b = [ [ 1 2 ];
            [ 3 4 ];
            [ 5 6] ]
    #Test single stream complete match
    LKL_result = [ i for i in LockandKeyLookup( eachrow( a ), eachrow.( [ b ] ), X -> X[:,1], X -> X[:,1] ) ]
    @test all( LKL_result .== [ 1 => ( 1, 1 ) , 2 => ( 1, 2 ), 3 => ( 1, 3 ) ] )
end


text:
module CP

using Compat: eachrow, isnothing
using LinearAlgebra: det
using Parameters: @with_kw
using Setfield: get, set, @lens, @set
using StaticArrays: SVector, SMatrix, FieldVector
using Unitful
using UnitfulAtomic

using ..Inputs:
    Namelist,
    QuantumESPRESSOInput,
    Card,
    optionof,
    optionpool,
    allnamelists,
    allcards,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards

import AbInitioSoftwareBase.Inputs: asstring, groupname
import Crystallography: Bravais, Lattice
# import Pseudopotentials: pseudoformat
import ..Inputs:
    optionpool,
    allnamelists,
    allcards,
    optionof,
    required_namelists,
    optional_namelists,
    required_cards,
    optional_cards
using ..Formatter: delimiter, newline, indent, floatfmt, intfmt

export ControlNamelist,
    SystemNamelist,
    ElectronsNamelist,
    IonsNamelist,
    CellNamelist,
    PressAiNamelist,
    WannierNamelist,
    AtomicSpecies,
    AtomicSpeciesCard,
    AtomicPosition,
    AtomicPositionsCard,
    CellParametersCard,
    RefCellParametersCard,
    AtomicVelocity,
    AtomicVelocitiesCard,
    AtomicForce,
    AtomicForcesCard,
    CPInput,
    optconvert

include("namelists.jl")
include("cards.jl")

"""
    CPInput <: QuantumESPRESSOInput
    CPInput(control, system, electrons, ions, cell, press_ai, wannier, atomic_species, atomic_positions, atomic_velocities, cell_parameters, ref_cell_parameters, constraints, occupations, atomic_forces, plot_wannier, autopilot)

Construct a `CPInput` which represents the input of program `cp.x`.

# Arguments
- `control::ControlNamelist=ControlNamelist()`: the `CONTROL` namelist of the input. Optional.
- `system::SystemNamelist=SystemNamelist()`: the `SYSTEM` namelist of the input. Optional.
- `electrons::ElectronsNamelist=ElectronsNamelist()`: the `ELECTRONS` namelist of the input. Optional.
- `ions::IonsNamelist=IonsNamelist()`: the `IONS` namelist of the input. Optional.
- `cell::CellNamelist=CellNamelist()`: the `CELL` namelist of the input. Optional.
- `press_ai::PressAiNamelist=PressAiNamelist()`: the `PRESS_AI` namelist of the input. Optional.
- `wannier::WannierNamelist=WannierNamelist()`: the `WANNIER` namelist of the input. Optional.
- `atomic_species::AtomicSpeciesCard`: the `ATOMIC_SPECIES` card of the input. Must be provided explicitly.
- `atomic_positions::AtomicPositionsCard`: the `ATOMIC_POSITIONS` card of the input. Must be provided explicitly.
- `atomic_velocities::AtomicVelocitiesCard`: the `ATOMIC_VELOCITIES` card of the input. Must be provided explicitly.
- `cell_parameters::Union{Nothing,CellParametersCard}`: the `CELL_PARAMETERS` card of the input. Must be either `nothing` or a `CellParametersCard`.
- `ref_cell_parameters::Union{Nothing,RefCellParametersCard}`: the `REF_CELL_PARAMETERS` card of the input. Must be either `nothing` or a `CellParametersCard`.
"""
@with_kw struct CPInput <: QuantumESPRESSOInput
    control::ControlNamelist = ControlNamelist()
    system::SystemNamelist = SystemNamelist()
    electrons::ElectronsNamelist = ElectronsNamelist()
    ions::IonsNamelist = IonsNamelist()
    cell::CellNamelist = CellNamelist()
    press_ai::PressAiNamelist = PressAiNamelist()
    wannier::WannierNamelist = WannierNamelist()
    atomic_species::AtomicSpeciesCard
    atomic_positions::AtomicPositionsCard
    atomic_velocities::AtomicVelocitiesCard
    cell_parameters::Union{Nothing,CellParametersCard} = nothing
    ref_cell_parameters::Union{Nothing,RefCellParametersCard} = nothing
    constraints::Union{Nothing,Float64} = nothing
    occupations::Union{Nothing,Float64} = nothing
    atomic_forces::Union{Nothing,AtomicForcesCard} = nothing
    plot_wannier::Union{Nothing,Float64} = nothing
    autopilot::Union{Nothing,Float64} = nothing
    @assert(
        !(isnothing(cell_parameters) && system.ibrav == 0),
        "Cannot specify `ibrav = 0` with an empty `cell_parameters`!"
    )
end # struct CPInput

"""
    optconvert(new_option::AbstractString, card::AbstractCellParametersCard)

Convert the option of an `AbstractCellParametersCard` from "bohr" to "angstrom", or its reverse.

!!! warning
    It does not support conversion between `"alat"` and the others.
"""
function optconvert(new_option::AbstractString, card::AbstractCellParametersCard)
    old_option = optionof(card)
    if new_option == old_option
        return card  # No conversion is needed
    else
        return typeof(card)(
            if (old_option => new_option) == ("bohr" => "angstrom")
                @. ustrip(u"angstrom", card.data * u"bohr")
            elseif (old_option => new_option) == ("angstrom" => "bohr")
                @. ustrip(u"bohr", card.data * u"angstrom")
            else
                error("unknown conversion rule $(old_option => new_option)!")
            end,
        )
    end
end # function optconvert

groupname(::Type{ControlNamelist}) = "CONTROL"
groupname(::Type{SystemNamelist}) = "SYSTEM"
groupname(::Type{ElectronsNamelist}) = "ELECTRONS"
groupname(::Type{IonsNamelist}) = "IONS"
groupname(::Type{CellNamelist}) = "CELL"
groupname(::Type{AtomicSpeciesCard}) = "ATOMIC_SPECIES"
groupname(::Type{AtomicPositionsCard}) = "ATOMIC_POSITIONS"
groupname(::Type{<:CellParametersCard}) = "CELL_PARAMETERS"

"""
    Bravais(nml::SystemNamelist)

Return a `Bravais` from a `SystemNamelist`.
"""
Bravais(nml::SystemNamelist) = Bravais(nml.ibrav)

"""
    Lattice(nml::SystemNamelist)

Return a `Lattice` from a `SystemNamelist`.
"""
Lattice(nml::SystemNamelist) = Lattice(Bravais(nml), nml.celldm)

"""
    cellvolume(nml::SystemNamelist)

Return the volume of the cell based on the information given in a `SystemNamelist`, in atomic unit.
"""
cellvolume(nml::SystemNamelist) = cellvolume(Lattice(nml))
"""
    cellvolume(card)

Return the cell volume of a `CellParametersCard` or `RefCellParametersCard`, in atomic unit.

!!! warning
    It will throw an error if the option is `"alat"`.
"""
function cellvolume(card::AbstractCellParametersCard)
    option = optionof(card)
    if option == "bohr"
        abs(det(card.data))
    elseif option == "angstrom"
        ustrip(u"bohr^3", abs(det(card.data)) * u"angstrom^3")
    else  # option == "alat"
        error("information not enough! Parameter `celldm[1]` needed!")
    end
end # function cellvolume

function asstring(data::AtomicSpecies)
    return join(
        (sprintf1("%3s", data.atom), sprintf1(floatfmt(data), data.mass), data.pseudopot),
        delimiter(data),
    )
end
function asstring(card::AtomicSpeciesCard)
    # Using generator expressions in `join` is faster than using `Vector`s.
    return "ATOMIC_SPECIES" *
           newline(card) *
           join((indent(card) * asstring(x) for x in unique(card.data)), newline(card))
end
function asstring(data::AtomicPosition)
    f(x) = x ? "" : "0"
    return join(
        [
            sprintf1("%3s", data.atom)
            map(x -> sprintf1(floatfmt(data), x), data.pos)
            map(f, data.if_pos)
        ],
        delimiter(data),
    )
end
function asstring(card::AtomicPositionsCard)
    return "ATOMIC_POSITIONS { $(optionof(card)) }" *
           newline(card) *
           join((indent(card) * asstring(x) for x in card.data), newline(card))
end
function asstring(card::CellParametersCard)
    it = (
        indent * join((sprintf1(floatfmt(card), x) for x in row), delimiter(card)) for
        row in eachrow(card.data)
    )
    return "CELL_PARAMETERS { $(optionof(card)) }" * newline(card) * join(it, newline)
end

optionof(::AtomicVelocitiesCard) = "a.u"

optionpool(::Type{<:AtomicPositionsCard}) =
    ("alat", "bohr", "angstrom", "crystal", "crystal_sg")
optionpool(::Type{<:CellParametersCard}) = ("alat", "bohr", "angstrom")
optionpool(::Type{<:AtomicVelocity}) = ("a.u",)
optionpool(::Type{<:RefCellParametersCard}) = ("bohr", "angstrom")

allnamelists(::Type{CPInput}) =
    (:control, :system, :electrons, :ions, :cell, :press_ai, :wannier)
allnamelists(x::CPInput) = (getfield(x, f) for f in allnamelists(typeof(x)))

allcards(::Type{CPInput}) = (
    :atomic_species,
    :atomic_positions,
    :atomic_velocities,
    :cell_parameters,
    :ref_cell_parameters,
    :constraints,
    :occupations,
    :atomic_forces,
    :plot_wannier,
    :autopilot,
)
allcards(x::CPInput) = (getfield(x, f) for f in allcards(typeof(x)))

required_namelists(::Type{CPInput}) = (:control, :system, :electrons)
required_namelists(x::CPInput) = (getfield(x, f) for f in required_namelists(typeof(x)))

optional_namelists(::Type{CPInput}) = (:ions, :cell, :press_ai, :wannier)
optional_namelists(x::CPInput) = (getfield(x, f) for f in optional_namelists(typeof(x)))

required_cards(::Type{CPInput}) = (:atomic_species, :atomic_positions)
required_cards(x::CPInput) = (getfield(x, f) for f in required_cards(typeof(x)))

optional_cards(::Type{CPInput}) = (
    :atomic_velocities,
    :cell_parameters,
    :ref_cell_parameters,
    :constraints,
    :occupations,
    :atomic_forces,
    :plot_wannier,
    :autopilot,
)
optional_cards(x::CPInput) = (getfield(x, f) for f in optional_cards(typeof(x)))

end


text:
module CommunalHelperMelvin

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/Melvin" Melvin(
    x::Integer,
    y::Integer,
    width::Integer=24,
    height::Integer=24,
    weakTop::Bool=false,
    weakRight::Bool=false,
    weakBottom::Bool=false,
    weakLeft::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Melvin (All Strong) (Communal Helper)" => Ahorn.EntityPlacement(
        Melvin,
        "rectangle",
    ),
    "Melvin (All Weak) (Communal Helper)" => Ahorn.EntityPlacement(
        Melvin,
        "rectangle",
        Dict{String,Any}(
            "weakTop" => true,
            "weakRight" => true,
            "weakBottom" => true,
            "weakLeft" => true,
        ),
    ),
    "Melvin (Horizontally Weak) (Communal Helper)" => Ahorn.EntityPlacement(
        Melvin,
        "rectangle",
        Dict{String,Any}(
            "weakRight" => true,
            "weakLeft" => true,
        ),
    ),
    "Melvin (Vertically Weak) (Communal Helper)" => Ahorn.EntityPlacement(
        Melvin,
        "rectangle",
        Dict{String,Any}(
            "weakTop" => true,
            "weakBottom" => true,
        ),
    ),
)

Ahorn.minimumSize(entity::Melvin) = 24, 24
Ahorn.resizable(entity::Melvin) = true, true

Ahorn.selection(entity::Melvin) = Ahorn.getEntityRectangle(entity)

const strongTiles = "objects/CommunalHelper/melvin/block_strong"
const weakTiles = "objects/CommunalHelper/melvin/block_weak"
const weakCornersHTiles = "objects/CommunalHelper/melvin/corners_weak_h"
const weakCornersVTiles = "objects/CommunalHelper/melvin/corners_weak_v"

const insideTiles = "objects/CommunalHelper/melvin/inside"
const eye = "objects/CommunalHelper/melvin/eye/idle_small00"

const kevinColor = (98, 34, 43) ./ 255

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Melvin, room::Maple.Room)
    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))

    eyeSprite = Ahorn.getSprite(eye, "Gameplay")

    tileRands = Ahorn.getDrawableRoom(Ahorn.loadedState.map, room).fgTileStates.rands

    weakTop = Bool(get(entity.data, "weakTop", false))
    weakBottom = Bool(get(entity.data, "weakBottom", false))
    weakRight = Bool(get(entity.data, "weakRight", false))
    weakLeft = Bool(get(entity.data, "weakLeft", false))

    x, y = Ahorn.position(entity)
    rtx = floor(Int, x / 8)
    rty = floor(Int, y / 8)

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, kevinColor)

    # edges and interior
    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 24 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 24 : 8)
        isEdge = (i == 1 || i == tileWidth || j == 1 || j == tileHeight)
        isCorner = ((i == 1 || i == tileWidth) && (j == 1 || j == tileHeight))
        tileChoice = mod(get(tileRands, (rty + j, rtx + i), 0), 4)

        if tx == 8
            if tileChoice == 1 || tileChoice == 3
                tx += 8
            end
        end
        if ty == 8
            if tileChoice == 2 || tileChoice == 3
                ty += 8
            end
        end
        if !isEdge
            tx -= 8
            ty -= 8
        end
        if !isEdge
            Ahorn.drawImage(ctx, insideTiles, (i - 1) * 8, (j - 1) * 8, tx, ty, 8, 8)
        elseif !isCorner
            block = strongTiles
            if i == 1
                block = weakLeft ? weakTiles : strongTiles
            elseif i == tileWidth
                block = weakRight ? weakTiles : strongTiles
            elseif j == 1
                block = weakTop ? weakTiles : strongTiles
            elseif j == tileHeight
                block = weakBottom ? weakTiles : strongTiles
            end
            Ahorn.drawImage(ctx, block, (i - 1) * 8, (j - 1) * 8, tx, ty, 8, 8)
        end
    end

    topleft = strongTiles
    topright = strongTiles
    bottomleft = strongTiles
    bottomright = strongTiles

    # topleft corner
    if weakLeft && weakTop
        topleft = weakTiles
    elseif weakLeft && !weakTop
        topleft = weakCornersHTiles
    elseif !weakLeft && weakTop
        topleft = weakCornersVTiles
    end
    # topright corner
    if weakRight && weakTop
        topright = weakTiles
    elseif weakRight && !weakTop
        topright = weakCornersHTiles
    elseif !weakRight && weakTop
        topright = weakCornersVTiles
    end
    # bottomleft corner
    if weakLeft && weakBottom
        bottomleft = weakTiles
    elseif weakLeft && !weakBottom
        bottomleft = weakCornersHTiles
    elseif !weakLeft && weakBottom
        bottomleft = weakCornersVTiles
    end
    # bottomright corner
    if weakRight && weakBottom
        bottomright = weakTiles
    elseif weakRight && !weakBottom
        bottomright = weakCornersHTiles
    elseif !weakRight && weakBottom
        bottomright = weakCornersVTiles
    end

    Ahorn.drawImage(ctx, topleft, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, topright, width - 8, 0, 24, 0, 8, 8)
    Ahorn.drawImage(ctx, bottomleft, 0, height - 8, 0, 24, 8, 8)
    Ahorn.drawImage(ctx, bottomright, width - 8, height - 8, 24, 24, 8, 8)

    Ahorn.drawImage(ctx, eyeSprite, div(width - eyeSprite.width, 2), div(height - eyeSprite.height, 2))
end

end


text:
# Queue
# =====

mutable struct Queue{T}
    data::Vector{T}
    offset::Int
    first::Int
    last::Int
end

function Queue{T}(bufsize::Integer=2^4) where {T}
    if bufsize ≤ 0
        throw(ArgumentError("buffer size must be positive"))
    elseif !ispow2(bufsize)
        throw(ArgumentError("buffer size must be a power of two"))
    end
    return Queue(Vector{T}(undef, bufsize), 0, 1, 0)
end

function Base.eltype(::Type{Queue{T}}) where T
    return T
end

function Base.length(queue::Queue)
    return queue.last - queue.first + 1
end

function Base.isempty(queue::Queue)
    return queue.first > queue.last
end

function Base.push!(queue::Queue{T}, elm::T) where T
    if length(queue.data) < length(queue) + 1
        index_first = dataindex(queue, queue.first)
        index_last = dataindex(queue, queue.last)
        index_end = lastindex(queue.data)
        # NOTE: resize factor must be a power of two
        resize!(queue.data, 2 * length(queue.data))
        @assert ispow2(length(queue.data))
        if !isempty(queue) && index_last < index_first
            # make the circular data linear
            copyto!(queue.data, index_end + 1, queue.data, 1, index_last)
        end
        copyto!(queue.data, 1, queue.data, index_first, length(queue))
        queue.offset = 0
    end
    queue.data[dataindex(queue, queue.last + 1)] = elm
    queue.last += 1
    return queue
end

function Base.popfirst!(queue::Queue)
    if isempty(queue)
        throw(ArgumentError("empty"))
    end
    elm = queue.data[dataindex(queue, queue.first)]
    queue.first += 1
    queue.offset += 1
    return elm
end

function Base.iterate(queue::Queue)
    return iterate(queue, queue.first)
end

function Base.iterate(queue::Queue, i)
    if !(queue.first <= i <= queue.last)
        return nothing
    else
        queue.data[dataindex(queue, i)], i + 1
    end
end

function dataindex(queue::Queue, i::Integer)
    return ((i - queue.first + queue.offset) & (length(queue.data) - 1)) + 1
end

function Base.firstindex(queue::Queue)
    return queue.first
end

function Base.lastindex(queue::Queue)
    return queue.last
end

function Base.getindex(queue::Queue, i::Integer)
    if !(queue.first <= i <= queue.last)
        throw(BoundsError())
    end

    return queue.data[dataindex(queue, i)]
end


text:
# Implement random turnover of network components
# Alex Tam, 12/10/2020

"Random turnover of actin filaments"
function actin_turnover(s, af, mm, parN, parA, Lxx, Lyy)
    N = length(af); j = 0; # Indexing pre-allocation
    Pt = 1 - exp(-parA.k_off*parN.dt); # Probability that a filament will turn over in the time step
    # Remove actin filaments
    for i = 1:N
        j += 1; # Index (dynamic) of current filament
        random = rand(); # Random number [0, 1]
        if random <= Pt
            # Remove data from state and filament list
            deleteat!(s.an, j)
            deleteat!(af, j)
            # Update motors
            for m in mm
                # Mark motors for removal if filament turns over
                if m.f1.index == j
                    m.to_remove = true;
                end
                if m.f2.index == j
                    m.to_remove = true;
                end
                # Update index of filaments associated with motors
                if m.f1.index > j
                    m.f1.index -= 1; # Reduce index by 1 if necessary
                end
                if m.f2.index > j
                    m.f2.index -= 1; # Reduce index by 1 if necessary
                end
            end
            # Update indices in actin filament list
            for a = 1:length(af)
                af[a].index = a;
            end
            j -= 1; # Ensure correct index for next filament
        end
    end
    # Add new actin filaments
    nNew = parN.nA - length(af); # Immediately replace filaments
    # nNew = rand(Poisson(parN.nA*parA.k_off*parN.dt)); # Draw number of filaments to add from Poisson distribution
    nOld = length(af); # Initial number of filaments for indexing
    for i = 1:nNew
        push!(af, Actin_Filament(s, nOld+i, parA, Lxx, Lyy)); # Create new filament and state data
    end
    return af, mm, s
end

"Random turnover of myosin motors"
function myosin_turnover(s, xl, mm, parN, parM, Lxx, Lxy, Lyx, Lyy)
    N = length(mm); j = 0;
    for i = 1:N
        j += 1; # Index (dynamic) of current motor
        random = 1; Pt = 0; # Pre-allocate
        if mm[j].to_remove == false
            random = rand(); # Generate random number [0, 1] for current motor
            x1, y1, x2, y2 = get_motor_pos(mm[j], s, Lxx, Lxy, Lyx, Lyy); # Extract dimensional motor positions
            Lm = sqrt((x1-x2)^2 + (y1-y2)^2); # Compute motor length
            off_rate = parM.k_off*exp(parM.k*Lm/parM.F_ref); # Compute off rate of current motor using Bell's Law
            Pt = 1 - exp(-off_rate*parN.dt); # Probability that random turnover will occur at least once in the time step
        end
        # Replace motors
        if any(s.mp[j] .>= 1) || (random <= Pt) || (mm[j].to_remove == true)
            xl_index = rand(1:length(xl));
            # Delete data from State and motor list
            deleteat!(s.mp, j);
            deleteat!(mm, j);
            # Generate replacement motor
            push!(mm, Myosin_Motor(s, xl[xl_index], j, Lxx, Lxy, Lyx, Lyy));
            deleteat!(xl, xl_index);
            # Update order of motor indices
            for m = 1:N
                mm[m].index = m;
            end
            j -= 1; # Ensure correct index for next motor
        end
    end
    return mm, xl, s
end


text:
function equivalencenot_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        CP.EquivalenceNot{MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x2]),
        CP.EquivalenceNot(MOI.EqualTo(1), MOI.EqualTo(1)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) != 1
    end
end

function equivalencenot_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.EquivalenceNot{MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOIU.vectorize([x1, x2]),
        CP.EquivalenceNot(MOI.EqualTo(1), MOI.EqualTo(1)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) != 1
    end
end

const equivalencenottests = Dict(
    "equivalencenot_singlevariable" => equivalencenot_singlevariable_test,
    "equivalencenot_scalaraffinefunction" =>
        equivalencenot_scalaraffinefunction_test,
)

MOIT.@moitestset equivalencenot


text:
using Combinatorics
using DelimitedFiles

function read_data_unweighted(dataset::String, order::Int64, exact_match::Bool=false)
    if order < 3 || order > 5; error("Only support for order-3,4,5 tensors"); end

    read(filename::String) = convert(Vector{Int64}, readdlm(filename, Int64)[:, 1])
    name = match(r"\\([^\\]+)$",dataset).captures[1]
    simplices = read("$dataset\\$name-simplices.txt")
    nverts = read("$dataset\\$name-nverts.txt")

    curr_ind = 1
    S = Set{NTuple{order,Int64}}()
    for nvert in nverts[1:10000]
        simp = sort(simplices[curr_ind:(curr_ind + nvert - 1)])
        curr_ind += nvert
        if (exact_match && nvert != order); continue; end
        for hedge in combinations(simp, order)
            push!(S, NTuple{order,Int64}(hedge))
        end
    end

    I = Vector{Vector{Int64}}(undef, order)
    for j in 1:order; I[j] = Vector{Int64}(); end
    for hedge in S
        for (j, v) in enumerate(hedge)
            push!(I[j], v)
        end
    end
    dim = maximum([maximum(I[j]) for j in 1:order])
    V = ones(Float64, length(I[1])) / length(I[1])
    if order == 3; return SymTensor3(I[1], I[2], I[3], V, dim); end
    if order == 4; return SymTensor4(I[1], I[2], I[3], I[4], V, dim); end
    if order == 5; return SymTensor5(I[1], I[2], I[3], I[4], I[5], V, dim); end
end

function read_node_labels(dataset::String)
    labels = String[]
    open("$folder\\$name-node-labels.txt") do f
        for line in eachline(f)
            data = split(strip(line))
            push!(labels, join(data[2:end], " "))
        end
    end
    return labels
end
;

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
@testset "4.6.3.1 (a+b csc)^m (d csc)^n (A+B csc)" begin
    (A, a, c, d, x, ) = @variables A a c d x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*csc(e+f*x))^m*(d*csc(e+f*x))^n*(A+B*csc(e+f*x))=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+a*csc(e+f*x))^m*(d*csc(e+f*x))^n*(A+B*csc(e+f*x))=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A+B*csc(c+d*x))*(a+a*csc(c+d*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A+B*csc(c+d*x))*(a+a*csc(c+d*x))^(n/2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A+A*csc(c+d*x))*(a+a*csc(c+d*x))^n=#

    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^3, x, 7, -((7*a*A*atanh(cos(c + d*x)))/(8*d)) - (2*a*A*cot(c + d*x))/d - (2*a*A*cot(c + d*x)^3)/(3*d) - (7*a*A*cot(c + d*x)*csc(c + d*x))/(8*d) - (a*A*cot(c + d*x)*csc(c + d*x)^3)/(4*d)]
    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^2, x, 7, -((a*A*atanh(cos(c + d*x)))/d) - (5*a*A*cot(c + d*x))/(3*d) - (a*A*cot(c + d*x)*csc(c + d*x))/d - (a*A*cot(c + d*x)*csc(c + d*x)^2)/(3*d)]
    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^1, x, 6, -((3*a*A*atanh(cos(c + d*x)))/(2*d)) - (2*a*A*cot(c + d*x))/d - (a*A*cot(c + d*x)*csc(c + d*x))/(2*d)]
    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*sin(c + d*x)^1, x, 5, 2*a*A*x - (a*A*atanh(cos(c + d*x)))/d - (a*A*cos(c + d*x))/d]
    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*sin(c + d*x)^2, x, 5, (3*a*A*x)/2 - (2*a*A*cos(c + d*x))/d - (a*A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [(A + A*csc(c + d*x))*(a + a*csc(c + d*x))^1*sin(c + d*x)^3, x, 7, a*A*x - (2*a*A*cos(c + d*x))/d + (a*A*cos(c + d*x)^3)/(3*d) - (a*A*cos(c + d*x)*sin(c + d*x))/d]

    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^3, x, 4, -((a*A*atanh(cos(c + d*x)))/(8*d)) - (a*A*cot(c + d*x)*csc(c + d*x))/(8*d) + (a*A*cot(c + d*x)*csc(c + d*x)^3)/(4*d)]
    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^2, x, 3, (a*A*cot(c + d*x)^3)/(3*d)]
    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^1, x, 3, -((a*A*atanh(cos(c + d*x)))/(2*d)) + (a*A*cot(c + d*x)*csc(c + d*x))/(2*d)]
    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*sin(c + d*x)^1, x, 4, (a*A*atanh(cos(c + d*x)))/d - (a*A*cos(c + d*x))/d]
    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*sin(c + d*x)^2, x, 3, (-(1/2))*a*A*x - (a*A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [(A + A*csc(c + d*x))*(a - a*csc(c + d*x))^1*sin(c + d*x)^3, x, 3, (a*A*cos(c + d*x)^3)/(3*d)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A-A*csc(c+d*x))*(a+a*csc(c+d*x))^n=#

    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^3, x, 4, -((a*A*atanh(cos(c + d*x)))/(8*d)) - (a*A*cot(c + d*x)*csc(c + d*x))/(8*d) + (a*A*cot(c + d*x)*csc(c + d*x)^3)/(4*d)]
    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^2, x, 3, (a*A*cot(c + d*x)^3)/(3*d)]
    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1*csc(c + d*x)^1, x, 3, -((a*A*atanh(cos(c + d*x)))/(2*d)) + (a*A*cot(c + d*x)*csc(c + d*x))/(2*d)]
    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1/csc(c + d*x)^1, x, 4, (a*A*atanh(cos(c + d*x)))/d - (a*A*cos(c + d*x))/d]
    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1/csc(c + d*x)^2, x, 3, (-(1/2))*a*A*x - (a*A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [(A - A*csc(c + d*x))*(a + a*csc(c + d*x))^1/csc(c + d*x)^3, x, 3, (a*A*cos(c + d*x)^3)/(3*d)]

    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^3, x, 7, -((7*a*A*atanh(cos(c + d*x)))/(8*d)) + (2*a*A*cot(c + d*x))/d + (2*a*A*cot(c + d*x)^3)/(3*d) - (7*a*A*cot(c + d*x)*csc(c + d*x))/(8*d) - (a*A*cot(c + d*x)*csc(c + d*x)^3)/(4*d)]
    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^2, x, 7, (a*A*atanh(cos(c + d*x)))/d - (5*a*A*cot(c + d*x))/(3*d) + (a*A*cot(c + d*x)*csc(c + d*x))/d - (a*A*cot(c + d*x)*csc(c + d*x)^2)/(3*d)]
    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1*csc(c + d*x)^1, x, 6, -((3*a*A*atanh(cos(c + d*x)))/(2*d)) + (2*a*A*cot(c + d*x))/d - (a*A*cot(c + d*x)*csc(c + d*x))/(2*d)]
    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1/csc(c + d*x)^1, x, 5, -2*a*A*x - (a*A*atanh(cos(c + d*x)))/d - (a*A*cos(c + d*x))/d]
    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1/csc(c + d*x)^2, x, 5, (3*a*A*x)/2 + (2*a*A*cos(c + d*x))/d - (a*A*cos(c + d*x)*sin(c + d*x))/(2*d)]
    @test_int [(A - A*csc(c + d*x))*(a - a*csc(c + d*x))^1/csc(c + d*x)^3, x, 7, (-a)*A*x - (2*a*A*cos(c + d*x))/d + (a*A*cos(c + d*x)^3)/(3*d) + (a*A*cos(c + d*x)*sin(c + d*x))/d]

    #= ::Section:: =#
    #=Integrands*of*the*form*(a+a*csc(e+f*x))^m*(d*sin(e+f*x))^n*(A+B*csc(e+f*x))=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*csc(e+f*x))^m*(d*csc(e+f*x))^n*(A+B*csc(e+f*x))=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A+B*csc(c+d*x))*(a+b*csc(c+d*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(c+d*x)^m*(A+B*csc(c+d*x))*(a+b*csc(c+d*x))^(n/2)=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(a+b*csc(e+f*x))^m*(d*sin(e+f*x))^n*(A+B*csc(e+f*x))=#
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.3.1 (a+b csc)^m (d csc)^n (A+B csc).jl


text:
using DelimitedFiles

"""Load union-find test data."""
function loadufdata(datapath)
    rawdata = readdlm(datapath)
    nv = Int(rawdata[1, 1])
    edges = Tuple{Int, Int}[]
    for (d1, d2) in eachrow(rawdata[2:end, :])
        push!(edges, (Int(d1), Int(d2)).+1)
    end
    (nv=nv, edges=edges)
end

"""Load edge-weighted graph data."""
function loadewg(datapath)
    rawdata = readdlm(datapath)
    nv = Int(rawdata[1, 1])
    ne = Int(rawdata[2, 1])
    edges = Tuple{Int, Int, Float64}[]
    for (d1, d2, d3) in eachrow(rawdata[3:end, :])
        push!(edges, (Int(d1)+1, Int(d2)+1, d3))
    end
    (nv=nv, ne=ne, edges=edges)
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
# from grid/grid.jl
__precompile__()

module GridArrays

include("ModCartesianIndices.jl")
using ..ModCartesianIndicesBase

using DomainSets, StaticArrays, RecipesBase, Test, FastGaussQuadrature, GaussQuadrature, FillArrays

using DomainSets: endpoints

import Base: *, size, length, @propagate_inbounds, step, ndims, unsafe_getindex,
    checkbounds, IndexStyle, ==, ≈, getindex, eachindex, convert, in, ^, string, axes, convert
import Base.Broadcast: broadcast

import DomainSets: cartesianproduct, element, elements, numelements, ×, cross, minimum, maximum, dimension

export AbstractGrid, AbstractGrid1d, AbstractGrid3d,
        AbstractEquispacedGrid, EquispacedGrid, PeriodicEquispacedGrid,
        FourierGrid, MidpointEquispacedGrid, RandomEquispacedGrid,
        AbstractIntervalGrid, eachelement, ScatteredGrid, ×, cartesianproduct,
        TensorSubGrid, instantiate, support, float_type, isperiodic,
        boundary, subgrid, mask, randomgrid, boundingbox, TensorSubGrid,
        iscomposite, dimension, prectype, ChebyshevTNodes, ChebyshevUNodes,
        LaguerreNodes, HermiteNodes, LegendreNodes, JacobiNodes
export ChebyshevNodes, ChebyshevGrid, ChebyshevPoints, ChebyshevExtremae, ×

# from grid/productgrid.jl
export ProductGrid

# from grid/subgrid.jl
export subindices, supergrid, issubindex, similar_subgrid

# from grid/mappedgrid.jl
export MappedGrid, mapped_grid, apply_map

import Base: isapprox

# TODO move to domainsets
isapprox(d1::DomainSets.ProductDomain,d2::DomainSets.ProductDomain) =
    reduce(&, map(isapprox,DomainSets.elements(d1), DomainSets.elements(d2)))

"Assign a floating point type to a domain element type T."
float_type(::Type{T}) where {T <: Real} = T
float_type(::Type{Complex{T}}) where {T <: Real} = Complex{T}
float_type(::Type{SVector{N,T}}) where {N,T} = T
float_type(::Type{NTuple{N,T}}) where {N,T} = T

# Fallback: we return Float64
float_type(::Type{T}) where {T} = Float64

include("grid.jl")
include("productgrid.jl")
include("intervalgrids.jl")
include("mappedgrid.jl")
include("scattered_grid.jl")

include("domains/boundingbox.jl")
include("domains/broadcast.jl")

include("randomgrid.jl")

include("subgrid/AbstractSubGrids.jl")

include("recipes.jl")

export test_generic_grid, test_interval_grid
export grid_iterator1, grid_iterator2
include("test/test_grids.jl")

≈(d1::DomainSets.Interval, d2::DomainSets.Interval) =
    1≈1+abs(DomainSets.infimum(d1)-DomainSets.infimum(d2))+abs(DomainSets.supremum(d1)-DomainSets.supremum(d2))

end # module


text:
module Scratchpad

sp = zeros(UInt8, 256)

function store(sX, ss) #Should store value sX, into ss. #Both sX, and ss are passed as integer value
    sp[(0xff & ss) + 1] = sX
end

function fetch(ss) #Should return the value at address ss #ss is passed by value
    return sp[ss + 1]
end

end

max_stars_repo_name: kinoroy/KCPSM-Simulator

max_issues_repo_name: kinoroy/KCPSM-Simulator

max_forks_repo_name: kinoroy/KCPSM-Simulator


text:
module Kepler

export rv_model

kepler_f(M, E, e) = E - e*sin(E) - M
kepler_fp(E, e) = one(E) - e*cos(E)
kepler_fpp(E,e) = e*sin(E)
kepler_fppp(E,e) = e*cos(E)

function kepler_solve_ea(n, e, t)
    M = mod(n*t, 2.0*pi)

    if M < pi
        E = M + 0.85*e
    else
        E = M - 0.85*e
    end

    f = kepler_f(M, E, e)

    while abs(f) > 1e-8
        fp = kepler_fp(E,e)
        disc = sqrt(abs(16.0*fp*fp - 20.0*f*kepler_fpp(E,e)))
        if fp > 0.0
            d = -5.0*f/(fp + disc)
        else
            d = -5.0*f/(fp - disc)
        end

        E += d
        f = kepler_f(M, E, e)
    end

    E
end

function kepler_solve_ta(n, e, t)
    E = kepler_solve_ea(n, e, t)

    f = 2.0*atan(sqrt((one(e)+e)/(one(e)-e))*tan(E/2.0))

    if f < 0.0
        f += 2.0*pi
    end

    f
end

"""
`rv_model(ts, K, e, omega, chi, P)`

Returns the radial velocity at times `ts` for the system with
semi-amplitude `K`, eccentricity `e`, argument of pericentre `omega`,
phase at a `t=0` of `2*pi*chi`, and period `P`.  These parameter
arguments can also be vectorised, in which case the returned array has
size `(ntimes, nplanets)`.
"""
function rv_model(ts::Array{Float64, 1},
                  K::Float64,
                  e::Float64,
                  omega::Float64,
                  chi::Float64,
                  P::Float64)

    nts = size(ts, 1)

    rvs = zeros(nts)

    t0 = -chi*P
    n = 2.0*pi/P

    ecw = e*cos(omega)

    for i in 1:nts
        t = ts[i]
        f = kepler_solve_ea(n, e, t-t0)

        rvs[i] = K*(cos(f + omega) + ecw)
    end

    rvs
end

function rv_model(ts::Array{Float64, 1},
                  Ks::Array{Float64, 1},
                  es::Array{Float64, 1},
                  omegas::Array{Float64, 1},
                  chis::Array{Float64, 1},
                  Ps::Array{Float64, 1})
    npl = size(Ks, 1)
    nts = size(ts, 1)

    rvs = zeros(nts, npl)

    for j in 1:npl
        rvs[:,j] = rv_model(ts, Ks[j], es[j], omegas[j], chis[j], Ps[j])
    end

    rvs
end

end

max_stars_repo_name: farr/RVChallenge

max_issues_repo_name: farr/RVChallenge

max_forks_repo_name: farr/RVChallenge


text:
using MultilinearOpt, JuMP, Gurobi

N = 1:11
S = 1:5
I = 6:7
T = 8:11

K = 1:4

#    1   2   3   4   5   6   7   8   9   10   11
c = [0   0   0   0   0   7   0   0   0    0    0  # 1
     0   0   0   0   0   3   0   0   0    0    0  # 2
     0   0   0   0   0   0   2   0   0    0    0  # 3
     0   0   0   0   0   0  10   0   0    0    0  # 4
     0   0   0   0   0   0   5   0   0    0    0  # 5
     0   0   0   0   0   0   0 -16 -25  -15  -10  # 6
     0   0   0   0   0   0   0 -16 -25  -15  -10] # 7

# A = zeros(11,11);
# A[1,6] = A[2,6] = A[3,7] = A[4,7] = A[5,7] = A[6,8] = A[6,9] = A[6,10] = A[6,11] = A[7,8]= A[7,9]= A[7,10]= A[7,11] = 1

A = [(1,6),(2,6),(3,7),(4,7),(5,7),(6,8),(6,9),(6,10),(6,11),(7,8),(7,9),(7,10),(7,11)]

N⁺(i) = [j for j in N if (i,j) in A]
N⁻(i) = [j for j in N if (j,i) in A]

q = [1.00    6.00    4.00    0.50
     4.00    1.00    3.00    2.00
     4.00    5.50    3.00    0.90
     3.00    3.00    3.00    1.00
     1.00    2.70    4.00    1.60
     0.00    0.00    0.00    0.00
     0.00    0.00    0.00    0.00
     3.00    3.00    3.25    0.75
     4.00    2.50    3.50    1.50
     1.50    5.50    3.90    0.80
     3.00    4.00    4.00    1.80]

b = [75,75,75,75,75,75,75,10,25,30,10]

m = Model(with_optimizer(Gurobi.Optimizer))
@variable(m, 0 ≤ f[i in N, j in N; (i,j) in A] ≤ maximum(b))
for i in N
    if i in S
        @constraint(m, sum(f[i,j] for j in N⁺(i)) ≤ b[i])
    else
        @constraint(m, sum(f[j,i] for j in N⁻(i)) ≤ b[i])
    end
    @constraint(m, sum(f[i,j] for j in N⁺(i)) == sum(f[j,i] for j in N⁻(i)))
end
@variable(m, 0 ≤ w[i=N,k=K] ≤ q[i,k])
for i in N, k in K
    if i in S
        @constraint(m, w[i,k] == q[i,k])
    else
        @constraint(m, w[i,k]*sum(f[j,i] for j in N⁻(i)) == sum(w[j,k]*f[j,i] for j in N⁻(i)))
    end
end

@objective(m, Min, sum(c[i,j]*f[i,j] for (i,j) in A))

relaxbilinear!(m, method=:Logarithmic2D)
optimize!(m)
stat = JuMP.termination_status(m)


text:
function rankstandard{T<:Real}(a::Array{T,1})
    r = collect(1:length(a))
    1 < r[end] || return r
    for i in ties(a)
        r[a.==i] = r[a.==i][1]
    end
    return r
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module ScrapeSEC

    # Dependencies
    using DataFrames
    import Dates
    import CSV
    import HTTP
    using ProgressMeter
    import ZipFile

    # source files
    include("download_metadata.jl")
    include("main_index.jl")
    include("download_filings.jl")

    export
    # metadata functions
    download_metadata_files, create_main_index,

    # filing downloaders
    get_quarterly_filings, download_filing

end


text:
module Profiles

using ..Types: CONFIG_DIR
using Base: UUID

using Pkg.TOML
import Base.Dict

struct Profile
    uuid::UUID
    name::Union{AbstractString,Nothing}
    date::Union{AbstractString,Nothing}
    about::Union{AbstractString,Nothing}
    homepage::Union{AbstractString,Nothing}
    email::Union{AbstractString,Nothing}
    facebook::Union{AbstractString,Nothing}
    twitter::Union{AbstractString,Nothing}
    github::Union{AbstractString,Nothing}
end

function Dict(profile::Profile)
    dict = Dict()

    dict["uuid"] = string(profile.uuid)
    isnothing(profile.name) || (dict["name"] = profile.name)
    isnothing(profile.date) || (dict["date"] = profile.date)
    isnothing(profile.about) || (dict["about"] = profile.about)
    isnothing(profile.homepage) || (dict["homepage"] = profile.homepage)
    isnothing(profile.email) || (dict["email"] = profile.email)
    isnothing(profile.facebook) || (dict["facebook"] = profile.facebook)
    isnothing(profile.twitter) || (dict["twitter"] = profile.twitter)
    isnothing(profile.github) || (dict["github"] = profile.github)

    return dict
end

### I could use this for the tests with something like Profile(Dict("name"=>"account1"))
function Profile(dict::Dict)
    uuid = UUID(dict["uuid"])
    haskey(dict,"name") ? (name=dict["name"]) : (name=nothing)
    haskey(dict,"date") ? (date=dict["date"]) : (date=nothing)
    haskey(dict,"about") ? (about=dict["about"]) : (about=nothing)
    haskey(dict,"homepage") ? (homepage=dict["hompage"]) : (homepage=nothing)
    haskey(dict,"email") ? (email=dict["email"]) : (email=nothing)
    haskey(dict,"facebook") ? (facebook=dict["facebook"]) : (facebook=nothing)
    haskey(dict,"twitter") ? (twitter=dict["twitter"]) : (twitter=nothing)
    haskey(dict,"github") ? (twitter=dict["github"]) : (github=nothing)

    # One may also add membership of each deme as part of the profile additionally with signatures in a envelope form.

    return Profile(uuid,name,date,about,homepage,email,facebook,twitter,github)
end

function save(fname::AbstractString,profile::Profile)
    mkpath(dirname(fname))
    dict = Dict(profile)
    open(fname, "w") do io
        TOML.print(io, dict)
    end
end

save(profile::Profile) = save(CONFIG_DIR * "/Profile.toml",profile)

function Profile(fname::AbstractString)
    @assert isfile(fname) "No profile file found"
    dict = TOML.parsefile(fname)
    return Profile(dict)
end

Profile() = Profile(CONFIG_DIR * "/Profile.toml")

end


text:
include("../src/Rhapsodie.jl")
#include("test_data_simulation.jl")

tau=0.1;

par=Rhapsodie.readdlm("test_results/Parameters.txt")
DSIZE=Int64(par[1]);
NTOT=Int64(par[2]);
Nframe=Int64(par[3]);
Nrot=Int64(par[4]);
Nangle=NTOT÷(Nframe*4)
Center=par[5:6];
DerotAng=Vector{Float64}();
Epsilon=Vector{Tuple{Array{Float64,1},Array{Float64,1}}}();

for iter=1:NTOT
    ind=div(iter-1, NTOT/4)+1
    push!(Epsilon,([0. ,0. ],par[end-1:end]));
    push!(DerotAng, deg2rad(-par[7]));
end

psf_center=Rhapsodie.readdlm("../data/PSF_centers_Airy.txt");

Rhapsodie.load_parameters((DSIZE, 2*DSIZE, NTOT), Nframe, Nrot, Nangle, Center, (psf_center[1:2], psf_center[3:4]), Epsilon, DerotAng)

Rhapsodie.Load_Data("test_results/DATA_$tau-$DSIZE.fits",
                    "test_results/WEIGHT_$tau-$DSIZE.fits")

PSF=Rhapsodie.read(Rhapsodie.FitsArray, "../data/PSF_parametered_Airy.fits");
const A=Rhapsodie.set_fft_op(PSF[1:end÷2,:]'[:,:],psf_center[1:2]);

X0 = Rhapsodie.read("stokes", "test_results/Results_Separable_Linear_$tau-$DSIZE.fits");

regularisation_parameters = [0.5 , 0. , -1., -3.]; #(in log10)

@time x = Rhapsodie.apply_rhapsodie(X0, A, Rhapsodie.dataset, regularisation_parameters,
                          maxeval=3000, maxiter=3000);

Rhapsodie.write(x, "test_results/RHAPSODIE_nonlinearresults_$tau-$DSIZE.fits")


text:
using SafeTestsets

@time begin
    @time @safetestset "BLAS" begin include("gemm_tests.jl") end
end


text:
function ipopt_style_inertia_correction!(newt::abstract_newton_direction, vars::class_variables, settings::class_settings)
	try
    delta_style! = update_newton_diag!;

		MAX_IT = 25;
		j = 1;
		new_delta = 0.0;

		# try delta = 0.0
    old_delta = newt.delta;
    newt.delta = new_delta;
    update_newton!(newt, vars, settings);
    inertia = delta_style!(newt, vars, settings)

		if inertia == 1
			new_delta = 0.0; # no need to modify newton system
		else
			newt.delta = old_delta;

      if newt.delta <= settings.delta_min;
            newt.delta = settings.delta_start;
      end

			for j = 2:MAX_IT
				inertia = delta_style!(newt, vars, settings)

				if inertia == 1
            new_delta = newt.delta * settings.delta_decrease
            break
				elseif inertia == 0
            newt.delta = newt.delta * settings.delta_increase;
				elseif inertia == -1
					  error("numerical stability issues when computing KKT system !!!")
				else
					  error("inertia_corection")
				end
			end
		end

		if j == MAX_IT
			error("maximum iterations for inertia_corection reached")
		else
      form_woodbury!(newt, vars)
			return new_delta, j
		end

	catch e
		println("ERROR in ipopt_style_inertia_correction")
		throw(e)
	end
end

function iterative_trust_region!(newt::abstract_newton_direction, vars::class_variables, settings::class_settings)
    # Goal: choose δ such that Δ = |log(x) - log(x + d_x)|_∞ = 0.5
    # or |d_x/x|_∞ <= 0.5
    # or |d_x^2/x^2| <= 1.0
    #
    #
    # compute d'
    try
      newt.delta = 0.0;
      inertia = update_newton!(newt, vars, settings);

      if inertia
          best_δ = 0.0;
      else
          target = 2.0;
          best_val = Inf;
          δ_vals = 2.^linspace(0.0,4.0,20)

          for δ = δ_vals
              newt.delta = δ
              inertia = update_newton_diag!(newt, vars, settings);
              if inertia
                  form_woodbury!(newt, vars)

									try
		                  compute_newton_direction!(newt, vars, class_theta(0.99,0.5,0.5));
		                  alpha_max = maximum_step(vars, newt.direction);

		                  val = abs(log(alpha_max) - log(target));
		                  if best_val > val
		                      best_val = val;
		                      best_δ = δ;
		                  end
									catch e
											warn(e.msg)
											warn("computation skipped")
									end
              end
          end

          @assert(best_val < Inf)
      end

      form_woodbury!(newt, vars)
      return best_δ, 1
    catch e
        println("ERROR in delta.iterative_trust_region!")
        throw(e)
	  end
end

#vars = class_variables(1,1);
#dir = deepcopy(vars);
#x(dir,[-0.5])
#alpha_max = maximum_step(vars,dir)

max_forks_repo_name: ohinder/Everest


text:
"""
    MOKernel

Abstract type for kernels with multiple outpus.
"""
abstract type MOKernel <: Kernel end


text:
# This is test of the code for computing a transit model integrated over a time step.
# Computes derivatives over the timestep.
using PyPlot

include("../../../src/integrate_transit_simpson.jl")
#include("../../../src/integrate_transit_cubature.jl")

# Test it out:

include("../../../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 1000; dt = (t2-t1)/nt*100
t = zeros(nt)
t .= linearspace(t1,t2,nt)
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(nt,5+nu)
favg1 = zeros(nt,5+nu)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

@time for i=1:nt
        trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
        favg0[i,1] =  transit_poly!(trans)
        favg0[i,2] =  trans.dfdrb[1]
        favg0[i,3] =  trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t[i])
        favg0[i,4] =  trans.dfdrb[2]*param[2]/trans.b*(t[i]-param[1])^2
        favg0[i,5] =  trans.dfdrb[2]*param[3]/trans.b
        favg0[i,6:5+nu] = trans.dfdu
      end

function integrate_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},dt::T,favg1::Array{T,2},nt::Int64,tol::T,maxdepth::Int64) where {T <: Real}
dtinv = inv(dt)
for i=1:nt
  ftmp = integrate_timestep_gradient(param,trans,t[i]-0.5*dt,t[i]+0.5*dt,tol*trans.r^2,maxdepth)*dtinv
  favg1[i,1:5]=ftmp[1:5]
  # Convert from d_n to u_n derivatives:
  favg1[i,6:5+nu]=BLAS.gemv('T',1.0,trans.dgdu,ftmp[6:6+nu])
end
return
end

@time integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-5,8)

#@time for i=1:nt
#        ftmp = integrate_timestep_gradient(param,trans,t[i],dt,1e-6*r^2,32)*dtinv
#        favg1[i,1:5]=ftmp[1:5]
#        # Convert from d_n to u_n derivatives:
#        favg1[i,6:5+nu]=BLAS.gemv('T',1.0,trans.dgdu,ftmp[6:6+nu])
#      end

# Now plot the results:
fig,axes = subplots(2,3,figsize=(12,8))
ax = axes[1]
ax[:plot](t,favg0[:,1],linestyle="--",linewidth=2,label=L"$\mathcal{F}$")
ax[:plot](t,favg1[:,1],linewidth=2,label=L"$\overline{\mathcal{F}}$")
ax[:plot]([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
ax[:annotate](L"$\Delta t$",xy=[-0.08,0.991])
ax[:legend](loc="upper center")
ax[:axis]([-1.5,1.5,0.988,1.002])
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux")
ax = axes[2]
ax[:plot](t,favg0[:,2]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial r}$")
ax[:plot](t,favg1[:,2]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial r}$")
ax[:legend](loc="upper center")
ax[:axis]([-1.5,1.5,-30,2])
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax = axes[3]
ax[:plot](t,favg0[:,3]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial t_0}$")
ax[:plot](t,favg1[:,3]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial t_0}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
#ax[:axis]([-1.5,1.5,-1,6])
ax[:legend](loc="lower left")
ax = axes[4]
ax[:plot](t,favg0[:,4]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial v}$")
ax[:plot](t,favg1[:,4]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial v}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax[:axis]([-1.5,1.5,-1,6])
ax[:legend](loc="upper center")
ax = axes[5]
ax[:plot](t,favg0[:,5]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial b_0}$")
ax[:plot](t,favg1[:,5]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial b_0}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax[:legend](loc="upper center")
ax = axes[6]
ax[:plot](t,favg0[:,6]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial u_1}$")
ax[:plot](t,favg1[:,6]*1000.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial u_1}$")
ax[:plot](t,favg0[:,7]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial u_2}$")
ax[:plot](t,favg1[:,7]*1000.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial u_2}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [ppt]")
ax[:legend](loc="upper center")
#b = sqrt.(param[3]^2+(param[2]*(t-param[1])).^2)
#plot(dt^2*(d2fdb2.*(param[2]*(t-param[1])./b).^2+dfdb.*param[2]^2./b)/24.0)
savefig("integrate_transit_gradient.pdf", bbox_inches="tight")


text:
mf = (x -> 100.0, x -> -100.0)
rsampler = GaussianRewards(mf)
@test rsampler(rand(1, 1), [1])[1, 1] > 90
@test rsampler(rand(1, 1), [2])[1, 1] < -90


text:
#  Unit "QnLogLike.jl" of the Diagonalization.jl Package for Julia language
#
#  MIT License
#  Copyright (c) 2020-2021,
#  Marco Congedo°, Ronald Phlypo, CNRS, UGA, Grenoble-INP, France
#  Alexandre Gramfort¨, INRIA, U. Paris Saclay, France
#  ° https://sites.google.com/site/marcocongedo/
#  ¨ http://alexandre.gramfort.net/

# ? CONTENTS :
#  Quasi-Newton Approximate Joint Diagonalization (AJD) algorithm of:
#  P. Ablin, J.F. Cardoso and A. Gramfort. Beyond Pham's algorithm
#  for joint diagonalization. Proc. ESANN 2019.
#  https://hal.archives-ouvertes.fr/hal-01936887v1
#
#  Code adapted in Julia from Python code provided from the authors at:
#  https://github.com/pierreablin/qndiag/blob/master/qndiag/qndiag.py
#
#  NOTATION
#  lower-case letter: a scalar, e.g.: a
#  Upper case letter: a matrix, e.g., A
#  Bold lover case letter: a vector of matrices, e.g., 𝐀
#
#  ADVANCED QnLogLike algorithm (no PRIMITIVE version is available):
#  The algorithm takes as input a vector of k real symmetric matrices
#  𝐂 and find a non-singular matrix B such that the congruences
#  B'*𝐂_κ*B are as diagonal as possible for all κ=1:k.
#
#  `w` is an optional vector of k positive weights for each matrix in 𝐂.
#  if `w` is different from `nothing` (default), the input matrices are
#  weighted with these weights.
#  A function can be passed as the `w` argument, in which case the kth weight
#  is found as the output of the function applied to the kth matrix in 𝐂.
#  A good choice in general is the `nonD` function declared in tools.jl unit.
#
#  if `whitening` = true is passed, the Jeffrey mean of the matrices in 𝐂 is
#  computed (using the PosDefManifold.jl package) and the matrices in 𝐂
#  are pre-transformed using the whitening matrix of the mean.
#  Dimensionality reduction can be obtained at this stage using optional
#  arguments `eVar` and `eVarMeth` (see documentation of the AJD constructors).
#
#  if `sort`=true (default) the column vectors of the B matrix are normalized
#  to unit norm and permuted so as to sort in descending order the mean over
#  κ=1:k of the diagonal elements of B'*𝐂_κ*B.
#  Note that if `whitening` is true the output B will not have unit norm
#  columns, as it is multiplied by the whitener after being scaled and sorted
#  and before being returned.
#
#  if  `whitening` = false (default), a matrix can be provided with the `init`
#  argument in order to initialize B. In this case the actual AJD
#  will be given by init*B, where B is the output of the algorithm.
#
#  `tol` is the convergence to be attained. It defaults to 1e-6.
#
#  `maxiter` is the maximum number of iterations allowed. It defaults to 1000.
#
#  `𝜆min` is used to reguarize Hessian coefficients; all coefficients smaller
#  than 𝜆min  will be set to 𝜆min. It defauts tp 1e-4.
#
#  `lsmax` is the maximum number of steps in the line search. It defaults to 10.
#
#  if `verbose`=true, the convergence attained at each iteration and other
#  information will be printed.
#
#  if `threaded` is true, the algorithm runs in multi-threaded mode.
#  By default this is activated if the number of threads Julia is instructed
#  to use is higher than 1 and if the number of matrices in 𝐂 is higher than
#  twice the number of threads julia is instructed to use. For small matrices
#  and/or small number of matrices, setting threaded to false may result
#  in better performances in term of speed.
#
#  return a 5-tuple holding: B, its pseudo-inverse,
# 		the mean diagonal elements of B'*mean(𝐂)*B,
#       the number of iterations and the convergence attained.
#
#  Note on the implementation:
#  GRADIENT:
#  D_1,...,D_n are the matrices in 𝐃;
#  the jth column in D_i is divided by D_i[j, j],
#  the mean of these matrices is then taken and finally
#  the identity is subtracted
#
#  HESSIAN COEFFICIENTS
#  for each COLUMN vector dg_1,...,dg_n
#  (the diagonal part of the n matrices in 𝐃)
#  we form a m·m matrix stacking vertically m copies of these vectors
#  transposed and dividing each of them by the mth element, as
#   _                   _
#  |  (dg_i)'/dg_i[1]   |
#  |  (dg_i)'/dg_i[...] |
#  |  (dg_i)'/dg_i[m]   |
#  _                   _
#  finally we take the mean of all the n matrices created in this way.

# function to get the weights from argment `w`
function _qnlogLikeWeights!(w, 𝐂)
	if w isa Function w=[w(C) for C∈𝐂] end
	return w./mean(w)
end

function qnLogLike( 𝐂::Union{Vector{Hermitian}, Vector{Symmetric}};
				w           :: 	Twf   = ○,
				preWhite    :: 	Bool = false,
				sort        :: 	Bool = true,
				init        :: 	Union{Matrix, Nothing} = ○,
				tol         :: 	Real = 1e-6,
				maxiter     :: 	Int  = 1000,
				𝜆min        ::	 Real = 1e-4,
				lsmax       :: 	Int  = 10,
				verbose     :: 	Bool = false,
				threaded	:: 	Bool =
							begin
								thr = Threads.nthreads()
								length(𝐂) ≥ 2*thr && thr>1
							end,
			eVar     :: TeVaro = ○,
			eVarMeth :: Function = searchsortedfirst)

	# # # # # # # # #  internal functions
	# half sum of mean of log of Diagonals
	hsmld(𝐀) = 0.5*sum(mean(log, [𝔻(A) for A ∈ 𝐀]))
	hsmld(𝐀, 𝛎) = 0.5*sum(mean(log, [𝔻(A)*ν for (A, ν) ∈ zip(𝐀, 𝛎)]))

	# minus log abs det
	mlad(A) = -logabsdet(A)[1]

	# line search: update B and 𝐃
	@inline function _linesearch!(B, B₊, M, 𝐃, 𝐃₊, 𝒟, 𝐯, loss, lsmax)
	    for i ∈ 1:lsmax
	        M[:] = (1.0/i * 𝒟) + I
			B₊[:] = B * M
			### 𝐃₊ = [Hermitian(M'*D*M) for D in 𝐃]
			@threads for j ∈ eachindex(𝐃) 𝐃₊[j] = Hermitian(M'*𝐃[j]*M) end
			### iter>2 && (loss₊ = w===○ ? mlad(B₊)+hsmld(𝐃₊) : mlad(B₊)+hsmld(𝐃₊, 𝐯))
			loss₊ = w===○ ? mlad(B₊)+hsmld(𝐃₊) : mlad(B₊)+hsmld(𝐃₊, 𝐯)
	        loss₊ < loss && break
	    end
		B[:] = B₊
		𝐃[:] = 𝐃₊
	    return loss₊
	end
	# # # # # # # # #

	# pre-whiten or initialize or nothing
	W, 𝐃 = _preWhiteOrInit(𝐂, preWhite, Jeffrey, eVar, eVarMeth, init, :Hvector)
	𝐯 = w===○ ? ○ : _qnlogLikeWeights!(w, 𝐂) # if w is `nonD` function, apply it to the original input 𝐂

	# set variables
	iter, conv, 😋, sqrtn, loss₊ = 1, Inf, false, √size(𝐃[1], 1), Inf
	B, ⩫ = Matrix{eltype(𝐃[1])}(I, size(𝐃[1])), similar
	B₊, 𝒟, M, ∇, ℌ, 𝐃₊ = ⩫(B), ⩫(B), ⩫(B), ⩫(B), ⩫(B), ⩫(𝐃)
	𝕯 = [zeros(eltype(𝐃[1]), size(𝐃[1], 1)) for i = 1:length(𝐃)]
	loss = w===○ ? 	hsmld(𝐃) : hsmld(𝐃, 𝐯) #loss = Inf

	verbose && println("Iterating quasi-Newton LogLike algorithm...")
	while true
		for j ∈ eachindex(𝐃) 𝕯[j] = diag(𝐃[j]) end
	    #𝕯 = [diag(D) for D ∈ 𝐃]

	    # gradient
		w===○ ? ∇[:] = mean(D./𝔡 for (D, 𝔡) ∈ zip(𝐃, 𝕯)) - I :
				∇[:] = mean(v.*(D./𝔡) for (v, D, 𝔡) ∈ zip(𝐯, 𝐃, 𝕯)) - I
	    conv = norm(∇)/sqrtn # relative norm of ∇ with respect to the identity : ||∇-I||/||I||

	    verbose && println("iteration: ", iter, "; convergence: ", conv)
	    (overRun = iter > maxiter) && @warn("qnLogLike: reached the max number of iterations before convergence:", iter-1)
	    (😋 = conv <= tol) || overRun==true ? break : iter += 1

		# Hessian Coefficients
		w===○ ? ℌ[:] = mean(𝔡'./𝔡 for 𝔡 ∈ 𝕯) :
				ℌ[:] = mean(v.*(𝔡'./𝔡) for (v, 𝔡) ∈ zip(𝐯, 𝕯))

		# Quasi-Newton Direction 𝒟
	    𝒟[:] = -(∇' .* ℌ - ∇)./replace(x -> x<𝜆min ? 𝜆min : x, @. (ℌ'*ℌ) - 1.)

	    loss = _linesearch!(B, B₊, M, 𝐃, 𝐃₊, 𝒟, 𝐯, loss, lsmax) # Line Search
	end
	verbose && @info("Convergence has "*(😋 ? "" : "not ")*"been attained.\n\n")

	# scale and permute the vectors of B
	λ = sort ? _permute!(_scale!(B, mean(𝔻(D) for D ∈ 𝐃), size(𝐃[1], 1))...) :
	            diag(mean(𝔻(D) for D ∈ 𝐃))

	return preWhite ? (W.F*B, pinv(B)*W.iF, λ, iter, conv) :
	                  (B, pinv(B), λ, iter, conv)
end