text:
module ShellModule
export Shell, expandShell
using ..BaseModule
using ..BasisFunctionsModule
using ..BasisSetModule
using ..GeometryModule
import Base.display

"""
    Shell(center::Position, lqn::LQuantumNumber, exponents::[Float], coefficients::[Float64]; renorm=true)
A shell of basis functions is the collection of basis functions that share a center and radial definition (e.g.
the px, py and pz orbital that are defined together by their center, exponents and coefficients. Contrary to other codes we also consider a shell to have a unique LQuantumNumber not allowing for shared-sp tricks. Working with shells
instead of the individual basis functions allows for some performance gains in integral evaluation.
If renorm is set to false, then the coefficients are read in "as is". With renorm set to true the coefficients are expected to be given as with a basis set file and normalized accordingly. This is the default in accordance with LibInt2Shell, where this is necessary due to performance reasons. (As such the normalization scaling is taken from the renorm() function in libint2/shell.h.
"""
type Shell
  center::Position
  lqn::LQuantumNumber
  exponents::Vector{Float64}
  coefficients::Vector{Float64}
  function Shell(center,lqn,exponents,coefficients;renorm=true)
    L = lqn.exponent
    if renorm
      nprim = length(exponents)
      # renorm primitive-by-primitive
      for p in 1:nprim
	coefficients[p] *= sqrt( (2^L*(2*exponents[p])^(L+3/2)) / (π^(3/2)*doublefactorial(2*L-1)) )
      end
      # normalize total shell
      norm = 0.
      for p in 1:nprim, q in 1:nprim
        norm += (doublefactorial(2*L-1)*π^(3/2)*coefficients[p]*coefficients[q])/
        	((2^L)*(exponents[p]+exponents[q])^(L+3/2))
      end
      for p in 1:nprim
        coefficients[p] /= sqrt(norm)
      end
    end
    new(center,lqn,exponents,coefficients)
  end
end

function display(sh::Shell)
  println("$(sh.lqn.symbol)-Shell   @    $(sh.center)")
  print("  Exponents:    ")
  for exp in sh.exponents
    @printf("  %f",exp)
  end
  println("")
  print("  Coefficients: ")
  for coeff in sh.coefficients
    @printf("  %f",coeff)
  end
  println("")
end

"""
Compute a list of ContractedGaussianBasisFunctions which are described within the input shell.
"""
function expandShell(sh::Shell)
  result = ContractedGaussianBasisFunction[]
  for mqn in MQuantumNumbers(sh.lqn)
    cgbf = ContractedGaussianBasisFunction(sh.coefficients,[])
    for exp in sh.exponents
      push!(cgbf.primitiveBFs,PrimitiveGaussianBasisFunction(sh.center,exp,mqn))
    end
    push!(result,cgbf)
  end
  return result
end

"""
compute the number of basisfunctions (cartesian) described by the shell
"""
function nbf(shell::Shell)
  div((shell.lqn.exponent+1)^2+(shell.lqn.exponent+1),2)
end
end # module


text:
using PDESolver
using NavierStokesMod
using Base.Test
using ODLCommonTools
import ODLCommonTools.sview

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))  # get TestSystem

#------------------------------------------------------------------------------
# define tests and tags

using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const NSTests = TestList()

# include test files
include("test_viscous.jl")

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Euler tests -----" begin

  runTestSystem(NSTests, solvePDE, ARGS)
end

ext: jl


text:
module NNlibCUDA

using NNlib
using CUDA
using Random, Statistics

const IntOrIntTuple = Union{Integer, NTuple{N,<:Integer} where N}

include("upsample.jl")
include("activations.jl")
include("batchedmul.jl")
include("scatter.jl")
include("gather.jl")
include("cudnn/cudnn.jl")
include("cudnn/conv.jl")
include("cudnn/pooling.jl")
include("cudnn/softmax.jl")
include("cudnn/activations.jl")
include("cudnn/batchnorm.jl")

end # module


text:
struct TimeSeriesCacheKey
    component_uuid::Base.UUID
    time_series_type::Type{<:IS.TimeSeriesData}
    name::String
    multiplier_id::Int
end

# TODO: Marge all structs (ModelInternal, ModelStoreParams and SimulationInfo) to a single Internal Struct

mutable struct SimulationInfo
    number::Int
    requires_rebuild::Bool
    sequence_uuid::Base.UUID
end

mutable struct ModelInternal
    container::OptimizationContainer
    ic_model_container::Union{Nothing, OptimizationContainer}
    status::BuildStatus
    run_status::RunStatus
    base_conversion::Bool
    executions::Int
    execution_count::Int
    output_dir::Union{Nothing, String}
    simulation_info::Union{Nothing, SimulationInfo}
    time_series_cache::Dict{TimeSeriesCacheKey, <:IS.TimeSeriesCache}
    recorders::Vector{Symbol}
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    store_parameters::Union{Nothing, ModelStoreParams}
    ext::Dict{String, Any}
end

function ModelInternal(
    container::OptimizationContainer;
    ext=Dict{String, Any}(),
    recorders=[],
)
    return ModelInternal(
        container,
        nothing,
        BuildStatus.EMPTY,
        RunStatus.READY,
        true,
        1, #Default executions is 1. The model will be run at least once
        0,
        nothing,
        nothing,
        Dict{TimeSeriesCacheKey, IS.TimeSeriesCache}(),
        recorders,
        Logging.Warn,
        Logging.Info,
        nothing,
        ext,
    )
end

function add_recorder!(internal::ModelInternal, recorder::Symbol)
    push!(internal.recorders, recorder)
    return
end

get_recorders(internal::ModelInternal) = internal.recorders

function configure_logging(internal::ModelInternal, file_mode)
    return IS.configure_logging(
        console=true,
        console_stream=stderr,
        console_level=internal.console_level,
        file=true,
        filename=joinpath(internal.output_dir, PROBLEM_LOG_FILENAME),
        file_level=internal.file_level,
        file_mode=file_mode,
        tracker=nothing,
        set_global=false,
    )
end


text:
export BrickletAnalogOutV2Identity
struct BrickletAnalogOutV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAnalogOutV2
"""
Generates configurable DC voltage between 0V and 12V
"""
mutable struct BrickletAnalogOutV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAnalogOutV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_analog_out_v2")
        deviceInternal = package.BrickletAnalogOutV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_output_voltage
"""
    $(SIGNATURES)

Sets the voltage.
"""
function set_output_voltage(device::BrickletAnalogOutV2, voltage)
    device.deviceInternal.set_output_voltage(voltage)
end

export get_output_voltage
"""
    $(SIGNATURES)

Returns the voltage as set by :func:`Set Output Voltage`.
"""
function get_output_voltage(device::BrickletAnalogOutV2)
    return device.deviceInternal.get_output_voltage()
end

export get_input_voltage
"""
    $(SIGNATURES)

Returns the input voltage.
"""
function get_input_voltage(device::BrickletAnalogOutV2)
    return device.deviceInternal.get_input_voltage()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAnalogOutV2)
    return device.deviceInternal.get_identity()
end


text:
# This file was generated, do not modify it. # hide
@load BayesianLDA pkg=MultivariateStats

