text:
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Tensor Manipulation: Psi4 and NumPy manipulation routines
# Contracting tensors together forms the core of the Psi4Julia project. First let us consider the popluar [Einstein Summation Notation](https://en.wikipedia.org/wiki/Einstein_notation) which allows for very succinct descriptions of a given tensor contraction.
#
# For example, let us consider a [inner (dot) product](https://en.wikipedia.org/wiki/Dot_product):
# $$c = \sum_{ij} A_{ij} * B_{ij}$$
#
# With the Einstein convention, all indices that are repeated are considered summed over, and the explicit summation symbol is dropped:
# $$c = A_{ij} * B_{ij}$$
#
# This can be extended to [matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication):
# \begin{align}
# \rm{Conventional}\;\;\;  C_{ik} &= \sum_{j} A_{ij} * B_{jk} \\
# \rm{Einstein}\;\;\;  C &= A_{ij} * B_{jk} \\
# \end{align}
#
# Where the $C$ matrix has *implied* indices of $C_{ik}$ as the only repeated index is $j$.
#
# However, there are many cases where this notation fails. Thus we often use the generalized Einstein convention. To demonstrate let us examine a [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)):
# $$C_{ij} = \sum_{ij} A_{ij} * B_{ij}$$
#
#
# This operation is nearly identical to the dot product above, and is not able to be written in pure Einstein convention. The generalized convention allows for the use of indices on the left hand side of the equation:
# $$C_{ij} = A_{ij} * B_{ij}$$
#
# Usually it should be apparent within the context the exact meaning of a given expression.
#
# Finally we also make use of Matrix notation:
# \begin{align}
# {\rm Matrix}\;\;\;  \bf{D} &= \bf{A B C} \\
# {\rm Einstein}\;\;\;  D_{il} &= A_{ij} B_{jk} C_{kl}
# \end{align}
#
# Note that this notation is signified by the use of bold characters to denote matrices and consecutive matrices next to each other imply a chain of matrix multiplications!

# ## Tensor Operations
#
# To perform most operations we turn to tensor packages (here we use [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl)). Those allow Einstein convention as an input. In addition to being much easier to read, manipulate, and change, it has (usually) optimal performance.
# First let us import our normal suite of modules:

using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy")
using TensorOperations: @tensor

# We can then use conventional Julia loops or `@tensor` to perform the same task.

using BenchmarkTools: @btime, @belapsed

# With `@btime`/`@belapsed` we average time over several executions to have more reliable timings than `@time`/`@elapsed` (single execution).

# **WARNING: We are using Julia's global variables, and those are known to be less efficient than local variables. It is better to wrap code inside function.**

# To begin let us consider the construction of the following tensor (which you may recognize):
# $$G_{pq} = 2.0 * I_{pqrs} D_{rs} - 1.0 * I_{prqs} D_{rs}$$
#
# Keep size relatively small as these 4-index tensors grow very quickly in size.

# +
dims = 20

@assert dims <= 30 "Size must be smaller than 30."
D = rand(dims, dims)
I = rand(dims, dims, dims, dims)

# Build the Fock matrix using loops, while keeping track of time
println("Time for loop G build:")
Gloop = @btime begin
   Gloop = np.zeros((dims, dims))
   @inbounds for ind in CartesianIndices(I)
       p, q, r, s = Tuple(ind)
       Gloop[p, q] += 2I[p, q, r, s] * D[r, s]
       Gloop[p, q] -=  I[p, r, q, s] * D[r, s]
   end
   Gloop
end

# Build the Fock matrix using einsum, while keeping track of time
println("Time for @tensor G build:")
G = @btime @tensor G[p,q] := 2I[p,q,r,s] * D[r,s] - I[p,r,q,s] * D[r,s]

# Make sure the correct answer is obtained
println("Loop and einsum builds of the Fock matrix match?    ", np.allclose(G, Gloop))
println()
# Print out relative times for explicit loop vs einsum Fock builds
#println("G builds with einsum are $(g_loop_time/einsum_time) times faster than Julia loops!")
# -

# As you can see, the `@tensor` macro can be considerably faster than plain Julia loops.

# ## Matrix multiplication chain/train
#
# Now let us turn our attention to a more canonical matrix multiplication example such as:
# $$D_{il} = A_{ij} B_{jk} C_{kl}$$
#
# Matrix multiplication is an extremely common operation in all branches of linear algebra. Thus, these functions have been optimized to be extremely efficient. `@tensor` uses it. The matrix product will explicitly compute the following operation:
# $$C_{ij} = A_{ij} * B_{ij}$$
#
#
# This is Julia's matrix multiplication method `*` for matrices.

# +
dims = 200
A = rand(dims, dims)
B = rand(dims, dims)
C = rand(dims, dims)

# First compute the pair product
tmp_dot = A * B
@tensor tmp_tensor[i,k] := A[i,j] * B[j,k]
println("Pair product allclose? ", np.allclose(tmp_dot, tmp_tensor))
# -

# Now that we have proved exactly what `*` product does, let us consider the full chain and do a timing comparison:

D_dot = A * B * C
@tensor D_tensor[i,l] := A[i,j] * B[j,k] * C[k,l]
println("Chain multiplication allclose? ", np.allclose(D_dot, D_tensor))

# +
println()
println("* time:")
@btime A * B * C

println()
println("@tensor time:")
@btime @tensor D_tensor[i,l] := A[i,j] * B[j,k] * C[k,l];
# -

# Both have similar timings, and both call [Basic Linear Algebra Subprograms (BLAS)](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms). The BLAS routines are highly optimized and threaded versions of the code.
#  - The `@tensor` code will factorize the operation by default; Thus, the overall cost is not ${\cal O}(N^4)$ (as there are four indices) rather it is the factored $(\bf{A B}) \bf{C}$ which runs ${\cal O}(N^3)$.
#
# Therefore you do not need to factorize the expression yourself (sometimes you might need):

println("@tensor factorized time:")
@btime @tensor begin
   tmp[i,k]  := A[i,j] * B[j,k]
   tmp2[i,l] := tmp[i,k] * C[k,l]
end
nothing

# On most machines the three have similar timings. The BLAS usage is usually recommended. Thankfully, in Julia its syntax is very clear. The Psi4Julia project tends to lean toward usage of tensor packages but if Julia's built-in matrix multiplication is significantly cleaner/faster we would use it. The real value of tensor packages will become tangible for more complicated expressions.

# ## Complicated tensor manipulations
# Let us consider a popular index transformation example:
# $$M_{pqrs} = C_{pi} C_{qj} I_{ijkl} C_{rk} C_{sl}$$
#
# Here, a naive loop implementation would scale like $\mathcal{O}(N^8)$ which translates to an extremely costly computation for all but the smallest $N$. A smarter implementation (factorizing the whole expression) would scale
# like $\mathcal{O}(N^5)$.

# **WARNING: First execution is slow because of compilation time. Successive are more honest to the running time.**

# +
# Grab orbitals
dims = 15
@assert dims <= 15 || "Size must be smaller than 15."

C = rand(dims, dims)
I = rand(dims, dims, dims, dims)

# @tensor full transformation.
print("\nStarting @tensor full transformation...")
n8_time = @elapsed @tensor MO_n8[I,J,K,L] := C[p,I] * C[q,J] * I[p,q,r,s] * C[r,K] * C[s,L]
print("complete in $n8_time s\n")

# @tensor factorized N^5 transformation.
print("\nStarting @tensor factorized N^5 transformation with einsum ... ")
n5_time = @elapsed @tensor begin
   MO_n5[A,q,r,s] := C[p,A] * I[p,q,r,s]
   MO_n5[A,B,r,s] := C[q,B] * MO_n5[A,q,r,s]
   MO_n5[A,B,C,s] := C[r,C] * MO_n5[A,B,r,s]
   MO_n5[A,B,C,D] := C[s,D] * MO_n5[A,B,C,s]
end
print("complete in $n5_time s \n")
println("  @tensor factorized is $(n8_time/n5_time) faster than full @tensor algorithm!")
println("  Allclose? ", np.allclose(MO_n8, MO_n5))

# Julia's GEMM N^5 transformation.
# Try to figure this one out!
print("\nStarting Julia's factorized transformation with * ... ")
dgemm_time = @elapsed begin
   MO = C' * reshape(I, dims, :)
   MO = reshape(MO, :, dims) * C
   MO = permutedims(reshape(MO, dims, dims, dims, dims), (2, 1, 4, 3))

   MO = C' * reshape(MO, dims, :)
   MO = reshape(MO, :, dims) * C
   MO = permutedims(reshape(MO, dims, dims, dims, dims),(2, 1, 4, 3))
end
print("complete in $dgemm_time s \n")
println("  * factorized is $(n8_time/dgemm_time) faster than full @tensor algorithm!")
println("  Allclose? ", np.allclose(MO_n8, MO))

# There are still several possibilities to explore:
# @inbounds, @simd, LinearAlgebra.LAPACK calls, Einsum.jl, Tullio.jl, ...
# -
# None of the above algorithms is $\mathcal{O}(N^8)$. `@tensor` factorizes the expression to achieve better performance. There is a small edge in doing the factorization manually. Factorized algorithms have similar timings, although it is clear that with `@tensor` is easier than with Julia's built-in `*`. To use the usual matrix multiplication with tensors we have to reshape and permute their dimensions, subtracting appeal to the simple `*` syntax.


text:
# code generation routines
module Codegen
using ..kRPC.krpc.schema
using ..kRPC
using ProtoBuf

abstract type ProcedureType end
struct Standard <: ProcedureType
    name::String
end
struct ServiceGetter <: ProcedureType
    name::String
end
struct ServiceSetter <: ProcedureType
    name::String
end
struct ClassGetter <: ProcedureType
    clss::String
    name::String
end
struct ClassSetter <: ProcedureType
    clss::String
    name::String
end
struct StaticMember <: ProcedureType
    clss::String
    name::String
end
struct Member <: ProcedureType
    clss::String
    name::String
end

function sanitizeName(kName::Symbol)
    if kName == :String
        return :kString
    elseif kName == :Bool
        return :kBool
    end
    return kName
end

function makeparam(kName::String)
    kName = Symbol(kName)
    if kName == :type
        return :ktype
    elseif kName == :begin
        return :kbegin
    elseif kName == :end
        return :kend
    elseif kName == :function
        return :kfunction
    end
    return kName
end

struct SymbolContext
    current::Symbol # the service we're currently in
    needed::Base.Set{Symbol} # services required to be imported by the current service
    neededEnums::Base.Set{Tuple{Symbol, Symbol}} # service/enum combinations needed
end

