text:
@muladd function perform_step!(integrator, cache::SKenCarpConstantCache)
  @unpack t,dt,uprev,u,g,p,f = integrator
  @unpack nlsolver = cache
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,ebtilde1,ebtilde2,ebtilde3,ebtilde4 = cache.tab
  @unpack nb021,nb043 = cache.tab
  alg = unwrap_alg(integrator, true)
  OrdinaryDiffEq.markfirststage!(nlsolver)

  sqrt3 = sqrt(3one(eltype(integrator.W.dW)))
  chi2 = (integrator.W.dW + integrator.W.dZ/sqrt3)/2 #I_(1,0)/h

  if typeof(integrator.f) <: SplitSDEFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  # calculate W
  repeat_step = false

  z₁ = dt*f( uprev,p,t)
  nlsolver.c = 2γ

  g1 = g(uprev,p,t)
  tmp = uprev + γ*z₁ + nb021*chi2.*g1
  if typeof(integrator.f) <: SplitSDEFunction
    # This assumes the implicit part is cheaper than the explicit part
    k1 = dt*f2(uprev,p,t)
    tmp += ea21*k1
  end
  nlsolver.tmp = tmp

  if alg.extrapolant == :min_correct
    z₂ = zero(z₁)
  elseif alg.extrapolant == :trivial
    z₂ = z₁
  end
  nlsolver.z = z₂

  z₂ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  nlsolver.c = c3
  if typeof(integrator.f) <: SplitSDEFunction
    u = tmp + γ*z₂
    k2 = dt*f2(u,p,t + 2γ*dt)
    tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
  else
    tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.tmp = tmp

  if alg.extrapolant == :min_correct
    z₃ = zero(z₂)
  elseif alg.extrapolant == :trivial
    z₃ = z₂
  end
  nlsolver.z = z₃

  z₃ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  nlsolver.c = one(nlsolver.c)

  # Note: Can use g1 since c13 = 0 => g3 == g1

  if typeof(integrator.f) <: SplitSDEFunction
    u = tmp + γ*z₃
    k3 = dt*f2(u,p,t + c3*dt)
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3 + nb043*chi2.*g1
  else
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + nb043*chi2.*g1
  end
  nlsolver.tmp = tmp

  if alg.extrapolant == :min_correct
    z₄ = zero(z₂)
  elseif alg.extrapolant == :trivial
    z₄ = z₂
  end
  nlsolver.z = z₄

  z₄ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  u = tmp + γ*z₄
  g4 = g(uprev,p,t+dt)

  E₂ = chi2.*(g1-g4)

  if typeof(integrator.f) <: SplitSDEFunction
    k4 = dt*f2(u,p,t+dt)
    u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + eb1*k1 + eb2*k2 + eb3*k3 + eb4*k4 + integrator.W.dW.*g4 + E₂
  else
    u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + integrator.W.dW.*g4 + E₂
  end

  ################################### Finalize

  if integrator.opts.adaptive

    if alg.ode_error_est
      if typeof(integrator.f) <: SplitSDEFunction
        tmp = btilde1*z₁  + btilde2*z₂  + btilde3*z₃ + btilde4*z₄ + ebtilde1*k1 + ebtilde2*k2 + ebtilde3*k3 + ebtilde4*k4 + chi2*(g1-g4)
      else
        tmp = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + chi2*(g1-g4)
      end
      if alg.smooth_est # From Shampine
        E₁ = OrdinaryDiffEq._reshape(OrdinaryDiffEq.get_W(nlsolver) \OrdinaryDiffEq._vec(tmp), axes(tmp))
      else
        E₁ = tmp
      end
    else
      E₁ = z₁ + z₂ + z₃ + z₄
    end

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(resids,t)
  end

  integrator.u = u
end