clf = machine(BayesianLDA(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
import MutableArithmetics
const MA = MutableArithmetics

function all_tests(a, b, c, d, e, f, g)
    a_copy = deepcopy(a)
    b_copy = deepcopy(b)
    c_copy = deepcopy(c)
    d_copy = deepcopy(d)
    e_copy = deepcopy(e)
    f_copy = deepcopy(f)
    g_copy = deepcopy(g)
    @testset "Scalar" begin
        MA.Test.scalar_test(a)
        MA.Test.scalar_test(b)
        MA.Test.scalar_test(c)
        MA.Test.scalar_test(d)
        MA.Test.scalar_test(e)
        MA.Test.scalar_test(f)
        MA.Test.scalar_test(g)
    end
    @test isequal(a, a_copy)
    @test isequal(b, b_copy)
    @test isequal(c, c_copy)
    @test isequal(d, d_copy)
    @test isequal(e, e_copy)
    @testset "Quadratic" begin
        MA.Test.quadratic_test(a, b, c, d)
        MA.Test.quadratic_test(b, c, d, e)
        MA.Test.quadratic_test(c, d, e, f)
        MA.Test.quadratic_test(d, e, f, g)
        MA.Test.quadratic_test(e, f, g, a)
        MA.Test.quadratic_test(f, g, a, b)
        MA.Test.quadratic_test(g, a, b, c)
    end
    @test isequal(a, a_copy)
    @test isequal(b, b_copy)
    @test isequal(c, c_copy)
    @test isequal(d, d_copy)
    @test isequal(e, e_copy)
    @testset "Sparse" begin
        MA.Test.sparse_test(a, b, [a b c; b c a; a b a])
    end
    @test isequal(a, a_copy)
    @test isequal(b, b_copy)
    @test isequal(c, c_copy)
    @test isequal(d, d_copy)
    @test isequal(e, e_copy)
    @testset "Vector" begin
        MA.Test.array_test([a, b, c])
        MA.Test.array_test([b, c, d])
        MA.Test.array_test([c, d, e])
        MA.Test.array_test([d, e, a])
        MA.Test.array_test([e, a, b])
    end
    @test isequal(a, a_copy)
    @test isequal(b, b_copy)
    @test isequal(c, c_copy)
    @test isequal(d, d_copy)
    @test isequal(e, e_copy)
    @testset "Matrix" begin
        MA.Test.array_test([a b; c d])
        MA.Test.array_test([c e; e d])
        MA.Test.array_test([a b c; b c a; a b a])
        MA.Test.array_test([d b c; d c e; e b a])
    end
    @test isequal(a, a_copy)
    @test isequal(b, b_copy)
    @test isequal(c, c_copy)
    @test isequal(d, d_copy)
    @test isequal(e, e_copy)
end

Mod.@polyvar x y

@testset "MutableArithmetics with variables" begin
    # Creating 7 different variables here gives a long compile time for TypedPolynomials
    all_tests(x, y, x, y, x, y, x)
end

@testset "MutableArithmetics with monomials" begin
    a = x^2
    b = y^2
    c = x
    d = x^2 * y
    e = x^3
    f = x^5
    g = x^4
    all_tests(a, b, c, d, e, f, g)
end

@testset "MutableArithmetics with terms in $T" for T in [Int, BigInt]
    if MA.mutability(T) isa MA.IsMutable && MA.mutability(typeof(x * y)) isa MA.IsMutable
        @testset "Int" begin
            MA.Test.int_test(termtype(x, T), exclude = ["int_add", "int_add_mul", "int_zero"])
        end
    end
    a = T(2) * x^2
    b = T(4) * y^2
    c = T(3) * x
    d = T(5) * x^2 * y
    e = T(5) * x^3
    f = T(1) * x^5
    g = T(2) * x^4
    all_tests(a, b, c, d, e, f, g)
end

@testset "MutableArithmetics with polynomials in $T" for T in [Int, BigInt]
    if MA.mutability(T) isa MA.IsMutable
        @testset "Int" begin
            MA.Test.int_test(polynomialtype(x, T))
        end
    end
    a = T(2) * x^2 + T(3) * x * y + T(4) * y
    b = T(4) * y^2 - T(1) * x * y + T(4) * x
    c = T(1) * x^2 + T(3) * x - T(4)
    d = T(5) * x^2 * y - T(3) * x + T(4) * y
    e = T(5) * x^3 - T(3) * x + T(4)
    f = x^5 + x^4 + x^2 + x + T(1)
    g = T(2)x^4 + x^3
    all_tests(a, b, c, d, e, f, g)
end


text:
#include("orbit/orbit.jl")
include("external_model/external_model.jl")
#include("dynamic_model/dynamic_model.jl")
#include("dynamics/dynamics.jl")
#include("satellite/satellite.jl")
#include("plot/plot_plots.jl")
include("plot/plot_makie.jl")
include("coordinates.jl")

using SatelliteToolbox
using Dates
using LinearAlgebra

using HDF5

function main()
	#=
	設定パラメータ
	=#
	DataNum =100 #シミュレータ反復回数
	dt = 5 ##シミュレータの計算間隔 [s]
	start_time = DateTime(2019, 12, 19, 3, 27, 10)	#シミュレート開始時刻
	TLEFileName = "./orbit/ISS_TLE.txt"
	# 進行方向とSCOFx軸とのずれ（内積）
	dotvs = zeros(DataNum)

	# 軌道計算については先に行い、全時間分を配列に保存する
	# JD_log, x_ecef_log, v_ecef_log, x_geod_log = orbit_cal(DataNum,dt,start_time,TLEFileName)
	JD_log = h5read("./orbit/orbit_data.h5", "JD_log")
	x_ecef_log = h5read("./orbit/orbit_data.h5", "x_ecef_log")
	v_ecef_log = h5read("./orbit/orbit_data.h5", "v_ecef_log")
	x_geod_log = h5read("./orbit/orbit_data.h5", "x_geod_log")

	# SCOF上での衛星の進行方向
	direct_on_SCOFs = zeros(DataNum,3)

	# 衛星外環境モデル用変数
	mag_vecs = zeros(DataNum,3)
	sun_vecs = zeros(DataNum,3)
	atoms_denses = zeros(DataNum)

  	# 衛星内環境モデル用変数
  	torqe = zeros(DataNum,3)

	# plot_makie
	# mp_earth()

	for i=1:DataNum

		# 時刻表示
		current_time = start_time+Second(dt) * i
		println("")
		println("Current DateTime:",current_time)
		println("              JD:",JD_log[i])
		println("      x_ecef_log:",x_ecef_log[i,:])
		println("      v_ecef_log:",v_ecef_log[i,:])
		println("      x_geod_log:",x_geod_log[i,:])
		dotv = dot(x_ecef_log[i,:] , v_ecef_log[i,:])
		dotvs[i] = dotv
		println("            dotv:",dotv)

		# 衛星外環境モデル計算
		mag_vec, sun_vec, atoms_dens = external_model(current_time,x_ecef_log[i,:],x_geod_log[i,:])
		println("         mag_vec:",mag_vec)
		mag_vecs[i,:] = mag_vec
		println("         sun_vec:",sun_vec)
		sun_vecs[i,:] = sun_vec
		println("      atoms_dens:",atoms_dens)
		atoms_denses[i] = atoms_dens

		direct_on_SCOFs[i,:] = ecef_to_DCM(x_ecef_log[i,:],v_ecef_log[i,:],true) * v_ecef_log[i,:]

		#mp_sat(my_scene,x_ecef_log[i,:])
		sleep(1/24)

	end

	#plot_2scalar(JD_log,dotvs,"dotvs")
	#plot_vec(direct_on_SCOFs,"direct_on_SCOFs")
	#
	#plot_2scalar(JD_log,atoms_denses,"atoms_dens")
	#plot_2scalar(JD_log,x_geod_log[:,3],"height")
	#plot_vec(mag_vecs,"mag_vec")
	#plot_vec(sun_vecs,"sun_vec")
	#plot_vec(x_ecef_log,"x_ecef_log")
	#plot_vec(v_ecef_log,"v_ecef_log")
	#plot_vec(x_geod_log,"x_geod_log")
	#plot_2scalar(x_geod_log[:,2],x_geod_log[:,1],"x_geod_log_2d")

end

my_scene = Scene()
mp_basic_coordinate(my_scene)

main()

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
"""
	solve_particle_scattering(k0, kin, P, sp::ScatteringProblem, u::Einc; get_inner = true, verbose = true) -> beta, inner

Solve the scattering problem `sp` with outer wavenumber `k0`, inner wavenumber
`kin`, `2P+1` cylindrical harmonics per inclusion and incident TM field
`u`. Solves multiple-scattering equation directly.
Returns the cylindrical harmonics basis `beta` along with potential
densities (in case of arbitrary inclusion) or inner cylindrical coefficients (in
case of circular). By default, incident wave propagates left->right.

Inner coefficients are only calculated if `get_inner` is true, and timing is
printed if `verbose` is true.
"""
function solve_particle_scattering(k0, kin, P, sp::ScatteringProblem, u::Einc;
								get_inner = true, verbose = true)
	#TODO: incorporate matrix-lexx rotation rotateMultipole!
	shapes = sp.shapes;	ids = sp.ids; centers = sp.centers; φs = sp.φs
	Ns = size(sp)
	#first solve for single scatterer densities
	dt1 = @elapsed begin
		scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids)
	end
	dt2 = @elapsed begin
		T = Array{Complex{Float64}}(undef, Ns*(2*P+1), Ns*(2*P+1))
		a = u2α(k0, u, centers, P)
		Btemp = Array{Complex{Float64}}(undef, 2*P+1, 2*P+1)
		for ic1 = 1:Ns
			rng1 = (ic1-1)*(2*P+1) .+ (1:2*P+1)
			if φs[ic1] == 0.0 || typeof(shapes[ids[ic1]]) == CircleParams
				RotScatMat = scatteringMatrices[ids[ic1]]
			else
				Rot = sparse(Diagonal(Complex{Float64}[exp(-1.0im*φs[ic1]*l) for l=-P:P])) #rotation matrix
				RotScatMat = Rot*(scatteringMatrices[ids[ic1]]*conj(Rot))
			end
			for ic2 = 1:Ns
				rng2 = (ic2-1)*(2*P+1) .+ (1:2*P+1)
				if ic1 == ic2
					T[rng1,rng2] = Matrix{Complex{Float64}}(I, 2*P+1, 2*P+1)
				else
					M2Lmatrix!(Btemp, k0, P, centers[ic1,:] - centers[ic2,:])
					T[rng1,rng2] = -RotScatMat*Btemp
				end
			end
			a[rng1] = RotScatMat*a[rng1]
		end
	end
	dt3 = @elapsed begin
		beta = T\a
	end
	#recover full incoming expansion - in sigma_mu terms for parametrized shape,
	#in multipole expansion for circle
	if get_inner
		#find LU factorization once for each shape
		dt4 = @elapsed begin
	        scatteringLU = [lu(scatteringMatrices[iid]) for iid = 1:length(shapes)]
			inner = Array{Vector{Complex{Float64}}}(undef, Ns)
			α_c = Array{Complex{Float64}}(undef, 2*P+1)
			for ic = 1:Ns
				rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
				if typeof(shapes[ids[ic]]) == ShapeParams
					if φs[ic] == 0.0
						α_c[:] = scatteringLU[ids[ic]]\beta[rng]
					else
						Rot = sparse(Diagonal(Complex{Float64}[exp(-1.0im*φs[ic]*l) for l=-P:P])) #rotation matrix
						α_c[:] = scatteringLU[ids[ic]]\(conj(Rot)*beta[rng])
					end
					inner[ic] = innerExpansions[ids[ic]]*α_c
				else
					inner[ic] = innerExpansions[ids[ic]]*beta[rng]
				end
			end
		end
	end
	verbose && begin
		println("Direct solution timing:")
		println("Scattering matrix solution: $dt1 s")
		println("Matrix construction: $dt2 s")
		println("Matrix solution: $dt3 s")
        get_inner && println("Retrieving inner coefficients: $dt4 s")
    end
	get_inner ? (return beta, inner) : (return beta)
end

function M2Lmatrix!(T, k, P, d)
	#builds off-diagonal translation matrix from particle 2 to particle 1 (d=1-2)
	#only calculates 2*(2*P+1) values
	kd = k*sqrt(sum(abs2,d))
	td = atan(d[2],d[1])
	bess = besselh.(0:2*P,1,kd)
	for ix = 1:2*P #lower diagonals
		rng = ix+1:1+(2*P+1):(2*P+1)^2-(2*P+1)*ix
		T[rng] .= exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1]
	end
	for ix = 0:2*P #central and upper diagonals
		rng = ix*(2*P+1)+1:1+(2*P+1):(2*P+1)^2-ix
		T[rng] .= exp(1im*td*ix)*bess[ix+1]
	end
end

function scattered_field_multipole(k0, beta, centers::Array{Float64,2}, points::Array{Float64,2}; recurrence = false)
	Ns = size(centers,1)
	P = div(div(length(beta),Ns)-1,2)
	len = size(points,1)
	Esc = zeros(Complex{Float64}, len)

	if recurrence
		scattered_field_multipole_recurrence!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len)
	else
		scattered_field_multipole!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len)
	end
	return Esc
end

function scattered_field_multipole!(Esc::Array{Complex{Float64},1}, k0, beta, P, centers::Array{Float64,2}, ind_centers, points::Array{Float64,2}, ind_points)
	for ic in ind_centers
		ind = (ic-1)*(2*P+1) + P + 1
		for ip in ind_points
			points_moved1 = points[ip,1] - centers[ic,1]
			points_moved2 = points[ip,2] - centers[ic,2]

			rs_moved = hypot(points_moved1, points_moved2)
			ts_moved = atan(points_moved2, points_moved1)
			if rs_moved == 0
				error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip")
			end
			Esc[ip] += beta[ind]*besselh(0, 1, k0*rs_moved)
			for p = 1:P
				Esc[ip] += besselh(p, 1, k0*rs_moved)*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved))
			end
		end
	end
end

function scattered_field_multipole_recurrence!(Esc::Array{Complex{Float64},1}, k0, beta, P, centers::Array{Float64,2}, ind_centers, points::Array{Float64,2}, ind_points)
	for ic in ind_centers
		ind = (ic-1)*(2*P+1) + P + 1
		for ip in ind_points
			points_moved1 = points[ip,1] - centers[ic,1]
			points_moved2 = points[ip,2] - centers[ic,2]

			rs_moved = hypot(points_moved1, points_moved2)
			ts_moved = atan(points_moved2, points_moved1)
			if rs_moved == 0
				error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip")
			end
			Hₚ₋₂ = besselh(-1, 1, k0*rs_moved)
			Hₚ₋₁ = besselh(0, 1, k0*rs_moved)
			Esc[ip] += beta[ind]*Hₚ₋₁
			for p = 1:P
				Hₚ = (2*(p-1)/(k0*rs_moved))*Hₚ₋₁ - Hₚ₋₂
				Esc[ip] += Hₚ*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved))
				Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ
			end
		end
	end
end

function circleScatteringMatrix(kout, kin, R, P; gamma = false)
    #non-vectorized, reuses bessel
    S = Array{Complex{Float64}}(undef, 2*P+1)
    gamma && (G = Array{Complex{Float64}}(undef, 2*P+1))

    pre_J0 = besselj(-P-1,kout*R)
    pre_J1 = besselj(-P-1,kin*R)
    pre_H = besselh(-P-1,kout*R)
    for p = -P:0
        J0 = besselj(p,kout*R)
        J1 = besselj(p,kin*R)
        H = besselh(p,kout*R)

        dJ0 = kout*(pre_J0 - (p/kout/R)*J0)
        dJ1 = kin*(pre_J1 - (p/kin/R)*J1)
        dH = kout*(pre_H - (p/kout/R)*H)

        S[p+P+1] = -(dJ0*J1 - J0*dJ1)/(dH*J1 - H*dJ1)
        p != 0 && (S[P+1-p] = S[p+P+1])

        if gamma #also returns gamma_n as a function of beta_n
            G[p+P+1] = (J0*dH - dJ0*H)/(J0*dJ1 - dJ0*J1)
            p != 0 && (G[P+1-p] = G[p+P+1])
        end
        pre_J0 = J0
        pre_J1 = J1
        pre_H = H
    end
    ScatMat = sparse(Diagonal(S))
    if gamma
        InnerMat = sparse(Diagonal(G))
        return ScatMat,InnerMat
    else
        return ScatMat
    end
end

function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Float64,2})
    #non-vectorized, reuses bessel
    P = div(length(gamma)-1,2)

	len = size(points,1)
	points_moved = similar(points)
    points_moved[:,1] = points[:,1] .- center[1]
	points_moved[:,2] = points[:,2] .- center[2]

	rs_moved = sqrt.(sum(abs2,points_moved, dims=2))
	ts_moved = atan.(points_moved[:,2], points_moved[:,1])

	bess = [besselj(p,kin*rs_moved[ii]) for ii=1:len, p=0:P]
	E = gamma[P + 1]*bess[:,1]
	for p = 1:P
		E += bess[:,p+1].*(gamma[p + P + 1]*exp.(1.0im*p*ts_moved) +
							(-1)^p*gamma[-p + P + 1]*exp.(-1.0im*p*ts_moved))
	end
	return E
end

function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Float64,1})
    #non-vectorized, reuses bessel
    P = div(length(gamma)-1,2)

    points_moved = points - center
	rs_moved = sqrt(sum(abs2,points_moved))
	ts_moved = atan(points_moved[2], points_moved[1])

	bess = [besselj(p,kin*rs_moved) for p=0:P]
	E = gamma[P + 1]*bess[1]
	for p = 1:P
		E += bess[p+1]*(gamma[p + P + 1]*exp(1.0im*p*ts_moved) + (-1)^p*gamma[-p + P + 1]*exp(-1.0im*p*ts_moved))
	end
	return E
end

function particleExpansion(k0, kin, shapes, P, ids)
    scatteringMatrices = Array{Any}(undef, 0)
	innerExpansions = Array{Any}(undef, 0)
	for i = 1:length(shapes)
        #no use in computing matrices if shape doesn't actually show up! push garbage to maintain order
		#TODO: don't even try to compute LU for these
		if all(ids .!= i)
			push!(scatteringMatrices, [NaN NaN])
			push!(innerExpansions, [NaN NaN])
			continue
		end
		if typeof(shapes[i]) == ShapeParams
            AB = shapeMultipoleExpansion(k0, shapes[i].t, shapes[i].ft, shapes[i].dft, P)
    		sigma_mu_mult = get_potential(k0, kin, P, shapes[i].t, shapes[i].ft, shapes[i].dft)
    		ScatMat = AB*sigma_mu_mult
            push!(scatteringMatrices, ScatMat)
            push!(innerExpansions, sigma_mu_mult)
		else
			ScatMat, inner = circleScatteringMatrix(k0, kin, shapes[i].R, P, gamma = true)
			push!(scatteringMatrices, ScatMat)
			push!(innerExpansions, inner)
		end
	end
    return scatteringMatrices,innerExpansions
end

function rotateMultipole!(arr,phi,P)
    #applies e^{-i \phi l}, p=-P...P
    #if slow, change from view to direct indexing (should be negligible despite higher memory use)
    for p = -P:P
        arr[P + 1 + p] *= exp(-1.0im*phi*p)
    end
end