function kRPCToJuliaType(kType::_Type,ctx::SymbolContext,delayed::Bool)
    kcode = kType.code
    types = kType.types
    if kcode == Type_TypeCode.NONE # ?
        return Nothing()
    elseif kcode == Type_TypeCode.DOUBLE return :Float64
    elseif kcode == Type_TypeCode.FLOAT return :Float32
    elseif kcode == Type_TypeCode.SINT32 return :Int32
    elseif kcode == Type_TypeCode.SINT64 return :Int64
    elseif kcode == Type_TypeCode.UINT32 return :UInt32
    elseif kcode == Type_TypeCode.UINT64 return :UInt64
    elseif kcode == Type_TypeCode.BOOL return :Bool
    elseif kcode == Type_TypeCode.STRING return :String
    elseif kcode == Type_TypeCode.BYTES return :(Array{UInt8, 1})
    elseif kcode == Type_TypeCode.CLASS
        sname = Symbol(kType.service)
        cname = sanitizeName(Symbol(kType.name))
        if sname == ctx.current && !delayed
            return cname
        end
        push!(ctx.needed, sname)
        return :($sname.$cname)
    elseif kcode == Type_TypeCode.ENUMERATION
        sname = Symbol(kType.service)
        emname = Symbol("Enum_" * kType.name)
        ename = Symbol("E"*kType.name)
        push!(ctx.neededEnums, (sname, emname))
        if sname == ctx.current
            return :($emname.$ename)
        end
        push!(ctx.needed, sname)
        return :($sname.$emname.$ename)
    elseif kcode == Type_TypeCode.EVENT
        return :kEvent
    elseif kcode == Type_TypeCode.PROCEDURE_CALL
        return :(kRPC.krpc.schema.ProcedureCall)
    elseif kcode == Type_TypeCode.STREAM
        return :kStream
    elseif kcode == Type_TypeCode.STATUS
        return :(Service)
    elseif kcode == Type_TypeCode.SERVICES
        return :(Services)
    elseif kcode == Type_TypeCode.TUPLE
        return quote Tuple{$(map(x-> kRPCToJuliaType(x,ctx,delayed), types)...)} end
    elseif kcode == Type_TypeCode.LIST
        return quote Array{$(kRPCToJuliaType(types[1],ctx,delayed)),1} end
    elseif kcode == Type_TypeCode.SET
        return quote Base.Set{$(kRPCToJuliaType(types[1],ctx,delayed))} end
    elseif kcode == Type_TypeCode.DICTIONARY
        kty = kRPCToJuliaType(types[1],ctx,delayed)
        vty = kRPCToJuliaType(types[2],ctx,delayed)
        return quote Dict{$kty, $vty} end
    end
end

function getArgument(idx::Int, arg)
    return Argument(position=idx, value=arg)
end

function generateRemoteCallGenerator(service::AbstractString, proc::Procedure, jlArgs, rtype::Union{Symbol, Expr, Nothing})
    return :(ProcedureCall(service=$service, procedure=$(proc.name), arguments=$jlArgs))
end

function generateRemoteCall(service::AbstractString, proc::Procedure, jlArgs, rtype::Union{Symbol, Expr, Nothing})
    postamble = Array{Expr,1}()
    if typeof(rtype) != Nothing
        push!(postamble, quote return kRPC.getJuliaValue(result.results[1].value,$rtype) end)
    end
    return Expr(:block,
        Expr(:line,1,Symbol(service * ":" * proc.name)),
        quote request = [$(generateRemoteCallGenerator(service, proc, jlArgs, rtype))] end,
        #quote println(request) end,
        quote result = kRPC.SendBiMessage(kRPC.Remote.conn, Request(calls=request)) end,
        quote kRPC.handle_potential_error(result) end,
        postamble...)
end

function generateRemoteCallDelayed(service::AbstractString, proc::Procedure, jlArgs, rtype::Union{Symbol, Expr, Nothing}, darg)
    return quote
        return kRPC.kPC($(generateRemoteCallGenerator(service, proc, jlArgs, rtype)), $rtype)
    end
end

twopartTypeR = r"^([^_]+)_([^_]+)$"
threepartTypeR = r"^([^_]+)_([^_]+)_([^_]+)$"

function parseKind(name::String) :: ProcedureType
    if occursin(twopartTypeR, name)
        captures = match(twopartTypeR, name).captures
        if captures[1] == "get"
            return ServiceGetter(captures[2])
        elseif captures[1] == "set"
            return ServiceSetter(captures[2])
        else
            return Member(captures[1], captures[2])
        end
    elseif occursin(threepartTypeR, name)
        captures = match(threepartTypeR, name).captures
        if captures[2] == "static"
            return StaticMember(captures[1], captures[3])
        elseif captures[2] == "get"
            return ClassGetter(captures[1], captures[3])
        elseif captures[2] == "set"
            return ClassSetter(captures[1], captures[3])
        else
            throw("Unexpected method name $name encountered")
        end
    else
        return Standard(name)
    end
end

function makeParameter(i::Int, param)
    if length(param.default_value) > 0
        quote
            if !isa($(Symbol(param.name)), Nothing)
                push!(args, Argument(position=$(i-1), value=kRPC.getWireValue($(makeparam(param.name)))))
            end
        end
    else
        quote
            push!(args, Argument(position=$(i-1), value=kRPC.getWireValue($(makeparam(param.name)))))
        end
    end
end

function generateProcedureCall(service::String, proc::Procedure, ctx::SymbolContext)
    kind = parseKind(proc.name)
    if isa(kind, Standard) || isa(kind, ServiceGetter) || isa(kind, ClassGetter) || isa(kind, StaticMember) || isa(kind, Member)
        name = kind.name
    elseif isa(kind, ServiceSetter) || isa(kind, ClassSetter)
        name = kind.name * "!"
    end
    name = sanitizeName(Symbol(name))

    exprargs = Any[ :($(makeparam(param.name)) :: $(kRPCToJuliaType(param._type,ctx,false))) for param in filter(p -> length(p.default_value) == 0, proc.parameters)]
    dexprargs = Any[ :($(makeparam(param.name)) :: $(kRPCToJuliaType(param._type,ctx,true))) for param in filter(p -> length(p.default_value) == 0, proc.parameters)]

    optargs = Any[ Expr(:kw,Expr(:(::),makeparam(param.name),:(Union{Nothing,$(kRPCToJuliaType(param._type,ctx,false))})),Nothing()) for param in filter(p -> length(p.default_value) > 0, proc.parameters)]
    doptargs = Any[ Expr(:kw,Expr(:(::),makeparam(param.name),:(Union{Nothing,$(kRPCToJuliaType(param._type,ctx,true))})),Nothing()) for param in filter(p -> length(p.default_value) > 0, proc.parameters)]
    if length(optargs) > 0
        pushfirst!(exprargs, Expr(:parameters, optargs...))
        pushfirst!(dexprargs, Expr(:parameters, doptargs...))
    end

    if ProtoBuf.has_field(proc, :return_type)
        retty = kRPCToJuliaType(proc.return_type,ctx,false)
        if proc.return_is_nullable
            retty = :(Nullable{$retty})
        end
    else
        retty = Nothing()
    end

    headargs = copy(exprargs)
    dheadargs = copy(dexprargs)
    fnhead = Expr(:call, pushfirst!(headargs, name) ...)
    fnheadd = Expr(:call, pushfirst!(dheadargs, name) ...)
    paramgen = quote
        args = Array{kRPC.krpc.schema.Argument,1}();
        $([makeParameter(i, proc.parameters[i]) for i=1:length(proc.parameters)]...);
        args
    end
    fnexpr = :($(Expr(:function, fnhead, generateRemoteCall(service, proc, paramgen, retty))); @doc $(proc.documentation) $(fnhead))
    fndexpr = Expr(:function, fnheadd,
        generateRemoteCallDelayed(service, proc, paramgen, retty, :context))
    return (:($fnexpr; export $name), :($fndexpr; export $name))
end

function makeService(name::String, body::Array{Expr,1}, delayed::Array{Expr, 1})
    return quote
        module $(Symbol(name))
            println($(name))
            $(body...)
            module Delayed
                $(delayed...)
            end
        end;
        #@doc $(service.documentation) $(Symbol(name))
    end
end

function collapse_modules(expr)
    if !(expr isa Expr) return [expr] end
    if expr.head == :block
        return collect(Iterators.flatten(map(collapse_modules, expr.args)))
    elseif expr.head == :module
        internal = collapse_modules(expr.args[3])
        return [Expr(:module, expr.args[1], expr.args[2], Expr(:block, internal...))]
    end
    return [expr]
end

function generateHelpers(info::Services,conn::kRPC.kRPCConnection)
    addons = Array{Expr,1}()
    services = Dict{String, Array{Expr, 1}}()
    dservices = Dict{String, Array{Expr, 1}}()
    for service in info.services
        services[service.name] = Array{Expr,1}()
        dservices[service.name] = Array{Expr,1}()
        push!(services[service.name], :(using kRPC; using kRPC.krpc.schema; using ..Remote))
        push!(dservices[service.name], :(using kRPC; using kRPC.krpc.schema; using ..Remote))
    end
    for service in info.services
        for cls in service.classes
            clsname = sanitizeName(Symbol(cls.name))
            push!(services[service.name], :(struct $clsname <: kRPC.kRPCTypes.Class handle::Base.UInt64 end;
                                            @doc $(cls.documentation) $clsname;
                                            export $clsname))
        end
        for enu in service.enumerations
            enumarr = Array{Expr,1}()
            enum_mod_name = Symbol("Enum_" * enu.name)
            enum_name = Symbol("E" * enu.name)
            push!(enumarr,:(abstract type $enum_name <: kRPC.kRPCTypes.Enum end; export $enum_name))
            for val in enu.values
                val_name,val_value = (Symbol(val.name), val.value)

                push!(enumarr,:(struct $val_name <: $enum_name end; export $val_name))
                push!(addons,:(function getEnumVal(::$(Symbol(service.name)).$enum_mod_name.$val_name) return $val_value end))
            end
            push!(services[service.name],quote module $enum_mod_name using ..kRPC; $(enumarr...) end end)
        end
        symbolContext = SymbolContext(Symbol(service.name), Base.Set{Symbol}(), Base.Set{Tuple{Symbol, Symbol}}())
        for proc in service.procedures
            sproc,dproc = generateProcedureCall(service.name, proc, symbolContext)
            push!(services[service.name], sproc)
            push!(dservices[service.name], dproc)
        end
        for needed in symbolContext.needed
            pushfirst!(services[service.name], Expr(:using, Expr(:., :(.), :(.), needed)))
            pushfirst!(dservices[service.name], Expr(:using, Expr(:.,:(.), :(.), :(.), needed)))
        end
        for needed in symbolContext.neededEnums
            pushfirst!(dservices[service.name], Expr(:import, Expr(:.,:(.), :(.), needed[2])))
        end
        pushfirst!(dservices[service.name], Expr(:using, Expr(:.,:(.), :(.), Symbol(service.name))))
    end

    opexpr = Array{Expr,1}()
    push!(opexpr, makeService("SpaceCenter", services["SpaceCenter"], dservices["SpaceCenter"]))
    append!(opexpr, [if first(p) == "SpaceCenter"
            :(begin end)
        else
            makeService(first(p), p[2], dservices[first(p)])
        end for p in services])
    #println(opexpr)
    expr = :(module Remote
                using ..kRPC
                function change_connection(iconn::kRPCConnection)
                    global conn = iconn
                end
                $(opexpr...)
                $(addons...)
            end)
    return collapse_modules(expr)[1]
end
end

ext: jl


text:
#modifying TTVFaster functions for efficiency

#include("../TTVFaster-noah/Julia/compute_ttv.jl")
#include("alt_ttv_succinct.jl")
#using ForwardDiff
#import TTVFaster.Planet_plane_hk

# Computes transit-timing variations to linear order in
# eccentricity for non-resonant, plane-parallel planets.
# Input:
#   jmax:  Maximum j over which to sum the TTV calculation for both planets
#     p1:  Planet type for inner planet
#     p2:  Planet type for outer planet
#  time1:  Transit times for inner planet
#  time2:  Transit times for outer planet
# alpha0:  Initial alpha for Taylor expansion of coefficients
#     b0:  Laplace coefficients and derivatives for use in Taylor expansion
# Output:
#   ttv1: TTVs of the inner planet
#   ttv2: TTVs of the outer planet
#     f1: TTV coefficients for inner planet
#     f2: TTV coefficients for outer planet
#      b: Laplace coefficients (& derivatives) for outer planet