@muladd function perform_step!(integrator, cache::SKenCarpCache)
  @unpack t,dt,uprev,u,g,p,f = integrator
  @unpack z₁,z₂,z₃,z₄,k1,k2,k3,k4,atmp = cache
  @unpack g1,g4,chi2,nlsolver = cache
  @unpack z,tmp = nlsolver
  @unpack k,dz = nlsolver.cache # alias to reduce memory
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4 = cache.tab
  @unpack ebtilde1,ebtilde2,ebtilde3,ebtilde4 = cache.tab
  @unpack nb021,nb043 = cache.tab
  alg = unwrap_alg(integrator, true)
  OrdinaryDiffEq.markfirststage!(nlsolver)

  repeat_step = false

  # Some aliases

  E₁ = g4
  E₂ = dz

  if typeof(integrator.f) <: SplitSDEFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  sqrt3 = sqrt(3one(eltype(integrator.W.dW)))
  if typeof(integrator.W.dW) <: Union{SArray,Number}
    chi2 = (integrator.W.dW + integrator.W.dZ/sqrt3)/2 #I_(1,0)/h
  else
    @.. chi2 = (integrator.W.dW + integrator.W.dZ/sqrt3)/2 #I_(1,0)/h
  end

  # precalculations

  γdt = γ*dt

  if !repeat_step && !integrator.last_stepfail
    f(z₁, integrator.uprev, p, integrator.t)
    z₁ .*= dt
  end

  ##### Step 2

  # TODO: Add a cache so this isn't overwritten near the end, so it can not repeat on fail
  g(g1,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @.. z₄ = chi2*g1 # use z₄ as storage for the g1*chi2
  else
    mul!(z₄,g1,chi2) # use z₄ as storage for the g1*chi2
  end

  @.. tmp = uprev + γ*z₁ + nb021*z₄

  if alg.extrapolant == :min_correct
    @.. z₂ = zero(eltype(dz))
  elseif alg.extrapolant == :trivial
    @.. z₂ = z₁
  end

  if typeof(integrator.f) <: SplitSDEFunction
    # This assumes the implicit part is cheaper than the explicit part
    if !repeat_step && !integrator.last_stepfail
      f2(k1,integrator.uprev,integrator.p,integrator.t)
      k1 .*= dt
    end
    @.. tmp += ea21*k1
  end

  nlsolver.z = z₂
  nlsolver.c = 2γ

  z₂ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  if typeof(integrator.f) <: SplitSDEFunction
    @.. u = tmp + γ*z₂
    f2(k2,u,p,t + 2γ*dt); k2 .*= dt
    @.. tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
  else
    @.. tmp = uprev + a31*z₁ + a32*z₂
  end

  if alg.extrapolant == :min_correct
    @.. z₃ = zero(eltype(dz))
  elseif alg.extrapolant == :trivial
    @.. z₃ = z₂
  end
  nlsolver.z = z₃
  nlsolver.c = c3
  z₃ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitSDEFunction
    @.. u = tmp + γ*z₃
    f2(k3,u,p,t + c3*dt); k3 .*= dt
    # z₄ is storage for the g1*chi2 from earlier
    @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3 + nb043*z₄
  else
    @unpack α41,α42 = cache.tab
    # z₄ is storage for the g1*chi2
    @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + nb043*z₄
  end

  if alg.extrapolant == :min_correct
    @.. z₄ = zero(eltype(dz))
  elseif alg.extrapolant == :trivial
    @.. z₄ = z₂
  end

  nlsolver.z = z₄
  nlsolver.c = one(nlsolver.c)
  z₄ = OrdinaryDiffEq.nlsolve!(nlsolver, integrator, cache, repeat_step)
  OrdinaryDiffEq.nlsolvefail(nlsolver) && return

  g(g4,u,p,t+dt)

  if typeof(integrator.f) <: SplitSDEFunction
    @.. u = tmp + γ*z₄
    f2(k4,u,p,t+dt); k4 .*= dt
    if is_diagonal_noise(integrator.sol.prob)
      @.. E₂ = chi2*(g1-g4)
      @.. u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + eb1*k1 + eb2*k2 + eb3*k3 + eb4*k4 + integrator.W.dW*g4 + E₂
    else
      g1 .-= g4
      mul!(E₂,g1,chi2)
      mul!(tmp,g4,integrator.W.dW)
      @.. u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + eb1*k1 + eb2*k2 + eb3*k3 + eb4*k4 + tmp + E₂
    end
  else
    if is_diagonal_noise(integrator.sol.prob)
      @.. E₂ = chi2*(g1-g4)
      @.. u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + integrator.W.dW*g4 + E₂
    else
      g1 .-= g4
      mul!(E₂,g1,chi2)
      mul!(tmp,g4,integrator.W.dW)
      @.. u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + tmp + E₂
    end
  end

  ################################### Finalize

  if integrator.opts.adaptive

    if alg.ode_error_est
      if typeof(integrator.f) <: SplitSDEFunction
        @.. g1 = btilde1*z₁  + btilde2*z₂  + btilde3*z₃ + btilde4*z₄ + ebtilde1*k1 + ebtilde2*k2 + ebtilde3*k3 + ebtilde4*k4
      else
        @.. g1 = btilde1*z₁  + btilde2*z₂  + btilde3*z₃ + btilde4*z₄
      end
      if alg.smooth_est # From Shampine
        nlsolver.cache.linsolve(vec(E₁),OrdinaryDiffEq.get_W(nlsolver),vec(g1),false)
      else
        E₁ .= dz
      end
    else
      @.. E₁ = z₁ + z₂ + z₃ + z₄
    end

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end


text:
# Julia wrapper for header: /usr/include/expat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Nothing XML_SetElementDeclHandler (XML_Parser, XML_ElementDeclHandler) libexpat
@c Nothing XML_SetAttlistDeclHandler (XML_Parser, XML_AttlistDeclHandler) libexpat
@c Nothing XML_SetXmlDeclHandler (XML_Parser, XML_XmlDeclHandler) libexpat
@c XML_Parser XML_ParserCreate (Ptr{XML_Char},) libexpat
@c XML_Parser XML_ParserCreateNS (Ptr{XML_Char}, XML_Char) libexpat
@c XML_Parser XML_ParserCreate_MM (Ptr{XML_Char}, Ptr{XML_Memory_Handling_Suite}, Ptr{XML_Char}) libexpat
@c XML_Bool XML_ParserReset (XML_Parser, Ptr{XML_Char}) libexpat
@c Nothing XML_SetEntityDeclHandler (XML_Parser, XML_EntityDeclHandler) libexpat
@c Nothing XML_SetElementHandler (XML_Parser, XML_StartElementHandler, XML_EndElementHandler) libexpat
@c Nothing XML_SetStartElementHandler (XML_Parser, XML_StartElementHandler) libexpat
@c Nothing XML_SetEndElementHandler (XML_Parser, XML_EndElementHandler) libexpat
@c Nothing XML_SetCharacterDataHandler (XML_Parser, XML_CharacterDataHandler) libexpat
@c Nothing XML_SetProcessingInstructionHandler (XML_Parser, XML_ProcessingInstructionHandler) libexpat
@c Nothing XML_SetCommentHandler (XML_Parser, XML_CommentHandler) libexpat
@c Nothing XML_SetCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler, XML_EndCdataSectionHandler) libexpat
@c Nothing XML_SetStartCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler) libexpat
@c Nothing XML_SetEndCdataSectionHandler (XML_Parser, XML_EndCdataSectionHandler) libexpat
@c Nothing XML_SetDefaultHandler (XML_Parser, XML_DefaultHandler) libexpat
@c Nothing XML_SetDefaultHandlerExpand (XML_Parser, XML_DefaultHandler) libexpat
@c Nothing XML_SetDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler, XML_EndDoctypeDeclHandler) libexpat
@c Nothing XML_SetStartDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler) libexpat
@c Nothing XML_SetEndDoctypeDeclHandler (XML_Parser, XML_EndDoctypeDeclHandler) libexpat
@c Nothing XML_SetUnparsedEntityDeclHandler (XML_Parser, XML_UnparsedEntityDeclHandler) libexpat
@c Nothing XML_SetNotationDeclHandler (XML_Parser, XML_NotationDeclHandler) libexpat
@c Nothing XML_SetNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler, XML_EndNamespaceDeclHandler) libexpat
@c Nothing XML_SetStartNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler) libexpat
@c Nothing XML_SetEndNamespaceDeclHandler (XML_Parser, XML_EndNamespaceDeclHandler) libexpat
@c Nothing XML_SetNotStandaloneHandler (XML_Parser, XML_NotStandaloneHandler) libexpat
@c Nothing XML_SetExternalEntityRefHandler (XML_Parser, XML_ExternalEntityRefHandler) libexpat
@c Nothing XML_SetExternalEntityRefHandlerArg (XML_Parser, Ptr{Nothing}) libexpat
@c Nothing XML_SetSkippedEntityHandler (XML_Parser, XML_SkippedEntityHandler) libexpat
@c Nothing XML_SetUnknownEncodingHandler (XML_Parser, XML_UnknownEncodingHandler, Ptr{Nothing}) libexpat
@c Nothing XML_DefaultCurrent (XML_Parser,) libexpat
@c Nothing XML_SetReturnNSTriplet (XML_Parser, Int32) libexpat
@c Nothing XML_SetUserData (XML_Parser, Ptr{Nothing}) libexpat
@c Int32 XML_SetEncoding (XML_Parser, Ptr{XML_Char}) libexpat
@c Nothing XML_UseParserAsHandlerArg (XML_Parser,) libexpat
@c Int32 XML_UseForeignDTD (XML_Parser, XML_Bool) libexpat
@c Int32 XML_SetBase (XML_Parser, Ptr{XML_Char}) libexpat
@c Ptr{XML_Char} XML_GetBase (XML_Parser,) libexpat
@c Int32 XML_GetSpecifiedAttributeCount (XML_Parser,) libexpat
@c Int32 XML_GetIdAttributeIndex (XML_Parser,) libexpat
@c Int32 XML_Parse (XML_Parser, Ptr{UInt8}, Int32, Int32) libexpat
@c Ptr{Nothing} XML_GetBuffer (XML_Parser, Int32) libexpat
@c Int32 XML_ParseBuffer (XML_Parser, Int32, Int32) libexpat
@c Int32 XML_StopParser (XML_Parser, XML_Bool) libexpat
@c Int32 XML_ResumeParser (XML_Parser,) libexpat
@c Nothing XML_GetParsingStatus (XML_Parser, Ptr{XML_ParsingStatus}) libexpat
@c XML_Parser XML_ExternalEntityParserCreate (XML_Parser, Ptr{XML_Char}, Ptr{XML_Char}) libexpat
@c Int32 XML_SetParamEntityParsing (XML_Parser, Nothing) libexpat
@c Int32 XML_SetHashSalt (XML_Parser, UInt32) libexpat
@c Int32 XML_GetErrorCode (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentLineNumber (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentColumnNumber (XML_Parser,) libexpat
@c XML_Index XML_GetCurrentByteIndex (XML_Parser,) libexpat
@c Int32 XML_GetCurrentByteCount (XML_Parser,) libexpat
@c Ptr{UInt8} XML_GetInputContext (XML_Parser, Ptr{Int32}, Ptr{Int32}) libexpat
@c Nothing XML_FreeContentModel (XML_Parser, Ptr{XML_Content}) libexpat
@c Ptr{Nothing} XML_MemMalloc (XML_Parser, Csize_t) libexpat
@c Ptr{Nothing} XML_MemRealloc (XML_Parser, Ptr{Nothing}, Csize_t) libexpat
@c Nothing XML_MemFree (XML_Parser, Ptr{Nothing}) libexpat
@c Nothing XML_ParserFree (XML_Parser,) libexpat
@c Ptr{XML_LChar} XML_ErrorString (Cint,) libexpat
@c Ptr{XML_LChar} XML_ExpatVersion () libexpat
@c XML_Expat_Version XML_ExpatVersionInfo () libexpat
@c Ptr{XML_Feature} XML_GetFeatureList () libexpat

ext: jl


text:
########## precompiling section

@assert precompile(Proto.meta, (Type{Proto.Type_Res}, ))
@assert precompile(Proto.meta, (Type{Proto.Type_Req}, ))
@assert precompile(Proto.meta, (Type{Proto.Transaction_Res}, ))
@assert precompile(Proto.meta, (Type{Proto.Transaction_ResPart}, ))

@assert precompile(Proto.which_oneof, (Proto.Transaction_ResPart, Symbol))
@assert precompile(Proto.which_oneof, (Proto.Transaction_Res, Symbol))
@assert precompile(_process_request, (BidirectionalStream, Proto.Transaction_Req, Bool))
@assert precompile(_open_result_channel, (BidirectionalStream, Proto.Transaction_Req, Bool))
@assert precompile(collect_result, (Channel{Union{Proto.Transaction_Res,Proto.Transaction_ResPart}}, BidirectionalStream))
@assert precompile(_is_stream_respart_done, (Proto.Transaction_Res, BidirectionalStream))
@assert precompile(_is_stream_respart_done, (Proto.Transaction_ResPart, BidirectionalStream))


text:
module KRPC
using Sockets
using ProtoBuf
using LightXML
import MacroTools

include("proto/krpc.jl")

abstract type Request{S, P, R} end

"""
KRPC stream listener. See `add_stream` for usage information.
"""
struct Listener{T}
	streams::Dict{UInt64, Pair{Request, Int}}
	connection::T
	current_value
	channel::Channel
end

"""
A connection to Kerbal Space Program.
"""
mutable struct KRPCConnection
    conn::TCPSocket
    stream_conn::TCPSocket
    identifier::Array{UInt8, 1}

    str_listener::Union{Nothing, Task}
    listeners::Dict{UInt64, Listener}
    active::Channel
end

function Base.show(io::IO, conn::KRPCConnection)
	print(io, "KRPC Connection(main fd: $(conn.conn), stream fd: $(conn.stream_conn))")
end

include("raw.jl")
include("connection.jl")
include("types.jl")
include("message.jl")
include("streams.jl")
include("codegen.jl")
if isfile("generated.jl")
	include("generated.jl")
end

struct AddStream_Phantom <: Request{:KRPC, :AddStream, kRPCTypes.kStream}
    call::KRPC.Request
    start::Bool
end
struct RemoveStream_Phantom <: Request{:KRPC, :RemoveStream, Nothing}
    id::UInt64
end
struct GetStatus_Phantom <: Request{:KRPC, :GetStatus, kRPCTypes.Service} end

export kerbal_connect, kerbal, add_stream, generate_stubs

end


text:
#=
kspace_spiral.jl
Jing Dong

Based on mri_kspace_spiral.m that was based on m-files from Valur Olafsson
that he got from Brad Sutton who got them from Doug Noll...
=#

export mri_kspace_spiral

using Interpolations

"""
    kspace, omega, gradxy = mri_kspace_spiral( [options] )

Make k-space spiral trajectory based on GE 3T scanner constraints

Option:
- `N` dimention of reconstructed image
- `Nt` # of time points
- `fov` field of view in cm
- `dt` time sampling interval out; default `5e-6` sec
- `gamp::Real` design gradient amplitude max, G/cm; default 2.2
- `gslew::Int` design slew rate, mT/m/ms; default 180

Out:
- `kspace [Nt,2]` kspace trajectory `[kx ky]` in cycles/cm, NO: cycles/FOV
- `omega [Nt,2]` "" in radians
- `gradxy [Nt 2]` gradient waveforms in (units?)
"""
function mri_kspace_spiral( ; fov::Real = 22,
		N::Int = 64, Nt::Int = -1,
		dt::Float64 = 5e-6, nl::Int = 1,
		gamp::Real = 2.2, gslew::Int = 180,
		warn_nk::Bool = true)

	if Nt == -1
		if fov == 20
			Nt = 4026
		elseif fov == 22
			Nt = 3770
		else
			Nt = 0 # let algorithm choose
		end
	end

	# generate spiral k-space trajectory
	kx, ky, gx, gy = genkspace(fov, N, Nt, nl, gamp, gslew, dt, warn_nk)
	if (nl == 1)
		kspace = [kx ky]
		gradxy = [gx gy]
	else
		kspace = permutedims(cat(dims = 3, kx, ky), [1,3,2])
		gradxy = permutedims(cat(dims = 3, gx, gy), [1,3,2])
	end
	omega = 2π * [kx ky] / N

	maximum(omega) > π && throw("bad spiral")

	return kspace, omega, gradxy
end

"""
    genkspace

Generate the proper length of k-space trajectory.

It linearly interpolates the output of `genspiral` to the correct `length()`
& takes care of the rotations for the interleaves.
* `ld` is the length of the data
* `nint` is the number of interleaves
Brad Sutton; University of Michigan
"""
function genkspace(
    FOV, N, ld, nint, gamp, gslew, tsamp, warn_nk ;
    rotamount::Int = 0,
)
	nk = ld/nint

	flag = (nk == 0) # auto determine number of k-space points
	# just input ld = 0.3

	#dts = 4e-6 # 5e-6 [sec]

	Gx, Gy, kxi, kyi, sx, sy, dts =
		genspi(FOV, N, nl=nint, gamp=gamp, gslew=gslew)

	ik = 0:(length(kxi)-1)
	tk = 0:(dts/tsamp*length(kxi)-1)
	tk = tk * tsamp
	kxt = interp1(ik*dts, kxi, tk)
	kyt = interp1(ik*dts, kyi, tk)

	ig = 0:(length(Gx)-1)
	tg = 0:(dts/tsamp*length(Gx)-1)
	tg = tg * tsamp
	gxt = interp1(ig*dts, Gx, tg)
	gyt = interp1(ig*dts, Gy, tg)

	if flag == 1
		nk = length(kxt) - 2
	end
	nk = Int64.(nk)

	if (nk > length(kxt))
		warn_nk && @warn "reduce nk from $nk to $(length(kxt))"
		nk = min(nk, length(kxt))
	end
	kxo = kxt[1:nk]
	kyo = kyt[1:nk]

	gxo = gxt[1:nk]
	gyo = gyt[1:nk]

	# rotate matrix for proper orientation
	phir = -rotamount*pi/2
	kxop = kxt*cos(phir) - kyt*sin(phir)
	kyop = kyt*cos(phir) + kxt*sin(phir)
	gxop = gxt*cos(phir) - gyt*sin(phir)
	gyop = gyt*cos(phir) + gxt*sin(phir)

	kx = zeros(nk, nint)
	ky = zeros(nk, nint)
	gx = zeros(nk, nint)
	gy = zeros(nk, nint)

	if (length(kxop) > length(nk))
		kx = zeros(length(kxop), nint)
		ky = zeros(length(kyop), nint)
		gx = zeros(length(gxop), nint)
		gy = zeros(length(gyop), nint)
	end

	kx[:,1] = kxop
	ky[:,1] = kyop
	gx[:,1] = gxop
	gy[:,1] = gyop
	phi = 2*π/nint
	if nint > 1
		phi = 2*π/nint
		for ii = 1:(nint-1)
			kx[:,ii+1] = kxop*cos(ii*phi) - kyop*sin(ii*phi)
			ky[:,ii+1] = kyop*cos(ii*phi) + kxop*sin(ii*phi)
			gx[:,ii+1] = gxop*cos(ii*phi) - gyop*sin(ii*phi)
			gy[:,ii+1] = gyop*cos(ii*phi) + gxop*sin(ii*phi)
		end
	end

	return kx, ky, gx, gy
end

"""
	Gx, Gy, kx, ky, sx, sy, gts = genspi(...)

This is translation of C code from scanner: exactly what is played out
to gradients at 4us.

multi-shot spiral design
uses Duyn's approximate slewrate limited design
augmented with archimedian `gmax` limit
in [args]
* `D` = FOV; cm
* `N` = matrix size()

* `Tmax` = longest acquisition allowed; s
* `dts` = output sample spacing; s
* `gtype` = trajectory type()

option [CVs]
* `nl` = number of interleaves
* `gamp` = design grad max; G/cm
* `gslew` = design slew rate; mT/m/ms
* `nramp` = number of rampdown points; default 0

out
* `Gx; Gy`

time is in sec()
* rev 0 12/26/98	original
* rev 1 4/15/99	little better calc of ts

Borrowed from Doug Noll; Univ. of Michigan.
Modified to take more input cv's.
"""
function genspi(D, N ;
	nl::Int = 1, gamp::Real = 202, gslew::Int = 180, nramp::Int = 0)

	########## Predefined variables
	GRESMAX = 21000
	# nramp=100
	MAX_PG_WAMP = 32766
	gts = 4e-6 # [sec]

	Tmax = GRESMAX * gts
	dts = gts
	opfov = D

	#################################
	gamma = 2*π*4.257e3
	gambar = gamma / (2*π)

	gx = zeros(2*GRESMAX)
	gy = zeros(2*GRESMAX)

	q = 5
	S0 = gslew*100
	dt = dts*.5

	# slew-rate limited approximation

	Ts = 0.666667 / nl*sqrt(((pi*N)^3)/(gamma*D*S0))

	a2 = N*π/(nl*(Ts^(0.666667)))
	a1 = 1.5*S0/a2
	beta = S0*gamma*D/nl
	Gmax = a1*(Ts^0.333333)
	gmax = 0

	t = 0:dt:Ts
	x = t .^ 1.333333
	theta = (t.^2) .* (.5 * beta ./ (q .+ .5 * beta ./ a2 .* x))
	y = q .+ .5 .* beta./a2 .* x
	dthdt = t .* (beta .* (q .+ 0.166667 * beta ./ a2 .* x) ./ (y .* y))
	c = cos.(theta)
	s = sin.(theta)
	gx = (nl/(D * gamma)) .* dthdt .* (c - theta .* s)
	gy = (nl/(D * gamma)) .* dthdt .* (s + theta .* c)
	gabs = abs.(gx + 1im*gy)
	# cut short if over peak
	gmax = abs.(gamp ./ (theta .+ eps()).+ 1im * gamp)
	l1 = length(t) - sum(gabs .> gmax)
	ts = t[l1]
	thetas = theta[l1]

	# gmax limited approximation
	l3 = 0
	T=ts
	if Gmax > gamp
		T = ((pi*N/nl)*(pi*N/nl) - thetas*thetas)/(2*gamma*gamp*D/nl) + ts
		t = (ts+dt):dt:T
		theta = sqrt.(thetas * thetas .+ (2 * gamma * gamp * D) .* (t .- ts) / nl)
		c = cos.(theta)
		s = sin.(theta)
		ind2 = l1 .+ (1:length(t))
		if (maximum(ind2) > length(gx))
			gx = [gx; zeros(maximum(ind2) - length(gx))]
		end
		if (maximum(ind2) > length(gy))
			gy = [gy; zeros(maximum(ind2) - length(gy))]
		end
		gx[floor.(Int,ind2)] = gamp.*(c./theta - s)
		gy[floor.(Int,ind2)] = gamp.*(s./theta + c)
		l3 = length(t)
	end

	l2 = l1 + l3
	Gx = gx[1:2:l2] # | gx[1:2:l2]*MAX_PG_WAMP/gamp
	Gy = gy[1:2:l2] # | gy[1:2:l2]*MAX_PG_WAMP/gamp
	g = Gx + 1im *Gy # slew rate vector
	s = diff(g)./(gts*1000)  # grad vector
	Kx = cumsum(Gx)*gts*opfov*gambar
	Ky = cumsum(Gy)*gts*opfov*gambar
	k = Kx + 1im*Ky  # kspace vector
	t = collect(0:gts:T) # time vector
	matrix = maximum(abs.(k))*2
	maxg = maximum(abs.(g))
	maxs = maximum(abs.(s))
	maxt = maximum(t)*1000

	kx = real(k)
	ky = imag(k)
	sx = real(s)
	sy = imag(s)

	return Gx, Gy, kx, ky, sx, sy, gts
end


text:
"""
This file initializes the IscaTools module by defining the directories relevant to the particular Isca experiment being resorted and analysed by IscaTools.
"""

function irootdict()

    idict = Dict{AbstractString,AbstractString}()
    idict["root"] = ""; idict["raw"] = ""; idict["ana"] = "";
    idict["experiment"] = ""; idict["configuration"] = "";
    idict["spinup"] = ""; idict["control"] = "";
    idict["fname"] = "";

    return idict

end

function iscaspin(iroot::AbstractDict)

    efol = joinpath(iroot["root"],iroot["experiment"]);
    sfol = joinpath(iroot["root"],iroot["experiment"],"spinup");
    if isdir(sfol)
        @info "$(Dates.now()) - A spinup configuration folder has been identified in $(efol)."
        return true;
    else
        @info "$(Dates.now()) - No spinup configuration folder was identified in $(efol)."
        return false;
    end

end

function iscawelcome()

    ftext = joinpath(@__DIR__,"../extra/welcome.txt");
    lines = readlines(ftext); count = 0; nl = length(lines);
    for l in lines; count += 1;
       if any(count .== [1,2]); print(Crayon(bold=true),"$l\n");
       elseif count == nl;      print(Crayon(bold=false),"$l\n\n");
       else;                    print(Crayon(bold=false),"$l\n");
       end
    end

end

function iscaroot(
    experiment::AbstractString,
    config::AbstractString,
    prjpath::AbstractString
)

    iroot = irootdict(); iroot["root"] = prjpath;
    iroot["raw"] = joinpath(prjpath,experiment,config)

    if !isdir(iroot["raw"])
        error("$(Dates.now()) - The folder $(iroot["raw"]) does not exist.  Please ensure that you have entered the correct project PATH, EXPERIMENT (if applicable), and CONFIGURATION details.")
    end

    iroot["ana"] = joinpath(prjpath,experiment,config,"analysis")
    iroot["experiment"] = experiment; iroot["configuration"] = config;

    @info """$(Dates.now()) - $(BOLD("PROJECT DETAILS:"))
      $(BOLD("Project Directory:")) $prjpath
      $(BOLD("Raw Data Directory:")) $(iroot["raw"])
      $(BOLD("Analysis Directory:")) $(iroot["ana"])
      $(BOLD("Experiment | Configuration:")) $experiment | $config
    """

    if iscaspin(iroot)
        iroot["spinup"]  = replace(iroot["raw"],config=>"spinup")
        iroot["control"] = replace(iroot["raw"],config=>"control")
        @info """$(Dates.now()) - $(BOLD("SPINUP DIRECTORIES:"))
          $(BOLD("Spinup Directory:"))  $(iroot["spinup"])
          $(BOLD("Control Directory:")) $(iroot["control"])
        """
    end

    return iroot

end

function iscaroot(
    experiment::AbstractString,
    config::AbstractString,
    rawpath::AbstractString,
    anapath::AbstractString
)

    iroot = irootdict();
    iroot["raw"] = joinpath(rawpath,experiment,config)

    if !isdir(iroot["raw"])
        error("$(Dates.now()) - The folder $(iroot["raw"]) does not exist.  Please ensure that you have entered the correct project PATH, EXPERIMENT (if applicable), and CONFIGURATION details.")
    end

    iroot["ana"] = joinpath(anapath,experiment,config)
    if !isdir(iroot["ana"])
        @info "$(Dates.now()) - The folder $(iroot["ana"]) does not exist.  Creating now ..."
    end

    iroot["experiment"] = experiment; iroot["configuration"] = config;

    @info """$(Dates.now()) - $(BOLD("PROJECT DETAILS:"))
      $(BOLD("Raw Data Directory:")) $(iroot["raw"])
      $(BOLD("Analysis Directory:")) $(iroot["ana"])
      $(BOLD("Experiment | Configuration:")) $experiment | $config
    """

    if iscaspin(iroot)
        iroot["spinup"]  = replace(iroot["raw"],config=>"spinup")
        iroot["control"] = replace(iroot["raw"],config=>"control")
        @info """$(Dates.now()) - $(BOLD("SPINUP DIRECTORIES:"))
          $(BOLD("Spinup Directory:"))  $(iroot["spinup"])
          $(BOLD("Control Directory:")) $(iroot["control"])
        """
    end

    return iroot

end

function iscaroot(;
    prjpath::AbstractString="",
    rawpath::AbstractString="",
    anapath::AbstractString="",
    experiment::AbstractString="",
    config::AbstractString
)

    if prjpath == "" && !(rawpath == "" && anapath == "")
        prjpath = ENV["GFDL_DATA"];
    end

    if !(rawpath == "" && anapath == "")
          return iscaroot(experiment,config,rawpath,anapath)
    else; return iscaroot(experiment,config,prjpath)
    end

end

function iscastartup(;
    prjpath::AbstractString="",
    rawpath::AbstractString="",
    anapath::AbstractString="",
    experiment::AbstractString="",
    config::AbstractString,
    fname::AbstractString,
    welcome::Bool=true
)

    if welcome; iscawelcome() end
    iroot = iscaroot(prjpath=prjpath,rawpath=rawpath,anapath=anapath,
                     experiment=experiment,config=config)

    fnc = joinpath(iroot["raw"],"run0001","$fname.nc");
    if !isfile(fnc)
        error("$(Dates.now()) - The output file \"$fname.nc\" does not exist.  Please double-check the filename in which the Isca raw data was save into and try again.")
    else
        iroot["fname"] = "$fname.nc";
    end

    if isfile("$(iroot["raw"])/init.jld2")
        @load "$(iroot["raw"])/init.jld2" init
    else
        init = retrievetime(fnc); retrieveruns!(init,iroot)
        ds = Dataset(fnc);
        init["phalf"] = ds["phalf"][:]*100;
        init["pfull"] = ds["pfull"][:]*100;
        init["sealp"] = ds["phalf"][end]*100;;
        init["lon"]   = ds["lon"][:]*1;
        init["lat"]   = ds["lat"][:]*1;
        close(ds);
    end

    return init,iroot

end


text:
using Decimals
using Base.Test

tests = ["constructor",
         "decimal",
         "norm",
         "arithmetic",
         "equals",
         "round"]

println("Running tests:")

for t in tests
    tfile = string("test_", t, ".jl")
    println(" * $(t) ...")
    include(tfile)
end


text:
using BenchmarkTools
using AugmentedGaussianProcesses
using LinearAlgebra, CSV
using Random
const AGP = AugmentedGaussianProcesses

## Benchmark variational parameters update

## Benchmark hyperparameter optimization

## Benchmark full update

## Benchmark fit performance

## Benchmark

compat = Dict{String,Dict{String,Vector{String}}}()
likelihoodnames = ["Gaussian", "StudentT", "Logistic", "BayesianSVM", "LogisticSoftMax"]
inferencenames = ["AnalyticVI", "AnalyticSVI"]
modelnames = ["GP", "VGP", "SVGP"]
funcs = [
    "init",
    "elbo",
    "computematrices",
    "updatevariational",
    "updatehyperparam",
    "predic",
    "predicproba",
]
compat["GP"] = Dict{String,Vector{String}}("Gaussian" => ["AnalyticVI"])
compat["VGP"] = Dict{String,Vector{String}}()
compat["VGP"]["StudentT"] = ["AnalyticVI"]
compat["VGP"]["Logistic"] = ["AnalyticVI"]
compat["VGP"]["BayesianSVM"] = ["AnalyticVI"]
compat["VGP"]["LogisticSoftMax"] = ["AnalyticVI"]
compat["SVGP"] = Dict{String,Vector{String}}()
compat["SVGP"]["Gaussian"] = ["AnalyticVI", "AnalyticSVI"]
compat["SVGP"]["StudentT"] = ["AnalyticVI", "AnalyticSVI"]
compat["SVGP"]["Logistic"] = ["AnalyticVI", "AnalyticSVI"]
compat["SVGP"]["BayesianSVM"] = ["AnalyticVI", "AnalyticSVI"]
compat["SVGP"]["LogisticSoftMax"] = ["AnalyticVI", "AnalyticSVI"]

const SUITE = BenchmarkGroup(["Models"])
Random.seed!(1234)
D = 20;
N = 3000;
data = CSV.read("benchmarkdata.csv")
X = Matrix(data[:, 1:D])
y_key = Dict(
    "Gaussian" => :y_reg,
    "StudentT" => :y_reg,
    "BayesianSVM" => :y_class,
    "Logistic" => :y_class,
    "LogisticSoftMax" => :y_multi,
)
n_ind = 50;
batchsize = 50;
ν = 5.0;
function convertl(lname::String)
    return lname *
           (lname != "BayesianSVM" ? "Likelihood" : "") *
           "(" *
           (lname == "StudentT" ? "ν" : "") *
           ")"
end
converti(iname::String) = iname * "(" * (iname == "AnalyticSVI" ? "batchsize" : "") * ")"
add_ind(mname::String) = mname == "SVGP" ? ",n_ind" : ""
kernel = RBFKernel([2.0]; variance=1.0, dim=D)
models = Dict{String,Dict{String,Dict{String,AbstractGP}}}()
SUITE["Models"] = BenchmarkGroup(modelnames)
for model in String.(keys(compat))
    SUITE["Models"][model] = BenchmarkGroup(String.(keys(compat[model])))
    models[model] = Dict{String,Dict{String,AbstractGP}}()
    for likelihood in String.(keys(compat[model]))
        SUITE["Models"][model][likelihood] = BenchmarkGroup(compat[model][likelihood])
        models[model][likelihood] = Dict{String,AbstractGP}()
        for i in compat[model][likelihood]
            SUITE["Models"][model][likelihood][i] = BenchmarkGroup(funcs)
            if model == "GP"
                models[model][likelihood][i] = eval(
                    Meta.parse(
                        model *
                        "(X,Vector(data[:$((y_key[likelihood]))]),kernel,atfrequency=1)",
                    ),
                )
                SUITE["Models"][model][likelihood][i]["init"] = eval(
                    Meta.parse(
                        "@benchmarkable $model(\$X,y_train,\$kernel,atfrequency=1) setup=(y_train = Vector(\$D[:\$((y_key[likelihood]))])",
                    ),
                )
            else
                # println(Meta.parse(model*"(X,y[\"$likelihood\"],kernel,$(convertl(likelihood)) ,$(converti(i))$(add_ind(model)),atfrequency=1)"))
                models[model][likelihood][i] = eval(
                    Meta.parse(
                        model *
                        "(X,Vector(data[:$((y_key[likelihood]))]),kernel,$(convertl(likelihood)) ,$(converti(i))$(add_ind(model)),atfrequency=1)",
                    ),
                )
                SUITE["Models"][model][likelihood][i]["init"] = eval(
                    Meta.parse(
                        "@benchmarkable $model(\$X,y_train,\$kernel,$(convertl(likelihood)),$(converti(i)) $(add_ind(model)),atfrequency=1) setup=(y_train = Vector(\$data[:\$((y_key[likelihood]))]))",
                    ),
                )
            end
            SUITE["Models"][model][likelihood][i]["elbo"] = @benchmarkable ELBO(gpmodel) setup = (
                gpmodel = deepcopy($(models[model][likelihood][i]))
            )
            SUITE["Models"][model][likelihood][i]["computematrices"] = @benchmarkable AGP.computeMatrices!(
                gpmodel
            ) setup = (gpmodel = deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["updatevariational"] = @benchmarkable AGP.variational_updates!(
                gpmodel
            ) setup = (gpmodel = deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["updatehyperparam"] = @benchmarkable AGP.update_hyperparameters!(
                gpmodel
            ) setup = (gpmodel = deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["predic"] = @benchmarkable predict_y(
                gpmodel, $X
            ) setup = (gpmodel = deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["predicproba"] = @benchmarkable proba_y(
                gpmodel, $X
            ) setup = (gpmodel = deepcopy($(models[model][likelihood][i])))
        end
    end
end


text:
"""A scirntific module (calculate exp(x))"""
module SCI

abstract type AbstractProblem end
struct Problem{T} <: AbstractProblem x::T end

abstract type AbstractAlgorithm end
struct Builtin <: AbstractAlgorithm end
Base.@kwdef struct Taylor <: AbstractAlgorithm n::Int = 10 end
default_algorithm(prob::Problem) = Builtin()

struct Solution{R, P<:AbstractProblem, A<:AbstractAlgorithm} result::R; prob::P; alg::A end
solve(prob::AbstractProblem) = solve(prob, default_algorithm(prob))
solve(prob::AbstractProblem, alg::Builtin) = Solution(exp(prob.x), prob, alg)
solve(prob::AbstractProblem, alg::Taylor) = Solution(sum(prob.x^k/factorial(k) for k in 0:alg.n), prob, alg)

using COM

COM.prettify(sol::Solution{R, P, A}) where {R, P<:AbstractProblem, A<:Builtin} = """
Problem:   x = $(sol.prob.x)
Algorithm: builtin exp(x)
Result:    $(sol.result)
"""

COM.prettify(sol::Solution{R, P, A}) where {R, P<:AbstractProblem, A<:Taylor} = """
Problem:   x = $(sol.prob.x)
Algorithm: Taylor series of exp(x) upto degree $(sol.alg.n)
Result:    $(sol.result)
"""

end

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
mutable struct ScanningStrategy{T<:AbstractFloat, I<:Int, AA<:AbstractArray{T}, AS<:AbstractString}
    nside::I
    duration::I
    sampling_rate::I
    alpha::T
    beta::T
    prec_rpm::T
    spin_rpm::T
    hwp_rpm::T
    FP_theta::AA
    FP_phi::AA
    start_point::AS
    start_angle::T
end

rpm2angfreq(rpm) = (2.0π / 60.0) * rpm
period2rpm(period) = 1.0 / period

"""
    gen_ScanningStrategy(args***)

This function generate scanning strategy.
    # Arguments
    ...
    nside::Int
    duration::Int
    sampling_rate::Int
    alpha::Float
    beta::Float
    prec_rpm::Float
    spin_rpm::Float
    hwp_rpm::Float
    FP_theta::Array
    FP_phi::Array
    start_point::String
    start_angle::Float
    ...

    # Returns
    ...
    - `scanning_strategy_structure`::ScanningStrategy
    ...
"""
function gen_ScanningStrategy(;
        nside=128, duration=60*60*24*365, sampling_rate=1, alpha=45, beta=50,
        prec_rpm=period2rpm(192.348), spin_rpm=0.05, hwp_rpm=0,
        FP_theta=[0.0], FP_phi=[0.0], start_point="equator", start_angle=0.0)
    scanning_strategy_structure = ScanningStrategy(
        nside,
        duration,
        sampling_rate,
        Float64(alpha),
        Float64(beta),
        Float64(prec_rpm),
        Float64(spin_rpm),
        Float64(hwp_rpm),
        Float64.(FP_theta),
        Float64.(FP_phi),
        start_point,
        start_angle
    )
    return scanning_strategy_structure
end

@inline function initial_pointings_rot(ss::ScanningStrategy)
    if ss.start_point =="pole"
        #= Set the initial position of the boresight to near the zenith. =#
        initial_vec = @views @SVector [sind(ss.alpha-ss.beta), 0, cosd(ss.alpha-ss.beta)]
    end
    if ss.start_point == "equator"
        #= Set the initial position of the boresight to near the equator. =#
        initial_vec = @views @SVector [cosd(ss.alpha-ss.beta), 0 , sind(ss.alpha-ss.beta)]
    end
    z_axis = [0, 0, 1]
    bore_quat = Quaternion(0.0, initial_vec)

    q = quaternion_rotator(ss.start_angle, 1.0, z_axis)
    bore_quat_rot = q * bore_quat / q

    return @views @SVector [bore_quat_rot.q1, bore_quat_rot.q2, bore_quat_rot.q3]
end

@inline function initial_pointings(ss::ScanningStrategy)
    if ss.start_point == "pole"
        #= Set the initial position of the boresight to near the zenith. =#
        initial_vec = @views @SVector [sind(ss.alpha-ss.beta), 0, cosd(ss.alpha-ss.beta)]
    end
    if ss.start_point == "equator"
        #= Set the initial position of the boresight to near the equator. =#
        initial_vec = @views @SVector [cosd(ss.alpha-ss.beta), 0 , sind(ss.alpha-ss.beta)]
    end
    return initial_vec
end

"""
    get_pointings_tuple(SS::ScanningStrategy, start, stop)

This function will return pointing tod as tuple.
    # Arguments
    ...
    - `SS::ScanningStrategy`: the ScanningStrategy struct.
    - `start::Int`: the initial time for pointing calculation.
    - `stop::Int`: the finish time for pointing calculation.
    ...

    # Examples
    ```jldoctest
    julia> ss = gen_ScanningStrategy()
    julia> pointings = get_pointings_tuple(s, 0, 100)

    # Returns
    ...
    - `pointings[1]`: TOD of theta. shape:(duration*sampling_rate, numOfdet)
    - `pointings[2]`: TOD of phi. shape:(duration*sampling_rate, numOfdet)
    - `pointings[3]`: TOD of psi. shape:(duration*sampling_rate, numOfdet)
    - `pointings[4]`: Array of time. shape:(duration*sampling_rate)
    ...
"""
@inline function get_pointings_tuple(SS::ScanningStrategy, start, stop)
    resol = Resolution(SS.nside)
    omega_spin = rpm2angfreq(SS.spin_rpm)
    omega_prec = rpm2angfreq(SS.prec_rpm)
    omega_revol = (2π) / (60.0 * 60.0 * 24.0 * 365.25)

    #time_array = @views Vector(start:1/SS.sampling_rate:stop-1/SS.sampling_rate)
    time_array = start:1/SS.sampling_rate:stop-1/SS.sampling_rate
    loop_times = length(time_array)
    numof_det = length(SS.FP_theta)

    psi_tod = zeros(loop_times, numof_det)
    theta_tod = zeros(loop_times, numof_det)
    phi_tod = zeros(loop_times, numof_det)

    antisun_axis = @SVector [1.0, 0.0, 0.0]
    y_axis = @SVector [0.0, 1.0, 0.0]
    z_axis = @SVector [0.0, 0.0, 1.0]

    antisun_axis = vector_rotator(antisun_axis, SS.start_angle, z_axis)
    spin_axis = @SVector [cosd(SS.alpha), 0.0, sind(SS.alpha)]
    spin_axis = vector_rotator(spin_axis, SS.start_angle, z_axis)

    #bore_0 = initial_pointings(SS)
    bore_0 = initial_pointings_rot(SS)

    boresight_0 = Quaternion(0.0, bore_0)
    #detector_vec_0 = @SVector [-boresight_0.q3, 0.0, boresight_0.q1]
    detector_vec_0 = @SVector [-boresight_0.q3, boresight_0.q2, boresight_0.q1]
    detector_orientation_0 = Quaternion(0.0, detector_vec_0)
    travel_direction_vec_0 = bore_0 × detector_vec_0
    travel_direction_0 = Quaternion(0.0, travel_direction_vec_0)

    @views @inbounds @simd for j = eachindex(SS.FP_theta)
        q_theta_in_FP = quaternion_rotator(deg2rad(SS.FP_theta[j]), 1, y_axis)
        q_phi_in_FP = quaternion_rotator(deg2rad(SS.FP_phi[j]), 1, bore_0)
        q_for_FP = q_phi_in_FP * q_theta_in_FP
        pointing = q_for_FP * boresight_0 / q_for_FP

        @views @inbounds @threads for i = eachindex(time_array)
            t = time_array[i]
            q_revol = quaternion_rotator(omega_revol, t, z_axis)
            q_prec = quaternion_rotator(omega_prec, t, antisun_axis)
            q_spin = quaternion_rotator(omega_spin, t, spin_axis)
            Q = q_revol * q_prec * q_spin
            P = Q * pointing / Q
            q_det = Q * detector_orientation_0 / Q
            travel_direction = Q * travel_direction_0 / Q
            pointing_t = @SVector [P.q1, P.q2, P.q3]
            travel_direction_vec = @SVector [travel_direction.q1, travel_direction.q2, travel_direction.q3]

            longitude = pointing_t × (pointing_t × z_axis)
            ang_t = vec2ang_ver2(pointing_t[1], pointing_t[2], pointing_t[3])

            theta_tod[i, j] = ang_t[1]
            phi_tod[i, j] = ang_t[2]

            divergent_vec = longitude × travel_direction_vec
            cosK = dot(travel_direction_vec, longitude) / (norm(travel_direction_vec) * norm(longitude))

            cosK = ifelse(abs(cosK) > 1.0, sign(cosK), cosK)

            divergent_vec3 = sign(divergent_vec[3])
            divergent_vec3 = ifelse(divergent_vec3==0, -1, divergent_vec3)
            pointing_t3 = sign(pointing_t[3])
            pointing_t3 = ifelse(pointing_t3==0, 1, pointing_t3)

            psi_tod[i, j] = acos(cosK) * divergent_vec3 * pointing_t3
        end
    end
    return (theta_tod, phi_tod, psi_tod, time_array)
end

@inline function get_pointings_xyz_tuple(SS::ScanningStrategy, start, stop)
    resol = Resolution(SS.nside)
    omega_spin = rpm2angfreq(SS.spin_rpm)
    omega_prec = rpm2angfreq(SS.prec_rpm)
    omega_revol = (2π) / (60.0 * 60.0 * 24.0 * 365.25)

    #time_array = @views Vector(start:1/SS.sampling_rate:stop-1/SS.sampling_rate)
    time_array = start:1/SS.sampling_rate:stop-1/SS.sampling_rate
    loop_times = length(time_array)
    numof_det = length(SS.FP_theta)

    psi_tod = zeros(loop_times, numof_det)
    #theta_tod = zeros(loop_times, numof_det)
    #phi_tod = zeros(loop_times, numof_det)
    xyz_tod = zeros(3, loop_times, numof_det)

    antisun_axis = @SVector [1.0, 0.0, 0.0]
    y_axis = @SVector [0.0, 1.0, 0.0]
    z_axis = @SVector [0.0, 0.0, 1.0]

    antisun_axis = vector_rotator(antisun_axis, SS.start_angle, z_axis)
    spin_axis = @SVector [cosd(SS.alpha), 0.0, sind(SS.alpha)]
    spin_axis = vector_rotator(spin_axis, SS.start_angle, z_axis)

    #bore_0 = initial_pointings(SS)
    bore_0 = initial_pointings_rot(SS)

    boresight_0 = Quaternion(0.0, bore_0)
    #detector_vec_0 = @SVector [-boresight_0.q3, 0.0, boresight_0.q1]
    detector_vec_0 = @SVector [-boresight_0.q3, boresight_0.q2, boresight_0.q1]
    detector_orientation_0 = Quaternion(0.0, detector_vec_0)
    travel_direction_vec_0 = bore_0 × detector_vec_0
    travel_direction_0 = Quaternion(0.0, travel_direction_vec_0)

    @views @inbounds @simd for j = eachindex(SS.FP_theta)
        q_theta_in_FP = quaternion_rotator(deg2rad(SS.FP_theta[j]), 1, y_axis)
        q_phi_in_FP = quaternion_rotator(deg2rad(SS.FP_phi[j]), 1, bore_0)
        q_for_FP = q_phi_in_FP * q_theta_in_FP
        pointing = q_for_FP * boresight_0 / q_for_FP

        @views @inbounds @threads for i = eachindex(time_array)
            t = time_array[i]
            q_revol = quaternion_rotator(omega_revol, t, z_axis)
            q_prec = quaternion_rotator(omega_prec, t, antisun_axis)
            q_spin = quaternion_rotator(omega_spin, t, spin_axis)
            Q = q_revol * q_prec * q_spin
            P = Q * pointing / Q
            q_det = Q * detector_orientation_0 / Q
            travel_direction = Q * travel_direction_0 / Q
            pointing_t = @SVector [P.q1, P.q2, P.q3]
            travel_direction_vec = @SVector [travel_direction.q1, travel_direction.q2, travel_direction.q3]

            longitude = pointing_t × (pointing_t × z_axis)
            ang_t = vec2ang_ver2(pointing_t[1], pointing_t[2], pointing_t[3])

            #theta_tod[i, j] = ang_t[1]
            #phi_tod[i, j] = ang_t[2]
            xyz_tod[:, i, j] = pointing_t

            divergent_vec = longitude × travel_direction_vec
            cosK = dot(travel_direction_vec, longitude) / (norm(travel_direction_vec) * norm(longitude))

            cosK = ifelse(abs(cosK) > 1.0, sign(cosK), cosK)

            divergent_vec3 = sign(divergent_vec[3])
            divergent_vec3 = ifelse(divergent_vec3==0, -1, divergent_vec3)
            pointing_t3 = sign(pointing_t[3])
            pointing_t3 = ifelse(pointing_t3==0, 1, pointing_t3)

            psi_tod[i, j] = acos(cosK) * divergent_vec3 * pointing_t3
        end
    end
    return (xyz_tod, psi_tod, time_array)
end

"""
    get_pointing_pixels(SS::ScanningStrategy, start, stop)

This function will return pointing pixel tod as tuple.
    # Arguments
    ...
    - `SS::ScanningStrategy`: the ScanningStrategy struct.
    - `start::Int`: the initial time for pointing calculation.
    - `stop::Int`: the finish time for pointing calculation.
    ...

    # Examples
    ```jldoctest
    julia> ss = gen_ScanningStrategy()
    julia> pointings = get_pointing_pixels(s, 0, 100)

    # Returns
    ...
    - `pointings[1]`: TOD of pixel indicies. shape:(duration*sampling_rate, numOfdet)
    - `pointings[2]`: TOD of psi. shape:(duration*sampling_rate, numOfdet)
    - `pointings[3]`: Array of time. shape:(duration*sampling_rate)
    ...
"""
@inline function get_pointing_pixels(SS::ScanningStrategy, start, stop)
    theta_tod, phi_tod, psi_tod, time_array = get_pointings_tuple(SS, start, stop)
    loop_times = length(theta_tod[:,1])
    numOfdet = length(theta_tod[1,:])
    resol = Resolution(SS.nside)

    pix_tod = @views zeros(Int64, loop_times, numOfdet)
    for j in eachindex(theta_tod[1,:])
        for i in eachindex(theta_tod[:,1])
            pix_tod[i, j] = @views ang2pixRing(resol, theta_tod[i, j], phi_tod[i, j])
        end
    end
    return (pix_tod, psi_tod, time_array)
end

"""
    get_pointings(SS::ScanningStrategy, start, stop)

This function will return pointing tod as dictionary type.
    # Arguments
    ...
    - `SS::ScanningStrategy`: the ScanningStrategy struct.
    - `start::Int`: the initial time for pointing calculation.
    - `stop::Int`: the finish time for pointing calculation.
    ...

    # Examples
    ```jldoctest
    julia> ss = gen_ScanningStrategy()
    julia> pointings = get_pointings(s, 0, 100)

    # Returns
    ...
    - `pointings["theta"]`: TOD of theta. shape:(duration*sampling_rate, numOfdet)
    - `pointings["phi"]`: TOD of phi. shape:(duration*sampling_rate, numOfdet)
    - `pointings["psi"]`: TOD of psi. shape:(duration*sampling_rate, numOfdet)
    - `pointings["time"]`: Array of time. shape:(duration*sampling_rate)
    ...
"""
function get_pointings(SS::ScanningStrategy, start, stop)
    pointings = get_pointings_tuple(SS, start, stop)
    pointings = Dict{String,AbstractArray{Float64}}(
        "theta" => pointings[1],
        "phi" => pointings[2],
        "psi" => pointings[3],
        "time" => pointings[4]
    )
    return pointings
end


text:
# If you don't do the branch at j == 9 or so, you will totally blow up
# compilation time. Really trying to keep that under control while still
# squeezing out all allocations.

open("uk_polys_hardcoded.jl", "w") do out
  uk_polys = BesselK.Uk_polynomials(20)
  redirect_stdout(out) do
    for (j,pj) in enumerate(uk_polys)
      c = pj.coeffs
      if j <= 10
        str = string("uk_", j-1, "_eval(x)= @evalpoly x", map(x->string(" ", x, " "), c)...)
        println(str)
      else
        stem = string("uk_", j-1)
        str1 = string("const ", stem, "_coef=[", map(x->string(x, ", "), c)..., "]")
        println(str1)
        cnm  = string(stem, "_coef")
        str2 = string(stem, "_eval(x) = evalpoly(x, $cnm)")
        println(str2)
      end
      println()
    end
    str = string(
    "const UK_POLYS = [uk_0_eval, uk_1_eval, uk_2_eval, uk_3_eval, uk_4_eval,
    uk_5_eval, uk_6_eval, uk_7_eval, uk_8_eval, uk_9_eval, uk_10_eval,
    uk_11_eval, uk_12_eval, uk_13_eval, uk_14_eval, uk_15_eval, uk_16_eval,
    uk_17_eval, uk_18_eval, uk_19_eval, uk_20_eval]")
    println(str)
  end
end


text:
module CVUI
include("CVProcessing.jl")
using Gtk, ImageView, .CVProcessing

export initUI

const image_filters = ("*.png", "*.jpg", "*.gif", "*.tiff")
# image drawing

function initUI()
    builder = GtkBuilder(
        filename = joinpath(
            @__DIR__,
            "../resource/gtk",
            "ComputerVisionWithJuliaBook.glade",
        ),
    )

    # windows and dialogs
    global winDraw = builder["main_window_draw"]
    winToolbar = builder["main_window_toolbar"]

    #resize
    global winEditResize = builder["edit_resize_window"]

    #resize
    global winEditCrop = builder["edit_crop_window"]

    #resize
    global winEditScale = builder["edit_scale_window"]

    #sliders
    global widthSlider = builder["slider_width_adjustment"]
    global heightSlider = builder["slider_height_adjustment"]
    global scaleSlider = builder["slider_scale_adjustment"]

    # crop sliders
    global topLeftXSlider = builder["crop_top_left_x"]
    global topLeftYSlider = builder["crop_top_left_y"]

    global bottomRightXSlider = builder["crop_bottom_right_x"]
    global bottomRightYSlider = builder["crop_bottom_right_y"]
    #  signals handling

    # open/save
    connect_widget(open_file, builder["open_image_btn"], "clicked", true)
    connect_widget(save_file, builder["save_image_btn"], "clicked", true)

    # Dimensions
    connect_widget(resize_image, builder["image_dimensions_resize_btn"], "clicked", true)
    connect_widget(crop_image, builder["image_dimensions_crop_btn"], "clicked", true)
    connect_widget(scale_image, builder["image_dimensions_scale_btn"], "clicked", true)

    #resize sliders
    connect_widget(widthset, builder["slider_width_adjustment"], "value-changed")
    connect_widget(heightset, builder["slider_height_adjustment"], "value-changed")

    # resize ok, cancel
    connect_widget(preview_resize, builder["resize_window_preview_btn"], "clicked")
    connect_widget(apply_resize, builder["resize_window_apply_btn"], "clicked")
    connect_widget(cancel_resize, builder["resize_window_cancel_btn"], "clicked")

    # crop ok, cancel
    connect_widget(preview_crop, builder["crop_window_preview_btn"], "clicked")
    connect_widget(apply_crop, builder["crop_window_apply_btn"], "clicked")
    connect_widget(cancel_crop, builder["crop_window_cancel_btn"], "clicked")

    # scale ok, cancel
    connect_widget(preview_scale, builder["scale_window_preview_btn"], "clicked")
    connect_widget(apply_scale, builder["scale_window_apply_btn"], "clicked")
    connect_widget(cancel_scale, builder["scale_window_cancel_btn"], "clicked")

    connect_widget(scaleset, builder["slider_scale_adjustment"], "value-changed")

    global originalImage = CVProcessing.defaultimage()
    global processedImage = copy(originalImage)
    @sync redrawImage(originalImage)

    showall.([winDraw, winToolbar])

    if !isinteractive()
        c = Condition()
        signal_connect(winDraw, :destroy) do widget
            notify(c)
        end
        wait(c)
    end
end

connect_widget(handler, w, event, after = true) = signal_connect(handler, w, event, after)

set_size!(w::Gtk.GtkWidget, width::Int, height::Int) = (@sigatom ccall(
    (:gtk_widget_set_size_request, Gtk.libgtk),
    Nothing,
    (Ptr{GObject}, Cint, Cint),
    w,
    width,
    height,
))

set_value!(w::Gtk.GObject, value::Float64) = (@sigatom ccall(
    (:gtk_adjustment_set_value, Gtk.libgtk),
    Nothing,
    (Ptr{GObject}, Cdouble),
    w,
    value,
))

set_upper_value!(w::Gtk.GObject, value::Float64) = (@sigatom ccall(
    (:gtk_adjustment_set_upper, Gtk.libgtk),
    Nothing,
    (Ptr{GObject}, Cdouble),
    w,
    value,
))

set_lower_value!(w::Gtk.GObject, value::Float64) = (@sigatom ccall(
    (:gtk_adjustment_set_lower, Gtk.libgtk),
    Nothing,
    (Ptr{GObject}, Cdouble),
    w,
    value,
))

get_value(w::Gtk.GObject) =
    Float64(ccall((:gtk_adjustment_get_value, Gtk.libgtk), Cdouble, (Ptr{GObject},), w))

set_sensitive!(w::Gtk.GtkWidget, state::Bool) = (@sigatom ccall(
    (:gtk_widget_set_sensitive, Gtk.libgtk),
    Nothing,
    (Ptr{GObject}, Cint),
    w,
    state,
))

enable(w) = @sync set_sensitive!(w, true)

disable(w) = @sync set_sensitive!(w, false)

function redrawImage(img)
    global currentHeight, currentWidth = size(img)

    drawAreaFrame, canvas = ImageView.frame_canvas(:none)
    set_size!(winDraw, currentWidth, currentHeight)

    empty!(winDraw)
    push!(winDraw, drawAreaFrame)
    imshow!(canvas, img)
    showall(winDraw)
end

function open_file(w)
    path = open_dialog("Open image file", GtkNullContainer(), image_filters)
    global originalImage = CVProcessing.open_file(path)
    global processedImage = copy(originalImage)
    @sync redrawImage(originalImage)
end

function save_file(w)
    path = save_dialog("Save image file", GtkNullContainer(), image_filters)
    CVProcessing.save_file(path, processedImage)
end

function resize_image(w)
    set_value!(widthSlider, convert(Float64, currentWidth))
    set_value!(heightSlider, convert(Float64, currentHeight))

    show(winEditResize)
end

widthset(widget) = global currentWidth = Gtk.GAccessor.value(widget)

function preview_resize(w)
    disable(w)
    @sync redrawImage(CVProcessing.resize_image(
        processedImage,
        currentWidth,
        currentHeight,
    ))
    enable(w)
end

function apply_resize(widget)
    hide(winEditResize)
    global processedImage =
        CVProcessing.resize_image(processedImage, currentWidth, currentHeight)
    @sync redrawImage(processedImage)
end

heightset(widget) = global currentHeight = Gtk.GAccessor.value(widget)

function cancel_resize(w)
    hide(winEditResize)
    @sync redrawImage(processedImage)
end

function crop_image(w)
    set_value!(topLeftXSlider, 1.0)
    set_upper_value!(topLeftXSlider, currentWidth / 2.0)

    set_value!(topLeftYSlider, 1.0)
    set_upper_value!(topLeftYSlider, currentHeight / 2.0)

    set_lower_value!(bottomRightXSlider, currentWidth / 2.0)
    set_upper_value!(bottomRightXSlider, convert(Float64, currentWidth))
    set_value!(bottomRightXSlider, convert(Float64, currentWidth))

    set_lower_value!(bottomRightYSlider, currentHeight / 2.0)
    set_upper_value!(bottomRightYSlider, convert(Float64, currentHeight))
    set_value!(bottomRightYSlider, convert(Float64, currentHeight))

    show(winEditCrop)
end

function preview_crop(w)
    disable(w)
    @sync redrawImage(CVProcessing.crop_image(
        processedImage,
        Gtk.GAccessor.value(topLeftXSlider),
        Gtk.GAccessor.value(topLeftYSlider),
        Gtk.GAccessor.value(bottomRightXSlider),
        Gtk.GAccessor.value(bottomRightYSlider),
    ))
    enable(w)
end

function apply_crop(w)
    hide(winEditCrop)
    global processedImage = CVProcessing.crop_image(
        processedImage,
        Gtk.GAccessor.value(topLeftXSlider),
        Gtk.GAccessor.value(topLeftYSlider),
        Gtk.GAccessor.value(bottomRightXSlider),
        Gtk.GAccessor.value(bottomRightYSlider),
    )
    @sync redrawImage(processedImage)
end

function cancel_crop(widget)
    hide(winEditCrop)
    @sync redrawImage(processedImage)
end

function scale_image(w)
    global currentScale = 0
    set_value!(widthSlider, convert(Float64, currentWidth))
    set_value!(heightSlider, convert(Float64, currentHeight))
    set_value!(scaleSlider, convert(Float64, currentScale))
    show(winEditScale)
end

function preview_scale(w)
    disable(w)

    @sync redrawImage(CVProcessing.scale_image(
        processedImage,
        currentWidth,
        currentHeight,
        currentScale == 0 ? 0 : currentScale / 100,
    ))

    enable(w)
end

function apply_scale(widget)
    hide(winEditScale)
    global processedImage = CVProcessing.scale_image(
        processedImage,
        currentWidth,
        currentHeight,
        currentScale == 0 ? 0 : currentScale / 100,
    )
    @sync redrawImage(processedImage)
end

function cancel_scale(widget)
    hide(winEditScale)
    @sync redrawImage(processedImage)
end

scaleset(widget) = global currentScale = convert(Int, Gtk.GAccessor.value(widget))
end

max_stars_repo_name: bemcho/ComputerVisionWithJuliaBook

max_issues_repo_name: bemcho/ComputerVisionWithJuliaBook

max_forks_repo_name: bemcho/ComputerVisionWithJuliaBook


text:
# Translocation of reserve rejected from the SU during catabolism

abstract type PassiveTranslocation end

"""
    passive_translocation!(source, dest)

Realocate state rejected from synthesizing units.
"""
function passive_translocation! end
passive_translocation!(source, dest) =
    passive_translocation!(passivetrans_pars(source), source, dest)
passive_translocation!(f::Nothing, source, dest) = nothing

"""
    DissipativeRejection(y_EC_ECT, y_EN_ENT)

Substrate rejected from synthesizing units during catabolism is returned to
reserve, but with some fraction of loss specified by yield parameters.

$(FIELDDOCTABLE)
"""
@columns struct DissipativePassiveTranslocation{MoMo} <: PassiveTranslocation
#   Field          | Default | Unit       | Bounds     | Log | Description
    y_EC_ECT::MoMo | 1.0     | mol*mol^-1 | (0.0, 1.0) | _   | "yield of translocated C-reserve"
    y_EN_ENT::MoMo | 1.0     | mol*mol^-1 | (0.0, 1.0) | _   | "yield of Translocated N-reserve"
end

passive_translocation!(f::DissipativePassiveTranslocation, source, dest) = begin
    Js, Jd = flux(source), flux(dest)
    Jd[:C,:tra] += -Js[:C,:rej] * f.y_EC_ECT
    Jd[:N,:tra] += -Js[:N,:rej] * f.y_EN_ENT
    nothing
end

"""
    LosslessPassiveTranslocation()

Parameterless rejection where substrate rejected from synthesizing units
during catabolism is returned to reserve without loss.

$(FIELDDOCTABLE)
"""
struct LosslessPassiveTranslocation <: PassiveTranslocation end

passive_translocation!(rejected::LosslessPassiveTranslocation, source, dest) = begin
    Js, Jd = flux(source), flux(dest)
    Jd[:C,:tra] += -Js[:C,:rej]
    Jd[:N,:tra] += -Js[:N,:rej]
    nothing
end

# Active translocation. Not required in practice, as
# translocation of rejected reserves gives resonable behaviour.

abstract type ActiveTranslocation end

@mix @flattenable @columns struct Active{F}
#   Field          | Flat | Default  | Unit       | Bounds    | Log | Description
    κtra::F        | true | 0.6      | _          | (0.0,1.0) | _   | "Reserve flux allocated to translocation"
end

κtra(activetrans_pars::ActiveTranslocation) = activetrans_pars.κtra

"""
    active_translocation!(o1, o2)

Translocation that actively moves a fraction of catabolised
reserve between organs.
"""
function active_translocation! end
active_translocation!(o1, o2) =
    active_translocation!(activetrans_pars(o1), o1, o2)
active_translocation!(p::Nothing, o1, o2) = nothing

"""
    DissipativeTranslocation(κtra, y_E_ET)

Translocation with dissipative losses to the environment.

$(FIELDDOCTABLE)
"""
@Active struct DissipativeActiveTranslocation{MoMo} <: ActiveTranslocation
    y_E_ET::MoMo   | true | 0.8      | mol*mol^-1 | (0.0,1.0) | _   | "yield of translocated reserve:"
end

active_translocation!(p::DissipativeActiveTranslocation, o1, o2) = begin
    # outgoing translocation
    trans = κtra(o1) * E_ctb(o1) / (1 - κsoma(o1))
    reserve_drain!(o1, :tra, trans)

    # incoming translocation
    transx = κtra(o2) * o2.J1[:E,:ctb] / (1 - κsoma(o2))
    o1.J[:E,:tra] += transx * y_E_ET(o2)

    nothing
end

"""
    DissipativeTranslocation(κtra, y_E_ET)

Perfect translocation between structures.

$(FIELDDOCTABLE)
"""
@Active struct LosslessActiveTranslocation{} <: ActiveTranslocation end

active_translocation!(p::LosslessActiveTranslocation, o1, o2) = begin
    # outgoing translocation
    trans = κtra(o1) * o1.J1[:E,:ctb] / (1 - κsoma(o1))
    reserve_drain!(o1, :tra, trans)

    # incoming translocation
    transx = κtra(o2) * o2.J1[:E,:ctb] / (1 - κsoma(o1))

    o1.J[:E,:tra] += transx
    nothing
end

"""
    translocation!(organs::Tuple)

Translocation occurs between adjacent organs in
both directions.

Both active and passive translocation are applied,
although `nothing` valued formulations do not translocate.
"""
function translocation! end
translocation!(organs::Tuple{T1,T2}) where {T1,T2} = begin
    passive_translocation!(organs[1], organs[2])
    passive_translocation!(organs[2], organs[1])
    active_translocation!(organs[1], organs[2])
    active_translocation!(organs[2], organs[1])
end
# translocation!(organs::Tuple) = translocation!(organs...)
translocation!(organs::Tuple{T}) where T = nothing
translocation!(organs::Tuple{}) = nothing


text:
@testset "CoreOperations Tests" begin

    @testset "Sub Tensor Tests" begin
        indices = [1 1 2 2; 1 2 2 2; 1 3 3 3;2 2 2 3] #fully connected
        vals = [1.0,2.0,3.0,4.0]

        DictA = SSSTensor([(indices[i,:],vals[i]) for i in 1:4],3)
        COOA  = COOTen(indices,vals,3)

        sub_DictA = get_sub_tensor(DictA,[1,2],false)
        sub_COOA  = get_sub_tensor(COOA,[1,2],false)

        @test sub_COOA == COOTen(indices[1:2,:],vals[1:2],2)
    end

    @testset "save/load tests" begin
        _, indices, vals = set_up(n,ord,nnz)

        A = COOTen(indices,vals,n)
        save(A, tempath)
        B = load(tempath;enforceFormatting = true,type = "COOTen")
        @test A == B

        A = SSSTensor([(indices[i,:],vals[i]) for i in 1:nnz],n)
        save(A, tempath)
        B = load(tempath;enforceFormatting = true,type ="DICTen")
        @test A == B
        #check for .ssten file extension functionality

        rm(tempath)
    end

    @testset  "Data Representation" begin

        @testset "flatten" begin
            x = rand(n)
            A = SSSTensor([([2,3,3],1.0),([1,2,4],-1.0),([3,4,5],1.0)])
            flattened_A = flatten(A)
            kron_x = reduce(kron,[x for _ in 1:(order(A)-1)])
            kron_mul = flattened_A*kron_x
            @test LinearAlgebra.norm(SparseSymmetricTensors.contract(A,x,2) -kron_mul)/LinearAlgebra.norm(kron_mul) < tol
        end

        @testset "dense" begin
            @test ones(2,2,2,2) == dense(SSSTensor(ones(2,2,2,2)))
        end
    end

end


text:
using Revise
using DataFrames, CSV, JDF
using WeakRefStrings

a[!, :stringarr] = StringArray(rand(["a", "a", "b"], size(a,1)))
a[!, :cate] = categorical(a[!, :stringarr])

@time a = CSV.read("c:/data/feature_matrix_cleaned.csv");

@time savejdf("c:/data/feature_matrix_cleaned.csv.jdf", a)

a = nothing
@time a = loadjdf("c:/data/feature_matrix_cleaned.csv.jdf")

type_compress!(a, compress_float=true)
@time savejdf("c:/data/feature_matrix_cleaned.csv.compressed.jdf", a)

using BenchmarkTools
@benchmark a = loadjdf("c:/data/feature_matrix_cleaned.csv.jdf")


text: print("Checking")

max_stars_repo_name: alexflorez/LearningJulia

max_issues_repo_name: alexflorez/LearningJulia

max_forks_repo_name: alexflorez/LearningJulia


text:
################################################################################
#
#    NfOrd/NfOrdFracIdl.jl : Fractional ideals of generic
#                               orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export frac_ideal

################################################################################
#
#  Consistency
#
################################################################################

function isconsistent(x::NfOrdFracIdl)
  return isconsistent(numerator(x))
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, A::FakeFmpqMat, A_in_hnf::Bool = false) -> NfOrdFracIdl

Creates the fractional ideal of $\mathcal O$ with basis matrix $A$. If A_in_hnf
is set, then it is assumed that the numerator of $A$ is already in lower left
HNF.
"""
function frac_ideal(O::NfOrd, x::FakeFmpqMat, x_in_hnf::Bool = false)
  !x_in_hnf ? x = hnf(x) : nothing
  z = NfOrdFracIdl(O, x)
  return z
end

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, A::fmpz_mat, b::fmpz, A_in_hnf::Bool = false) -> NfOrdFracIdl

Creates the fractional ideal of $\mathcal O$ with basis matrix $A/b$. If
A_in_hnf is set, then it is assumed that $A$ is already in lower left HNF.
"""
function frac_ideal(O::NfOrd, x::fmpz_mat, y::fmpz=fmpz(1), x_in_hnf::Bool = false)
  !x_in_hnf ? x = _hnf(x, :lowerleft) : nothing
  y = FakeFmpqMat(x, y)
  z = NfOrdFracIdl(O, y)
  return z
end

frac_ideal(O::NfOrd, x::fmpz_mat, y::Integer) = frac_ideal(O, x, fmpz(y))

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, I::NfOrdIdl) -> NfOrdFracIdl

Turns the ideal $I$ into a fractional ideal of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdIdl)
  z = NfOrdFracIdl(O, x, fmpz(1))
  return z
end

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, I::NfOrdIdl, b::fmpz) -> NfOrdFracIdl

Creates the fractional ideal $I/b$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdIdl, y::fmpz)
  z = NfOrdFracIdl(O, x, deepcopy(y)) # deepcopy x?
  return z
end

frac_ideal(O::NfOrd, x::NfOrdIdl, y::Integer) = frac_ideal(O, x, fmpz(y))

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, a::nf_elem) -> NfOrdFracIdl

Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::nf_elem)
  z = NfOrdFracIdl(O, deepcopy(x))
  return z