function rotateMultipole!(dest_arr,source_arr,phi,P)
    #applies e^{-i \phi l}, p=-P...P
    #if slow, change from view to direct indexing (should be negligible despite higher memory use)
    for p = -P:P
        dest_arr[P + 1 + p] = source_arr[P + 1 + p]*exp(-1.0im*phi*p)
    end
end


text:
# # chapter1
#

using Joe


text:
using ConfigEnv; dotenv()

TOKEN = ENV["TOKEN"]
MASTER_ID = parse(Int, ENV["MASTER_ID"])
API_URL = ENV["API_URL"]

COMMANDS = ["lucky", "move", "shuffle", "top_playlist", "fork", "merge", "group", "news"]

START_MESSAGE = """
Hi! I am a spotify-utils bot made by iwsylit living at https://github.com/iwsylit/spotify-utils-bot.\n
I can:
    /lucky - create a playlist containing all your songs in random order
    /move - move n songs from the bottom to the top of a playlist
    /shuffle - shuffle songs in a playlist
    /top_playlist - create a playlist containing your top songs
    /fork - steal someone's playlist
    /merge - merge two playlists
    /group - group playlists
    /news - create a playlist containing all your songs ordered by date added
"""

STRANGER_MESSAGE = """
I won't respond you in any way, you are not my master.
If you want to become one, go to https://github.com/iwsylit/spotify-utils-bot and try to tame me.
One day my master will make me available for you, too, but I am not sure.
"""

max_stars_repo_name: iwsylit/spotify-utils-bot

max_issues_repo_name: iwsylit/spotify-utils-bot

max_forks_repo_name: iwsylit/spotify-utils-bot


text:
InputFolder = './Images/NucleiNoisy/';
OutputFolder = './Results/Images/NucleiNoisy/';

@iA = '*.tif';

@fxg_gDenoiseBM3 [iA] > [D, An];
params.AddNoiseVar = 0;
/endf

/show iA >;
/show D >;
/keep D > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
# Testing the Small Demo
# Copyright (c) 2015 Harold Soh

# based on MATLAB code by Mark Schmidt
# For more details, see: http://www.cs.ubc.ca/~schmidtm/Software/UGM/small.html

include("../src/UGM.jl")

using UGM
using Base.Test

include("testUtils.jl")

# ======================================
# Simple Independent Node Test
# ======================================
srand(0)
nodepot = [ 1 3;
            9 1;
            1 3;
            9 1 ]

nsamples = 100;
samples = zeros(4, nsamples)
for i = 1:4
  for s = 1:nsamples
    samples[i,s] = sampleDiscrete(nodepot[i,:]/sum(nodepot[i,:]));
  end
end