function compute_inner_ttv!{T1<:Number}(jmax::Integer, p1::Planet_plane_hk{T1}, p2::Planet_plane_hk{T1}, time1::Vector, ttv1::Vector,f1::Array,f2::Array,b::Array,alpha0::Number,b0::Array)
#works when tested
# Compute the semi-major axis ratio of the planets:
const alpha = abs(p1.period/p2.period)^(2//3)  # Julia supports rational numbers!
# Number of times:
const ntime1 = length(time1)
#const ntime2 = length(time2)
# Compute the coefficients:
#if typeof(alpha)==ForwardDiff.Dual
#  alpha=alpha.value
#end
ttv_succinct!(jmax+1,alpha,f1,f2,b,alpha0,b0)  # I need to compute coefficients one higher than jmax
# Compute TTVs for inner planet (equation 33):
# Compute since of \pomegas:
e1 = sqrt(p1.esinw*p1.esinw+p1.ecosw*p1.ecosw)
e2 = sqrt(p2.esinw*p2.esinw+p2.ecosw*p2.ecosw)
sin1om=p1.esinw/e1
sin2om=p2.esinw/e2
cos1om=p1.ecosw/e1
cos2om=p2.ecosw/e2
# Compute mean motions:
n1=2pi/p1.period
n2=2pi/p2.period
# Compute initial longitudes:
lam10=-n1*p1.trans0 + 2*p1.esinw # 2*p1.eccen*sin1om
lam20=-n2*p2.trans0 + 2*p2.esinw # 2*p2.eccen*sin2om
@inbounds for i=1:ntime1
# Compute the longitudes of the planets at times of transit of planet 1 (equation 49):
  lam11 = n1*time1[i]+lam10
  lam21 = n2*time1[i]+lam20
  psi1  = lam11-lam21 # Compute difference in longitudes at times of transit of planet 1
  sinpsi1=sin(psi1)
  cospsi1=cos(psi1)
  sinlam11 = sin(lam11)
  coslam11 = cos(lam11)
  sinlam1om1=sinlam11*cos1om-coslam11*sin1om
  coslam1om1=coslam11*cos1om+sinlam11*sin1om
  sinlam1om2=sinlam11*cos2om-coslam11*sin2om
  coslam1om2=coslam11*cos2om+sinlam11*sin2om
  ttv1[i]=zero(p1.period) #0.0
  sinjm1psi1=zero(p1.period) #0.0
  cosjm1psi1=one(p1.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi1=sinjm1psi1*cospsi1+cosjm1psi1*sinpsi1
    cosjpsi1=cosjm1psi1*cospsi1-sinjm1psi1*sinpsi1
    ttv1[i] += f1[j+1,1]*sinjpsi1
    ttv1[i] += f1[j+1,2]*e1*(sinjpsi1*coslam1om1-cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j+1,3]*e1*(sinjpsi1*coslam1om1+cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j  ,4]*e2*(sinjpsi1*coslam1om2-cosjpsi1*sinlam1om2)
    ttv1[i] += f1[j+2,5]*e2*(sinjpsi1*coslam1om2+cosjpsi1*sinlam1om2)
    sinjm1psi1=sinjpsi1
    cosjm1psi1=cosjpsi1
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv1[i] = ttv1[i]*p1.period*p2.mass_ratio/(2pi)
end
# Finished!
return
end

function compute_outer_ttv!{T1<:Number}(jmax::Integer,p1::Planet_plane_hk{T1},p2::Planet_plane_hk{T1},time2::Vector,ttv2::Vector,f1::Array,f2::Array,b::Array,alpha0::Number,b0::Array)
#works when tested
# Compute the semi-major axis ratio of the planets:
# println(p1.period,p2.period)
const alpha = abs(p1.period/p2.period)^(2//3)  # Julia supports rational numbers!
# Number of times:
#const ntime1 = length(time1)
const ntime2 = length(time2)
# Compute the coefficients:
ttv_succinct!(jmax+1,alpha,f1,f2,b,alpha0,b0)  # I need to compute coefficients one higher than jmax
# Compute TTVs for inner planet (equation 33):
# Compute since of \pomegas:
e1 = sqrt(p1.esinw*p1.esinw+p1.ecosw*p1.ecosw)
e2 = sqrt(p2.esinw*p2.esinw+p2.ecosw*p2.ecosw)
sin1om=p1.esinw/e1
sin2om=p2.esinw/e2
cos1om=p1.ecosw/e1
cos2om=p2.ecosw/e2
# Compute mean motions:
n1=2pi/p1.period
n2=2pi/p2.period
# Compute initial longitudes:
lam10=-n1*p1.trans0 + 2*p1.esinw # 2*p1.eccen*sin1om
lam20=-n2*p2.trans0 + 2*p2.esinw # 2*p2.eccen*sin2om

# Compute TTVs for outer planet (equation 33):
@inbounds for i=1:ntime2
# Compute the longitudes of the planets at times of transit of planet 2:
  lam12 = n1*time2[i]+lam10
  lam22 = n2*time2[i]+lam20
  sinlam22 = sin(lam22)
  coslam22 = cos(lam22)
  psi2  = lam12-lam22 # Compute difference in longitudes at times of transit of planet 2
  sinpsi2=sin(psi2)
  cospsi2=cos(psi2)
  sinlam2om1=sinlam22*cos1om-coslam22*sin1om
  coslam2om1=coslam22*cos1om+sinlam22*sin1om
  sinlam2om2=sinlam22*cos2om-coslam22*sin2om
  coslam2om2=coslam22*cos2om+sinlam22*sin2om
  ttv2[i]=zero(p2.period) #0.0
  sinjm1psi2=zero(p2.period) #0.0
  cosjm1psi2=one(p2.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi2=sinjm1psi2*cospsi2+cosjm1psi2*sinpsi2
    cosjpsi2=cosjm1psi2*cospsi2-sinjm1psi2*sinpsi2
    ttv2[i] += f2[j+1,1]*sinjpsi2
    ttv2[i] += f2[j+1,2]*e2*(sinjpsi2*coslam2om2-cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+1,3]*e2*(sinjpsi2*coslam2om2+cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+2,4]*e1*(sinjpsi2*coslam2om1-cosjpsi2*sinlam2om1)
    ttv2[i] += f2[j  ,5]*e1*(sinjpsi2*coslam2om1+cosjpsi2*sinlam2om1)
    sinjm1psi2=sinjpsi2
    cosjm1psi2=cosjpsi2
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv2[i] = ttv2[i]*p2.period*p1.mass_ratio/(2pi)
end
# Finished!
return
end

max_stars_repo_name: nwtuchow/TTVFaster

max_issues_repo_name: nwtuchow/TTVFaster

max_forks_repo_name: nwtuchow/TTVFaster


text:
using LightGraphs: AbstractGraph

## Linear algebra API for AbstractGraph

function degrees(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = adjacency_matrix(g, T; dir=dir)
    degrees(adj, T; dir=dir)
end

function degree_matrix(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = adjacency_matrix(g, T; dir=dir)
    degree_matrix(adj, T; dir=dir)
end

function inv_sqrt_degree_matrix(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = adjacency_matrix(g, T; dir=dir)
    inv_sqrt_degree_matrix(adj, T; dir=dir)
end

function laplacian_matrix(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = adjacency_matrix(g, T; dir=dir)
    laplacian_matrix(adj, T; dir=dir)
end

function normalized_laplacian(g::AbstractGraph, T::DataType=eltype(g); selfloop::Bool=false)
    adj = adjacency_matrix(g, T)
    selfloop && (adj += I)
    normalized_laplacian(adj, T)
end

function scaled_laplacian(g::AbstractGraph, T::DataType=eltype(g))
    adj = adjacency_matrix(g, T)
    scaled_laplacian(adj, T)
end


text:
using Flux, Zygote, LinearAlgebra, Statistics
println("NNPDEHAN_tests")
using Test, StochasticDiffEq
println("Starting Soon!")
using NeuralPDE

using Random
Random.seed!(100)

# one-dimensional heat equation
x0 = [11.0f0]  # initial points
tspan = (0.0f0,5.0f0)
dt = 0.5   # time step
time_steps = div(tspan[2]-tspan[1],dt)
d = 1      # number of dimensions
m = 10     # number of trajectories (batch size)
batch_size = 8192
train_steps = 400
lr_boundaries = [100, 200, 300]
lr_values = [1e-1,1e-2,1e-3,1e-4]

g(X) = sum(X.^2)   # terminal condition
f(y,z,v_y,v_z,dv_y,dv_z,t) = sin(v_y) - v_z * π^(d/2)*σ_sampling^2 # function from solved equation
μ_f(X,p,t) = 0.0
σ_f(X,p,t) = 1.0
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = d + 50 #hidden layer size
#neural network approximating solutions at the desired point
# with batch normalisation
# u0 = Flux.Chain(Dense(d,hls),
#                 BatchNorm(hls,tanh),
#                 Dense(hls,hls,tanh),
#                 BatchNorm(hls,tanh),
#                 Dense(hls,1))
# without batch normalisation
u0 = Flux.Chain(Dense(d,hls,tanh),
                Dense(hls,hls,tanh),
                Dense(hls,1))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random, LinearAlgebra

@testset "middle" begin
    @test middle(3) === 3.0
    @test middle(2, 3) === 2.5
    let x = ((realmax(1.0)/4)*3)
        @test middle(x, x) === x
    end
    @test middle(1:8) === 4.5
    @test middle([1:8;]) === 4.5

    # ensure type-correctness
    for T in [Bool,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,Float16,Float32,Float64]
        @test middle(one(T)) === middle(one(T), one(T))
    end
end

@testset "median" begin
    @test median([1.]) === 1.
    @test median([1.,3]) === 2.
    @test median([1.,3,2]) === 2.

    @test median([1,3,2]) === 2.0
    @test median([1,3,2,4]) === 2.5

    @test median([0.0,Inf]) == Inf
    @test median([0.0,-Inf]) == -Inf
    @test median([0.,Inf,-Inf]) == 0.0
    @test median([1.,-1.,Inf,-Inf]) == 0.0
    @test isnan(median([-Inf,Inf]))

    X = [2 3 1 -1; 7 4 5 -4]
    @test all(median(X, dims=2) .== [1.5, 4.5])
    @test all(median(X, dims=1) .== [4.5 3.5 3.0 -2.5])
    @test X == [2 3 1 -1; 7 4 5 -4] # issue #17153

    @test_throws ArgumentError median([])
    @test isnan(median([NaN]))
    @test isnan(median([0.0,NaN]))
    @test isnan(median([NaN,0.0]))
    @test isequal(median([NaN 0.0; 1.2 4.5], dims=2), reshape([NaN; 2.85], 2, 1))

    @test median!([1 2 3 4]) == 2.5
    @test median!([1 2; 3 4]) == 2.5

    @test invoke(median, Tuple{AbstractVector}, 1:10) == median(1:10) == 5.5
end

@testset "mean" begin
    @test_throws ArgumentError mean(())
    @test mean((1,2,3)) === 2.
    @test mean([0]) === 0.
    @test mean([1.]) === 1.
    @test mean([1.,3]) == 2.
    @test mean([1,2,3]) == 2.
    @test mean([0 1 2; 4 5 6], dims=1) == [2.  3.  4.]
    @test mean([1 2 3; 4 5 6], dims=1) == [2.5 3.5 4.5]
    @test mean(i->i+1, 0:2) === 2.
    @test mean(isodd, [3]) === 1.
    @test mean(x->3x, (1,1)) === 3.

    @test isnan(mean([NaN]))
    @test isnan(mean([0.0,NaN]))
    @test isnan(mean([NaN,0.0]))

    @test isnan(mean([0.,Inf,-Inf]))
    @test isnan(mean([1.,-1.,Inf,-Inf]))
    @test isnan(mean([-Inf,Inf]))
    @test isequal(mean([NaN 0.0; 1.2 4.5], dims=2), reshape([NaN; 2.85], 2, 1))

    # Check that small types are accumulated using wider type
    for T in (Int8, UInt8)
        x = [typemax(T) typemax(T)]
        g = (v for v in x)
        @test mean(x) == mean(g) == typemax(T)
        @test mean(identity, x) == mean(identity, g) == typemax(T)
        @test mean(x, dims=2) == [typemax(T)]'
    end
end

@testset "quantile" begin
    @test quantile([1,2,3,4],0.5) == 2.5
    @test quantile([1,2,3,4],[0.5]) == [2.5]
    @test quantile([1., 3],[.25,.5,.75])[2] == median([1., 3])
    @test quantile(100.0:-1.0:0.0, 0.0:0.1:1.0) == 0.0:10.0:100.0
    @test quantile(0.0:100.0, 0.0:0.1:1.0, sorted=true) == 0.0:10.0:100.0
    @test quantile(100f0:-1f0:0.0, 0.0:0.1:1.0) == 0f0:10f0:100f0
    @test quantile([Inf,Inf],0.5) == Inf
    @test quantile([-Inf,1],0.5) == -Inf
    @test quantile([0,1],1e-18) == 1e-18
    @test quantile([1, 2, 3, 4],[]) == []
    @test quantile([1, 2, 3, 4], (0.5,)) == (2.5,)
    @test quantile([4, 9, 1, 5, 7, 8, 2, 3, 5, 17, 11], (0.1, 0.2, 0.4, 0.9)) == (2.0, 3.0, 5.0, 11.0)
    @test quantile([1, 2, 3, 4], ()) == ()
end

# StatsBase issue 164
let y = [0.40003674665581906, 0.4085630862624367, 0.41662034698690303, 0.41662034698690303, 0.42189053966652057, 0.42189053966652057, 0.42553514344518345, 0.43985732442991354]
    @test issorted(quantile(y, range(0.01, stop=0.99, length=17)))
end

@testset "Issue #17153 and PR #17154" begin
    a = rand(10,10)
    b = copy(a)
    x = median(a, dims=1)
    @test b == a
    x = median(a, dims=2)
    @test b == a
    x = mean(a, dims=1)
    @test b == a
    x = mean(a, dims=2)
    @test b == a
end

# dimensional correctness
isdefined(Main, :TestHelpers) || @eval Main include("TestHelpers.jl")
using .Main.TestHelpers: Furlong
@testset "Unitful elements" begin
    r = Furlong(1):Furlong(1):Furlong(2)
    a = Vector(r)
    @test sum(r) == sum(a) == Furlong(3)
    @test cumsum(r) == Furlong.([1,3])
    @test mean(r) == mean(a) == median(a) == median(r) == Furlong(1.5)

    # Issue #21786
    A = [Furlong{1}(rand(-5:5)) for i in 1:2, j in 1:2]
    @test mean(mean(A, dims=1), dims=2)[1] === mean(A)
end

# Issue #22901
@testset "quantile of Any arrays" begin
    x = Any[1, 2, 4, 10]
    y = Any[1, 2, 4, 10//1]
    @test quantile(x, 0.5)  === 3.0
    @test quantile(x, 1//2) === 3//1
end

@testset "Mean along dimension of empty array" begin
    a0  = zeros(0)
    a00 = zeros(0, 0)
    a01 = zeros(0, 1)
    a10 = zeros(1, 0)
    @test isequal(mean(a0, dims=1)      , fill(NaN, 1))
    @test isequal(mean(a00, dims=(1, 2)), fill(NaN, 1, 1))
    @test isequal(mean(a01, dims=1)     , fill(NaN, 1, 1))
    @test isequal(mean(a10, dims=2)     , fill(NaN, 1, 1))
end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
export Lazy, eachdim′, eachcol′, eachrow′, eachslice′

### import many funtions
import Base.Broadcast: throwdm, AbstractArrayStyle, Unknown, combine_eltypes,
    preprocess, Broadcasted, DefaultArrayStyle, ischunkedbroadcast, chunkedcopyto!, bitcache_size,
    dumpbitcache, bitcache_chunks, materialize!, BroadcastStyle, combine_styles, instantiate,
    broadcastable, broadcast_unalias, Style, _broadcast_getindex, broadcasted
import Base: size, axes, setindex!, unalias, mightalias, unaliascopy, IndexStyle, parent,
    unsafe_setindex!, copyto!, IndexStyle, tail

# getsame
@noinline _getsamethrow(@nospecialize(f)) =
    throw(ArgumentError("Inputs with different $f"))
@inline getsame(f::F, x) where {F} = f(x)
@inline getsame(f::F, x, y, zs...) where {F} = begin
    fx = getsame(f, y, zs...)
    f(x) == fx || _getsamethrow(f)
    fx
end

const FilledBC = Broadcasted{<:AbstractArrayStyle{0}}
module LazyCollect
import Base: @propagate_inbounds, getindex
import Base.Broadcast: BroadcastStyle, broadcastable, extrude, newindex, broadcasted,
    instantiate, Broadcasted, DefaultArrayStyle, combine_styles
import Base.Iterators: ProductIterator

export Lazy

ndims(x::Tuple) = 1
size(x::Tuple) = (length(x),)
ndims(x) = Base.ndims(x)
size(x) = Base.size(x)
# Like reshape, but allow tuple inputs.
struct FakeDim{N,S,T,D} <: AbstractArray{T,N}
    data::D
    FakeDim{N,S}(data::T) where {N,S,T} = new{N,S,eltype(data),T}(data)
end
FakeDim{N,S}(data::FakeDim) where {N,S} = FakeDim{N,S}(data.data)
Base.axes(id::FakeDim{N,S}) where {N,S} = (ntuple(Returns(Base.OneTo(1)), Val(S - 1))..., axes(id.data)...)
Base.size(id::FakeDim{N,S}) where {N,S} = (ntuple(Returns(1), Val(S - 1))..., size(id.data)...)
extrude(x::FakeDim) = x
struct FakeIndex{N}
    inds::NTuple{N,Int}
    FakeIndex(inds::NTuple{N,Int}) where {N} = new{N}(inds)
end
@inline pick_inds(I::NTuple{N,Int}, ::Val{S}, ::Val{E}) where {N,S,E} = ntuple(i -> I[S + i - 1], Val(E - S + 1))
@inline newindex(::FakeDim{N,S}, I::CartesianIndex) where {N,S} = FakeIndex(pick_inds(I.I, Val(S), Val(N)))
@propagate_inbounds getindex(id::FakeDim, I::FakeIndex{N}) where {N} = id.data[I.inds...]
@propagate_inbounds getindex(id::FakeDim{N,S}, I::Vararg{Int, N}) where {N,S} = id.data[I[S:N]...]
@inline BroadcastStyle(::Type{ID}) where {ID<:FakeDim} =
    ID.parameters[4] <: Tuple ? DefaultArrayStyle{ID.parameters[1]}() :
                                BroadcastStyle(ID.parameters[4])

@inline fakedim(pre_ax::Tuple, data) = begin
    data′ = broadcastable(data)
    ax = axes(data)
    m, n = length(pre_ax), length(ax)
    data″ = m == 0 || n == 0 ? data′ : FakeDim{m + n, m + 1}(data′)
    data″, (pre_ax..., ax...)
end
@inline fakedim(pre_ax::Tuple, g::Base.Generator) = begin
    iter, ax = fakedim(pre_ax, g.iter)
    Broadcasted(g.f, (iter,), ax), ax
end
@inline fakedim(pre_ax::Tuple, i::Iterators.ProductIterator) = begin
    iter, ax = fakedims(pre_ax, i.iterators...)
    Broadcasted(tuple, iter, ax), ax
end
@inline fakedims(pre_ax::Tuple) = (), pre_ax
@inline fakedims(pre_ax::Tuple, x, y...) = begin
    x′, pre_ax′ = fakedim(pre_ax, x)
    y′, ax = fakedims(pre_ax′, y...)
    (x′, y′...), ax
end
# Lazy is used to wrap Generator to avoid collect before broadcast.
Lazy(x) = fakedim((), x)[1]

end
using .LazyCollect
## better each
@inline unsafe_view(A, I...) = Base.unsafe_view(A, to_indices(A, I)...)
eachcol′(A::AbstractVecOrMat) = Lazy(unsafe_view(A, :, i) for i in axes(A, 2))
eachrow′(A::AbstractVecOrMat) = Lazy(unsafe_view(A, i, :) for i in axes(A, 1))
function eachdim′(A::AbstractArray; dim::Val{D}) where {D}
    D <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    axes_all = ntuple(d -> d == D ? Ref(:) : axes(A, d), ndims(A))
    Lazy(unsafe_view(A, i...) for i in Iterators.product(axes_all...))
end
function eachslice′(A::AbstractArray; dim::Val{D}) where {D}
    D <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    inds_before = ntuple(Returns(:), D - 1)
    inds_after = ntuple(Returns(:), ndims(A) - D)
    Lazy(unsafe_view(A, inds_before..., i, inds_after...) for i in axes(A, D))
end

## faster Base.read!
function Base.read!(s::IO, a::StridedArray{T}) where {T}
    _style(a) = ndims(a) < 2 ? IndexLinear() : IndexStyle(a)
    if isbitstype(T) && _style(a) isa Union{IndexLinear, Base.IndexSCartesian2} && isone(stride(a,1))
        GC.@preserve a unsafe_read(s, pointer(a), sizeof(a))
        return a
    end
    Base.@invoke read!(s::IO, a::AbstractArray)
end

function Base.copyto_unaliased!(::IndexLinear, dest::AbstractArray, ::IndexLinear, src::AbstractArray)
    isempty(src) && return dest
    length(dest) < length(src) && throw(BoundsError(dest, LinearIndices(src)))
    Δi = firstindex(dest) - firstindex(src)
    for i in eachindex(IndexLinear(), src)
        @inbounds dest[i + Δi] = src[i]
    end
    return dest
end

splitloop(iter::CartesianIndices) = begin
    indices = iter.indices
    CartesianIndices(tail(indices)), first(indices)
end

function Base.copyto_unaliased!(::IndexLinear, dest::AbstractArray, ::IndexCartesian, src::AbstractArray)
    isempty(src) && return dest
    length(dest) < length(src) && throw(BoundsError(dest, LinearIndices(src)))
    iter, j = CartesianIndices(src), firstindex(dest) - 1
    inner_len = size(iter, 1)
    if inner_len >= 16
        # manually expand the inner loop similar to @simd
        outer, inner = splitloop(iter)
        @inbounds for II in outer
            for i in inner
                dest[j += 1] = src[i, II]
            end
        end
    else
        @inbounds for I in iter
            dest[j += 1] = src[I]
        end
    end
    return dest
end

function Base.copyto_unaliased!(::IndexCartesian, dest::AbstractArray, ::IndexLinear, src::AbstractArray)
    isempty(src) && return dest
    length(dest) < length(src) && throw(BoundsError(dest, LinearIndices(src)))
    iter, i = CartesianIndices(dest), firstindex(src) - 1
    inner_len = size(iter, 1)
    if inner_len >= 16
        # manually expand the inner loop similar to @simd
        outer, inner = splitloop(iter)
        if length(dest) == length(src)
            @inbounds for II in outer
                for j in inner
                    dest[j, II] = src[i += 1]
                end
            end
        else
            final = lastindex(src)
            @inbounds for II in outer
                if i + inner_len >= final
                    j = first(inner)
                    while i < final
                        dest[j, II] = src[i += 1]
                        j += step(inner)
                    end
                    break
                end
                for j in inner
                    dest[j, II] = src[i += 1]
                end
            end
        end
    elseif length(dest) == length(src)
        @inbounds for I in iter
            dest[I] = src[i += 1]
        end
    else
        # use zip based interator
        @inbounds for (I, J) in zip(eachindex(src), iter)
            dest[J] = src[I]
        end
    end
    return dest
end

function Base.copyto_unaliased!(::IndexCartesian, dest::AbstractArray, ::IndexCartesian, src::AbstractArray)
    isempty(src) && return dest
    length(dest) < length(src) && throw(BoundsError(dest, LinearIndices(src)))
    iterdest, itersrc = CartesianIndices(dest), CartesianIndices(src)
    if iterdest == itersrc
        iter = itersrc
        inner_len = size(iter, 1)
        if inner_len >= 16
            # manually expand the inner loop similar to @simd
            outer, inner = splitloop(iter)
            @inbounds for II in outer
                for i in inner
                    dest[i, II] = src[i, II]
                end
            end
        else
            @inbounds for I in iter
                dest[I] = src[I]
            end
        end
    else
        @inbounds for (J, I) in zip(iterdest, itersrc)
            dest[J] = src[I]
        end
    end
    return dest
end


text:
using Zygote, LinearAlgebra
using KrylovKit: eigsolve

test_eigmin(A::AbstractMatrix) = eigsolve(A, 1, :SR)[1][1]

A = rand(4, 4); A = A + A';

@show eigmin(A)
@show test_eigmin(A)
@show Zygote.gradient(test_eigmin, A)

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
"""
	varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)

VARIMAX perform varimax (or quartimax, equamax, parsimax) rotation to the column vectors of the input matrix.

# Input Arguments
- `A::Matrix{Float64}`: input matrix, whose column vectors are to be rotated. d, m = size(A).
- `gamma::Float64`: default is 1. gamma = 0, 1, m/2, and d(m - 1)/(d + m - 2), corresponding to quartimax, varimax, equamax, and parsimax.
- `minit::Int`: default is 20. Minimum number of iterations, in case of the stopping criteria fails initially.
- `maxit::Int`: default is 1000. Maximum number of iterations.
- `reltol::Float64`: default is 1e-12. Relative tolerance for stopping criteria.

# Output Argument
- `B::Matrix{Float64}`: output matrix, whose columns are already been rotated.

Implemented by Haotian Li, Aug. 20, 2019
"""
function varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)
	# Get the sizes of input matrix
	d, m = size(A)

	# If there is only one vector, then do nothing.
	if m == 1
		return A
	end

	if d == m && rank(A) == d
		return Matrix{Float64}(I, d, m)
	end

	# Warm up step: start with a good initial orthogonal matrix T by SVD and QR
	T = Matrix{Float64}(I, m, m)
	B = A * T
	L,_,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
	T = L * M'
	if norm(T-Matrix{Float64}(I, m, m)) < reltol
		T,_ = qr(randn(m,m)).Q
		B = A * T
	end

	# Iteration step: get better T to maximize the objective (as described in Factor Analysis book)
	D = 0
	for k in 1:maxit
		Dold = D
		L,s,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
		T = L * M'
		D = sum(s)
		B = A * T
		if (abs(D - Dold)/D < reltol) && k >= minit
			break
		end
	end

	# Adjust the sign of each rotated vector such that the maximum absolute value is positive.
	for i in 1:m
		if abs(maximum(B[:,i])) < abs(minimum(B[:,i]))
			B[:,i] .= - B[:,i]
		end
	end

	return B
end


text:
include("src/A2.jl")
using .A2

using GRUtils

a = 0
b = 1
n = 1000

y, z = solve(
    Problem(
        (A) -> 1 / 3,
        (A) -> -A,
        (A) -> -2,
        (A) -> A,
        (x, z) -> sqrt(z^2 / 2),
        (x, z) -> - z / (z^2 + 1),
        a,
        b,
        -5,
        -4,
        4,
    ),
    Options(n),
)

open("result", "w") do io
    println(io, "y z")
    for i in 1:n+1
        println(io, y[i], " ", z[i])
    end
end

h = (b - a) / n
x = [a + i * h for i in 0:n]

plot(x, y)
plot(x, z, hold=true)
legend("y", "z")
savefig("figures/result.pdf")


text:
@testset "serialize" begin
    @testset "MVHistory" begin
        expected = MVHistory()
        push!(expected, :a, 1, 1.2)
        push!(expected, :a, 2, 2.2)
        push!(expected, :b, 1, 3)
        push!(expected, :b, 2, 4)
        push!(expected, :c, 1, [:inlier])
        push!(expected, :c, 2, [:outlier])
        push!(expected, :d, 1, ones(2, 1))
        push!(expected, :d, 2, zeros(2, 1))
        push!(expected, :e, 1, [1])
        push!(expected, :e, 2, [2])
        push!(expected, :f, 1, OneClassActiveLearning.ConfusionMatrix(1,1,1,1))
        push!(expected, :f, 2, OneClassActiveLearning.ConfusionMatrix(2,2,2,2))

        actual = Unmarshal.unmarshal(MVHistory, JSON.parse(JSON.json(expected)))

        for k in keys(expected)
            @test expected[k] == actual[k]
        end
    end

    @testset "al summary" begin
        expected = Dict(:m1 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]),
                           :m2 => Dict(:a => 1.0, :b => 1, :c => [1.0], :d => [1, 2]))

        actual = OneClassActiveLearning.unmarshal_al_summary(JSON.parse(JSON.json(expected)))

        for (k_metrics, v_metrics) in expected
            for (k_summary, v_summary) in v_metrics
                @test expected[k_metrics][k_summary] == actual[k_metrics][k_summary]
            end
        end
    end

    @testset "Result" begin
        id = 42
        experiment = Dict{Symbol, Any}(
                :hash => id,
                :data_file => TEST_DATA_FILE,
                :output_file => TEST_OUTPUT_FILE,
                :model => Dict( :type => :(SVDD.RandomOCClassifier),
                                :init_strategy => :(SVDD.FixedParameterInitialization(GaussianKernel(2), 0.5))),
                :query_strategy => Dict(:type => :(OneClassActiveLearning.QueryStrategies.RandomPQs), :param => Dict()),
                :split_strategy => OneClassActiveLearning.DataSplits(trues(123), OneClassActiveLearning.FullSplitStrat()),
                :param => Dict(:num_al_iterations => 5,
                               :solver => Dict(:type => TEST_SOLVER.optimizer_constructor,
                                               :flags => Dict(TEST_SOLVER.params)),
                               :initial_pools => fill(:U, 123)))

       al_history = MVHistory()
       push!(al_history, :a, 1, 1.2)
       push!(al_history, :a, 2, 2.2)
       push!(al_history, :b, 1, 3)
       push!(al_history, :b, 2, 4)
       push!(al_history, :query_labels, 1, [:inlier])
       push!(al_history, :query_labels, 2, [:outlier])
       push!(al_history, :query_history, 1, ones(2, 1))
       push!(al_history, :query_history, 2, zeros(2, 1))
       push!(al_history, :c, 1, OneClassActiveLearning.ConfusionMatrix(1,1,1,1))
       push!(al_history, :c, 2, OneClassActiveLearning.ConfusionMatrix(2,2,2,2))

        al_summary = Dict(:m1 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]),
                   :m2 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]))

        worker_info = OneClassActiveLearning.get_worker_info()
        expected = OneClassActiveLearning.Result(id, experiment, worker_info, DataStats(1, 1, 0.1, 0.1, [1,2,3,4], [5,6,7]), al_history, al_summary)

        OneClassActiveLearning.write_result_to_file(TEST_OUTPUT_FILE, expected)
        actual = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parsefile(TEST_OUTPUT_FILE))
        @test actual.id == 42
        @test length(actual.experiment[:param][:initial_pools]) == 123
        @test haskey(actual.worker_info, :hostname)
        @test haskey(actual.al_history, :c)
        @test haskey(actual.al_summary, :m1)
        @test haskey(actual.status, :exit_code)
    end

    @testset "Result-SubSVDD" begin
        id = 43
        subspaces = [[1,2], [3,4,5]]
        experiment = Dict{Symbol, Any}(
                :hash => id,
                :data_file => "$(@__DIR__)/../data/input/dummy.csv",
                :output_file => "OneClassActiveLearning.jl/data/output/scenarioA/data_qs_model_id.tmp",
                :model => Dict(:type => :SubSVDD,
                               :param => Dict{Symbol, Any}(:subspaces => subspaces,
                                                           :weight_update_strategy => SVDD.FixedWeightStrategy(10.0, 0.01)),
                               :init_strategy => SimpleSubspaceStrategy(FixedGammaStrategy([MLKernels.GaussianKernel(5), MLKernels.GaussianKernel(8)]),
                                                 FixedCStrategy(0.1),
                                                 gamma_scope=Val(:Subspace))),
                :query_strategy => Dict(:type => SubspaceQs{RandomPQs},
                                        :param => Dict{Symbol, Any}(:scale_fct => min_max_normalize,
                                                                    :combination_fct => max,
                                                                    :subspaces => subspaces)),
                :split_strategy => OneClassActiveLearning.DataSplits(trues(123), OneClassActiveLearning.FullSplitStrat()),
                :param => Dict(:num_al_iterations => 5,
                               :solver => Dict(:type => TEST_SOLVER.optimizer_constructor,
                                               :flags => Dict(TEST_SOLVER.params)),
                               :initial_pools => fill(:U, 124)))

       al_history = MVHistory()
       push!(al_history, :a, 1, 1.2)
       push!(al_history, :a, 2, 2.2)
       push!(al_history, :query_labels, 1, [:inlier])
       push!(al_history, :query_labels, 2, [:outlier])
       push!(al_history, :query_history, 1, [4])
       push!(al_history, :query_history, 2, [7])
       push!(al_history, :c, 1, OneClassActiveLearning.ConfusionMatrix(1,1,1,1))
       push!(al_history, :c, 2, OneClassActiveLearning.ConfusionMatrix(2,2,2,2))

        al_summary = Dict(:m1 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]),
                   :m2 => Dict(:a => 1.0, :b => 1, :c => [1.0, 2.0], :d => [1, 2]))

        worker_info = OneClassActiveLearning.get_worker_info()
        expected = OneClassActiveLearning.Result(id, experiment, worker_info, DataStats(1, 1, 0.1, 0.1, [1,2,3,4], [5,6,7]), al_history, al_summary)

        actual = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parse(JSON.json(expected)))
        @test actual.id == 43
        @test length(actual.experiment[:param][:initial_pools]) == 124
        @test haskey(actual.worker_info, :hostname)
        @test haskey(actual.al_history, :c)
        @test haskey(actual.al_summary, :m1)
        @test haskey(actual.status, :exit_code)
    end