end

@doc Markdown.doc"""
    frac_ideal(O::NfOrd, a::NfOrdElem) -> NfOrdFracIdl

Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdElem)
  z = NfOrdFracIdl(O, elem_in_nf(x))
  return z
end

function frac_ideal_from_z_gens(O::NfAbsOrd{S, T}, b::Vector{T}) where {S, T}
  d = degree(O)
  den = lcm([ denominator(bb, O) for bb in b ])
  num = ideal_from_z_gens(O, [ O(den*bb) for bb in b ])
  return frac_ideal(O, num, den)
end

################################################################################
#
#  String I/O
#
################################################################################

export basis_matrix, norm, inv, ==, *, integral_split

export parent, order, basis_matrix, basis_mat_inv, basis, norm,
       ring_of_multipliers, ==

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfOrdFracIdl) = NfOrdFracIdlSet(order(a), false)

function FracIdealSet(O::NfAbsOrd)
   return NfOrdFracIdlSet(O, false)
end

function Base.hash(a::NfOrdFracIdl, h::UInt)
  b = simplify(a)
  return hash(numerator(b, copy = false), hash(denominator(b, copy = false), h))
end

elem_type(::Type{NfOrdFracIdlSet}) = NfOrdFracIdl

elem_type(::NfOrdFracIdlSet) = NfOrdFracIdl

parent_type(::Type{NfOrdFracIdl}) = NfOrdFracIdlSet

################################################################################
#
#  Order
#
################################################################################

order(a::NfOrdFracIdlSet) = a.order

@doc Markdown.doc"""
    order(a::NfOrdFracIdl) -> NfOrd