sum_samples = sum(samples,2)
println("Sum Samples: ", sum_samples')
gold_result = [172; 107; 168; 110] # at seed 0
@test all(sum_samples .== gold_result)
print("Independent Nodes Test: ")
printPassed()

# ======================================
# Set up Dependent Nodes
# ======================================
nnodes = 4
nstates = 2

edgelist =  [ 1 2;
              2 1;
              2 3;
              3 2;
              3 4;
              4 3]

# Make Edge Structure
es = EdgeStruct(edgelist = edgelist, nstates = nstates, maxiter = 25)

# edgestruct tests
@test es.edgeends == [1 2; 2 3; 3 4]
@test es.edgedict[1] == Set([2])
@test es.edgedict[2] == Set([1,3])
@test es.edgedict[3] == Set([2,4])
@test es.edgedict[4] == Set([3])
@test es.nnodes == 4
@test es.nedges == 3
@test all(es.nstates .== [2.0; 2.0; 2.0; 2.0])
@test es.maxiter == 25
print("Edge Structure Test: ")
printPassed()

# Node potentials
nodepot = [1 3;
          9 1;
          1 3;
          9 1]

# Edge potentials
maxstate = maximum(es.nstates);
edgepot = zeros( int(maxstate), int(maxstate), int(es.nedges));
for e = 1:es.nedges
   edgepot[:,:,e] = [2 1 ; 1 2];
end

# ======================================
# Decoding
# ======================================
optdecoding = decodeExact(nodepot,edgepot,es)
println("Optimal Decoding: ", optdecoding')
@test all(optdecoding .== [2.0; 1.0; 1.0; 1.0])
print("Exact Decoding Test: ")
printPassed()

# ======================================
# Inference
# ======================================
nodebel,edgebel,logZ,Z = inferExact(nodepot,edgepot,es)

@test nodebel == [0.3596306068601583 0.6403693931398416
                   0.8430079155672823 0.15699208443271767
                   0.4862796833773087 0.5137203166226912
                   0.8810026385224274 0.11899736147757256]

@test edgebel[:,:,1] == [0.33720316622691293 0.022427440633245383
                         0.5058047493403693 0.1345646437994723]

@test edgebel[:,:,2] ==  [0.45118733509234826 0.391820580474934
                          0.03509234828496042 0.12189973614775726]

@test edgebel[:,:,3] == [0.46068601583113455 0.02559366754617414
                         0.4203166226912929 0.09340369393139841]

@test_approx_eq_eps(logZ, 8.240121298, 1e-8)
print("Exact Inference Test: ")
printPassed()

# ======================================
# Sampling Test
# ======================================
es.maxiter = 100
samples = sampleExact(nodepot,edgepot,es)
@test all(samples[1,:] .== [2.0 1.0 2.0 2.0])
@test all(samples[4,:] .== [2.0 2.0 2.0 1.0])
@test all(samples[100,:] .== [2.0 1.0 2.0 1.0])
print("Exact Sampling Test: ")
printPassed()

# println(samples)
# using Gadfly
# spy(samples)


text:
module WebBasedWorkspace

using HTTP
using Pkg
using TOML
using URIs

include("utils.jl")
include("package_specs.jl")
include("github.jl")

end

max_stars_repo_name: MarkNahabedian/WebBasedWorkspace

max_issues_repo_name: MarkNahabedian/WebBasedWorkspace

max_forks_repo_name: MarkNahabedian/WebBasedWorkspace


text:
module FindNeighbours

abstract type PrevNext end
struct Prev <: PrevNext end
struct Next <: PrevNext end

ind(x, ::Prev) = x.i_prev
ind(x, ::Next) = x.i_next
ind(x::AbstractVector, prevnext::Vector) = Tuple(ind.(x, prevnext))

wgt(x, ::Prev) = x.λ
wgt(x, ::Next) = 1 - x.λ
wgt(x::AbstractVector, prevnext::Vector) = prod(wgt.(x, prevnext))

# """
#   Given a ordered vector [x_1 < x_2 < ... < x_N] and y ∈ [x_1, x_N] find i and λ such that y = λ x_i + (1-λ) x_{i+1}
# """
function findneighbours_line(vector, point)
  i_next = searchsortedfirst(vector, point)

  if i_next == 1
    i_prev = 1
    λ = 1.0
  else
    i_prev = i_next - 1
    prev = vector[i_prev]
    next = vector[i_next]

    Δ = next - prev
    λ = (next - point) / Δ
  end

  (i_prev=i_prev, i_next=i_next, λ=λ)
end

function findneighbours(grids, point)

  ngb_vec = [findneighbours_line(grids[i], point[i]) for i in 1:length(point)]

  prev_next_combis = collect(Iterators.product(repeat([[Prev(), Next()]], length(grids))...))

  ngb = map(prev_next_combis) do prev_next
    (i=ind(ngb_vec, collect(prev_next)), weight=wgt(ngb_vec, collect(prev_next)))
  end
  vec(ngb)
end

export findneighbours

using Test, Parameters

@testset "findneighbours_line" begin
  a_grid = LinRange(1, 10, 100)
  p = 7.23456
  @unpack i_prev, i_next, λ = findneighbours_line(a_grid, p)

  a_prev, a_next = a_grid[[i_prev; i_next]]

  @test λ * a_prev + (1-λ) * a_next ≈ p
  @test a_prev < p < a_next
end

@testset "findneighbours" begin
  grids = (a = LinRange(0.1,1,10), b=LinRange(1.1,2,10))
  point = (a=0.15, b=1.15)

  ngb = findneighbours(grids, point)

  a_check = sum(grids.a[ngb[i].i[1]] .* ngb[i].weight for i in 1:length(ngb))
  b_check = sum(grids.b[ngb[i].i[2]] .* ngb[i].weight for i in 1:length(ngb))

  @test a_check ≈ point.a
  @test b_check ≈ point.b

  point = (a=0.11111, b=1.56788)
  ngb = findneighbours(grids, point)

  a_check = sum(grids.a[ngb[i].i[1]] .* ngb[i].weight for i in 1:length(ngb))
  b_check = sum(grids.b[ngb[i].i[2]] .* ngb[i].weight for i in 1:length(ngb))
  @test a_check ≈ point.a
  @test b_check ≈ point.b

  grids = (a = LinRange(0.1,1,32), b=LinRange(1.1,2,71), c=LinRange(0,10,100))
  point = (a=0.15, b=1.15, c=π)

  ngb = findneighbours(grids, point)

  a_check = sum(grids.a[ngb[i].i[1]] .* ngb[i].weight for i in 1:length(ngb))
  b_check = sum(grids.b[ngb[i].i[2]] .* ngb[i].weight for i in 1:length(ngb))
  c_check = sum(grids.c[ngb[i].i[3]] .* ngb[i].weight for i in 1:length(ngb))

  @test a_check ≈ point.a
  @test b_check ≈ point.b
  @test c_check ≈ point.c
end

end # module

using .FindNeighbours

## Build up transition matrix
function controlled_markov_chain(endo, exo, policy)
  lin_ind = LinearIndices((length(endo), length(exo)))

  len_exo = length(exo)
  len = length(endo) * len_exo

  n_endo_neighbours = 2 ^ dimension(endo) # 2 neighbours on line, 4 on plane, 8 in ℝ³

  len_sparse = len * len_exo * n_endo_neighbours

  I = zeros(Int, len_sparse)
  J = zeros(Int, len_sparse)
  V = zeros(len_sparse)

  controlled_markov_chain!(I, J, V, lin_ind, endo, exo, policy)

  (I=I, J=J, V=V, n=len)
end

function controlled_markov_chain!(I, J, V, lin_ind, endo, exo, policy)

  len_exo = length(exo)
  len_endo = length(endo)

  j = 0

  for i_exo in 1:len_exo
    # extract jump probabilities
    π = exo.mc.p

    for i_endo in 1:len_endo
      p = policy[i_endo, i_exo]
      ijump_mass_vec = findneighbours(endo.grids, p)
      for ijump_mass in ijump_mass_vec
        @unpack i, weight = ijump_mass
        i_jump = linear_indices(endo)[i...]
        mass = weight
        for i_exo_next in 1:len_exo
          j += 1
          I[j] = lin_ind[i_endo, i_exo]
          J[j] = lin_ind[i_jump, i_exo_next]
          V[j] = mass * π[i_exo, i_exo_next]
        end
      end
    end
  end
end

## Stationary distribution

function stationary_distribution(endo, exo, policy)

  @unpack I, J, V, n = controlled_markov_chain(endo, exo, policy)

  transition = sparse(I, J, V, n, n)

  eigval, eigvec, _ = eigs(transpose(transition))
  is = findall(abs.(eigval) .≈ 1)
  @assert length(is) == 1
  @assert isreal(eigvec[:,is[1]])
  dist = real(eigvec[:,is[1]])

  dist ./= sum(dist)

  reshape(dist, size(policy))
end

# let value = value, policy=policy, z_mc = z_MC
#   lin_ind = LinearIndices(size(value))
#
#   ngp_exo = length(z_mc.state_values)
#   n = length(policy)
#   len = n * ngp_exo * 2
#
#   I = zeros(Int, len)
#   J = zeros(Int, len)
#   V = zeros(len)
#
#   @btime Aiyagari.controlled_markov_chain!($I, $J, $V, $lin_ind, $z_mc, $a_grid, $policy)
# end


text:
using NonlinearEigenproblems: nep_gallery, compute_Mder, contour_block_SS, contour_beyn
using FEASTSolver

using SuiteSparse

SuiteSparse.UMFPACK.umf_ctrl[8] = 0;

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-5]
        in_eig_conv = in_eig[in_res .<= 1e-5]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

nep = nep_gallery("nlevp_native_gun")
T(x) = compute_Mder(nep, x)

n = size(nep, 1)
C = complex(62500.0, 0.0)
R = 50000

e, v, res = nlfeast!(T, rand(ComplexF64,n,32), 2^3, 10, c=C, r=R, debug=true, ϵ=10e-10, store=true)
# e, v, res = nlfeast!(T, rand(ComplexF64,n,50), 2^4, 10, c=C, r=2*R, debug=true, ϵ=10e-16, store=false)
#
# e, v, res = nlfeast_moments!(T, rand(ComplexF64,n,4), 2^5, 0, moments=8, c=C, r=R, debug=true, ϵ=10e-16, store=false)
info(e,v,res,C,R)

# e, v = contour_beyn(nep, radius=R, σ=C, N=2^8, k=80)
# e, v = contour_block_SS(nep, radius=2*R, σ=C, N=2^8, k=80, K=4)
# print(e)


text:
"""
Given a form `value::type`, return `value` wrapped in a
`Nullable`. Errors of type `catch_errors` are caught and returned
as a `Nullable{type}()`, other errors are rethrown.

Example: @catch_as_Nullable Date(string)::Date
"""
macro catch_as_Nullable(expr, catch_errors=:ArgumentError)
    @assert expr.head ≡ :(::) "Use this macro with value::type."
    (expr_value, expr_type) = map(esc, expr.args)
    quote
        try
            Nullable{$expr_type}($expr_value)
        catch e
            if isa(e, $(esc(catch_errors)))
                Nullable{$expr_type}()
            else
                rethrow()
            end
        end
    end
end

"Test if `i` is below maximum number of iterations (when given)."
below_maxiter(i, maxiter::Nullable{Int}) = isnull(maxiter) || (i ≤ get(maxiter))
below_maxiter(i, maxiter::Int) = i ≤ maxiter

"""
When `key ∈ dict`, return the corresponding value as a `Nullable`.

Otherwise, with `probability`, add `key` with the value `default`, and
return that as a `Nullable`, or with `1-probability`, return a
`Nullable` without value.

Useful for keeping a random subsample of `key`s, with the given default.
"""
function random_push!{T,S}(dict::Dict{T,S}, key::T, probability::Real,
                           default::Function)
    if haskey(dict, key)
        Nullable(dict[key])
    elseif rand() ≤ probability
        value = default()
        dict[key] = value
        Nullable(value)
    else
        Nullable{S}()
    end
end

"""
Define using `key => value` pairs, define a dictionary in `dictvar`,
and an `enum` names `enumname` using `value`s.

Useful because one does not need to enumerate the values twice.
"""
macro enum_and_Dict(enumname, dictvar, pairs...)
    function expr_second(expr)
        @assert expr.head ≡ :(=>) "$expr is not a Pair."
        esc(expr.args[2])
    end
    quote
        @enum($(esc(enumname)), $(unique(map(expr_second, pairs))...))
        $(esc(dictvar)) = Dict($(map(esc, pairs)...))
    end
end

"""Open a zip compressed file, read (and discard) the first line when
`skip_first`, then call `f` with the stream, finally close it."""
function process_compressed_database(f, path; skip_first = true)
    io = open(path) |> Libz.ZlibInflateInputStream
    try
        skip_first && readline(io)  # skip header
        f(io)
    finally
        close(io)
    end
end

"""
Read lines and call function on the result.
"""
function Base.readlines(f::Function, io::IO, maxiter = Nullable{Int}())
    i = 1
    while !eof(io) && below_maxiter(i, maxiter)
        line = chomp(readline(io))
        f(line)
        i += 1
    end
end

"""
Read from the gzip-compressed `path`, and return the deserialied contents.
"""
deserialize_gz(path::String) = GZip.open(deserialize, path, "r")

# NOTE until the following issue is fixed:
# https://github.com/JuliaLang/DataStructures.jl/issues/259
function Base.normalize{T,V}(acc::Accumulator{T,V}, p=1)
    _norm = norm(collect(values(acc)), p)
    Dict(key => value/_norm for (key,value) in acc.map)
end


text:
function test_discrete_sets(T)
    n = 10
    s1 = DiscreteVectorDictionary{T}(n)
    @test domaintype(s1) == Int
    @test codomaintype(s1) == T
    @test length(s1) == n
    @test isdiscrete(s1)

    s2 = resize(s1, 2*n)
    @test domaintype(s2) == Int
    @test codomaintype(s2) == T
    @test length(s2) == 2*n

    # Support
    @test in_support(s1, 1, 4)
    @test !in_support(s1, 1, 0)
    @test in_support(s1, 1, length(s1))
    @test !in_support(s1, 1, length(s1)+1)
    @test !in_support(s1, 1, 2.5)

    @test support(s1) == ClosedInterval{Int}(1, length(s1))

    # Evaluation
    @test eval_element(s1, 1, 1) == one(T)
    @test eval_element(s1, 1, 2) == zero(T)
    @test eval_element(s1, n, n) == one(T)

    t1 = DiscreteArrayDictionary{T}((10,10))
    @test codomaintype(t1) == T
    @test domaintype(t1) == ProductIndex{2}
    @test size(t1) == (10,10)

    @test in_support(t1, (1,1), 5)
    @test !in_support(t1, (1,1), 101)
    @test in_support(t1, (1,1), (3,3))
    @test !in_support(t1, (1,1), (11,3))
    @test in_support(t1, (1,1), CartesianIndex(5,3))
    @test !in_support(t1, (1,1), CartesianIndex(5,11))

    t2 = resize(t1, (20,20))
    @test codomaintype(t2) == T
    @test domaintype(t2) == ProductIndex{2}
    @test size(t2) == (20,20)

    # Evaluation
    @test eval_element(t1, 1, 1) == one(T)
    @test eval_element(t1, 1, (1,1)) == one(T)
    @test eval_element(t1, 1, CartesianIndex(1,1)) == one(T)
    @test eval_element(t1, (1,1), 1) == one(T)
    @test eval_element(t1, CartesianIndex(1,1), 1) == one(T)
    @test eval_element(t1, CartesianIndex(1,1), (1,1)) == one(T)

    @test eval_element(t1, 2, 5) == zero(T)
    @test eval_element(t1, (2,1), 5) == zero(T)
    @test eval_element(t1, (2,1), CartesianIndex(3,3)) == zero(T)

end


text:
using AutoGrad

"""

    gradcheck(f, x...; kwargs...)

Numerically check the gradient of `f(x...)` and return a boolean
result.

Each argument can be a Number, Array, Tuple or Dict which in turn can contain other Arrays
etc.  Only 10 random entries in each large numeric array are checked by default.  If the
output of `f` is not a number, we check the gradient of `sum(f(x...))`. See also `gcheck`
for a different take on marking parameters.

# Keywords

* `args=:`: the argument indices to check gradients with respect to. Could be an array or
  range of indices or a single index. By default all arguments that have a `length` method
  are checked.

* `kw=()`: keyword arguments to be passed to `f`.

* `nsample=10`: number of random entries from each numeric array in gradient
  `dw=(grad(f))(w,x...;o...)` compared to their numerical estimates.

* `atol=0.01`: tolerance parameter.  See `isapprox` for an explanation.

* `rtol=0.05`: tolerance parameter.  See `isapprox` for an explanation.

* `delta=0.0001`: step size for numerical gradient calculation.

* `verbose=1`: 0 prints nothing, 1 shows failing tests, 2 shows all tests.

"""
function gradcheck(f, x...; kw=(), args=:, nsample=10, verbose=1, rtol=0.05, atol=0.01, delta=0.0001)
    args = isa(args, Colon) ? (1:length(x)) : args
    xrec  = Any[x...]
    for i in args; xrec[i] = Param(xrec[i]); end
    result = @diff gcsum(f, xrec...; kw...)
    f0 = value(result)
    xptr = Any[x...]
    gptr = Array{Any}(undef, length(x))
    for i in args; gptr[i] = full(grad(result, xrec[i])); end
    all(args) do i
        gcwalk(i, xptr, gptr, f0, f, xptr, kw, nsample, verbose, delta, rtol, atol)
    end
end

function gcsum(f,x...;o...)
    y = f(x...;o...)
    v = value(y)
    if isa(v,Number)
        return y
    elseif isempty(v)
        return 0
    else
        return sum(y)
    end
end

function gcwalk(i, xptr, gptr, f0, f, x, kw, nsample, verbose, delta, rtol, atol)
    if isa(value(xptr[i]), Number)
        if isa(xptr[i], Param)
            xi = xptr[i].value
            delta = delta > 0 ? oftype(float(xi),delta) : cbrt(eps(xi))
            xptr[i].value = xi >= 0 ? xi + delta : xi - delta
            f1 = gcsum(f, x...; kw...)
            nd = (f1 - f0) / (xptr[i] - xi)
            xptr[i].value = xi
        else
            xi = xptr[i]
            delta = delta > 0 ? oftype(float(xi),delta) : cbrt(eps(xi))
            xptr[i] = xi >= 0 ? xi + delta : xi - delta
            f1 = gcsum(f, x...; kw...)
            nd = (f1 - f0) / (xptr[i] - xi)
            xptr[i] = xi
        end
        ad = gcget(gptr,i,0)
        result = isapprox(nd, ad, rtol=rtol, atol=atol)
        if verbose >= 2 || (!result && verbose >= 1)
            #fn = (f==broadcast ? x[1] : f)
            pa = summary(xptr)
            @show pa,xi,f0,nd,ad
        end
        return result
    elseif !isempty(methods(length,(typeof(xptr[i]),)))
        n = length(xptr[i])
        if isa(xptr[i], Tuple)
            xptr[i] = Any[xptr[i]...]
            indices = 1:n
        elseif isa(xptr[i], Dict)
            indices = keys(xptr[i])
        elseif isbitstype(eltype(xptr[i]))
            indices = n <= nsample ? (1:n) : rand(1:n, nsample)
        else
            indices = 1:n
        end
        return all(indices) do j
            gcwalk(j, xptr[i], gcget(gptr,i,nothing), f0, f, x, kw, nsample, verbose, delta, rtol, atol)
        end
    else
        return true
    end
end

function gcget(g,i,default=nothing)
    if isa(g,Nothing)
        return default
    elseif isa(g,AbstractDict) && !haskey(g,i)
        return default
    elseif isa(g,AbstractArray) && !isassigned(g,i)
        return default
    elseif isa(g,Tuple) && !in(i,1:length(g))
        return default
    elseif isa(g[i],Nothing)
        return default
    end
    return g[i]
end

"Test a numeric function with Float32/64 randn scalars and randn arrays, possibly transforming the input to match the domain"
function randcheck(f,t1=identity,ts...; args=:, kw...)
    ts = [t1,ts...]
    if isa(args,Colon); args=1:length(ts); end
    x64 = map(t->t(randn()), ts)
    a64 = map(t->t.(randn(2)), ts)
    x32 = map(t->t(randn(Float32)), ts)
    a32 = map(t->t.(randn(Float32,2)), ts)
    gradcheck(f, x32...; args=args, kw...) &&
    gradcheck(broadcast, f, a32...; args=(args.+1), kw...) &&
    gradcheck(f, x64...; args=args, kw...) &&
    gradcheck(broadcast, f, a64...; args=(args.+1), kw...)
end

"""

    gcheck(f, x...; kw, o...)
    @gcheck f(x...; kw...) (opt1=val1,opt2=val2,...)

Numerically check the gradient of `f(x...; kw...)` and return a boolean result.

Example call: `gcheck(nll,model,x,y)` or `@gcheck nll(model,x,y)`. The parameters should be
marked as `Param` arrays in `f`, `x`, and/or `kw`.  Only 10 random entries in each large
numeric array are checked by default.  If the output of `f` is not a number, we check the
gradient of `sum(f(x...; kw...))`. Keyword arguments:

* `kw=()`: keyword arguments to be passed to `f`, i.e. `f(x...; kw...)`
* `nsample=10`: number of random entries from each param to check
* `atol=0.01,rtol=0.05`: tolerance parameters.  See `isapprox` for their meaning.
* `delta=0.0001`: step size for numerical gradient calculation.
* `verbose=1`: 0 prints nothing, 1 shows failing tests, 2 shows all tests.

"""
gcheck, :(@gcheck)

function gcheck(f, x...; kw=(), nsample=10, verbose=1, rtol=0.05, atol=0.01, delta=0.0001)
    y = @diff gcsum(f, x...; kw...)
    if !isa(y, Tape); @warn("Output independent of params"); return true; end
    f0 = value(y)
    ps = Param[ n.Value for n in y.list if isa(n.Value, Param) ]
    if isempty(ps); @error("Cannot find any params"); end
    #vs = value.(ps)
    gs = (p->full(grad(y,p))).(ps)
    all(1:length(ps)) do i
        #gcwalk(i, vs, gs, f0, f, x, kw, nsample, verbose, delta, rtol, atol)
        gcwalk(i, ps, gs, f0, f, x, kw, nsample, verbose, delta, rtol, atol)
    end
end

macro gcheck(fx,options=:(NamedTuple()))
    :(gcheck(()->$(esc(fx));$(esc(options))...))
end


text:
"""
Functions for estimating an intensity function Ω from observed data. Throughout, H is an object of type hypergraph and Z is a label vector of the same length as H.D. The entries of Z are integers between 1 and kmax, where kmax is the number of clusters. Z is assumed to contain at least one entry with value k for each k ∈ [kmax].
"""

function estimateΩEmpirically(H, Z; pseudocount=0.0, aggregator=identity, bigNums = true, two_groups = false)
    """
    Estimate an intensity function Ω directly from data, with no parameterization.
    The intensity function is assumed to be symmetric under permutations of node labels.
    Equivalently, it is a function of the partition vector p alone.

    # Arguments
    - H::hypergraph, a hypergraph
    - Z::Vector{Int4}, a label vector
    - min_val::Float64, the pseudocount used when no hyperedges with a given specified feature key are observed. Altering the default is not recommended.
    - aggregator::function, an aggregation function that maps a partition vector p to a feature vector summarizing p.
    - bigNums::bool, whether to use bigInt and bigFloat in computation. Altering the default is not recommended except on exceptionally small data instances and with great caution.

    # Returns

    - Ω::IntensityFunction, giving the maximum likelihood estimates of the intensity parameter for each possible value of the supplied aggregator.
    """

    ℓ = maximum(keys(H.E))  # size of largest hyperedge
    S = evalSums(Z,H,maximum(Z),true)[3]

    if two_groups
        P = keys(S)
    else
        P = partitionsUpTo(ℓ)
    end

    C = Dict(p => 0.0 for p in P)
    for k in keys(H.E)
        Ek = H.E[k]
        for e in keys(Ek)
            p = partitionize(Z[e])
            C[p] = get(C, p, 0) + Ek[e]
        end
    end

    if bigNums
        C_agg = Dict(aggregator(p) => big(0.0) for p in keys(S))
        S_agg = Dict(aggregator(p) => big(0.0) for p in keys(S))
    else
        C_agg = Dict(aggregator(p) => 0.0 for p in keys(S))
        S_agg = Dict(aggregator(p) => 0.0 for p in keys(S))
    end

    for p in keys(S)
        agg_key = aggregator(p)
        if bigNums
            C_agg[agg_key] += get(C, p, big(pseudocount))
            S_agg[agg_key] += get(S, p, big(pseudocount))
        else
            C_agg[agg_key] += get(C, p, pseudocount)
            S_agg[agg_key] += get(S, p, pseudocount)
        end
    end

    ω̂ = Dict(a => abs(C_agg[a] / S_agg[a]) for a in keys(C_agg))

    return empiricalIntensityFunction((p, α) -> ω̂[p], ℓ, aggregator)
end

function learnParameters(H, Z, Ω, α0; max_iters = 10, verbose = false, tol = 1e-2)
    """
    Estimate the parameter vector α for a parameterized IntensityFunction Ω from a hypergraph H and label vector Z.
    Currently uses a very simple coordinatewise method, and is quite slow for larger data. Tuning amin and amax may be necessary to achieve appropriate results.

    # Arguments

    - H::hypergraph, a hypergraph
    - Z::Vector{Int4}, a label vector
    - Ω::IntensityFunction, the intensity function parameterized by α
    - α0::Vector{Float64}, the initial guess for α.
    - n_iters::Int64, the number of out iterations in coordinate ascent to perform
    - amin::Float64 and amax::Float64, box constraints on the entries of α

    # Returns

    α: Vector{Float64}, the parameter vector obtained via repeated optimization of the likelihood objective.
    """
    kmax = length(α0) ÷ 2

    objs = formObjectives(H, Z, Ω)

    function G(a, α, k, pos)
        α_ = copy(α)
        α_[(pos == 1) ? k : (kmax + k)] = a[1]
        return objs[k](α_)
    end

    α = copy(α0);

    for i ∈ 1:max_iters
        old_val = sum(objs[k](α) for k in 1:kmax)
        for pos ∈ [2,1]
            for k ∈ 1:kmax
                ix = (pos == 1) ? k : (kmax + k)
                a0 = α[ix]
                res = Optim.optimize(a -> G(a, α, k, pos), [a0], Optim.LBFGS())
                α[ix] = Optim.minimizer(res)[1]
            end
        end
        new_val = sum(objs[k](α) for k in 1:kmax)
        if abs((new_val - old_val)) < tol
            return α
        end
        if verbose
            println("Q = $(Float64(modularity(H, Z, Ω;α = α)))")
        end
    end
    return(α)
end

function formObjectives(H, Z, Ω)
    ℓ       = maximum(Z)
    C       = evalCuts(H,Z,Ω)
    V, μ, S = evalSums(Z,H,ℓ,true)
    Ŝ       = aggregateSums(S,Ω)
    kmax    = maximum(keys(H.E))

    function obj(α, k)
        val = 0.0
        for p in keys(Ŝ)
            if sum(p) == k
                Op   = Ω.ω(p, α)
                val += get(C, p, 0)*log(Op) - Ŝ[p]*Op
            end
        end
        return - val
    end
    return [(α -> obj(α, k)) for k in 1:kmax]
end

function formObjective(H, Z, Ω)
    ℓ       = maximum(Z)
    C       = evalCuts(H,Z,Ω)
    V, μ, S = evalSums(Z,H,ℓ,true)
    Ŝ       = aggregateSums(S,Ω)

    function objective(α)
        obj = 0.0
        for p in keys(Ŝ)
            Op   = Ω.ω(p, α)
            obj += get(C, p, 0)*log(Op) - Ŝ[p]*Op
        end
        return -obj # sign is for minimization
    end
    return objective
end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
module ITensorUnicodePlots

using Graphs
using NetworkLayout
using Reexport
using Statistics
@reexport using ITensorVisualizationBase

using ITensorVisualizationBase:
  @Backend_str,
  default_vertex_labels,
  default_vertex_labels_prefix,
  default_vertex_size,
  default_vertex_textsize,
  default_edge_textsize,
  default_edge_widths,
  default_edge_labels,
  default_arrow_show,
  default_arrow_size,
  default_siteinds_direction,
  is_self_loop

using UnicodePlots: UnicodePlots

import ITensorVisualizationBase: visualize

function __init__()
  return ITensorVisualizationBase.set_backend!(Backend"UnicodePlots"())
end

function plot(::Backend"UnicodePlots"; xlim, ylim, width, height)
  plot = UnicodePlots.lineplot(
    [0.0],
    [0.0];
    border=:none,
    labels=false,
    grid=false,
    xlim=xlim,
    ylim=ylim,
    width=width,
    height=height,
  )
  return plot
end

function draw_edge!(b::Backend"UnicodePlots", plot, v1, v2; color)
  UnicodePlots.lineplot!(plot, [v1[1], v2[1]], [v1[2], v2[2]]; color=color)
  return plot
end

function annotate!(::Backend"UnicodePlots", plot, x, y, str)
  UnicodePlots.annotate!(plot, x, y, str)
  return plot
end

# Don't use new lines by default, it messes up the formatting
default_newlines(::Backend"UnicodePlots") = false

function visualize(
  b::Backend"UnicodePlots",
  g::AbstractGraph;
  interactive=false, # TODO: change to `default_interactive(b)`
  ndims=2, # TODO: change to `default_ndims(b)`
  layout=Spring(; dim=ndims), # TODO: change to `default_layout(b, ndims)`

  # vertex
  vertex_labels_prefix=default_vertex_labels_prefix(b, g),
  vertex_labels=default_vertex_labels(b, g, vertex_labels_prefix),
  vertex_size=default_vertex_size(b, g),
  vertex_textsize=default_vertex_textsize(b, g),

  # edge
  edge_textsize=default_edge_textsize(b),
  edge_widths=default_edge_widths(b, g),
  edge_labels=default_edge_labels(b, g),

  # arrow
  arrow_show=default_arrow_show(b, g),
  arrow_size=default_arrow_size(b, g),

  # siteinds direction
  siteinds_direction=default_siteinds_direction(b, g),
  width=50,
  height=20,
)
  edge_color = :blue # TODO: make into keyword argument

  edge_labels = ITensorVisualizationBase.edge_labels(b, edge_labels, g)

  node_pos = layout(g)
  edge_pos = [node_pos[src(edge)] => node_pos[dst(edge)] for edge in edges(g)]
  xmin = minimum(first.(node_pos))
  xmax = maximum(first.(node_pos))
  ymin = minimum(last.(node_pos))
  ymax = maximum(last.(node_pos))

  #vertex_size = vertex_size * (xmax - xmin)

  xscale = 0.1 * (xmax - xmin)
  yscale = max(0.3 * (ymax - ymin), 0.1 * xscale)
  xlim = [xmin - xscale, xmax + xscale]
  ylim = [ymin - yscale, ymax + yscale]

  site_vertex_shift = siteinds_direction

  #site_vertex_shift = -Point(0, 0.2 * abs(ylim[2] - ylim[1]))
  #site_vertex_shift = -Point(0, 0.001 * (xmax - xmin))

  # Initialize the plot
  plt = plot(b; xlim=xlim, ylim=ylim, width=width, height=height)

  # Add edges and nodes
  for (e_pos, e) in zip(edge_pos, edges(g))
    if is_self_loop(e)
      draw_edge!(b, plt, e_pos[1], e_pos[1] + site_vertex_shift; color=edge_color)
    else
      draw_edge!(b, plt, e_pos[1], e_pos[2]; color=edge_color)
    end
  end

  # Add edge labels and node labels
  for (n, e) in enumerate(edges(g))
    e_pos = edge_pos[n]
    edge_label = edge_labels[n]
    if is_self_loop(e)
      @assert e_pos[1] == e_pos[2]
      str_pos = e_pos[1] + 0.5 * site_vertex_shift
      annotate!(b, plt, str_pos..., edge_label)
    else
      annotate!(b, plt, mean(e_pos)..., edge_label)
    end
  end
  if length(vertex_labels) ≠ nv(g)
    throw(
      DimensionMismatch(
        "Number of vertex labels must equal the number of vertices. Vertex labels $(vertex_labels) of length $(length(vertex_labels)) does not equal the number of vertices $(nv(g)).",
      ),
    )
  end
  for v in vertices(g)
    x, y = node_pos[v]
    node_label = vertex_labels[v]
    annotate!(b, plt, x, y, node_label)
  end
  return plt
end

end


text:
@symbol_func function G_PU(cur_reactor::AbstractReactor)
  cur_left_part = K_RP(cur_reactor)

  cur_left_part *= cur_reactor.tau_FT

  cur_right_part = K_LP(cur_reactor)

  cur_right_part *= safe_symbol(cur_reactor, :R_0) ^ 2

  cur_right_part *= cur_reactor.T_bar ^ (3/2)

  cur_G = cur_left_part

  cur_G += cur_right_part

  cur_G
end


text:
#
# This is a package to convert to use zero origin array in a function.
# Usage:
#  @zeroorigin [single symbol or list/tuple for variables with zeroorigin] [definition of function]
#
# When we definie a function, the code `@zeroorigin (x,y)` is added to the head of definition.
# Then in the function, the ref of x and y, i.e., x[i] and y[i] will be changed to x[i+1] and y[i+1].
#
# Example:
#
# using Origin
# using Test
#
# @testset "Origin.jl" begin
#     @origin (a=>0, b=>2) function test()
#         a = collect(0:10)
#         b = collect(1:5)
#         @test a[0] == 0
#         @test a[10] == 10
#         @test b[2] == 1
#         @test b[4] == 3
#         @test b[end] == 5
#         @test b[2:end] == [1, 2, 3, 4, 5]
#         @test b[a[5]] == 4
#     end

#     test()
# end

export @origin

function peek(stack)
    return stack[end]
end

function _changeindex(expr, u)
    if Meta.isexpr(expr, :call) && expr.args[1] == :(:)
        if length(expr.args) == 3
            Expr(:call, :(:), Expr(:call, :+, expr.args[2], u), Expr(:call, :+, expr.args[3], u))
        elseif length(expr.args) == 4
            Expr(:call, :(:), Expr(:call, :+, expr.args[2], u), expr.args[3], Expr(:call, :+, expr.args[4], u))
        else
            error("The inner of index should be a simple expr or call :(:)")
        end
    else
        Expr(:call, :+, expr, u)
    end
end

function _replace(expr::Expr, symlist, stack)
    if Meta.isexpr(expr, :ref)
        push!(stack, expr.args[1])
        args = Any[_replace(x, symlist, stack) for x = expr.args]
        if haskey(symlist, peek(stack))
            m = symlist[peek(stack)]
            u = Expr(:call, :-, 1, m)
            for i = 2:length(args)
                args[i] = _changeindex(args[i], u)
            end
        end
        pop!(stack)
        Expr(expr.head, args...)
    elseif Meta.isexpr(expr, :macrocall) && expr.args[1] == Symbol("@origin")
        macroexpand(Origin, expr, recursive=false)
    else
        if Meta.isexpr(expr, :call) && expr.args[1] == :eachindex
            if _hassym(expr.args[2], symlist)
                println("Warning: This macro has not been applied to the block including 'eachindex' yet.")
            end
        end
        args = Any[_replace(x, symlist, stack) for x = expr.args]
        Expr(expr.head, args...)
    end
end

function _replace(expr::Symbol, symlist, stack)
    if expr == :end && haskey(symlist, peek(stack))
        m = symlist[peek(stack)]
        Expr(:call, :-, :end, Expr(:call, :-, 1, m))
    else
        expr
    end
end

function _replace(expr::Any, symlist, stack)
    expr
end

function _hassym(expr::Expr, symlist)
    any(Bool[_hassym(x, symlist) for x = expr.args])
end

function _hassym(expr::Symbol, symlist)
    haskey(symlist, expr)
end

function _add!(symlist, x)
    if typeof(x) <: Expr && x.head == :call && x.args[1] == :(=>)
        symlist[x.args[2]] = x.args[3]
    end
end

macro origin(vars, expr)
    symlist = Dict{Symbol,Any}()
    stack = []
    if typeof(vars) <: Expr && (vars.head == :vect || vars.head == :tuple)
        for x = vars.args
            _add!(symlist, x)
        end
    else
        _add!(symlist, vars)
    end
    ex = _replace(expr, symlist, stack)
    esc(ex)
end


text:
using CautiousEngine

@testset "adding and chanding GP data" begin
    N = 3
    x = rand(2, N)
    y = [sum(x[:, i]) for i=1:N]
    @test size(x, 2) == length(y)
    test_gp = CautiousEngine.train_gp_1dim(x, y)
    @test x == test_gp.x
    @test y == test_gp.y

    new_x = rand(2,1)
    new_y = sum(new_x)
    new_datapoint = (new_x, new_y)

    CautiousEngine.add_datapoint_to_gp(test_gp, new_datapoint)
    @test test_gp.nobs == N+1

    # Test Updating GP
    M = 5
    x_replace = rand(2,M)
    @test size(x_replace, 2) == M
    y_replace = [sum(x_replace[:, i]) for i=1:M]
    updated_gp = CautiousEngine.update_gp_data(test_gp, x_replace, y_replace)
    @test updated_gp.x === x_replace
    @test updated_gp.y === y_replace
    @test updated_gp.nobs == M

    # Just check if dimension is correct
    @test updated_gp.dim == 2

    # KNN Test
    P = 10
    x = rand(2,P)
    @info typeof(x)
    y = rand(1,P)
    center = [0.;0.]
    NN = 3
    x_sub, y_sub = CautiousEngine.get_local_data_knn(center, x, y; num_neighbors=NN)
    @test size(x_sub, 2) == NN
    @test length(y_sub) == NN
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 845202ce-23d9-11eb-085e-ad2c26fd531c
begin
	using TextAnalysis
	using CSV
	using DataFrames
end

# ╔═╡ cac12010-23d8-11eb-0321-3fc66858fc84
md"# Authorship Analysis using TextAnalysis.jl

In this notebook we will implement an authorship analysis technique that uses n-gram frequencies for identifying authors. This technique is described in [Grieve 2018](https://research.birmingham.ac.uk/portal/files/53402456/Bixby_PREPRINT.pdf) The technique counts token and character n-grams and compares them to a reference corpus. It counts the number of overlapping n-grams in the reference corpus and identifies the author based on which reference corpus contains more of the n-grams from the text."

# ╔═╡ f400d820-23d9-11eb-091a-898d693456ad
md"First off, the necessary imports are made."

# ╔═╡ c9d140dc-2463-11eb-06a2-2d0dcb8dbed2
md"Our data will be from [Spooky Author Identification](https://www.kaggle.com/c/spooky-author-identification/overview) Kaggle contest. The data comes in the form of a CSV file which contains sentences attributed to one of three authors. The data was manualy cleaned to remove rows that were corropted as well as making sure that each author had approximately equal word counts.

The data is read into a dataframe after which we iterate over the df and turn each sentence into a StringDocument which is then accumalted into a Corpus for each author."

# ╔═╡ 68255d02-23da-11eb-3232-613c0e1b0ee2
referenceData = CSV.read("train.csv", DataFrame);

# ╔═╡ c8b1e874-23e0-11eb-1aad-995c03f040ef
begin
	MWS = GenericDocument[]
	EAP = GenericDocument[]
	HPL = GenericDocument[]

	for line in eachrow(referenceData)
		if line.author == "MWS"
			push!(MWS, StringDocument(line.text))
		elseif line.author == "EAP"
			push!(EAP, StringDocument(line.text))
		else
			push!(HPL, StringDocument(line.text))
		end
	end

	MWS_corpus = Corpus(MWS)
	EAP_corpus = Corpus(EAP)
	HPL_corpus = Corpus(HPL);
end

# ╔═╡ 60390844-2472-11eb-010b-d155abe74997
md"The following code snippet allows one to retrieve the word count of a corpus. The word counts are as follows:
- MWS_corpus: 147914
- EAP_corpus: 147911
- HPL_corpus: 147916"

# ╔═╡ 7b9ce77e-23e4-11eb-39d4-03e26b33cad7
begin
	wordCount = 0
	for doc in HPL_corpus
		global wordCount += length(ngrams(doc, 1))
	end
end

# ╔═╡ 43a58408-246f-11eb-29ab-91cc00373c2a
wordCount

# ╔═╡ 96f5d60a-248b-11eb-2795-a51a9fbf3fc0
md"Now that we have our corpora it is time to extract the ngrams that are needed for the analysis. The original paper used both character ngrams and token ngrams. Due to TextAnalysis.jl's functionality I will be restricting this analysis to only using token ngrams."

# ╔═╡ 95691ad4-2472-11eb-346e-07a9f5dc8f6c
begin
	MWS_ngrams = []
	HPL_ngrams = []
	EAP_ngrams = []
	for i in 1:4
		push!(MWS_ngrams, map(x -> ngrams(x,i), MWS_corpus))
		push!(HPL_ngrams, map(x -> ngrams(x,i), HPL_corpus))
		push!(EAP_ngrams, map(x -> ngrams(x,i), EAP_corpus))
	end
end

# ╔═╡ e207949e-2526-11eb-3ae4-0712b7082681
md"We also retrieve the token ngrams from the files we want to analyze and classify."

# ╔═╡ 62443a4c-2525-11eb-0e70-393860da7b4f
begin
	# first we load our files from memory
	text1 = FileDocument("./MWS.txt")
	text2 = FileDocument("./HPL.txt")
	text3 = FileDocument("./EAP.txt")

	# this is the same as retrieveing ngrams for the reference corpora
	text1_ngrams = []
	text2_ngrams = []
	text3_ngrams = []
	for i in 1:4
		push!(text1_ngrams, ngrams(text1, i))
		push!(text2_ngrams, ngrams(text2, i))
		push!(text3_ngrams, ngrams(text3, i))
	end
end

# ╔═╡ 32909f1e-2527-11eb-0efd-6d81943df28e
md"Finally we will a script that calculates the intersection of the ngrams from the text and the three reference corpora. To try different files edit the second for loop."

# ╔═╡ 5d1320ee-2479-11eb-3664-c33f5ace1ba7
begin
	MWS_count = 0
	HPL_count = 0
	EAP_count = 0
	for i in 1:4
		for j in text1_ngrams[i]
			for k in MWS_ngrams[i]
				if haskey(k, j[1])
					global MWS_count += 1
					break
				end
			end
			for k in HPL_ngrams[i]
				if haskey(k, j[1])
					global HPL_count += 1
					break
				end
			end
			for k in EAP_ngrams[i]
				if haskey(k, j[1])
					global EAP_count += 1
					break
				end
			end
		end
	end
end

# ╔═╡ 16a5c100-252a-11eb-1258-b769d50904d3
md"The following block will output the author as identified by this script.
The correct authors should be as follows:
- text1 = MWS
- text2 = HPL
- text3 = EAP"

# ╔═╡ df84c6ee-247a-11eb-24f9-e92653c334bb
begin
	if MWS_count > HPL_count && MWS_count > EAP_count
		"MWS"
	elseif HPL_count > MWS_count && HPL_count > EAP_count
		"HPL"
	else
		"EAP"
	end
end

# ╔═╡ 34de486c-252b-11eb-2071-f12f0f445d9f
md"## Conclusion
While this is a very watered down analysis of what the paper presents it is a good starting point and shows that this technique is viable as it correctly predicts the authors of our random test data. One way to improve the current analysis is to add character ngrams as the paper does.

All in all this was a fun project to learn more about Julia and its NLP cabalities. This project only scratches the surface of what Julia has to offer given that the TextAnalysis.jl package offers things such as part-of-speech(POS) tagging, tf-idf measures, a sentiment analysis model, and much more."

# ╔═╡ Cell order:
# ╟─cac12010-23d8-11eb-0321-3fc66858fc84
# ╟─f400d820-23d9-11eb-091a-898d693456ad
# ╠═845202ce-23d9-11eb-085e-ad2c26fd531c
# ╟─c9d140dc-2463-11eb-06a2-2d0dcb8dbed2
# ╠═68255d02-23da-11eb-3232-613c0e1b0ee2
# ╠═c8b1e874-23e0-11eb-1aad-995c03f040ef
# ╟─60390844-2472-11eb-010b-d155abe74997
# ╠═7b9ce77e-23e4-11eb-39d4-03e26b33cad7
# ╠═43a58408-246f-11eb-29ab-91cc00373c2a
# ╟─96f5d60a-248b-11eb-2795-a51a9fbf3fc0
# ╠═95691ad4-2472-11eb-346e-07a9f5dc8f6c
# ╟─e207949e-2526-11eb-3ae4-0712b7082681
# ╠═62443a4c-2525-11eb-0e70-393860da7b4f
# ╟─32909f1e-2527-11eb-0efd-6d81943df28e
# ╠═5d1320ee-2479-11eb-3664-c33f5ace1ba7
# ╟─16a5c100-252a-11eb-1258-b769d50904d3
# ╠═df84c6ee-247a-11eb-24f9-e92653c334bb
# ╟─34de486c-252b-11eb-2071-f12f0f445d9f

max_stars_repo_name: Humans-of-Julia/Challenges

max_issues_repo_name: Humans-of-Julia/WeeklyContest

max_forks_repo_name: Humans-of-Julia/WeeklyContest


text:
mutable struct PKContext
    data::Ptr{Cvoid}

    function PKContext()
        ctx = new()
        ctx.data = Libc.malloc(32)

        ccall((:mbedtls_pk_init, libmbedcrypto), Cvoid, (Ptr{Cvoid},), ctx.data)

        finalizer(ctx->begin
            ccall((:mbedtls_pk_free, libmbedcrypto), Cvoid, (Ptr{Cvoid},), ctx.data)
            Libc.free(ctx.data)
        end, ctx)
        ctx
    end
end

function parse_keyfile!(ctx::PKContext, path, password="")
    @err_check ccall((:mbedtls_pk_parse_keyfile, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Cstring, Cstring),
        ctx.data, path, password)
end

function parse_keyfile(path, password="")
    ctx = PKContext()
    parse_keyfile!(ctx, path, password)
    ctx
end

function parse_public_keyfile!(ctx::PKContext, path)
    @err_check ccall((:mbedtls_pk_parse_public_keyfile, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Cstring),
        ctx.data, path)
end

function parse_public_keyfile(path)
    ctx = PKContext()
    parse_public_keyfile!(ctx, path)
    ctx
end

function parse_public_key!(ctx::PKContext, key)
    key_bs = String(key)
    @err_check ccall((:mbedtls_pk_parse_public_key, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cuchar}, Csize_t),
        ctx.data, key_bs, sizeof(key_bs) + 1)
end

function parse_key!(ctx::PKContext, key, maybe_pw = nothing)
    key_bs = String(key)
    if maybe_pw === nothing
        pw = C_NULL
        pw_size = 0
    else
        pw = String(maybe_pw)
        pw_size = sizeof(pw)  # Might be off-by-one
    end
    @err_check ccall((:mbedtls_pk_parse_key, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{Cuchar}, Csize_t, Ptr{Cuchar}, Csize_t),
        ctx.data, key_bs, sizeof(key_bs) + 1, pw, pw_size)
end

function bitlength(ctx::PKContext)
    sz = ccall((:mbedtls_pk_get_bitlen, libmbedcrypto), Csize_t,
        (Ptr{Cvoid},), ctx.data)
    sz >= 0 || mbed_err(sz)
    Int(sz)
end

function decrypt!(ctx::PKContext, input, output, rng)
    outlen_ref = Ref{Cint}(0)
    @err_check ccall((:mbedtls_pk_decrypt, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Csize_t, Ptr{Cvoid}, Ref{Cint}, Csize_t, Ptr{Cvoid}, Any),
        ctx.data, input, sizeof(input), output, outlen_ref, sizeof(output), c_rng[], rng)
    outlen = outlen_ref[]
    Int(outlen)
end

function encrypt!(ctx::PKContext, input, output, rng)
    outlen_ref = Ref{Cint}(0)
    @err_check ccall((:mbedtls_pk_encrypt, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Csize_t, Ptr{Cvoid}, Ref{Cint}, Csize_t, Ptr{Cvoid}, Any),
        ctx.data, input, sizeof(input), output, outlen_ref, sizeof(output), c_rng[], rng)
    outlen = outlen_ref[]
    Int(outlen)
end

function sign!(ctx::PKContext, hash_alg::MDKind, hash, output, rng)
    outlen_ref = Ref{Csize_t}(sizeof(output))
    @err_check ccall((:mbedtls_pk_sign, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Cint, Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ref{Csize_t}, Ptr{Cvoid}, Any),
        ctx.data, hash_alg, hash, sizeof(hash), output, outlen_ref, c_rng[], rng)
    outlen = outlen_ref[]
    Int(outlen)
end

function sign(ctx::PKContext, hash_alg::MDKind, hash, rng)
    n = Int64(ceil(bitlength(ctx) / 8))
    output = Vector{UInt8}(undef, n)
    @assert sign!(ctx, hash_alg, hash, output, rng) == n
    output
end

function verify(ctx::PKContext, hash_alg::MDKind, hash, signature)
    @err_check ccall((:mbedtls_pk_verify, libmbedcrypto), Cint,
        (Ptr{Cvoid}, Cint, Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
        ctx.data, hash_alg, hash, sizeof(hash), signature, sizeof(signature))
end

function get_name(ctx::PKContext)
    ptr = ccall((:mbedtls_pk_get_name, libmbedcrypto), Ptr{Cchar}, (Ptr{Cvoid},), ctx.data)
    unsafe_string(convert(Ptr{UInt8}, ptr))
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ErrLogEntry
    time::DateTime
    thread::Int
    process::Int
    error::Exception
end

const g_error_log = Ref{Union{Vector{ErrLogEntry},Nothing}}(nothing)

"""
    BAT.error_log()

*Experimental feature, not part of stable public API.*

Get a log of certain exceptions throws by BAT, e.g. density evaluation
errors.

The error log is disabled by default, use [`BAT.enable_error_log`](@ref)
to enable it.
"""
error_log() = g_error_log[]

"""
    BAT.enable_error_log(enable::Bool = true)

*Experimental feature, not part of stable public API.*

Enable/disable BAT's error (exception) log.

The error log is disabled by default.

See [`BAT.error_log`](@ref).
"""
function enable_error_log(enable::Bool = true)
    if enable
        if isnothing(g_error_log[])
            g_error_log[] = Vector{ErrLogEntry}[]
        end
    else
        g_error_log[] = nothing
    end
    nothing
end

function log_error(entry::ErrLogEntry)
    if !isnothing(error_log())
        push!(error_log(), entry)
    end
    if myid() != 1
        remotecall(() -> log_error(entry), 1)
    end
end

log_error(err::Exception) = log_error(ErrLogEntry(now(), threadid(), myid(), err))

macro throw_logged(expr)
    quote
        err = $(esc(expr))
        log_error(err)
        throw(err)
    end
end

macro rethrow_logged(expr)
    quote
        err = $(esc(expr))
        log_error(err)
        rethrow(err)
    end
end


text:
using GasChromatographyTools
using GasChromatographySimulator

# Options
opt = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "inlet", true)
GasChromatographySimulator.Options()
# System
L = 10.0
a_d = [0.25e-3]
a_df = [0.25e-6]
d(x) = GasChromatographySimulator.gf_const(x, a_d)
df(x) = GasChromatographySimulator.gf_const(x, a_df)
sp = "SLB5ms"
gas = "He"
sys = GasChromatographySimulator.constructor_System(L, a_d[1], a_df[1], sp, gas)

# Program
time_steps = [0.0, 60.0, 600.0, 300.0]
temp_steps = [40.0, 40.0, 300.0, 300.0]
pin_steps = [200.0, 200.0, 300.0, 300.0].*1000.0 .+ 101300
pout_steps = [101.3, 101.3, 101.3, 101.3].*1000.0
ΔT_steps = zeros(length(time_steps))
x₀_steps = zeros(length(time_steps))
L₀_steps = L.*ones(length(time_steps))
α_steps = zeros(length(time_steps))
a_gf = [ΔT_steps x₀_steps L₀_steps α_steps]
gf(x) = GasChromatographySimulator.gradient(x, a_gf)
T_itp = GasChromatographySimulator.temperature_interpolation(time_steps, temp_steps, gf, sys.L)
pin_itp = GasChromatographySimulator.pressure_interpolation(time_steps, pin_steps)
pout_itp = GasChromatographySimulator.pressure_interpolation(time_steps, pout_steps)
prog =  GasChromatographySimulator.constructor_Program(time_steps, temp_steps, pin_steps, pout_steps, ΔT_steps, x₀_steps, L₀_steps, α_steps, opt.Tcontrol, L)

# Substance
sub = Array{GasChromatographySimulator.Substance}(undef, 2)
sub[1] = GasChromatographySimulator.Substance("solute 1",
                                                "0-0-0",
                                                400.0,
                                                30.0,
                                                100.0,
                                                1e-3,
                                                "test",
                                                1e-4,
                                                0.0,
                                                0.0)
sub[2] = GasChromatographySimulator.Substance("solute 2",
                                                "1-0-0",
                                                420.0,
                                                31.0,
                                                100.0,
                                                1e-3,
                                                "test",
                                                1e-4,
                                                0.0,
                                                0.0)

# Parameters
par = GasChromatographySimulator.Parameters(sys, prog, sub, opt)

#----test RT-locking functions--------------------------------
sol_n1 = GasChromatographySimulator.solve_system_multithreads(par)
tR_lock_n1 = sol_n1[2].u[end][1]

stretched_prog = GasChromatographyTools.stretched_program(2.0, par)
sol_n2 = GasChromatographySimulator.solve_system_multithreads(stretched_prog)
tR_lock_n2 = sol_n2[2].u[end][1]

n_factor = GasChromatographyTools.RT_locking(par, 300.0, 1e-3, "solute 2")

#---test nb_solutes_RTlock.jl
opt_C30 = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "outlet", true)
sys_C30 = GasChromatographySimulator.constructor_System(4.0, 0.1e-3, 0.1e-6, "Rxi5MS", "He")
prog_C30 =  GasChromatographySimulator.constructor_Program([0.0, 10.0, 60.0, 120.0], [50.1518, 50.1518, 310.152, 310.152], [401300.0, 401300.0, 576300.0, 576300.0], [0.0, 0.0, 0.0, 0.0], [20.0, 20.0, 20.0, 20.0], [0.0, 0.0, 0.0, 0.0], [4.0, 4.0, 4.0, 4.0], [0.0, 0.0, 0.0, 0.0], opt_C30.Tcontrol, sys_C30.L)
sub_C30 = [GasChromatographySimulator.Substance("C30",
                                                "638-68-6",
                                                600.477,
                                                37.5137,
                                                181.004,
                                                1e-3,
                                                "Gaida.2021",
                                                5.17619e-5,
                                                0.0,
                                                0.0)]