end


text:
include("./shared.jl")
include("../scripts-part2/2D_SWE/shallow_water_2D_LF_xpu_mpi.jl")
include("../scripts-part2/2D_SWE/shallow_water_2D.jl")

# run the model for dam break in y-direction
H, xc = shallow_water_2D_xpu_mpi(; nx = n, ny = n, dam1D_x = false)
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(xc), 12)))
d = Dict(:X => xc[indsx], :H => H[indsx, indsy])

@testset "Ref-test shallow water 2D with dam break in y-direction" begin
      @test_reference "reftest-files/shallow_water_2D_y.bson" d by = comp
end

@reset_parallel_stencil()


text:
export Category, FreeCategory, Ob, Hom, dom, codom, id, compose, ⋅,
  Category2, FreeCategory2, Hom2, compose2

import Base: show

# Category
##########

""" Theory of *categories* (with no extra structure)

**Warning**: We compose functions from left to right, i.e., if f:A→B and g:B→C
then compose(f,g):A→C. Under this convention function are applied on the right,
e.g., if a∈A then af∈B.

We retain the usual meaning of the symbol ∘ (\\circ), i.e., g∘f = compose(f,g).
This usage is too entrenched to overturn, inconvenient though it may be.
We use symbol ⋅ (\\cdot) for diagrammatic composition: f⋅g = compose(f,g).
"""
@theory Category(Ob,Hom) begin
  # Unicode aliases.
  @op begin
    (→) := Hom
    (⋅) := compose
  end

  """ Object in a category """
  Ob::TYPE

  """ Morphism in a category """
  Hom(dom::Ob,codom::Ob)::TYPE

  id(A::Ob)::(A → A)
  compose(f::(A → B), g::(B → C))::(A → C) ⊣ (A::Ob, B::Ob, C::Ob)

  # Category axioms.
  ((f ⋅ g) ⋅ h == f ⋅ (g ⋅ h)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B), g::(B → C), h::(C → D)))
  f ⋅ id(B) == f ⊣ (A::Ob, B::Ob, f::(A → B))
  id(A) ⋅ f == f ⊣ (A::Ob, B::Ob, f::(A → B))