The order that was used to define the ideal $a$.
"""
order(a::NfOrdFracIdl) = a.order

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

@doc Markdown.doc"""
    basis_mat_inv(I::NfOrdFracIdl) -> FakeFmpqMat

Returns the inverse of the basis matrix of $I$.
"""
function basis_mat_inv(a::NfOrdFracIdl)
  if isdefined(a, :basis_mat_inv)
    return deepcopy(a.basis_mat_inv)
  else
    a.basis_mat_inv = inv(basis_matrix(a))
    return deepcopy(a.basis_mat_inv)
  end
end

################################################################################
#
#  Basis
#
################################################################################

@doc Markdown.doc"""
    basis(I::NfOrdFracIdl) -> Array{nf_elem, 1}

Returns the $\mathbf Z$-basis of $I$.
"""
function basis(a::NfOrdFracIdl)
  B = basis_matrix(a)
  d = degree(order(a))
  O = order(a)
  K = nf(O)
  Oba = O.basis_nf
  res = Array{nf_elem}(undef, d)

  for i in 1:d
    z = K()
    for j in 1:d
      z = z + B.num[i, j]*Oba[j]
    end
    z = divexact(z, B.den)
    res[i] = z
  end

  return res
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(x::NfOrdFracIdl, dict::IdDict)
  z = NfOrdFracIdl(numerator(x), denominator(x))
  if isdefined(x, :basis_matrix)
    z.basis_matrix = deepcopy(x.basis_matrix)
  end
  return z
end

################################################################################
#
#  Basis matrix
#
################################################################################

function assure_has_basis_matrix(a::NfOrdFracIdl)
  if isdefined(a, :basis_matrix)
    return nothing
  end
  if !isdefined(a, :num)
    error("Not a valid fractional ideal")
  end

  a.basis_matrix = FakeFmpqMat(basis_matrix(numerator(a, copy = false)), denominator(a))
  return nothing
end

@doc Markdown.doc"""
    basis_matrix(I::NfOrdFracIdl) -> FakeFmpqMat