par_C30 = GasChromatographySimulator.Parameters(sys_C30, prog_C30, sub_C30, opt_C30)

tR_lock_C30 = 60.0
tR_tol_C30 = 1e-3
n = GasChromatographyTools.RT_locking(par_C30, tR_lock_C30, tR_tol_C30, "C30"; opt_itp="linear")
n2 = GasChromatographyTools.RT_locking(par_C30, tR_lock_C30, tR_tol_C30, "C30"; opt_itp="spline")


text:
using Test
using LinearAlgebra
using NonNegLeastSquares
using PyCall

# wrapper functions for convienence
nnls(A,b) = nonneg_lsq(A,b;alg=:nnls)
pivot(A,b) = nonneg_lsq(A,b;alg=:pivot)

# Solve A*x = b for x, subject to x >=0
A = [ 0.53879488  0.65816267
      0.12873446  0.98669198
      0.24555042  0.00598804
      0.80491791  0.32793762 ]

b = [0.888,  0.562,  0.255,  0.077]

# Test that nnls produces the same solution as scipy
x = [0.15512102, 0.69328985] # approx solution from scipy
@test norm(pivot(A,b)-x) < 1e-5

## A second test case
A2 = [ -0.24  -0.82   1.35   0.36   0.35
       -0.53  -0.20  -0.76   0.98  -0.54
        0.22   1.25  -1.60  -1.37  -1.94
       -0.51  -0.56  -0.08   0.96   0.46
        0.48  -2.25   0.38   0.06  -1.29 ]