end

# Convenience constructors
compose(fs::Vector) = foldl(compose, fs)
compose(f, g, h, fs...) = compose([f, g, h, fs...])

@syntax FreeCategory(ObExpr,HomExpr) Category begin
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

function show_unicode(io::IO, expr::HomExpr{:compose}; kw...)
  Syntax.show_unicode_infix(io, expr, "⋅"; kw...)
end

function show_latex(io::IO, expr::HomExpr{:id}; kw...)
  Syntax.show_latex_script(io, expr, "\\mathrm{id}")
end
function show_latex(io::IO, expr::HomExpr{:compose}; paren::Bool=false, kw...)
  Syntax.show_latex_infix(io, expr, "\\cdot"; paren=paren)
end

function show(io::IO, ::MIME"text/latex", expr::GATExpr)
  print(io, "\$")
  show_latex(io, expr)
  print(io, "\$")
end

function show(io::IO, ::MIME"text/latex", expr::HomExpr)
  print(io, "\$")
  show_latex(io, expr)
  print(io, " : ")
  show_latex(io, dom(expr))
  print(io, " \\to ")
  show_latex(io, codom(expr))
  print(io, "\$")
end

# 2-category
############

""" Theory of (strict) *2-categories*
"""
@signature Category(Ob,Hom) => Category2(Ob,Hom,Hom2) begin
  """ 2-morphism in a 2-category """
  Hom2(dom::Hom(A,B), codom::Hom(A,B))::TYPE ⊣ (A::Ob, B::Ob)
  @op (⇒) := Hom2

  # Hom categories: Vertical composition
  id(f)::(f ⇒ f) ⊣ (A::Ob, B::Ob, f::(A ⇒ B))
  compose(α::(f ⇒ g), β::(g ⇒ h))::(f ⇒ h) ⊣
    (A::Ob, B::Ob, f::(A → B), g::(A → B), h::(A → B))

  # Horizontal compostion
  compose2(α::(f ⇒ g), β::(h ⇒ k))::((f ⋅ h) ⇒ (g ⋅ k)) ⊣
    (A::Ob, B::Ob, C::Ob, f::(A → B), g::(A → B), h::(B → C), k::(B → C))