Returns the basis matrix of $I$ with respect to the basis of the order.
"""
function basis_matrix(x::NfOrdFracIdl; copy::Bool = true)
  assure_has_basis_matrix(x)
  if copy
    return deepcopy(x.basis_matrix)
  else
    return x.basis_matrix
  end
end

################################################################################
#
#  Numerator and denominator
#
################################################################################

function assure_has_numerator_and_denominator(a::NfOrdFracIdl)
  if isdefined(a, :num) && isdefined(a, :den)
    return nothing
  end
  if !isdefined(a, :basis_matrix)
    error("Not a valid fractional ideal")
  end

  a.num = ideal(order(a), numerator(basis_matrix(a, copy = false)))
  a.den = denominator(basis_matrix(a, copy = false))
  return nothing
end

function numerator(x::NfOrdFracIdl; copy::Bool = true)
  assure_has_numerator_and_denominator(x)
  if copy
    return deepcopy(x.num)
  else
    return x.num
  end
end

function denominator(x::NfOrdFracIdl; copy::Bool = true)
  assure_has_numerator_and_denominator(x)
  if copy
    return deepcopy(x.den)
  else
    return x.den
  end
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, s::NfOrdFracIdlSet)
   print(io, "Set of fractional ideals of ")
   print(io, s.order)
end

function show(io::IO, id::NfOrdFracIdl)
  if isdefined(id, :num) && isdefined(id, :den)
    print(io, "1//", denominator(id, copy = false), " * ")
    print(io, numerator(id, copy = false))
  else
    print(io, "Fractional ideal with basis matrix\n")
    print(io, basis_matrix(id, copy = false))
  end
end

################################################################################
#
#  Norm
#
################################################################################

@doc Markdown.doc"""
    norm(I::NfOrdFracIdl) -> fmpq