b2 = [-1.6,  0.19,  0.17,  0.31, -1.27]
x2 = [2.2010416, 1.19009924, 0.0, 1.55001345, 0.0]
@test norm(pivot(A2,b2)-x2) < 1e-5

## Test a bunch of random cases against python

for i = 1:100
  	m,n = rand(1:10),rand(1:10)
  	A3 = randn(m,n)
  	b3 = randn(m)
  	x3,resid = pyopt.nnls(A3,b3)
    resid = norm(b3-A3*pivot(A3,b3))
    if resid > 1e-5
        @test norm(pivot(A3,b3)-x3) < 1e-5
    end
end

## Test a bunch of random cases against nnls
for i = 1:100
    m,n = rand(1:10),rand(1:10)
    A4 = randn(m,n)
    b4 = randn(m)
    x4 = nnls(A4,b4)
    resid = norm(b4-A4*pivot(A4,b4))
    if resid > 1e-5
        @test norm(pivot(A4,b4)-x4) < 1e-5
    end
end


text:
"""
This file defines `resolve()`, which ensures all dependencies are installed.
"""

function _resolve_top_env(load_path)
    top_env = ""
    for env in load_path
        proj = Base.env_project_file(env)
        is_condapkg = proj isa String && Base.project_file_name_uuid(proj, "").uuid == UUID
        depends_on_condapkg = Base.manifest_uuid_path(env, PKGID) !== nothing
        if is_condapkg || depends_on_condapkg
            top_env = proj isa String ? dirname(proj) : env
            break
        end
    end
    top_env == "" && error("no environment in the LOAD_PATH ($load_path) depends on CondaPkg")
    top_env