end

# Convenience constructors
compose2(αs::Vector) = foldl(compose2, αs)
compose2(α, β, γ, αs...) = compose2([α, β, γ, αs...])

""" Syntax for a 2-category.

Checks domains of morphisms but not 2-morphisms.
"""
@syntax FreeCategory2(ObExpr,HomExpr,Hom2Expr) Category2 begin
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
  compose(α::Hom2, β::Hom2) = associate(new(α,β))
  compose2(α::Hom2, β::Hom2) = associate(new(α,β))
end

function show_unicode(io::IO, expr::Hom2Expr{:compose}; kw...)
  Syntax.show_unicode_infix(io, expr, "⋅"; kw...)
end
function show_unicode(io::IO, expr::Hom2Expr{:compose2}; kw...)
  Syntax.show_unicode_infix(io, expr, "*"; kw...)
end

function show_latex(io::IO, expr::Hom2Expr{:compose}; kw...)
  Syntax.show_latex_infix(io, expr, "\\cdot"; kw...)
end
function show_latex(io::IO, expr::Hom2Expr{:compose2}; kw...)
  Syntax.show_latex_infix(io, expr, "*"; kw...)
end


text:
using Stheno, Plots, Random
using Stheno: @model

###########################  Define and inspect our model  ###########################

#=
We wish to perform inference in some latent process `f`, but it is corrupted by both some
iid noise, and some non-iid noise. Specifically, ω is some periodic behaviour in which we
are uninterested. Similarly, g is also drawn from a GP with an EQ kernel, but with a much
shorter length-scale. We wish to use observations of
`f + ω + g + GP(ZeroMean{Float64}(), Noise(0.1))`
to infer `f`. In this context, `ω + g + GP(ZeroMean{Float64}(), Noise(0.1))` is
regarded as our "noise process".

The described model is unidentiable. Conditioning on an observation of the periodic process
at some point, however, allows us to infer a particular `f`. Use the `model_identifiable`
function rather than `model` to see the results of this. Note the example of sequential
conditioning.
=#
@model function model()
    f = GP(EQ())
    ω = GP(periodic(EQ(), 1))
    return f, ω, ω + f, ω + f + GP(Noise(0.1))
end
@model function model_identifiable()
    f = GP(EQ())
    t = GP(periodic(EQ(), 1))
    ω = t | (t([0])←[0.0])
    return f, ω, ω + f, ω + f + GP(Noise(0.001))
end

# Select some input locations and sample from the prior.
rng, N, Nplot, S = MersenneTwister(123456), 100, 500, 100;
X, Xp = sort(rand(rng, N) * 10), range(-2.5, stop=12.5, length=Nplot);
f, ω, fpω, y = model();
fs, ωs, ys, ŷ = rand(rng, [f(Xp), ω(Xp), fpω(Xp), y(X)]);

# Compute posterior distribution over f′.
f′ = f | (y(X) ← ŷ);

# Sample from the posterior and write to file.
f′Xp = rand(rng, f′(Xp), S);

# Get posterior mean and marginals f′ and y′ and write them for plotting.
μf′, σf′ = marginals(f′(Xp));

###########################  Plot results - USE ONLY Julia-0.6!  ###########################

posterior_plot = plot();

# Prior over `f`.
plot!(posterior_plot, Xp, ys;
    linewidth=1.0,
    linecolor=:red,
    label="y");
plot!(posterior_plot, Xp, fs;
    linewidth=2.0,
    linecolor=:black,
    label="f");

# Posterior over `f`.
plot!(posterior_plot, Xp, μf′;
    linecolor=:blue,
    linewidth=2.0,
    label="f1");
plot!(posterior_plot, Xp, [μf′ μf′];
    linewidth=0.0,
    fillrange=[μf′ .- 3 .* σf′, μf′ .+ 3 * σf′],
    fillalpha=0.3,
    fillcolor=:blue,
    label="");
scatter!(posterior_plot, X, ŷ;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");
plot!(posterior_plot, Xp, f′Xp,
    linecolor=:blue,
    linealpha=0.2,
    label="");

display(posterior_plot);


text:
module ProxSDP

    using Arpack
    using KrylovKit
    using MathOptInterface
    using TimerOutputs

    using Printf
    using SparseArrays
    using LinearAlgebra
    import Random
    import Logging

    import LinearAlgebra: BlasInt

    include("structs.jl")
    include("util.jl")
    include("printing.jl")
    include("scaling.jl")
    include("equilibration.jl")
    include("pdhg.jl")
    include("residuals.jl")
    include("eigsolver.jl")
    include("prox_operators.jl")

    include("MOI_wrapper.jl")

    MOIU.@model _ProxSDPModelData () (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan) (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.PositiveSemidefiniteConeTriangle) () (MOI.SingleVariable,) (MOI.ScalarAffineFunction,) (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)

    Solver(;args...) = MOIU.CachingOptimizer(_ProxSDPModelData{Float64}(), ProxSDP.Optimizer(args))

end


text:
using Documenter
using MesoscaleML

theory = ["Introduction" => "equations.md",
    "Glossary" => "geoglossary.md",
    "2-Layer QG" => "2layerqg.md",
    "Mean Field 2-Layer QG" => "meanfield_2layerqg.md"
    ]

makedocs(
    sitename = "MesoscaleML",
    format = Documenter.HTML(),
    modules = [MesoscaleML],
    pages = [
        "Home" =>  "index.md",
        "Theory" => theory]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/upiterbarg/MesoscaleML.jl"
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RoutingPreference <: SwaggerModel
    routingChoice::Any # spec type: Union{ Nothing, String } # spec name: routingChoice
    publishMicrosoftEndpoints::Any # spec type: Union{ Nothing, Bool } # spec name: publishMicrosoftEndpoints
    publishInternetEndpoints::Any # spec type: Union{ Nothing, Bool } # spec name: publishInternetEndpoints

    function RoutingPreference(;routingChoice=nothing, publishMicrosoftEndpoints=nothing, publishInternetEndpoints=nothing)
        o = new()
        validate_property(RoutingPreference, Symbol("routingChoice"), routingChoice)
        setfield!(o, Symbol("routingChoice"), routingChoice)
        validate_property(RoutingPreference, Symbol("publishMicrosoftEndpoints"), publishMicrosoftEndpoints)
        setfield!(o, Symbol("publishMicrosoftEndpoints"), publishMicrosoftEndpoints)
        validate_property(RoutingPreference, Symbol("publishInternetEndpoints"), publishInternetEndpoints)
        setfield!(o, Symbol("publishInternetEndpoints"), publishInternetEndpoints)
        o
    end
end # type RoutingPreference

const _property_map_RoutingPreference = Dict{Symbol,Symbol}(Symbol("routingChoice")=>Symbol("routingChoice"), Symbol("publishMicrosoftEndpoints")=>Symbol("publishMicrosoftEndpoints"), Symbol("publishInternetEndpoints")=>Symbol("publishInternetEndpoints"))
const _property_types_RoutingPreference = Dict{Symbol,String}(Symbol("routingChoice")=>"String", Symbol("publishMicrosoftEndpoints")=>"Bool", Symbol("publishInternetEndpoints")=>"Bool")
Base.propertynames(::Type{ RoutingPreference }) = collect(keys(_property_map_RoutingPreference))
Swagger.property_type(::Type{ RoutingPreference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RoutingPreference[name]))}
Swagger.field_name(::Type{ RoutingPreference }, property_name::Symbol) =  _property_map_RoutingPreference[property_name]

const _allowed_RoutingPreference_routingChoice = ["MicrosoftRouting", "InternetRouting"]

function check_required(o::RoutingPreference)
    true
end

function validate_property(::Type{ RoutingPreference }, name::Symbol, val)
    if name === Symbol("routingChoice")
        Swagger.validate_param(name, "RoutingPreference", :enum, val, _allowed_RoutingPreference_routingChoice)
    end
end


text:
function TriMesh_PromoteL2Q(p,eConn)
#=
  TriMesh_PromoteL2Q
   Interpolate values at the midside nodes from the vertices in a
   triangular mesh.  Used to promote a linear to a quadratic representation
   of a finite element function.

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.
=#

  nElements = size(eConn,1)

  for nEl=1:nElements
    localNodes = eConn[nEl,:]
    p[ localNodes[4] ] = ( p[localNodes[1]] + p[localNodes[2]] )/2
    p[ localNodes[5] ] = ( p[localNodes[2]] + p[localNodes[3]] )/2
    p[ localNodes[6] ] = ( p[localNodes[3]] + p[localNodes[1]] )/2
  end

  return p
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export evaluate, selection_string

__ptr(selection::Selection) = __ptr(selection.__handle)
__const_ptr(selection::Selection) = __const_ptr(selection.__handle)

"""
Create a ``Selection`` from a selection string.
"""
function Selection(selection::AbstractString)
    ptr = @__check_ptr(lib.chfl_selection(pointer(selection)))
    return Selection(CxxPointer(ptr, is_const=false))
end

"""
Get the selection string used to create a given ``selection``.
"""
function selection_string(selection::Selection)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_selection_string(
            __const_ptr(selection), buffer, size
        ))
    )
end