Returns the norm of $I$
"""
function norm(A::NfOrdFracIdl)
  if isdefined(A, :norm)
    return deepcopy(A.norm)
  else
    A.norm = norm(numerator(A, copy = false))//denominator(A, copy = false)^degree(order(A))
    return deepcopy(A.norm)
  end
end

################################################################################
#
#  Minimum
#
################################################################################

function minimum(A::NfOrdFracIdl)
  return minimum(numerator(A, copy = false))//denominator(A, copy = false)
end

################################################################################
#
#  Inverse
#
################################################################################

@doc Markdown.doc"""
    inv(A::NfOrdFracIdl) -> NfOrdFracIdl

Returns the fractional ideal $B$ such that $AB = \mathcal O$.
"""
function inv(A::NfOrdFracIdl)
  B = inv(numerator(A, copy = false))
  g = gcd(denominator(B, copy = false), denominator(A, copy = false))
  B.den = divexact(denominator(B, copy = false), g)
  a = divexact(denominator(A, copy = false), g)
  return prod_by_int(B, a)
end

################################################################################
#
#  Simplification
#
################################################################################

function simplify_exact!(A::NfOrdFracIdl)
  assure_has_numerator_and_denominator(A)
  g = A.den

  A.den = fmpz(1)
  A.num = divexact(A.num, g)
end

function simplify(A::NfOrdFracIdl)
  assure_has_numerator_and_denominator(A)
  simplify(A.num)

  if has_2_elem(A.num)
    ZK = order(A)
    g = Base.reduce(gcd, coordinates(ZK(A.num.gen_two)))
    g = gcd(g, A.den)
    g = gcd(g, A.num.gen_one)
  else
    b = basis_matrix(A.num, copy = false)
    g = gcd(denominator(A), content(b))
  end

  if g != 1
    A.num = divexact(A.num, g)
    A.den = divexact(A.den, g)
  end

  return A
end

################################################################################
#
#  Primeness
#
################################################################################

isprime_known(A::NfOrdFracIdl) = isprime_known(numerator(A, copy = false))

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(x::NfOrdFracIdl, y::NfOrdFracIdl) -> Bool

Returns whether $x$ and $y$ are equal.
"""
function ==(A::NfOrdFracIdl, B::NfOrdFracIdl)
  #return B.den * basis_matrix(A.num) == A.den * basis_matrix(B.num)
  if !ismaximal_known(order(A)) || !ismaximal(order(A))
    return basis_matrix(A, copy = false) == basis_matrix(B, copy = false)
  end

  if !isdefined(A, :num) || !isdefined(B, :num)
    return basis_matrix(A, copy = false) == basis_matrix(B, copy = false)
  end
  D = inv(B)
  E = prod(A, D)
  C = simplify(E)
  return isone(denominator(C, copy = false)) && isone(norm(C))