end

function _resolve_can_skip_1(load_path, meta_file)
    meta = open(read_meta, meta_file)
    if meta !== nothing && meta.version == VERSION && meta.load_path == load_path
        timestamp = max(meta.timestamp, stat(meta_file).mtime)
        skip = true
        for env in [meta.load_path; meta.extra_path]
            dir = isfile(env) ? dirname(env) : isdir(env) ? env : continue
            if isdir(dir)
                if stat(dir).mtime > timestamp
                    skip = false
                    break
                else
                    fn = joinpath(dir, "CondaPkg.toml")
                    if isfile(fn) && stat(fn).mtime > timestamp
                        skip = false
                        break
                    end
                end
            end
        end
        return skip
    else
        return false
    end
end

_convert(::Type{T}, @nospecialize(x)) where {T} = convert(T, x)::T

function _resolve_find_dependencies(io, load_path)
    packages = Dict{String,Dict{String,PkgSpec}}() # name -> depsfile -> spec
    channels = ChannelSpec[]
    pip_packages = Dict{String,Dict{String,PipPkgSpec}}() # name -> depsfile -> spec
    extra_path = String[]
    for env in [load_path; [p.source for p in values(Pkg.dependencies())]]
        dir = isfile(env) ? dirname(env) : isdir(env) ? env : continue
        fn = joinpath(dir, "CondaPkg.toml")
        if isfile(fn)
            env in load_path || push!(extra_path, env)
            _log(io, "Found dependencies: $fn")
            pkgs, chans, pippkgs = parse_deps(TOML.parsefile(fn))
            for pkg in pkgs
                get!(Dict{String,PkgSpec}, packages, pkg.name)[fn] = pkg
            end
            if isempty(chans)
                push!(channels, ChannelSpec("conda-forge"))
            else
                append!(channels, chans)
            end
            for pkg in pippkgs
                get!(Dict{String,PipPkgSpec}, pip_packages, pkg.name)[fn] = pkg
            end
        end
    end
    (packages, channels, pip_packages, extra_path)
end

function _resolve_merge_packages(packages)
    specs = PkgSpec[]
    for (name, pkgs) in packages
        @assert length(pkgs) > 0
        for pkg in values(pkgs)
            @assert pkg.name == name
            push!(specs, pkg)
        end
    end
    sort!(unique!(specs), by=x->x.name)
end

function abspathurl(args...)
    path = abspath(args...)
    if Sys.iswindows()
        path = replace(path, '\\' => '/')
        @assert !startswith(path, "/") # TODO: handle \\machine\... paths
        path = "/$path"
    else
        @assert startswith(path, "/")
    end
    return "file://$path"
end

function _resolve_merge_pip_packages(packages)
    specs = PipPkgSpec[]
    for (name, pkgs) in packages
        @assert length(pkgs) > 0
        versions = String[]
        urls = String[]
        for (fn, pkg) in pkgs
            @assert pkg.name == name
            if startswith(pkg.version, "@")
                url = strip(pkg.version[2:end])
                if startswith(url, ".")
                    url = abspathurl(dirname(fn), url)
                end
                push!(urls, url)
            elseif pkg.version != ""
                push!(versions, pkg.version)
            end
        end
        sort!(unique!(urls))
        sort!(unique!(versions))
        if isempty(urls)
            version = join(versions, ",")
        elseif isempty(versions)
            length(urls) == 1 || error("multiple direct references ('@ ...') given for pip package '$name'")
            version = "@ $(urls[1])"
        else
            error("direct references ('@ ...') and version specifiers both given for pip package '$name'")
        end
        push!(specs, PipPkgSpec(name, version=version))
    end
    sort!(specs, by=x->x.name)
end

function _resolve_diff(old_specs, new_specs)
    # group by package name
    old_dict = Dict{String,Set{PkgSpec}}()
    for spec in old_specs
        push!(get!(Set{PkgSpec}, old_dict, spec.name), spec)
    end
    new_dict = Dict{String,Set{PkgSpec}}()
    for spec in new_specs
        push!(get!(Set{PkgSpec}, new_dict, spec.name), spec)
    end
    # find changes
    removed = collect(String, setdiff(keys(old_dict), keys(new_dict)))
    added = collect(String, setdiff(keys(new_dict), keys(old_dict)))
    changed = String[k for k in intersect(keys(old_dict), keys(new_dict)) if old_dict[k] != new_dict[k]]
    # don't remove pip, this avoids some flip-flopping when removing pip packages
    filter!(x->x!="pip", removed)
    return (removed, changed, added)
end

function _resolve_pip_diff(old_specs, new_specs)
    # make dicts
    old_dict = Dict{String,PipPkgSpec}(spec.name => spec for spec in old_specs)
    new_dict = Dict{String,PipPkgSpec}(spec.name => spec for spec in new_specs)
    # find changes
    removed = collect(String, setdiff(keys(old_dict), keys(new_dict)))
    added = collect(String, setdiff(keys(new_dict), keys(old_dict)))
    changed = String[k for k in intersect(keys(old_dict), keys(new_dict)) if old_dict[k] != new_dict[k]]
    return (removed, changed, added)
end

function _verbosity()
    parse(Int, get(ENV, "JULIA_CONDAPKG_VERBOSITY", "-1"))
end

function _verbosity_flags()
    n = _verbosity()
    n < 0 ? ["-"*"q"^(-n)] : n > 0 ? ["-"*"v"^n] : String[]