"""
Get the size of a ``selection``, *i.e.* the number of atoms we are selecting
together.
"""
function Base.size(selection::Selection)
    result = Ref{UInt64}(0)
    __check(lib.chfl_selection_size(__const_ptr(selection), result))
    return result[]
end

"""
Evaluate a ``selection`` on a given ``frame``. This function return a list of
indexes or tuples of indexes of atoms in the frame matching the selection.
"""
function evaluate(selection::Selection, frame::Frame)
    # Get the number of matching atoms
    count = Ref{UInt64}(0)
    __check(lib.chfl_selection_evaluate(__ptr(selection), __const_ptr(frame), count))
    count = count[]
    # Allocate memory and get matching atoms
    matches = Array{lib.chfl_match}(undef, count)
    __check(lib.chfl_selection_matches(__ptr(selection), pointer(matches), count))
    # Return a clearer array of tuples
    result = []
    selection_size = size(selection)
    for match in matches
        @assert match.size == selection_size
        if selection_size == 1
            push!(result, match.atoms_1)
        elseif selection_size == 2
            push!(result, (match.atoms_1, match.atoms_2))
        elseif selection_size == 3
            push!(result, (match.atoms_1, match.atoms_2, match.atoms_3))
        elseif selection_size == 4
            push!(result, (match.atoms_1, match.atoms_2, match.atoms_3, match.atoms_4))
        end
    end
    return result
end

"""
Make a deep copy of a ``selection``.
"""
function Base.deepcopy(selection::Selection)
    ptr = lib.chfl_selection_copy(__const_ptr(selection))
    return Selection(CxxPointer(ptr, is_const=false))
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to cell alignment.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Apply the column alignment obtained from regex to the data after conversion to
# string.
function _apply_alignment_anchor_regex!(
    data_str::Matrix{Vector{String}},
    alignment::Vector{Symbol},
    id_cols::Vector{Int},
    id_rows::Vector{Int},
    Δc::Int,
    # Configurations.
    alignment_anchor_fallback::Symbol,
    alignment_anchor_fallback_override::Dict{Int, Symbol},
    alignment_anchor_regex::Dict{Int, T} where T<:AbstractVector{Regex},
    cell_alignment_override::Dict{Tuple{Int, Int}, Symbol}
)

    num_printed_rows, num_printed_cols = size(data_str)

    # If we have a key `0`, then it will be used to align all the columns.
    if haskey(alignment_anchor_regex, 0)
        alignment_keys = collect(1:maximum(id_cols))
        global_regex   = true
    else
        alignment_keys = sort(collect(keys(alignment_anchor_regex)))
        global_regex   = false
    end

    @inbounds for jc in alignment_keys
        j = findfirst(x -> x == jc, id_cols)
        j === nothing && continue
        j += Δc

        # If `j` is larger than `num_printed_cols`, then we can stop the
        # processing, because the keys are ordered and anything larger than
        # `num_printed_cols` will not be printed.
        j > num_printed_cols && break

        regex = global_regex ? alignment_anchor_regex[0] : alignment_anchor_regex[jc]

        # Store in which column we must align the match.
        alignment_column = 0

        # We need to pass through the entire row searching for matches to
        # compute in which column we need to align the matches.
        for i = 1:num_printed_rows
            !isassigned(data_str, i, j) && continue
            haskey(cell_alignment_override, (id_rows[i], jc)) && continue

            for l = 1:length(data_str[i, j])
                line = data_str[i, j][l]

                m = nothing

                for r in regex
                    m_r = findfirst(r, line)
                    if m_r !== nothing
                        m = m_r
                        break
                    end
                end

                if m !== nothing
                    alignment_column_i = first(m)
                else
                    # If a match is not found, then the alignment column
                    # depends on the user selection.

                    fallback = haskey(alignment_anchor_fallback_override, jc) ?
                        alignment_anchor_fallback_override[jc] :
                        alignment_anchor_fallback

                    if fallback == :c
                        data_ijl_len = textwidth(line)
                        alignment_column_i = cld(data_ijl_len, 2)
                    elseif fallback == :r
                        data_ijl_len = textwidth(line)
                        alignment_column_i = data_ijl_len + 1
                    else
                        alignment_column_i = 0
                    end
                end

                alignment_column_i > alignment_column &&
                    (alignment_column = alignment_column_i)
            end
        end

        # Variable to store the largest width of a cell.
        largest_cell_width = 0

        # Now, we need to pass again applying the alignments.
        for i = 1:num_printed_rows
            !isassigned(data_str, i, j) && continue
            haskey(cell_alignment_override, (id_rows[i], jc)) && continue

            for l = 1:length(data_str[i, j])
                line = data_str[i, j][l]

                m = nothing

                for r in regex
                    m_r = findfirst(r, line)
                    if m_r !== nothing
                        m = m_r
                        break
                    end
                end

                if m !== nothing
                    match_column_k = first(m)
                    pad = alignment_column - match_column_k
                else
                    # If a match is not found, then the alignment column
                    # depends on the user selection.

                    fallback = haskey(alignment_anchor_fallback_override, jc) ?
                        alignment_anchor_fallback_override[jc] :
                        alignment_anchor_fallback

                    if fallback == :c
                        data_ijl_len = textwidth(line)
                        pad = alignment_column - cld(data_ijl_len, 2)
                    elseif fallback == :r
                        data_ijl_len = textwidth(line)
                        pad = alignment_column - data_ijl_len - 1
                    else
                        pad = alignment_column
                    end
                end

                # Make sure `pad` is positive.
                pad < 0 && (pad = 0)

                data_str[i, j][l]  = " "^pad * line
                data_ijl_len = textwidth(data_str[i, j][l])

                if data_ijl_len > largest_cell_width
                    largest_cell_width = data_ijl_len
                end
            end
        end

        # The third pass aligns the elements correctly. This is performed by
        # adding spaces to the right so that all the cells have the same width.
        for i = 1:num_printed_rows
            !isassigned(data_str, i, j) && continue
            haskey(cell_alignment_override, (id_rows[i], jc)) && continue

            for l = 1:length(data_str[i, j])
                pad = largest_cell_width - textwidth(data_str[i, j][l])
                pad < 0 && (pad = 0)
                data_str[i, j][l] = data_str[i, j][l] * " "^pad
            end
        end
    end

    return nothing
end

# Compute a list of cells in which the alignment is overridden by the user.
function _compute_cell_alignment_override(
    data::Any,
    id_cols::Vector{Int},
    id_rows::Vector{Int},
    Δc::Int,
    num_printed_cols::Int,
    num_printed_rows::Int,
    # Configurations.
    cell_alignment::Ref{Any}
)
    # Dictionary with the cells in which the alignment is overridden.
    cell_alignment_override = Dict{Tuple{Int, Int}, Symbol}()

    for j = (1 + Δc):num_printed_cols
        jc = id_cols[j - Δc]

        for i = 1:num_printed_rows
            ir = id_rows[i]

            for f in cell_alignment.x
                aux = f(_getdata(data), ir, jc)

                if (aux == :l) || (aux == :c) || (aux == :r) ||
                   (aux == :L) || (aux == :C) || (aux == :R)
                    cell_alignment_override[(ir, jc)] = aux
                end
            end
        end
    end

    return cell_alignment_override
end


text:
# Test plotting
module TestPlots

using OpenStreetMap
using Base.Test
import Winston

MAP_FILENAME = "tech_square.osm"

# Load and crop map to file bounds
nodesLLA, hwys, builds, feats = getOSMData(MAP_FILENAME)
boundsLLA = getBounds(parseMapXML(MAP_FILENAME))
cropMap!(nodesLLA, boundsLLA, highways=hwys, buildings=builds, features=feats, delete_nodes=true)

roads = roadways(hwys)
peds = walkways(hwys)
cycles = cycleways(hwys)
bldg_classes = classify(builds)
feat_classes = classify(feats)

p = plotMap(nodesLLA, highways=hwys, buildings=builds, features=feats, bounds=boundsLLA, width=500, feature_classes=feat_classes, building_classes=bldg_classes, cycleways=cycles, walkways=peds, roadways=roads)

@test typeof(p) == Winston.FramedPlot
@test Winston.getattr(p, "xlabel") == "Longitude (deg)"
@test Winston.getattr(p, "ylabel") == "Latitude (deg)"
@test Winston.getattr(p.x1, "draw_axis") == true
@test Winston.getattr(p.x1, "draw_grid") == false
@test Winston.getattr(p, "xrange") == (-71.0939, -71.0891)
@test Winston.getattr(p, "yrange") == (42.3626, 42.3659)

lla_ref = center(boundsLLA)
nodesENU = ENU(nodesLLA, lla_ref)
boundsENU = ENU(boundsLLA, lla_ref)

p2 = plotMap(nodesENU, highways=hwys, buildings=builds, features=feats, bounds=boundsENU, width=500, feature_classes=feat_classes, building_classes=bldg_classes, cycleways=cycles, walkways=peds, roadways=roads, km=true, fontsize=4)

@test typeof(p2) == Winston.FramedPlot
@test Winston.getattr(p2, "xlabel") == "East (km)"
@test Winston.getattr(p2, "ylabel") == "North (km)"
@test Winston.getattr(p2.x1, "draw_axis") == true
@test Winston.getattr(p2.x1, "draw_grid") == false
@test_approx_eq_eps Winston.getattr(p2, "xrange")[1] -0.19770898045741428 1e-6
@test_approx_eq_eps Winston.getattr(p2, "xrange")[2] 0.19770898045628862 1e-6
@test_approx_eq_eps Winston.getattr(p2, "yrange")[1] -0.18328257005093138 1e-6
@test_approx_eq_eps Winston.getattr(p2, "yrange")[2] 0.18328541309049148 1e-6

end # module TestPlots


text:
# This file holds the definition of the functions pertaining to the
#   Subpopulation type.

export  setSubpopulationName!,
        setSubpopulationSourceNode!,
        addSubpopulationCondition!,
        clearSubpopulationConditions!,
        setSubpopulationConditions!

"""
```
setSubpopulationName!(
    subpopulation::Subpopulation,
    name::String )
```
This function sets the name of the subpopulation `subpopulation` to `name`.

The function returns `true`, indicating the name is successfully set.
"""
function setSubpopulationName!( subpopulation::Subpopulation,
    name::String )::Bool

    subpopulation.name = name
    return true

end  # setSubpopulationName!( subpopulation, name )

"""
```
setSubpopulationSourceNode!(
    subpopulation::Subpopulation,
    sourceNode::String )
```
This function sets the source node of the subpopulation `subpopulation` to `node`. The source node can be any base node, compound node, or the entire population. The latter is entered by `""` or `"active"`.

The function returns `true`, indicating the source node is successfully set.
"""
function setSubpopulationSourceNode!( subpopulation::Subpopulation,
    sourceNode::String )::Bool

    subpopulation.sourceNode = lowercase( sourceNode ) ∈ ["active", ""] ?
        "active" : sourceNode
    return true

end  # setSubpopulationSourceNode!( subpopulation, sourceNode )

"""
```
addSubpopulationCondition!(
    subpopulation::Subpopulation,
    conditions::MPcondition... )
```
This function adds the conditions in `conditions` to the subpopulation `subpopulation`.

Several other conditions are also permitted, relating to time and/or the history of the personnel members:
* `"age"`: the age of the personnel members,
* `"tenure"`: the tenure of the personnel members,
* `"time in node"`: the time in the personnel members' current base node,
* `"had transition"`: the personnel member did a transition with the given name,
* `"started as"`: the personnel member entered the system in the given base node,
* `"was"`: the personnel member was in the given base node at any point.

Note that this function performs no checks on whether the conditions are feasible, contradict each other or contradict the source node of the subpopulation.

This function returns `true`, indicating the conditions are successfully added.
"""
function addSubpopulationCondition!( subpopulation::Subpopulation,
    conditions::MPcondition... )::Bool

    for cond in conditions
        tmpCond = MPcondition( lowercase( cond.attribute ), cond.operator,
            cond.value )  # MPcondition is immutable!

        if tmpCond.attribute ∈ timeAttributes
            push!( subpopulation.timeConds, tmpCond )
        elseif tmpCond.attribute ∈ histAttributes
            push!( subpopulation.historyConds, tmpCond )
        else
            push!( subpopulation.attributeConds, deepcopy( cond ) )
        end  # if tmpCond.attribute ∈ timeAttributes
    end  # for cond in conditions

    return true