end

################################################################################
#
#  Binary operations
#
################################################################################

function prod(a::NfOrdFracIdl, b::NfOrdFracIdl)
  A = numerator(a, copy = false)*numerator(b, copy = false)
  return NfOrdFracIdl(A, denominator(a, copy = false)*denominator(b, copy = false))
end

@doc Markdown.doc"""
    *(I::NfOrdFracIdl, J::NfOrdFracIdl) -> NfOrdFracIdl

Returns $IJ$.
"""
*(A::NfOrdFracIdl, B::NfOrdFracIdl) = prod(A, B)

################################################################################
#
#  Powering
#
################################################################################

function ^(A::NfOrdFracIdl, a::Int)
  if a == 0
    B = NfOrdFracIdl(ideal(order(A), 1), fmpz(1))
    return B
  end

  if a == 1
    return A # copy?
  end

  if a < 0
    return inv(A^(-a))
  end

  if a == 2
    return A*A
  end

  if mod(a, 2) == 0
    return (A^div(a, 2))^2
  else
    return A * A^(a - 1)
  end
end

function //(A::NfOrdFracIdl, B::NfOrdFracIdl)
  C = prod(A, inv(B))
  return C
end

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function prod_by_int(A::NfOrdFracIdl, a::fmpz)
  return NfOrdFracIdl(numerator(A, copy = false) * a, denominator(A))
end

*(A::NfOrdFracIdl, a::fmpz) = prod_by_int(A, a)
*(a::fmpz, A::NfOrdFracIdl) = prod_by_int(A, a)

function *(A::NfOrdIdl, B::NfOrdFracIdl)
  z = NfOrdFracIdl(A*numerator(B, copy = false), denominator(B))
  return z
end

*(A::NfOrdFracIdl, B::NfOrdIdl) = NfOrdFracIdl(numerator(A, copy = false)*B, denominator(A))

function *(A::NfOrdFracIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfOrdFracIdl) = A*a

function *(A::NfOrdIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfOrdIdl) = A*a

function //(A::NfOrdFracIdl, B::NfOrdIdl)
  C = prod(A, inv(B))
  return C
end

function //(A::NfOrdIdl, B::NfOrdIdl)
  return A*inv(B)
end

function //(A::NfOrdIdl, B::NfOrdFracIdl)
  return A*inv(B)
end

function //(A::NfOrdFracIdl, a::nf_elem)
  C = prod(A, Idl((order(A), inv(a))))
  return C
end

function //(A::NfOrdIdl, d::fmpz)
  return Hecke.NfOrdFracIdl(A, d)
end

function //(A::NfOrdIdl, d::Integer)
  return A//fmpz(d)
end

function +(A::NfOrdIdl, B::NfOrdFracIdl)
  return (A*denominator(B)+numerator(B))//denominator(B)
end

+(A::NfOrdFracIdl, B::NfOrdIdl) = B+A

function +(A::NfOrdFracIdl, B::Hecke.NfOrdFracIdl)
  d = lcm(denominator(A), denominator(B))
  ma = div(d, denominator(A))
  mb = div(d, denominator(B))
  return (numerator(A)*ma + numerator(B)*mb)//d
end

function *(x::nf_elem, y::NfOrd)
  b, z = _check_elem_in_order(denominator(x, y)*x, y)
  return NfOrdFracIdl(ideal(y, y(z)), denominator(x, y))
end

################################################################################
#
#  Ad hoc comparison
#
################################################################################

function ==(A::NfOrdIdl, B::NfOrdFracIdl)
  if order(A) !== order(B)
    return false
  end

  C = simplify(B)

  if denominator(C, copy = false) != 1 || numerator(C, copy = false) != A
    return false
  else
    return true
  end
end

==(A::NfOrdFracIdl, B::NfOrdIdl) = B == A

################################################################################
#
#  Conversion
#
################################################################################

function (ord::NfOrdIdlSet)(b::NfOrdFracIdl)
  denominator(b, copy = false) > 1 && error("not integral")
  return numerator(b, copy = false)
end

function ideal(O::NfOrd, a::nf_elem)
  return a*O
end

@doc Markdown.doc"""
    integral_split(A::NfOrdFracIdl) -> NfOrdIdl, NfOrdIdl
Computes the unique coprime integral ideals $N$ and $D$ s.th. $A = ND^{-1}$
"""
function integral_split(A::NfOrdFracIdl)
  d = simplify(inv(A + ideal(order(A), fmpz(1))))
  @assert denominator(d) == 1
  n = simplify(A*d)
  @assert denominator(n) == 1
  return numerator(n), numerator(d)
end

@doc Markdown.doc"""
    factor(I::NfOrdFracIdl) -> Dict{NfOrdIdl, Int}
The factorisation of $I$.
"""
function factor(I::NfOrdFracIdl)
  if iszero(norm(I))
    error("Cannot factor zero ideal")
  end
  n, d = integral_split(I)
  fn = factor(n)
  fd = factor(d)
  for (k, v) = fd
    if haskey(fn, k)
      fn[k] -= v
    else
      fn[k] = -v
    end
  end
  return fn
end

function one(A::NfOrdFracIdl)
  return NfOrdFracIdl(ideal(order(A), 1), fmpz(1))
end

@doc Markdown.doc"""
    valuation(A::NfOrdFracIdl, p::NfOrdIdl)
The valuation of $A$ at $p$.
"""
function valuation(A::NfOrdFracIdl, p::NfOrdIdl)
  return valuation(numerator(A, copy = false), p) - valuation(denominator(A, copy = false), p)
end

################################################################################
#
#  Coprime ideals
#
################################################################################

#Given I with v_p(I) = 0, returns element a \in I such that v_p(a) = 0
function coprime_to(I::NfOrdFracIdl, p::NfOrdIdl)
  pi = anti_uniformizer(p)
  a = basis(I)[1]
  l = valuation(a, p)
  @assert l >= 0
  if l > 0
    a = pi^l * a
  end
  @assert valuation(a, p) == 0
  @assert denominator(I)*a in numerator(I)
  return a
end

################################################################################
#
#  Colon
#
################################################################################

function colon(I::NfOrdFracIdl, J::NfOrdFracIdl)
  # Let I = a/c and J = b/d, a and b integral ideals, c, d \in Z, then
  # \{ x \in K | xJ \subseteq I \} = \{ x \in K | xcb \subseteq da \}

  II = numerator(I, copy = false)*denominator(J, copy = false)
  JJ = numerator(J, copy = false)*denominator(I, copy = false)
  return colon(II, JJ)
end

################################################################################
#
#  Move ideals to another order
#
################################################################################

function extend(I::NfOrdFracIdl, O::NfAbsOrd)
  J = extend(numerator(I, copy = false), O)
  return frac_ideal(O, J, denominator(I, copy = false))
end

*(I::NfOrdFracIdl, O::NfAbsOrd) = extend(I, O)
*(O::NfAbsOrd, I::NfOrdFracIdl) = extend(I, O)

function _as_frac_ideal_of_smaller_order(O::NfAbsOrd, I::NfAbsOrdIdl)
  M = basis_matrix(I, copy = false)
  M = M*basis_matrix(order(I), copy = false)*basis_mat_inv(O, copy = false)
  return frac_ideal(O, M)
end

function _as_frac_ideal_of_smaller_order(O::NfAbsOrd, I::NfOrdFracIdl)
  J = _as_frac_ideal_of_smaller_order(O, numerator(I, copy = false))
  return nf(O)(fmpq(1, denominator(I, copy = false)))*J
end


text:
module IPMeasures

using StatsBase
using LinearAlgebra
using Distances
using CUDA

using Distributions
using DistributionsAD
using ConditionalDists

const MetricOrFun = Union{PreMetric,Function}
const CuMatrix = CuArray{<:Float32,2}

"""
    samplecolumns(x::AbstractMatrix, n::Int)

Sample n columns from a matrix. Returns x if the matrix has less than n columns.
"""
function samplecolumns(x::AbstractMatrix, n::Int)
    (size(x,2) > n) ? x[:,sample(1:size(x,2), n, replace=false)] : x
end

"""
    split2(x::AbstractMatrix)

Splits a matrix into halves.
"""
function split2(x::AbstractMatrix)
	n = size(x,2)
	x[:,1:div(n,2)], x[:,div(n,2)+1:end]
end