end

function _resolve_conda_remove_all(io, conda_env)
    vrb = _verbosity_flags()
    cmd = conda_cmd(`remove $vrb -y -p $conda_env --all`, io=io)
    flags = append!(["-y", "--all"], vrb)
    _run(io, cmd, "Removing environment", flags=flags)
    nothing
end

function _resolve_conda_install(io, conda_env, specs, channels; create=false)
    args = String[]
    for spec in specs
        push!(args, specstr(spec))
    end
    for channel in channels
        push!(args, "-c", specstr(channel))
    end
    vrb = _verbosity_flags()
    cmd = conda_cmd(`$(create ? "create" : "install") $vrb -y -p $conda_env --override-channels --no-channel-priority $args`, io=io)
    flags = append!(["-y", "--override-channels", "--no-channel-priority"], vrb)
    _run(io, cmd, "Installing packages", flags=flags)
    nothing
end

function _resolve_conda_remove(io, conda_env, pkgs)
    vrb = _verbosity_flags()
    cmd = conda_cmd(`remove $vrb -y -p $conda_env $pkgs`, io=io)
    flags = append!(["-y"], vrb)
    _run(io, cmd, "Removing packages", flags=flags)
    nothing
end

function _resolve_pip_install(io, pip_specs, load_path)
    args = String[]
    for spec in pip_specs
        push!(args, specstr(spec))
    end
    vrb = _verbosity_flags()
    flags = vrb
    old_load_path = STATE.load_path
    try
        STATE.resolved = true
        STATE.load_path = load_path
        withenv() do
            pip = which("pip")
            cmd = `$pip install $vrb $args`
            _run(io, cmd, "Installing Pip packages", flags=flags)
        end
    finally
        STATE.resolved = false
        STATE.load_path = old_load_path
    end
    nothing
end

function _resolve_pip_remove(io, pkgs, load_path)
    vrb = _verbosity_flags()
    flags = append!(["-y"], vrb)
    old_load_path = STATE.load_path
    try
        STATE.resolved = true
        STATE.load_path = load_path
        withenv() do
            pip = which("pip")
            cmd = `$pip uninstall $vrb -y $pkgs`
            _run(io, cmd, "Removing Pip packages", flags=flags)
        end
    finally
        STATE.resolved = false
        STATE.load_path = old_load_path
    end
    nothing
end

function _log(printfunc::Function, io::IO, args...; label="CondaPkg", opts...)
    printstyled(io, lpad(label, 12), " ", color=:green, bold=true)
    printfunc(io, args...; opts...)
    println(io)
    flush(io)
end

function _log(io::IO, args...; opts...)
    _log(printstyled, io, args...; opts...)
end

function _cmdlines(cmd, flags)
    lines = String[]
    isarg = false
    for x in cmd.exec
        if isarg
            lines[end] *= " $x"
            isarg = false
        else
            push!(lines, x)
            isarg = length(lines) > 1 && startswith(x, "-") && x ∉ flags
        end
    end
    lines
end

function _run(io::IO, cmd::Cmd, args...; flags=String[])
    _log(io, args...)
    if _verbosity() ≥ 0
        lines = _cmdlines(cmd, flags)
        for (i, line) in enumerate(lines)
            pre = i==length(lines) ? "└ " : "│ "
            _log(io, label="") do io
                print(io, pre)
                printstyled(io, line, color=:light_black)
            end
        end
    end
    run(cmd)
end

function resolve(; force::Bool=false, io::IO=stderr, interactive::Bool=false, dry_run::Bool=false)
    # if frozen, do nothing
    if STATE.frozen
        return
    end
    # if backend is Null, assume resolved
    if backend() == :Null
        interactive && _log(io, "Using the Null backend, nothing to do")
        STATE.resolved = true
        return
    end
    # skip resolving if already resolved and LOAD_PATH unchanged
    # this is a very fast check which avoids touching the file system
    load_path = Base.load_path()
    if !force && STATE.resolved && STATE.load_path == load_path
        interactive && _log(io, "Dependencies already up to date")
        return
    end
    STATE.resolved = false
    STATE.load_path = load_path
    # find the topmost env in the load_path which depends on CondaPkg
    top_env = _resolve_top_env(load_path)
    STATE.meta_dir = meta_dir = joinpath(top_env, ".CondaPkg")
    meta_file = joinpath(meta_dir, "meta")
    conda_env = joinpath(meta_dir, "env")
    # skip resolving if nothing has changed since the metadata was updated
    if !force && isdir(conda_env) && isfile(meta_file) && _resolve_can_skip_1(load_path, meta_file)
        STATE.resolved = true
        interactive && _log(io, "Dependencies already up to date")
        return
    end
    # find all dependencies
    (packages, channels, pip_packages, extra_path) = _resolve_find_dependencies(io, load_path)
    # install pip if there are pip packages to install
    if !isempty(pip_packages) && !haskey(packages, "pip")
        get!(Dict{String,PkgSpec}, packages, "pip")["<internal>"] = PkgSpec("pip")
        if !any(c.name in ("conda-forge", "anaconda") for c in channels)
            push!(channels, ChannelSpec("conda-forge"))
        end
    end
    # sort channels
    # (in the future we might prioritise them)
    sort!(unique!(channels), by=c->c.name)
    # merge dependencies
    specs = _resolve_merge_packages(packages)
    # merge pip dependencies
    pip_specs = _resolve_merge_pip_packages(pip_packages)
    # find what has changed
    meta = isfile(meta_file) ? open(read_meta, meta_file) : nothing
    if meta === nothing
        removed_pkgs = String[]
        changed_pkgs = String[]
        added_pkgs = unique!(String[x.name for x in specs])
        removed_pip_pkgs = String[]
        changed_pip_pkgs = String[]
        added_pip_pkgs = unique!(String[x.name for x in pip_specs])
    else
        removed_pkgs, changed_pkgs, added_pkgs = _resolve_diff(meta.packages, specs)
        removed_pip_pkgs, changed_pip_pkgs, added_pip_pkgs = _resolve_pip_diff(meta.pip_packages, pip_specs)
    end
    changes = sort([
        (mod1(i,3), i>3 ? "$pkg (pip)" : pkg)
        for (i, pkgs) in enumerate([added_pkgs, changed_pkgs, removed_pkgs, added_pip_pkgs, changed_pip_pkgs, removed_pip_pkgs])
        for pkg in pkgs
    ])
    if !dry_run && !isempty(changes)
        _log(io, "Resolving changes")
        for (i, pkg) in changes
            char = i==1 ? "+" : i==2 ? "~" : "-"
            color = i==1 ? :green : i==2 ? :yellow : :red
            _log(io, char, " ", pkg, label="", color=color)
        end
    end
    # install/uninstall packages
    if !force && meta !== nothing && stat(conda_env).mtime < meta.timestamp && (isdir(conda_env) || (isempty(meta.packages) && isempty(meta.pip_packages)))
        # the state is sufficiently clean that we can modify the existing conda environment
        changed = false
        if !isempty(removed_pip_pkgs)
            dry_run && return
            changed = true
            _resolve_pip_remove(io, removed_pip_pkgs, load_path)
        end
        if !isempty(removed_pkgs)
            dry_run && return
            changed = true
            _resolve_conda_remove(io, conda_env, removed_pkgs)
        end
        if !isempty(specs) && (!isempty(added_pkgs) || !isempty(changed_pkgs) || (meta.channels != channels) || changed)
            dry_run && return
            changed = true
            _resolve_conda_install(io, conda_env, specs, channels)
        end
        if !isempty(pip_specs) && (!isempty(added_pip_pkgs) || !isempty(changed_pip_pkgs) || changed)
            dry_run && return
            changed = true
            _resolve_pip_install(io, pip_specs, load_path)
        end
        if !changed
            _log(io, "Dependencies already up to date")
        end
    else
        # the state is too dirty, recreate the conda environment from scratch
        dry_run && return
        # remove environment
        mkpath(meta_dir)
        if isdir(conda_env)
            _resolve_conda_remove_all(io, conda_env)
        end
        # create conda environment
        _resolve_conda_install(io, conda_env, specs, channels; create=true)
        # install pip packages
        if !isempty(pip_specs)
            _resolve_pip_install(io, pip_specs, load_path)
        end
    end
    # save metadata
    meta = Meta(
        timestamp = time(),
        load_path = load_path,
        extra_path = extra_path,
        version = VERSION,
        packages = specs,
        channels = channels,
        pip_packages = pip_specs,
    )
    open(io->write_meta(io, meta), meta_file, "w")
    # all done
    STATE.resolved = true
    return
end

function is_resolved()
    resolve(io=devnull, dry_run=true)
    STATE.resolved
end


text:
function get_dnl(::Type{UpperBounds}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::jint, arg3::jint)
    return jcall(UpperBounds, "getDnl", RealFieldElement, (RealFieldElement, RealFieldElement, jint, jint), arg0, arg1, arg2, arg3)
end

function get_dnl(::Type{UpperBounds}, arg0::jdouble, arg1::jdouble, arg2::jint, arg3::jint)
    return jcall(UpperBounds, "getDnl", jdouble, (jdouble, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function get_rnml(::Type{UpperBounds}, arg0::RealFieldElement, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint)
    return jcall(UpperBounds, "getRnml", RealFieldElement, (RealFieldElement, jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_rnml(::Type{UpperBounds}, arg0::jdouble, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint)
    return jcall(UpperBounds, "getRnml", jdouble, (jdouble, jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end


text:
# NNlib seems about 30x faster than using Tullio.

@doc raw"""
    conv_v2h(w, v)

Internal function used to compute inputs from a visible configurations `v` to the hidden
layer, where `w` are the convolutional RBM weights.

```math
I_\mu^{k_1,\dots,k_n,b} = \sum_{c,j_1,\dots,j_n} w_{c,j_1,\dots,j_n,\mu} v_{c,j_1+k_1-1,\dots,j_n+k_n-1,b}
```

Assumes that:

* `w` is of size `(C,J₁,...,Jₙ,M)`
* `v` is of size `(C,N₁,...,Nₙ,B)`

Here `C` is the flattened channel dimension,
`M` is the number of hidden units and `B` is the batch size.
Therefore the hidden and batch dimensions must be flattened before calling `conv_v2h`.
In the default case with `stride = 1`, `pad = 0`, `dilation = 1`, ...,
the output `I` is of size `(M, N₁ - J₁ + 1, ..., Nₙ - Jₙ + 1, B)`,

!!! warning
    This is an internal function and is not part of the public API.

!!! warning
    Only works for `n = 1, 2, 3` due to a technical limitations.
"""
function conv_v2h(
    w::AbstractTensor{N}, v::AbstractTensor{N};
    stride = 1, pad = 0, dilation = 1, groups = 1
) where {N}
    @assert size(w, 1) == size(v, 1) # channel size
    @assert all((size(w) .≤ size(v))[(begin + 1):(end - 1)]) # kernel dims
    @assert 3 ≤ N ≤ 5 # NNlib limitation
    w_ = nnlib_conv_permdims(w)
    v_ = nnlib_conv_permdims(v)
    I_ = NNlib.conv(v_, w_; stride, pad, dilation, groups, flipped=true)
    @assert ndims(I_) == N
    return nnlib_conv_invpermdims(I_)
end

@doc raw"""
    conv_h2v(w, v)

Internal function used to compute inputs from a hidden configurations `h` to the visible
layer, where `w` are the convolutional RBM weights.

```math
I_{i_1,\dots,i_n}^{k_1,\dots,k_n,b} = \sum_{c,j_1,\dots,j_n} w_{c,j_1,\dots,j_n,\mu} v_{c,j_1+k_1-1,\dots,j_n+k_n-1,b}
```

Assumes that:

* `w` is of size `(C,J₁,...,Jₙ,M)`
* `h` is of size `(M,K₁,...,Kₙ,B)`

Here `C` is the channel dimension, `M` is the number of hidden units and `B` is the batch size.
These three dimensions are flat, so that the hidden and batch dimensions must be
flattened before calling `conv_h2v`.
In the default case with `stride = 1`, `pad = 0`, `dilation = 1`, ...,
the output `I` is of size `(C, K₁ + J₁ - 1, ..., Kₙ + Jₙ - 1, B)`.

!!! warning
    This is an internal function, not part of the public API. It is subject to change.
"""
function conv_h2v(
    w::AbstractTensor{N}, h::AbstractTensor{N};
    stride = 1, pad = 0, dilation = 1, groups = 1
) where {N}
    @assert size(w)[end] == size(h, 1) # hidden units
    w_ = nnlib_conv_permdims(w)
    h_ = nnlib_conv_permdims(h)
    layer = Flux.ConvTranspose(w_, false; stride, pad, dilation, groups)
    cdims = NNlib.DenseConvDims(Flux.conv_transpose_dims(layer, h_); F=true)
    I_ = NNlib.∇conv_data(h_, w_, cdims)
    return nnlib_conv_invpermdims(I_)
end

function nnlib_conv_permdims(x::AbstractArray{T,N}) where {T,N}
    # NNlib expects dims in this order: spatial (convolved), channels, batch
    return permutedims(x, (ntuple(d -> d + 1, N - 2)..., 1, N))
end

function nnlib_conv_invpermdims(x::AbstractArray{T,N}) where {T,N}
    return permutedims(x, (N - 1, ntuple(identity, N - 2)..., N))
end