end  # addSubpopulationCondition!( subpopulation, conditions )

"""
```
clearSubpopulationConditions!( subpopulation::Subpopulation )
```
This function clears all conditions from the subpopulation `subpopulation`.

The function returns `true`, indicating the conditions are successfully cleared.
"""
function clearSubpopulationConditions!( subpopulation::Subpopulation )::Bool

    empty!( subpopulation.timeConds )
    empty!( subpopulation.historyConds )
    empty!( subpopulation.attributeConds )
    return true

end  # clearSubpopulationConditions!( subpopulation::Subpopulation )

"""
```
setSubpopulationConditions!(
    subpopulation::Subpopulation,
    conditions::Vector{MPcondition} )
```
This function sets the conditions in `conditions` to the subpopulation `subpopulation`.

Several other conditions are also permitted, relating to time and/or the history of the personnel members:
* `"age"`: the age of the personnel members,
* `"tenure"`: the tenure of the personnel members,
* `"time in node"`: the time in the personnel members' current base node,
* `"had transition"`: the personnel member did a transition with the given name,
* `"started as"`: the personnel member entered the system in the given base node,
* `"was"`: the personnel member was in the given base node at any point.

Note that this function performs no checks on whether the conditions are feasible, contradict each other or contradict the source node of the subpopulation.

This function returns `true`, indicating the conditions are successfully set.
"""
function setSubpopulationConditions!( subpopulation::Subpopulation,
    conditions::Vector{MPcondition} )::Bool

    clearSubpopulationConditions!( subpopulation )
    return addSubpopulationCondition!( subpopulation, conditions... )

end  # setSubpopulationConditions!( subpopulation, conditions )

function Base.show( io::IO, subpopulation::Subpopulation )

    print( io, "Subpopulation: ", subpopulation.name )
    print( io, "\n  Root node: ", subpopulation.sourceNode )

    if !all( [isempty( subpopulation.timeConds ),
        isempty( subpopulation.historyConds ),
        isempty( subpopulation.attributeConds )] )
        print( io, "\n  Conditions\n    ", join( vcat( subpopulation.timeConds,
            subpopulation.historyConds, subpopulation.attributeConds ),
            "\n    " ) )
    end  # if !all( [...] )

    return

end  # show( io, state )


text:
@testset "orient" begin

    @testset "orient_strands" begin
        seqs = ["GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                        "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATC",
                        "GTCGATCGACTAGCTGCATGACTGACATCGACATCGACGGAGCATGACTAGGACGACGAGCATCAAAAAAAAAAAAAAAAAAAAA"]
        refs = "GTCGAT"
        names = ["s1", "s2", "s3"]
        phreds = [[Int8(30) for i in 1:length(seq)] for seq in seqs]
        oriented, newphreds, newnames = orient_strands(seqs, phreds, names, refs)

        @test oriented == ["GATGCTCGTCGTCCTAGTCATGCTCCGTCGATGTCGATGTCAGTCATGCAGCTAGTCGATCGAC",
                           "GATGCTCGTCGTCCTAGTCATGCTCCGTCGATGTCGATGTCAGTCATGCAGCTAGTCGATCGAC",
                           "TTTTTTTTTTTTTTTTTTTTTGATGCTCGTCGTCCTAGTCATGCTCCGTCGATGTCGATGTCAGTCATGCAGCTAGTCGATCGAC"]
    end

end


text:
@testset "Topology" begin
  @testset "FullTopology" begin
    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    struc = FullTopology(elems)
    @test nvertices(struc) == 4
    @test nelements(struc) == 2

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    struc = FullTopology(elems)
    @test nvertices(struc) == 6
    @test nelements(struc) == 4

    # 1 triangle + 3 quadrangles + 1 triangle hole
    elems = connect.([(1,2,6,5),(2,4,7,6),(4,3,7),(3,1,5,7)])
    struc = FullTopology(elems)
    @test nvertices(struc) == 7
    @test nelements(struc) == 4
  end

  @testset "GridTopology" begin
    t = GridTopology(3)
    @test nvertices(t) == 4
    @test nelements(t) == 3
    @test nfacets(t) == 4
    @test size(t) == (3,)
    @test element(t, 1) == connect((1,2))
    @test element(t, 2) == connect((2,3))
    @test element(t, 3) == connect((3,4))
    @test faces(t, 1) == elements(t)

    t = GridTopology(3, 4)
    @test nvertices(t) == 20
    @test nelements(t) == 12
    @test nfacets(t) == 31
    @test size(t) == (3, 4)
    @test element(t, 1) == connect((1,2,6,5))
    @test element(t, 5) == connect((6,7,11,10))
    @test faces(t, 2) == elements(t)
    @test facet.(Ref(t), 1:24) == connect.([(1,2),(1,5),(2,3),(2,6),(3,4),(3,7),
                                            (5,6),(5,9),(6,7),(6,10),(7,8),(7,11),
                                            (9,10),(9,13),(10,11),(10,14),(11,12),(11,15),
                                            (13,14),(13,17),(14,15),(14,18),(15,16),(15,19)])
    @test facet.(Ref(t), 25:31) == connect.([(4,8),(17,18),(8,12),(18,19),(12,16),(19,20),(16,20)])

    t = GridTopology(3, 4, 2)
    @test nvertices(t) == 60
    @test nelements(t) == 24
    @test nfacets(t) == 3*24 + 3*4 + 4*2 + 3*2
    @test size(t) == (3, 4, 2)
    @test element(t, 1) == connect((1,2,6,5,21,22,26,25), Hexahedron)
    @test element(t, 5) == connect((6,7,11,10,26,27,31,30), Hexahedron)
    @test faces(t, 3) == elements(t)
  end

  @testset "HalfEdgeTopology" begin
    function test_halfedge(elems, structure)
      @test nelements(structure) == length(elems)
      for e in 1:nelements(structure)
        he = half4elem(e, structure)
        inds = indices(elems[e])
        @test he.elem == e
        @test he.head ∈ inds
      end
    end

    # 2 triangles as a list of half-edges
    h1  = HalfEdge(1, 1)
    h2  = HalfEdge(2, nothing)
    h3  = HalfEdge(2, 1)
    h4  = HalfEdge(3, 2)
    h5  = HalfEdge(3, 1)
    h6  = HalfEdge(1, nothing)
    h7  = HalfEdge(2, 2)
    h8  = HalfEdge(4, nothing)
    h9  = HalfEdge(4, 2)
    h10 = HalfEdge(3, nothing)
    h1.half =  h2; h2.half  = h1
    h3.half =  h4; h4.half  = h3
    h5.half =  h6; h6.half  = h5
    h7.half =  h8; h8.half  = h7
    h9.half = h10; h10.half = h9
    h1.prev = h5;  h1.next = h3
    h3.prev = h1;  h3.next = h5
    h4.prev = h9;  h4.next = h7
    h5.prev = h3;  h5.next = h1
    h7.prev = h4;  h7.next = h9
    h9.prev = h7;  h9.next = h4
    halves  = [(h1,h2),(h3,h4),(h5,h6),(h7,h8),(h9,h10)]
    struc = HalfEdgeTopology(halves)
    @test half4elem(1, struc) == h1
    @test half4elem(2, struc) == h4
    @test half4vert(1, struc) == h1
    @test half4vert(2, struc) == h3
    @test half4vert(3, struc) == h4
    @test half4vert(4, struc) == h9
    @test edge4pair((1,2), struc) == 1
    @test edge4pair((2,1), struc) == 1
    @test edge4pair((2,3), struc) == 2
    @test edge4pair((3,2), struc) == 2
    @test edge4pair((3,1), struc) == 3
    @test edge4pair((1,3), struc) == 3
    @test edge4pair((2,4), struc) == 4
    @test edge4pair((4,2), struc) == 4
    @test edge4pair((4,3), struc) == 5
    @test edge4pair((3,4), struc) == 5

    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    struc = HalfEdgeTopology(elems)
    @test nvertices(struc) == 4
    @test nelements(struc) == 2
    @test nfacets(struc) == 5
    test_halfedge(elems, struc)

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    struc = HalfEdgeTopology(elems)
    @test nvertices(struc) == 6
    @test nelements(struc) == 4
    @test nfacets(struc) == 9
    test_halfedge(elems, struc)

    # 1 triangle + 3 quadrangles + 1 triangle hole
    elems = connect.([(1,2,6,5),(2,4,7,6),(4,3,7),(3,1,5,7)])
    struc = HalfEdgeTopology(elems)
    @test nvertices(struc) == 7
    @test nelements(struc) == 4
    @test nfacets(struc) == 11
    test_halfedge(elems, struc)
  end
end


text:
## fem.jl ##
abstract type GaussianRandomFieldGenerator end

function GaussianRandomField(mean::Vector{T},cov::CovarianceFunction{d},method::GaussianRandomFieldGenerator,p::Matrix{T},t::Matrix{N};mode="nodes",kwargs...) where {d,T<:Real,N<:Int}
    size(p,2) == d || throw(DimensionMismatch("second dimension of points must be equal to $(d)"))
    size(t,2) == d+1 || throw(DimensionMismatch("second dimension of nodes must be equal to $(d+1)"))
	typeof(method) <: CirculantEmbedding && throw(ArgumentError("cannot use circulant embedding with a finite element mesh"))
    if mode == "center"
        length(mean) == size(t,1) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))
        pts = compute_centers(p,t)'
        tri = Matrix{N}(undef,0,0)
    elseif mode == "nodes"
        length(mean) == size(p,1) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))
        pts = p'
        tri = t'
    else
        throw(ArgumentError("unknown mode $(mode)"))
    end
    _GaussianRandomField(mean,cov,method,pts,tri;kwargs...)
end

GaussianRandomField(cov::CovarianceFunction,method::GaussianRandomFieldGenerator,p::Matrix{T},t::Matrix{N} where {N<:Int};kwargs...) where {T<:Real} = GaussianRandomField(0,cov,method,p,t;kwargs...)

function GaussianRandomField(mean::N where {N<:Real},cov::CovarianceFunction,method::GaussianRandomFieldGenerator,p::Matrix{T},t::Matrix{N} where {N<:Int};mode="nodes",kwargs...) where {T<:Real}
    if mode == "center"
        M = mean*ones(T,size(t,1))
    elseif mode == "nodes"
        M = mean*ones(T,size(p,1))
    else
        throw(ArgumentError("unknown mode $(mode)"))
    end
    GaussianRandomField(M,cov,method,p,t;mode=mode,kwargs...)
end

function compute_centers(p,t)
    n = size(t,1)
    d = size(p,2)
    pts = zeros(n,d)
    for i in 1:d
        x = p[t[:],i]
        x = reshape(x,size(t))
        pts[:,i] = mean(x,dims=2)
    end
    pts
end

shape(grf::GaussianRandomField{C,M,Tuple{T1,T2}}) where {C,M,T1<:AbstractMatrix,T2<:AbstractMatrix} = size(grf.pts[1],2)

show(io::IO,grf::GaussianRandomField{C,M,Tuple{T1,T2}}) where {C,M,T1<:AbstractMatrix,T2<:AbstractMatrix} = print(io, "Gaussian random field with $(grf.cov) on a mesh with $(size(grf.pts[1],2)) points and $(size(grf.pts[2],2)) elements, using a $(M())")