"""
	pairwisel2(x,y)
	pairwisel2(x)

Calculates pairwise squared euclidean distances of the columns of `x` and `y`
or `x` and `x`. The dispatches for CuArrays are necessary until
https://github.com/JuliaStats/Distances.jl/pull/142 is merged.
"""
pairwisel2(x::Matrix, y::Matrix) = pairwise(SqEuclidean(), x, y, dims=2)
pairwisel2(x::CuMatrix, y::CuMatrix) =
    -2 .* x' * y .+ sum(x.^2, dims=1)' .+ sum(y.^2,dims=1)
pairwisel2(x::AbstractMatrix) = pairwisel2(x,x)

"""
	pairwisel2(x,y)
	pairwisel2(x)

Calculates pairwise cosine distances of the columns of `x` and `y` or
`x` and `x`. The dispatches for CuArrays are necessary until
https://github.com/JuliaStats/Distances.jl/pull/142 is merged.
"""
pairwisecos(x::Matrix, y::Matrix) = pairwise(CosineDist(), x, y, dims=2)
pairwisecos(x::CuMatrix, y::CuMatrix) = 1 .- (x' * y / norm(x) / norm(y))
pairwisecos(x::AbstractMatrix) = pairwisecos(x,x)

include("kernels.jl")
include("mmd.jl")
include("kl_divergence.jl")

"""
	null_distribution(k::AbstractKernel, x, n, l = div(size(x,2),2))

Estimates the null distribution of samples `x` from `n` random draws of subsets of size `l`
"""
null_distribution(k::AbstractKernel, x, n, l = div(size(x,2),2)) =
    _null_distribution(k, pairwisel2(x,x), n, l )

function _null_distribution(k, d, n, l)
	ns = size(d,2)
	map(1:n) do _
		idx = sample(1:ns, l, replace = false)
		cidx = setdiff(1:ns, idx)
		cidx = sample(cidx, min(l, length(cidx)), replace = false)
		mmdfromdist(k, d, idx, cidx)
	end
end

end # module


text:
using Test, QuHamiltonian

@testset "constructor" begin
    @vertex i, j, k, l

    @subham Z[i] # => 1 instance
    @subham Z[i]Z[j] # => 2 instance
    @subham Z[i]Z[j]Z[k] # => 3 instances
    A = zeros(8, 8)
    B = zeros(4, 4)

    # multi-sites
    @subham A[i, j, k] # => 1 instance
    @subham A[i, j, k]B[j, k] # => 2 instance

    @subham(A[i, j, k]) == @subham(A[i, j, k]B[j, k])
end

@testset "vertex region" begin
    @vertex i, j

    LZ1 = @subham Z[i]
    LZ2 = @subham Z[j]
    LZ1 * LZ2 == @subham Z[i]Z[j] # i, j ∈ V1 × V2
    LZ1 * LZ1 == @subham Z[i]Z[i] # i, i ∈ V1^2
end

@testset "neighbor region" begin
    @nearest i, j, k

    LZ1 = @subham Z[i]
    LZ2 = @subham Z[j]
    LZ3 = @subham Z[k]

    @test LZ1 * LZ1 == @subham(Z[i]Z[i])
    # i, j ∈ E(V)
    @test LZ1 * LZ2 == @subham Z[i]Z[j]

    # error
    # neighbor relation should only support 2-local hamiltonians
    @test_throws ErrorException LZ1 * LZ2 * LZ3

    # this should be equivalent
    # name of place holders does not mean anything
    @test LZ1 * LZ2 == LZ2 * LZ3
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module Day7
export solvePart1, solvePart2

import DataStructures: PriorityQueue, enqueue!, dequeue!

function getRawInput()
    lines = []
    open("day7input.txt") do f
        for line in eachline(f)
            push!(lines, line)
        end
    end
    return lines
end

function parseInput(lines)
    parents = Dict()
    children = Dict()
    for s in lines
        step, required = s[37], s[6]
        parents[step] = union(get(parents, step, Set()), Set([required]))
        children[required] = union(get(children, required, Set()), Set([step]))
    end

    return (parents, children)
end

function getInput()
    parseInput(getRawInput())
end

function solvePart1()
    parents, children = getInput()

    queue = PriorityQueue{Char, Char}()
    roots = setdiff(Set(keys(children)), Set(keys(parents)))
    for root in roots
        enqueue!(queue, root, root)
    end

    order = []
    processed = Set()
    while length(queue) > 0
        c = dequeue!(queue)
        push!(processed, c)
        push!(order, c)

        childs = get(children, c, Set())
        for child in childs
            deps = get(parents, child, Set())
            if length(setdiff(parents[child], processed)) == 0
                enqueue!(queue, child, child)
            end
        end
    end
    join(order, "")
end

function timeToDo(task::Char)
    60 + task - 'A' + 1
end

function addToDictionary!(d, time, task)
    d[time] = [get(d, time, []); task]
end

function enqueueUnblockedChildren!(q, c, children, parents, processed)
    childs = get(children, c, Set())
    for child in childs
        deps = get(parents, child, Set())
        if length(setdiff(parents[child], processed)) == 0
            enqueue!(q, child, child)
        end
    end
end

function solvePart2()
    parents, children = getInput()
    freeWorkers = 5
    queue = PriorityQueue{Char, Char}()
    roots = setdiff(Set(keys(children)), Set(keys(parents)))
    for root in roots
        enqueue!(queue, root, root)
    end

    processed = Set()
    seconds = 0
    # a map with 15 -> ['A', 'B'], where 'A' and 'B' are tasks, and 15 is the
    # seconds count when they will finish.
    inProgress = Dict()
    while (length(queue) + length(keys(inProgress)) > 0)
        # Remove completed tasks
        if seconds in keys(inProgress)
            for completed in inProgress[seconds]
                push!(processed, completed)
                enqueueUnblockedChildren!(queue, completed, children, parents, processed)
                freeWorkers += 1
            end
            delete!(inProgress, seconds)
        end
        # To avoid an off-by-one error, break immediately when we finish
        if length(queue) + length(keys(inProgress)) == 0
            break
        end

        # Start as many workers as we have and have tasks for
        while freeWorkers > 0 && length(queue) > 0
            task = dequeue!(queue)
            addToDictionary!(inProgress, seconds + timeToDo(task), task)
            freeWorkers -= 1
        end
        seconds += 1
    end
    return seconds
end

end


text:
using FileIO: open, close, read, write
using Plots

function addToPlot!(plot, xlist, ylist)
    plot!(plot, xlist, ylist)
    empty!(xlist)
    empty!(ylist)
    return
end

xlist, ylist = [], []
p = plot([], [], xlim=(0, 1250), ylim=(0, 200))
open("rl001.fx", "r") do io
    let
        skip_flag = true
        skipchars(ch -> (skip_flag &= (ch != '\f')) || ch == '\f', io)
    end
    for line in eachline(io)
        if line == ""
            continue
        end

        xy = map(str -> tryparse(Float32, str), split(line))
        if length(xy) != 2 || xy[1] == Nothing || xy[2] == Nothing
            if line == "="
                addToPlot!(p, xlist, ylist)
            else
                println("Invalid format")
                exit(0)
            end
        else
            x, y = xy
            append!(xlist, x)
            append!(ylist, y)
        end
    end
end

if !isempty(xlist)
    addToPlot!(p, xlist, ylist)
end

display(p)


text:
immutable Kfold
    indices::Vector{Int}
    k::Int

    function Kfold( n::Int, k::Int )
        new( randperm(n), k )
    end

end

function partition( xval::Kfold, i::Int )
    part = xval.indices[ i : xval.k : size(xval.indices,1) ]
    return setdiff( xval.indices, part ), part
end

length( xval::Kfold ) = xval.k
Base.start( xval::Kfold ) = 1
Base.next( xval::Kfold, i::Int ) = partition( xval, i ), i + 1
Base.done( xval::Kfold, i::Int ) = i > xval.k

function cross_validate( train_f, error_f, partitions )
    scores = Float64[]
    for (train,test) in partitions
        model = train_f( train )
        push!( scores, error_f( model, test ) )
    end
    return scores
end

function kfold( train_f, error_f, n, k )
    return cross_validate( train_f, error_f, Kfold( n, k ) )
end

max_stars_repo_name: tbepler/julia-regression

max_issues_repo_name: tbepler/julia-regression

max_forks_repo_name: tbepler/julia-regression


text:
module RoguePkg

import Base: Pkg.Dir.path, Pkg.test, show

export @pkg_for_str, @pkg_at_str, @pkg_local_str

######################################################################
# general
######################################################################

abstract type PackageSpec end

"""
    ensure_dirpath(path)

If `path` is not a directory, append a `/`.
"""
ensure_dirpath(path) = isdirpath(path) ? path : path * "/"

"""
    pkg_root(package_spec)

Return the root directory for a package specification.  See the documentation
for ways to specify a package.
"""
function pkg_root end

"""
    pkg_name(package_spec)

Return the name of the package as a string. See the documentation for ways to
specify a package.
"""
pkg_name(pkg::PackageSpec) = splitdir(splitdir(pkg_root(pkg))[1])[2]

# strings are passed through
pkg_name(s::AbstractString) = s

######################################################################
# package for a module
######################################################################

"""
    PkgForModule(module_name)

Package for the module `module_name`, located dynamically using.
"""
struct PkgForModule{S <: AbstractString} <: PackageSpec
    module_name::S
end

show(io::IO, pkg::PkgForModule) =
    print(io, "package for module \"$(pkg.module_name)\" (located dynamically)")

"""
    pkg_for"Module"

Package for `"Module"`, located dynamically.
"""
macro pkg_for_str(module_name)
    PkgForModule(module_name)
end

if VERSION ≥ v"0.7-"
    _find_package(name) = Base.find_package(name)
else
    _find_package(name) = Base.find_in_path(name, nothing)
end

function pkg_root(pkg::PkgForModule)
    module_path = _find_package(pkg.module_name)
    module_path == nothing &&
        error("module \"$(pkg.module_name)\" could not be found in the load path")
    src_dir = splitdir(module_path)[1]
    normpath(src_dir, "..")
end

######################################################################
# package at a path
######################################################################

"""
    PkgAtPath(path)

Package at the given path. `~` is resolved.
"""
struct PkgAtPath{S <: AbstractString} <: PackageSpec
    package_path::S
end

function show(io::IO, pkg::PkgAtPath)
    path = pkg.package_path
    print(io, "package at path $(path)")
    '~' in path && print(io, "(~ expanded dynamically)")
    nothing
end

"""
    pkg_at"path"

Package at the given path. `~` is resolved.
"""
macro pkg_at_str(package_path)
    PkgAtPath(package_path)
end

pkg_root(pkg::PkgAtPath) = ensure_dirpath(expanduser(pkg.package_path))

######################################################################
# package in a (specified) local directory
######################################################################

const PKG_LOCAL_DIR_KEY = "JULIA_LOCAL_PACKAGES"

const ENV_PKG_LOCAL_DIR_KEY = "ENV[\"" * PKG_LOCAL_DIR_KEY * "\"]"

"""
    PkgLocal(directory)

Package in `ENV[$(PKG_LOCAL_DIR_KEY)]/directory`.
"""
struct PkgLocal{S <: AbstractString} <: PackageSpec
    package_name::S
end

function pkg_local_dir(expand = true)
    haskey(ENV, PKG_LOCAL_DIR_KEY) ||
        error("Specify local package directory with $(ENV_PKG_LOCAL_DIR_KEY)")
    dir = ENV[PKG_LOCAL_DIR_KEY]
    expand ? expanduser(dir) : dir
end

function show(io::IO, pkg::PkgLocal)
    print(io, "package \"$(pkg.package_name)\" in ")
    if haskey(ENV, PKG_LOCAL_DIR_KEY)
        print(io, pkg_local_dir(false))
    else
        print(io, ENV_PKG_LOCAL_DIR_KEY * " ")
        print_with_color(:red, io, "(not set)")
    end
end

macro pkg_local_str(package_name)
    PkgLocal(package_name)
end

pkg_root(pkg::PkgLocal) =
    ensure_dirpath(normpath(pkg_local_dir(), pkg_name(pkg)))

pkg_name(pkg::PkgLocal) = pkg.package_name

######################################################################
# provide methods for types above
######################################################################

# when pkg is the first argument, use its root
Pkg.Dir.path(pkg::PackageSpec, names::Union{PackageSpec, AbstractString}...) =
    normpath(pkg_root(pkg), pkg_name.(names)...)

# when pkg is not the first argument, use its name
Pkg.Dir.path(names::Union{PackageSpec, AbstractString}...) =
    Pkg.Dir.path(pkg_name.(names)...)

# just works for a single package
function Pkg.test(pkg::PackageSpec; coverage = false)
    errs = AbstractString[]
    nopkgs = AbstractString[]
    notests = AbstractString[]
    Base.cd(() -> Pkg.Entry.test!(pkg_root(pkg), errs, nopkgs, notests; coverage = coverage),
            Pkg.Dir.path())
    # FIXME check errs, nopks, notests
    nothing
end

Pkg.build(pkgs::PackageSpec...) =
    Pkg.cd(Pkg.Entry.build, [Pkg.dir.(pkgs)...])

end # module