text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
  	drillhole(collar::Collar, survey::Survey, intervals::Intervals)

Desurvey drill hole based on collar, survey and interval table(s) information.
The intervals can be passed as a single `Interval` or as an array of
`Interval`. Outputs a `DrillHole` object.
"""
function drillhole(collar::Collar,survey::Survey,intervals::Intervals)
	# pre process information
	pars  = getcolnames(survey,intervals)
	warns = validations(collar, survey, intervals)
	in("Error",warns[:,:TYPE]) && (return DrillHole(nothing,nothing,pars,warns))

	# create trace information
	trace = gettrace(collar, survey)
	fillxyz!(trace, pars)

	# merge interval tables
	table = mergetables(intervals, pars)
	fillxyz!(table, trace, pars)

	DrillHole(table,trace,pars,warns)
end

function getcolnames(s,i)
	f = i isa Interval ? i.from : i[1].from
	t = i isa Interval ? i.to   : i[1].to
	m = s.method == :tangential
	c = s.convention

	# get most common dip sign and assume it is downwards
	if c == :auto
		df = s.file isa String ? CSV.read(s.file, DataFrame) : s.file
		c  = sum(sign.(df[!,s.dip])) > 0 ? :positivedownwards : :negativedownwards
	end

	inv  = (c == :positivedownwards)
	pars = (holeid=s.holeid, at=s.at, azm=s.azm, dip=s.dip, from=f,
	       to=t, invdip=inv, tang=m)
end

function gettrace(c, s)
	collar = c.file isa String ? CSV.read(c.file, DataFrame) : c.file
	survey = s.file isa String ? CSV.read(s.file, DataFrame) : s.file

	# rename collar columns to match survey columns if necessary
	n1 = (c.x,c.y,c.z,c.holeid)
	n2 = ( :X, :Y, :Z,s.holeid)
	namepairs = [a=>b for (a,b) in zip(Symbol.(n1),Symbol.(n2)) if a!=b]
	length(namepairs) > 0 && rename!(collar, namepairs...)

	# force coordinates to floats if necessary
	for k in [:X,:Y,:Z]
		force_float = !(eltype(collar[!,k]) <: Float64)
		force_float && (collar[!, k] = convert.(Float64, collar[:, k]))
	end

	# merge collar coordinates to initial survey point
	collar[!,s.at] .= 0.0
	trace = leftjoin(survey,collar,on=[s.holeid,s.at])
	sort!(trace, [s.holeid,s.at])
end

# minimum curvature desurvey method
function mincurv(az1, dp1, az2, dp2, d12)
	dp1, dp2 = (90-dp1), (90-dp2)

    DL = acos(cosd(dp2-dp1)-sind(dp1)*sind(dp2)*(1-cosd(az2-az1)))
    RF = DL!=0.0 ? 2*tan(DL/2)/DL : 1

    dx = 0.5*d12*(sind(dp1)*sind(az1)+sind(dp2)*sind(az2))*RF
	dy = 0.5*d12*(sind(dp1)*cosd(az1)+sind(dp2)*cosd(az2))*RF
	dz = 0.5*d12*(cosd(dp1)+cosd(dp2))*RF
    dx, dy, dz
end

# tangential desurvey method
function tangential(az1, dp1, d12)
	dp1 = (90-dp1)
    dx  = d12*sind(dp1)*sind(az1)
	dy  = d12*sind(dp1)*cosd(az1)
	dz  = d12*cosd(dp1)
    dx, dy, dz
end

# find survey depths bounding given depth
function findbounds(depths::AbstractArray, at)
	# get closest survey
    nearid = findmin(abs.(depths.-at))[2]
    nearest = depths[nearid]

	# check if depth is after last interval
	nearid == length(depths) && nearest < at && return (nearid,nearid)

	# return (previous, next) survey ids for given depth
    nearest == at && return (nearid, nearid)
    nearest >  at && return (nearid-1, nearid)
    nearest <  at && return (nearid, nearid+1)
end

# convert survey angles to 3-D vector and vice versa
angs2vec(az,dp) = [sind(az)*cosd(-dp), cosd(az)*cosd(-dp), sind(-dp)]
vec2angs(i,j,k) = [atand(i,j), -asind(k)]

# average angle between two surveyed intervals
function weightedangs(angs1,angs2,d12,d1x)
	# angle to vectors
    v1 = angs2vec(angs1...)
    v2 = angs2vec(angs2...)

	# weight average vector according to distance to surveys
    p2   = d1x/d12
    p1   = 1-p2
    v12  = v1*p1 + v2*p2
    v12 /= sqrt(sum(abs2,v12))

	# convert average vector to survey angles and return it
    azm, dip = vec2angs(v12...)
    azm, dip
end

# fill xyz for dh trace files
function fillxyz!(trace, pars)
	# get column names
	at, az, dp, tang = pars.at, pars.azm, pars.dip, pars.tang
	f = pars.invdip ? -1 : 1

	# loop trace file
	for i in 1:size(trace,1)
		# pass depth 0 where collar coordinates are already available
        trace[i,at] == 0 && continue

		# get distances and angles; return increments dx,dy,dz
		d12      = trace[i,at] - trace[i-1,at]
        az1, dp1 = trace[i-1,az], f*trace[i-1,dp]
        az2, dp2 = trace[i,az], f*trace[i,dp]
        dx,dy,dz = tang ? tangential(az1,dp1,d12) : mincurv(az1,dp1,az2,dp2,d12)

		# add increments dx,dy,dz to previous coordinates
        trace[i,:X] = dx + trace[i-1,:X]
        trace[i,:Y] = dy + trace[i-1,:Y]
        trace[i,:Z] = dz + trace[i-1,:Z]
    end
end

# fill xyz for interval tables with from-to information
function fillxyz!(tab, trace, pars)
	# get column names
	bh, at, az, dp, tang = pars.holeid, pars.at, pars.azm, pars.dip, pars.tang
	from, to = pars.from, pars.to
	f = pars.invdip ? -1 : 1

	# initialize coordinate columns with float values
	tab[!,:X] .= -9999.9999
	tab[!,:Y] .= -9999.9999
	tab[!,:Z] .= -9999.9999

	# get first hole name and get trace of that hole
	lastbhid = tab[1,bh]
	dht = trace[(trace[!,bh] .== lastbhid),:]

	# loop all intervals
	for i in 1:size(tab,1)
		# get hole name and mid point depth
		bhid, atx = tab[i,bh], tab[i,from]+tab[i,:LENGTH]/2
		# update trace if hole name is different than previous one
		bhid != lastbhid && (dht = trace[(trace[!,bh] .== bhid),:])
		lastbhid = bhid
		# pass if no survey is available (WARN)
		size(dht, 1) == 0 && continue

		# get surveys bounding given depth
		b   = findbounds(dht[:,at],atx)
		d1x = atx-dht[b[1],at]

		if d1x == 0
			# if interval depth matches trace depth, get trace coordinates
			tab[i,:X] = dht[b[1],:X]
			tab[i,:Y] = dht[b[1],:Y]
			tab[i,:Z] = dht[b[1],:Z]
		else
			# if not, calculate coordinates increments dx,dy,dz
			d12 = dht[b[2],at]-dht[b[1],at]
			az1, dp1 = dht[b[1],az], f*dht[b[1],dp]
			az2, dp2 = dht[b[2],az], f*dht[b[2],dp]
			azx, dpx = b[1]==b[2] ? (az2, dp2) : weightedangs([az1,dp1],[az2,dp2],d12,d1x)
			dx,dy,dz = tang ? tangential(az1,dp1,d1x) : mincurv(az1,dp1,azx,dpx,d1x)

			# add increments dx,dy,dz to trace coordinates
			tab[i,:X] = dx + dht[b[1],:X]
			tab[i,:Y] = dy + dht[b[1],:Y]
			tab[i,:Z] = dz + dht[b[1],:Z]
		end
	end
	# check if some coordinate was not filled and return a warning if necessary
	filter!(row -> row.X != -9999.9999, tab)
end


text:
Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y

instantiate(x::ModelingToolkit.Variable{Real}) = 1.0
instantiate(x::ModelingToolkit.Variable) = oneunit(1*ModelingToolkit.vartype(x))
function instantiate(x::Num)
    x = value(x)
    if x.op isa Sym
        return instantiate(x.op)
    elseif x.op isa Differential
        instantiate(x.args[1])/instantiate(x.args[1].args[1])
    else
        x.op(instantiate.(x.args)...)
    end
end

function validate(eq::ModelingToolkit.Equation)
    try
        return typeof(instantiate(eq.lhs)) == typeof(instantiate(eq.rhs))
    catch
        return false
    end
end

function validate(sys::AbstractODESystem)
    all(validate.(equations(sys)))
end

ext: jl


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

@doc raw"""CustomResourceDefinitionSpec describes how a user wants their resource to appear

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec(;
        conversion=nothing,
        group=nothing,
        names=nothing,
        preserveUnknownFields=nothing,
        scope=nothing,
        versions=nothing,
    )

    - conversion::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceConversion : conversion defines conversion settings for the CRD.
    - group::String : group is the API group of the defined custom resource. The custom resources are served under `/apis/<group>/...`. Must match the name of the CustomResourceDefinition (in the form `<names.plural>.<group>`).
    - names::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames : names specify the resource and kind names for the custom resource.
    - preserveUnknownFields::Bool : preserveUnknownFields indicates that object fields which are not specified in the OpenAPI schema should be preserved when persisting to storage. apiVersion, kind, metadata and known fields inside metadata are always preserved. This field is deprecated in favor of setting `x-preserve-unknown-fields` to true in `spec.versions[*].schema.openAPIV3Schema`. See https://kubernetes.io/docs/tasks/access-kubernetes-api/custom-resources/custom-resource-definitions/#pruning-versus-preserving-unknown-fields for details.
    - scope::String : scope indicates whether the defined custom resource is cluster- or namespace-scoped. Allowed values are `Cluster` and `Namespaced`.
    - versions::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion} : versions is the list of all API versions of the defined custom resource. Version names are used to compute the order in which served versions are listed in API discovery. If the version string is \"kube-like\", it will sort above non \"kube-like\" version strings, which are ordered lexicographically. \"Kube-like\" versions start with a \"v\", then are followed by a number (the major version), then optionally the string \"alpha\" or \"beta\" and another number (the minor version). These are sorted first by GA > beta > alpha (where GA is a version with no suffix such as beta or alpha), and then by comparing major version, then minor version. An example sorted list of versions: v10, v2, v1, v11beta2, v10beta3, v3beta1, v12alpha1, v11alpha2, foo1, foo10.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec <: SwaggerModel
    conversion::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceConversion } # spec name: conversion
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    names::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames } # spec name: names
    preserveUnknownFields::Any # spec type: Union{ Nothing, Bool } # spec name: preserveUnknownFields
    scope::Any # spec type: Union{ Nothing, String } # spec name: scope
    versions::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion} } # spec name: versions

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec(;conversion=nothing, group=nothing, names=nothing, preserveUnknownFields=nothing, scope=nothing, versions=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("conversion"), conversion)
        setfield!(o, Symbol("conversion"), conversion)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("names"), names)
        setfield!(o, Symbol("names"), names)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("preserveUnknownFields"), preserveUnknownFields)
        setfield!(o, Symbol("preserveUnknownFields"), preserveUnknownFields)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("scope"), scope)
        setfield!(o, Symbol("scope"), scope)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("versions"), versions)
        setfield!(o, Symbol("versions"), versions)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec = Dict{Symbol,Symbol}(Symbol("conversion")=>Symbol("conversion"), Symbol("group")=>Symbol("group"), Symbol("names")=>Symbol("names"), Symbol("preserveUnknownFields")=>Symbol("preserveUnknownFields"), Symbol("scope")=>Symbol("scope"), Symbol("versions")=>Symbol("versions"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec = Dict{Symbol,String}(Symbol("conversion")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceConversion", Symbol("group")=>"String", Symbol("names")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames", Symbol("preserveUnknownFields")=>"Bool", Symbol("scope")=>"String", Symbol("versions")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion}")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec)
    (getproperty(o, Symbol("group")) === nothing) && (return false)
    (getproperty(o, Symbol("names")) === nothing) && (return false)
    (getproperty(o, Symbol("scope")) === nothing) && (return false)
    (getproperty(o, Symbol("versions")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }, name::Symbol, val)
end


text:
"""
This example shows how to create Tree visualizations.

The easiest starting point is a `Dict` object with the info we want
to show in the tree
"""

import Term.tree: Tree

data = Dict(
    "a" => 1,
    "b" => Int64,
    "c" => [1, 2, 3],
)

print(Tree(data))

"""
`Tree` has several parameters that can be used to change the appearance
"""

print("\n\n")
for guides_type in (:standardtree, :boldtree, :asciitree)
    print(Tree(data; title=string(guides_type), guides_type=guides_type))
end

# set the colors!
print("\n\n")
print(
    Tree(
        data,
        title="my colors",
        title_style="bold red",
        node_style="blue underline",
        leaf_style="green",
        guides_style="red dim"
    )
)

"""
Tree can handle nested data too!
"""
data = Dict(
    "lvl1" => Dict(
        "lvl2" => Dict(
            "a" => 1,
            "b" => 3,
        ),
        "l2a" => "nested dicts rule!",
        "l2b" => "wohoo"
    )
)

print(Tree(data))


text:
function right_angle_check(x1::Int, y1::Int, x2::Int, y2::Int)
    if x1 == 0 && x2 == 0
        return 0
    end

    if x1 == x2
        if y1 == 0
            return 1
        end
        return 0
    end

    if x1 == 0
        if y2 == y1 || y2 == 0 || y2//x2 * (y2 - y1)//(x2 - x1) == -1
            return 1
        end
        return 0
    end

    if x2 == 0
        if y1 == y2 || y1 == 0 || y1//x1 * (y2 - y1)//(x2 - x1) == -1
            return 1
        end
        return 0
    end

    g_1 = y1//x1
    g_2 = y2//x2
    g_3 = (y2 - y1)//(x2 - x1)

    if g_1 * g_2 == -1 || g_1 * g_3 == -1 || g_2 * g_3 == -1
        return 1
    end

    return 0
end

function count_triangles(dim::Int)
    count = 0
    for p in 1:(dim ^ 2 - 2)
        for q in (p + 1):(dim ^ 2 - 1)
            count += right_angle_check(p % dim, Int(floor(p/dim)), q % dim, Int(floor(q/dim)))
        end
    end
    println(count)
end

count_triangles(51)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
#=
---------------------------------------------------------------------
This function computes the Volume integral contribution to the RHS vector
for a unified CGDG method for the 2D Wave Equation.

Written by F.X. Giraldo on July 12, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function create_rhs_volume(q,u,De_x,De_y,intma,periodicity,DFloat)

    #Get lengths of arrays
    (Np,Np,Ne)=size(intma)
    Npoin=length(q)

    #Store local arrays
    rhs=zeros(DFloat,Npoin)

    #build differentiation matrix contribution
    for e=1:Ne
        for i=1:Np, j=1:Np
            I=intma[i,j,e]
            for k=1:Np, l=1:Np
                K=intma[k,l,e]
                rhs[periodicity[I]]+=De_x[i,j,k,l,e]*q[K]*u[1,K] + De_y[i,j,k,l,e]*q[K]*u[2,K]
            end #k,l
        end #i,j
    end #e

    #Return arrays
    return(rhs)
end


text:
using ClusterTrees, BEAST, CompScienceMeshes

a, h = 1.0, 0.25
Γ1 = meshsphere(a, h)
Γ2 = CompScienceMeshes.translate(Γ1, point(2.5,0,0))

tfs = raviartthomas(Γ1)
bfs = raviartthomas(Γ2)

p = positions(tfs)
q = positions(bfs)

p, tp, permp = clustertree(p)
q, tq, permq = clustertree(q)

# Collect all boxes at depth 5
B = Array{typeof(tp)}(0)
depthfirst(tp) do t, level
  level == 5 && push!(B,t)
end

using MATLAB

# Plot points in the boxes at depth 5 in different colours
@matlab figure()
@matlab hold("on")
for b in B
  ps = p[b[1].begin_idx : b[1].end_idx-1]
  ps = [p[i] for p in ps, i in 1:3]
  @mput ps
  #@matlab plot3(ps[:,1],ps[:,2],ps[:,3],'.')
  mat"plot3(ps(:,1),ps(:,2),ps(:,3),'.')"
end

function adm(b)
    nmin = 20
    η = 1.5
    I = b[1][1].begin_idx : b[1][1].end_idx-1
    J = b[2][1].begin_idx : b[2][1].end_idx-1
    length(I) < nmin && return true
    length(J) < nmin && return true
    ll1, ur1 = ClusterTrees.boundingbox(p[I]); c1 = (ll1+ur1)/2;
    ll2, ur2 = ClusterTrees.boundingbox(q[J]); c2 = (ll2+ur2)/2;
    diam1 = norm(ur1-c1)
    diam2 = norm(ur2-c2)
    dist12 = norm(c2-c1)
    @show (diam1, diam2, dist12)
    return dist12 >= η*max(diam1, diam2)
end

# Plot the pair of clusters making up block n
tb = ClusterTrees.admissable_partition((tp,tq),adm)
n = 1
τ, σ = tb[n][1][1], tb[n][2][1]
@matlab figure()
@matlab hold("on")
ps = [p[i] for p in p[τ.begin_idx:τ.end_idx-1], i in 1:3]
qs = [q[i] for q in q[σ.begin_idx:σ.end_idx-1], i in 1:3]
@mput ps qs
mat"""
plot3(ps(:,1),ps(:,2),ps(:,3),'.')
plot3(qs(:,1),qs(:,2),qs(:,3),'.')
axis('equal','square')
"""
# @matlab begin
#   plot3(ps[:,1],ps[:,2],ps[:,3],'.')
#   plot3(qs[:,1],qs[:,2],qs[:,3],'.')
#   axis("equal","square")
# end

# Visualise the blocktree by filling in the matrix with random entries,
# one per block.
A = zeros(numfunctions(tfs), numfunctions(bfs))
for b in tb
    τ, σ = b[1][1], b[2][1]
    I = τ.begin_idx : τ.end_idx-1
    J = σ.begin_idx : σ.end_idx-1
    A[I,J] += rand(1:10)
end

using Plots
plotlyjs()
heatmap(A)


text:
__precompile__()

"""

# GeoEfficiency Package
introduce a fast and flexible tool to calculate in batch or individually the `geometrical efficiency`
for a set of common radiation detectors shapes (cylindrical,Bore-hole, Well-type) as seen form
a source. The source can be a point, a disc or even a cylinder.

# Quick Usage
*  geoEff()	: Calculate the geometrical efficiency for one geometrical setup return only the value of the geometrical efficiency.\n
*  calc() 	: Calculate the geometrical efficiency for one geometrical setup and display full information on the console.\n
*  calcN()	: Calculate the geometrical efficiency for geometrical setup(s) and display full information on the console until the user quit.\n
*  batch()	: Calculate the geometrical efficiency using data in the **``$dataDir``** folder in batch mode.

**for more information and updates refer to the repository at [``GitHub.com``](https://github.com/DrKrar/GeoEfficiency.jl/)**

"""
module GeoEfficiency

export # Config
    about,

 # Input_Interface
	getDetectors,
	setSrcToPoint,
	typeofSrc,

 # Physics_Model
	Point,
	source,
	Detector,
	CylDetector,
	BoreDetector,
	WellDetector,

 # Calculations
	geoEff,

 # Output_Interface
	calc,
	calcN,
	batch,
	batchInfo

include("Config.jl")
include("Physics_Model.jl")
include("Input_Interface.jl")
include("Calculations.jl")
include("Output_Interface.jl")

about()
printstyled("""
\n\tBatch mode
\t-  read files from directory `$dataDir`
\t-  save results to directory `$resultdir`
\n\tfor more information see `batch`, `prepare_batch`.
\n"""
,  color=:green)

end #module


text:
# ----------- TUPLES -------------
#
# Equality between tuples of float64 is defined at 10^-6 rouding error
≈(a::Tuple{Float64,Float64},b::Tuple{Float64,Float64};atol::Real=10^-6)=≈(a[1],b[1],atol=atol)&&≈(a[2],b[2],atol=atol)
≈(a::Array{Tuple{Float64,Float64}},b::Array{Tuple{Float64,Float64}};atol::Real=10^-6)=length(a)==length(b) ? all([x[1]≈x[2] for x in zip(a,b)]) : false

mutable struct Node <: StreamObject
    name::AbstractString
    presence::Intervals
end

# ----------- LINK DEFINITIONS -------------
#
mutable struct Link <: StreamObject
    name::AbstractString
    presence::Intervals
    from::AbstractString
    to::AbstractString
    weight::Float64
end

# --- Operations on StreamObjects ---
⊆(i::Intervals,o::StreamObject)=i ⊆ o.presence
⊆(i::Tuple{Float64,Float64},o::StreamObject)=i ⊆ o.presence
⊆(a::StreamObject,b::StreamObject)=a.presence ⊆ b.presence
∈(t::Float64,o::StreamObject)=t ∈ o.presence
∩(a::StreamObject,b::StreamObject)=a.presence ∩ b.presence
∪(a::StreamObject,b::StreamObject)=a.presence ∪ b.presence

# --- Operations on Nodes ---
function merge(n1::Node,n2::Node)
    if n1.name!=n2.name
        throw("Cannot merge nodes with different names...")
    end
    return Node(n1.name, n1 ∪ n2)
end

function merge!(n1::Node,n2::Node)
    if n1.name!=n2.name
        throw("Cannot merge nodes with different names...")
    end
    n1.presence = n1 ∪ n2
end

# --- Operations on Links ---
from_match(l1::Link,l2::Link)=(l1.from==l2.from)
from_match(name::AbstractString,l::Link)=(name==l.from)
to_match(l1::Link,l2::Link)=(l1.to==l2.to)
to_match(name::AbstractString,l::Link)=(name==l.to)
match(l1::Link,l2::Link)=from_match(l1,l2)&to_match(l1,l2)
match(from::AbstractString,to::AbstractString,l::Link)=from_match(from,l)&to_match(to,l)

function merge(l1::Link,l2::Link)
    if !match(l1,l2)
        throw("Cannot merge links with different end points...")
    end
    return Link(l1.name, l1 ∪ l2, l1.from, l1.to, l1.weight)
end

function merge!(l1::Link,l2::Link)
    if !match(l1,l2)
        throw("Cannot merge links with different end points...")
    end
    l1.presence = l1 ∪ l2
end

# --- Operations on Vectors of Nodes ---
get_idx(n::Node,a::Vector{Node})=findall(i->i==n,a)
get_idx(name::AbstractString,a::Vector{Node})=findall(i->i.name==name,a)
get_idx(t::Float64,a::Vector{Node})=findall(i->t ∈ i,a)

function ∪(a::Vector{Node},b::Vector{Node})
    c = Vector(Node[])
    for aa in a
        idx = get_idx(aa.name,b)
        if length(idx)==0
            new = Node(aa.name, aa.presence)
        elseif length(idx)==1
            new = Node(aa.name, aa ∪ b[idx][1])
        else
            throw("More than one node named $aa.name in array.")
        end
        push!(c,new)
    end
    for bb in b
        idx = get_idx(bb,a)
        if length(idx)==0
            new = Node(bb.name, bb.presence)
            push!(c,new)
        elseif length(idx)!=1
            throw("More than one node named $aa.name in array.")
        end
    end
    return c
end

function ∩(a::Vector{Node},b::Vector{Node})
    c = Vector(Node[])
    for aa in a
        idx = get_idx(aa.name,b)
        if length(idx)==1
            new = Node(aa.name, aa ∩ b[idx][1])
            push!(c,new)
        elseif length(idx)!=0
            throw("More than one node named $aa.name in array.")
        end
    end
    return c
end

# --- Operations on Vectors of Links ---
from_match(l1::Link,l::Vector{Link})=findall(x->from_match(x,l1),l)
from_match(name::AbstractString,l::Vector{Link})=findall(x->from_match(name,x),l)
to_match(l1::Link,l::Vector{Link})=findall(x->to_match(x,l1),l)
to_match(name::AbstractString,l::Vector{Link})=findall(x->to_match(name,x),l)
match(l1::Link,l::Vector{Link})=findall(x->match(x,l1),l)
match(from::AbstractString,to::AbstractString,l::Vector{Link})=findall(x->match(from,to,x),l)

get_idx(l::Link,a::Vector{Link})=findall(i->i==l,a)
get_idx(name::AbstractString,a::Vector{Link})=findall(i->i.name==name,a)
get_idx(t::Float64,a::Vector{Link})=findall(i->t ∈ i,a)

function ∪(a::Vector{Link},b::Vector{Link})
    c = Vector(Link[])
    for aa in a
        idx = match(aa,b)
        if length(idx)==0
            new = Link(aa.name, aa.presence, aa.from, aa.to, aa.weight)
        elseif length(idx)==1
            new = Link(aa.name, aa ∪ b[idx][1], aa.from, aa.to, aa.weight)
        else
            throw("More than one link matching from $aa.from and to $aa.to in array.")
        end
        push!(c,new)
    end
    for bb in b
        idx = match(bb,a)
        if length(idx)==0
            new = Link(bb.name, bb.presence, bb.from, bb.to, bb.weight)
            push!(c,new)
        elseif length(idx)!=1
            throw("More than one link matching from $aa.from and to $aa.to in array.")
        end
    end
    return c
end

function ∩(a::Vector{Link},b::Vector{Link})
    c = Vector(Link[])
    for aa in a
        idx = match(aa,b)
        if length(idx)==1
            new = Link(aa.name, aa ∩ b[idx][1], aa.from, aa.to, aa.weight)
            push!(c,new)
        elseif length(idx)!=0
            throw("More than one link matching from $aa.from and to $aa.to in array.")
        end
    end
    return c
end

# ----------- STREAM DEFINITIONS -------------
#
struct LinkStream <: AbstractUndirectedStream
    name::AbstractString
    T::Intervals
    V::Set{AbstractString}
    E::Dict{AbstractString,Dict{AbstractString,Link}}
end

struct DirectedLinkStream <: AbstractDirectedStream
    name::AbstractString
    T::Intervals
    V::Set{AbstractString}
    E::Dict{AbstractString,Dict{AbstractString,Link}}
end

#Base.getindex(ls::LinkStream, name::AbstractString)=haskey(ls.E,name) ? ls.E[name] : []
#Base.getindex(ls::LinkStream, t::Float64)=[l for (k,v) in ls.E for (kk,l) in v if t ∈ l]

struct StreamGraph <: AbstractUndirectedStream
    name::AbstractString
    T::Intervals
    V::Set{AbstractString}
    W::Dict{AbstractString,Node}
    E::Dict{AbstractString,Dict{AbstractString,Link}}
end

struct DirectedStreamGraph <: AbstractDirectedStream
    name::AbstractString
    T::Intervals
    V::Set{AbstractString}
    W::Dict{AbstractString,Node}
    E::Dict{AbstractString,Dict{AbstractString,Link}}
end

LinkStream(name::AbstractString) = LinkStream(name, Intervals([]), Set(), Dict())
LinkStream(name::AbstractString,T::Intervals) = LinkStream(name, T, Set(), Dict())

DirectedLinkStream(name::AbstractString) = DirectedLinkStream(name, Intervals([]), Set(), Dict())
DirectedLinkStream(name::AbstractString,T::Intervals) = DirectedLinkStream(name, T, Set(), Dict())

StreamGraph(name::AbstractString) = StreamGraph(name, Intervals([]), Set(), Dict(), Dict())
StreamGraph(name::AbstractString,T::Intervals) = StreamGraph(name, T, Set(), Dict(), Dict())

function neighborhood(s::AbstractUndirectedStream, node::AbstractString)
    N=Dict{AbstractString,Intervals}()
    for l in links(s,node)
        if l.from==node
            if haskey(N,l.to)
                N[l.to]=N[l.to] ∪ (times(s,l.to) ∩ l.presence)
            else
                N[l.to]=times(s,l.to) ∩ l.presence
            end
        elseif l.to==node
            if haskey(N,l.from)
                N[l.from]=N[l.from] ∪ (times(s,l.from) ∩ l.presence)
            else
                N[l.from]=times(s,l.from) ∩ l.presence
            end
        end
    end
    Dict{AbstractString,Node}([k=>Node(k,v) for (k,v) in N])
end

function neighborhood(s::AbstractUndirectedStream, node::AbstractString, t::Float64)
    neighbors=AbstractString[]
    for l in links(s,node)
        if t ∈ l
            if l.from==node
                push!(neighbors,l.to)
            elseif l.to==node
                push!(neighbors,l.from)
            end
        end
    end
    neighbors
end

#Base.getindex(s::StreamGraph, name::AbstractString)=ls.W[name]
#Base.getindex(s::AbstractDirectedStream, from::AbstractString, to::AbstractString)=s.E[from][to]
#Base.getindex(s::AbstractUndirectedStream, from::AbstractString, to::AbstractString)=(haskey(s.E,from)&haskey(s.E[from],to)) ? s.E[from][to] : s.E[to][from]

# ----------- OPERATIONS ON STREAMS -------------
#
∈(t::Float64,s::AbstractStream)=t ∈ s.T
∈(n::Node,s::StreamGraph)=haskey(s.W,n.name) & s.W[n.name] == n
∈(n::Node,s::DirectedStreamGraph)=haskey(s.W,n.name) & s.W[n.name] == n
∈(l::Link,s::AbstractStream)=haskey(s.E,l.from) & haskey(s.E[l.from],l.to) & s.E[l.from][l.to] == l
∈(n::AbstractString,ls::Union{LinkStream,DirectedLinkStream})=n ∈ ls.V
∈(n::AbstractString,s::Union{StreamGraph,DirectedStreamGraph})=(n ∈ s.V) & haskey(s.W,n)

⊆(t::Tuple{Float64,Float64},s::AbstractStream)=t ⊆ s.T
# TO UPDATE -----
⊆(ls1::LinkStream,ls2::LinkStream)=(ls1.T ⊆ ls2.T)&(ls1.V ⊆ ls2.V)&(ls1.E ⊆ ls2.E)
⊆(s1::StreamGraph,s2::StreamGraph)=(s1.T ⊆ s2.T)&(s1.V ⊆ s2.V)&(s1.W ⊆ s2.W)&(s1.E ⊆ s2.E)

∩(ls1::LinkStream,ls2::LinkStream)=LinkStream("$ls1.name n $ls2.name", ls1.T ∩ ls2.T, ls1.V ∩ ls2.V, ls1.E ∩ ls2.E)
∩(s1::StreamGraph,s2::StreamGraph)=StreamGraph("$s1.name n $s2.name", s1.T ∩ s2.T, s1.V ∩ s2.V, s1.W ∩ s2.W, s1.E ∩ s2.E)

∪(ls1::LinkStream,ls2::LinkStream)=LinkStream("$ls1.name u $ls2.name", ls1.T ∪ ls2.T, ls1.V ∪ ls2.V, ls1.E ∪ ls2.E)
∪(s1::StreamGraph,s2::StreamGraph)=StreamGraph("$s1.name u $s2.name", s1.T ∪ s2.T, s1.V ∪ s2.V, s1.W ∪ s2.W, s1.E ∪ s2.E)
# -----

struct NodeEvent <: Event
    t::Float64
    arrive::Bool
    object::AbstractString
end

struct LinkEvent <: Event
    t::Float64
    arrive::Bool
    object::Tuple{AbstractString,AbstractString}
end

==(e1::Event,e2::Event)=(e1.t==e2.t)&(e1.arrive==e2.arrive)&(e1.object==e2.object)

# ----------- STATE -------------
#
mutable struct State
    t0::Float64
    t1::Float64
    nodes::Set{AbstractString}
    links::Set{Tuple{AbstractString,AbstractString}}
end

# ----------- TRANSITION -------------
#
mutable struct Transition
    tprev::Float64
    t::Float64
    tnxt::Float64
    node_arrivals::Set{AbstractString}
    node_departures::Set{AbstractString}
    link_arrivals::Set{Tuple{AbstractString,AbstractString}}
    link_departures::Set{Tuple{AbstractString,AbstractString}}
end

Δnodes(t::Transition)=length(t.node_arrivals)-length(t.node_departures)
Δlinks(t::Transition)=length(t.link_arrivals)-length(t.link_departures)

function apply!(s::State,τ::Transition)
    # Moving backward
    if s.t0 == τ.t
        s.t0=τ.tprev
        s.t1=τ.t
        length(τ.node_arrivals - s.nodes)>0 && throw("<Backward Move> Absent nodes cannot leave the state.")
        s.nodes=s.nodes - τ.node_arrivals
        length(s.nodes ∩ τ.node_departures)>0 && throw("<Backward Move> Nodes already present in state cannot arrive in Transition.")
        s.nodes=s.nodes ∪ τ.node_departures
        length(τ.link_arrivals - s.links)>0 && throw("<Backward Move> Absent links cannot leave the state.")
        s.links=s.links - τ.link_arrivals
        length(s.links ∩ τ.link_arrivals)>0 && throw("<Backward Move> Links already present in state cannot arrive in Transition.")
        s.links=s.links ∪ τ.link_departures
    # Moving forward
    elseif s.t1 == τ.t
        s.t0=τ.t
        s.t1=τ.tnxt
        length(τ.node_departures - s.nodes)>0 && throw("<Forward Move> Absent nodes cannot leave the state.")
        s.nodes=s.nodes - τ.node_departures
        length(s.nodes ∩ τ.node_arrivals)>0 && throw("<Forward Move> Nodes already present in state cannot arrive in Transition.")
        s.nodes=s.nodes ∪ τ.node_arrivals
        length(τ.link_departures - s.links)>0 && throw("<Forward Move> Absent links cannot leave the state.")
        s.links=s.links - τ.link_departures
        length(s.links ∩ τ.link_arrivals)>0 && throw("<Forward Move> Links already present in state cannot arrive in Transition.")
        s.links=s.links ∪ τ.link_arrivals
    else
        throw("Cannot apply transition to given state.")
    end
end

# ----------- TIME CURSOR -------------
#
mutable struct TimeCursor
    S::State
    T::Dict{Float64,Transition}
end

## Move the cursor in the Stream ##
next_transition(tc::TimeCursor)=haskey(tc.T,tc.S.t1) ? tc.T[tc.S.t1] : throw("No transition in TimeMap at t=$(tc.S.t1)")
previous_transition(tc::TimeCursor)=haskey(tc.T,tc.S.t0) ? tc.T[tc.S.t0] : throw("No transition in TimeMap at t=$(tc.S.t0)")
next!(tc::TimeCursor)=haskey(tc.T,tc.S.t1) && apply!(tc.S,next_transition(tc))
previous!(tc::TimeCursor)=haskey(tc.T,tc.S.t0) && apply!(tc.S,previous_transition(tc))

function goto!(tc::TimeCursor,t::Float64)
    while !(tc.S.t0 <= t < tc.S.t1)
        t >= tc.S.t1 ? next!(tc) : previous!(tc)
    end
end

function start!(tc::TimeCursor)
    while haskey(tc.T,tc.S.t0)
        previous!(tc)
    end
    next!(tc)
end

function end!(tc::TimeCursor)
    while haskey(tc.T,tc.S.t1)
        next!(tc)
    end
    previous!(tc)
end

function load!(tc::TimeCursor,events::Vector{Event})
    tprev::Float64=events[1].t-0.1
    t::Float64=events[1].t
    tc.S=State(tprev,t,Set(),Set())
    idx::Int64=1
    node_arrivals=Set{AbstractString}()
    node_departures=Set{AbstractString}()
    link_arrivals=Set{Tuple{AbstractString,AbstractString}}()
    link_departures=Set{Tuple{AbstractString,AbstractString}}()
    transitions=Transition[]
    while idx <= length(events)
        while (idx <= length(events)) && (events[idx].t==t)
            if typeof(events[idx])==NodeEvent
                if events[idx].arrive
                    push!(node_arrivals,events[idx].object)
                else
                    push!(node_departures,events[idx].object)
                end
            else
                if events[idx].arrive
                    push!(link_arrivals,events[idx].object)
                else
                    push!(link_departures,events[idx].object)
                end
            end
            idx += 1
        end
        if length(transitions)>0
            transitions[end].tnxt=t
        end
        push!(transitions,Transition(tprev,t,-1,deepcopy(node_arrivals),
                                                deepcopy(node_departures),
                                                deepcopy(link_arrivals),
                                                deepcopy(link_departures)))
        empty!(node_arrivals)
        empty!(node_departures)
        empty!(link_arrivals)
        empty!(link_departures)
        tprev=t
        if idx<=length(events)
            t=events[idx].t
        end
    end
    transitions[end].tnxt=t
    for τ in transitions
        tc.T[τ.t]=τ
    end
end

# ----------- JUMPS -------------
#
struct Jump
    t::Float64
    from::AbstractString
    to::AbstractString
end

struct DurationJump
    t::Float64
    from::AbstractString
    to::AbstractString
    δ::Float64
end

# ----------- PATHS -------------
#
mutable struct Path <: AbstractPath
    jumps::Vector{Jump}
end

Path()=Path([])

mutable struct DurationPath <: AbstractPath
    jumps::Vector{DurationJump}
end

DurationPath()=DurationPath([])

start(p::AbstractPath)= length(p.jumps) > 0 ? p.jumps[1].t : 0
finish(p::Path)=length(p.jumps) > 0 ? p.jumps[end].t : 0
finish(p::DurationPath)=length(p.jumps) > 0 ? p.jumps[end].t + p.jumps[end].δ : 0

function is_valid(p::Path)
    if length(p)<=1
        return true
    end
    for cpt in zip(p.jumps[1:end-1],p.jumps[2:end])
        if cpt[1].t > cpt[2].t || cpt[1].to != cpt[2].from
            return false
        end
    end
    return true
end

function is_valid(p::DurationPath)
    if length(p)<=1
        return true
    end
    for cpt in zip(p.jumps[1:end-1],p.jumps[2:end])
        if cpt[1].t + cpt[1].δ > cpt[2].t || cpt[1].to != cpt[2].from
            return false
        end
    end
    return true
end

function +(p1::T,p2::T) where T <: AbstractPath
    if !is_valid(T([p1.jumps[end],p2.jumps[1]]))
        throw("Cannot concat paths.")
    end
    p = deepcopy(p1)
    for jump in p2.jumps
        push!(p.jumps,jump)
    end
    return p
end

function push(p::AbstractPath,jump::Jump)
    p2 = deepcopy(p)
    push!(p2.jumps,jump)
    if is_valid(p2)
        return p2
    else
        throw("Cannot add jump $jump to path because this will break the path feature.")
    end
end

function ⊆(p1::T,p2::T) where T <: AbstractPath
    lp1 = length(p1)
    lp2 = length(p2)
    if lp1 > lp2
        return false
    end
    for i in range(1,lp2-lp1)
        if p2.jumps[i:i+lp1-1] == p1.jumps
            return true
        end
    end
    return false
end


text:
#=
cat.jl
Concatenation support for LinearMapAX
2018-01-19, Jeff Fessler, University of Michigan
=#

export lmaa_hcat, lmaa_vcat, lmaa_hvcat

using LinearAlgebra: UniformScaling, I

#=
function warndim(A::LinearMapAX)
    ((length(A._idim) > 1) || (length(A._odim) > 1)) && @warn("dim ignored")
    nothing
end
=#

# cat (hcat, vcat, hvcat) are tricky for avoiding type piracy
# It is especially hard to handle AbstractMatrix,
# so typically force the user to wrap it in LinearMap(AX) first.
#LMcat{T} = Union{LinearMapAM{T}, LinearMap{T}, UniformScaling{T},AbstractMatrix{T}}
LMcat{T} = Union{LinearMapAX{T}, LinearMap{T}, UniformScaling{T}} # settle
#LMelse{T} = Union{LinearMap{T},UniformScaling{T},AbstractMatrix{T}} # non AM

# convert to something suitable for LinearMap.*cat
function lm_promote(A::LMcat)
#    @show typeof(A)
    A isa LinearMapAX ? A._lmap :
    A isa UniformScaling ? A : # leave unchanged - ok for LinearMaps.*cat
    A # otherwise it is this
#    throw("bug") # should only be only of LMcat types
#    A isa AbstractMatrix ? LinearMap(A) :
#    A isa LinearMap ?
end

# single-letter codes for cat objects, e.g., [A I A] becomes "AIA"
lm_code(::LinearMapAM) = "A"
lm_code(::LinearMapAO) = "O"
lm_code(::UniformScaling) = "I"
lm_code(::LinearMap) = "L"
#lm_code(::AbstractMatrix) = "M"
#lm_code(::Any) = "?"

# concatenate the single-letter codes, e.g., [A I A] becomes "AIA"
lm_name = As -> *(lm_code.(As)...)
# lm_name = As -> nothing

# these rely on LinearMap.*cat methods
"`B = lmaa_hcat(A1, A2, ...)` `hcat` of multiple objects"
lmaa_hcat(As::LMcat...) =
    LinearMapAA(hcat(lm_promote.(As)...), (hcat=lm_name(As),))
"`B = lmaa_vcat(A1, A2, ...)` `vcat` of multiple objects"
lmaa_vcat(As::LMcat...) =
    LinearMapAA(vcat(lm_promote.(As)...), (vcat=lm_name(As),))
"`B = lmaa_hvcat(rows, A1, A2, ...)` `hvcat` of multiple objects"
lmaa_hvcat(rows::NTuple{nr,Int} where {nr}, As::LMcat...) =
    LinearMapAA(hvcat(rows, lm_promote.(As)...), (hvcat=lm_name(As),))

# a single leading LinearMapAM followed by others is clear
Base.hcat(A1::LinearMapAX, As::LMcat...) =
    lmaa_hcat(A1, As...)
Base.vcat(A1::LinearMapAX, As::LMcat...) =
    lmaa_vcat(A1, As...)
Base.hvcat(rows::NTuple{nr,Int} where {nr}, A1::LinearMapAX, As::LMcat...) =
    lmaa_hvcat(rows, A1, As...)
# or in 2nd position, beyond that, user can use lmaa_*
#Base.hcat(A1::LMelse, A2::LinearMapAM, As::LMcat...) = # fails!?
Base.hcat(A1::UniformScaling, A2::LinearMapAX, As::LMcat...) =
    lmaa_hcat(A1, A2, As...)
Base.vcat(A1::UniformScaling, A2::LinearMapAX, As::LMcat...) =
    lmaa_vcat(A1, A2, As...)
Base.hvcat(rows::NTuple{nr,Int} where nr,
        A1::UniformScaling, A2::LinearMapAX, As::LMcat...) =
    lmaa_hvcat(rows, A1, A2, As...)

# special handling for solely AO collections

function _hcat(tryop::Bool, As::LinearMapAO...)
    B = LinearMaps.hcat(map(A -> A._lmap, As)...)
    nblock = length(As)
    prop = (nblock = nblock,)
    if (tryop && nblock > 1 && same_idim(As...) && same_odim(As...)) # operator version
        return LinearMapAA(B ; prop=prop,
            idim = (As[1]._idim..., nblock),
            odim = As[1]._odim,
        )
    end
    return LinearMapAA(B ; prop=prop)
end

function _vcat(tryop::Bool, As::LinearMapAO...)
    B = LinearMaps.vcat(map(A -> A._lmap, As)...)
    nblock = length(As)
    prop = (nblock = nblock,)
    if (tryop && nblock > 1 && same_idim(As...) && same_odim(As...)) # operator version
        return LinearMapAA(B ; prop=prop,
            idim = As[1]._idim,
            odim = (As[1]._odim..., nblock),
        )
    end
    return LinearMapAA(B ; prop=prop)
end

"""
    hcat(As::LinearMapAO... ; tryop::Bool=true)
`hcat` with (by default) attempt to append `nblock` to `idim` if consistent blocks.
"""
Base.hcat(A1::LinearMapAO, As::LinearMapAO... ; tryop::Bool=true) =
    _hcat(tryop, A1, As...)

"""
    vcat(As::LinearMapAO... ; tryop::Bool=true)
`vcat` with (by default) attempt to append `nblock` to `odim` if consistent blocks.
"""
Base.vcat(A1::LinearMapAO, As::LinearMapAO... ; tryop::Bool=true) =
    _vcat(tryop, A1, As...)

"""
    hvcat(rows, As::LinearMapAO...)
`hvcat` that discards special `idim` and `odim` information (too hard!) # todo?
"""
Base.hvcat(rows::NTuple{nr,Int} where nr,
        A1::LinearMapAO, As::LinearMapAO...) =
    lmaa_hvcat(rows, undim(A1), undim.(As)...)


text:
module internal
    using Logging, DataStructures

    include("../structs/DataSet.jl")
    include("../functions/Data.jl")

    #set sok and parse data type
    function set_sok(instance, dat_sori, hol_sori)
        key_list = dat_sori ; sok_dict = Dict(dat => [] for dat in dat_sori) ; sok_type = Dict()
        for hol in hol_sori
            for i in eachindex(key_list)
                append!(sok_dict[key_list[i]], hol[i])
            end
        end
        for key in key_list
            sok_type[key], sok_dict[key] = Data.parse_data_list(sok_dict[key])["type"], Data.parse_data_list(sok_dict[key])["data"]
        end

        new_hol = []
        for hol in hol_sori
            tmp_hol = []
            for i in eachindex(key_list)
                if sok_type[key_list[i]] == String
                    append!(tmp_hol, string(hol[i]))
                else
                    append!(tmp_hol, parse(sok_type[key_list[i]], string(hol[i])))
                end
            end

            append!(new_hol, [tmp_hol])
        end
        instance.hol_sori = new_hol

        return [
            DataStructures.OrderedDict(key_list[i]=> DataSet.ColumnSet(sok_dict[key_list[i]], i, key_list[i], instance) for i in eachindex(key_list)),
            [DataSet.RowSet(DataStructures.OrderedDict(key => sok_dict[key][i] for key in key_list), i, instance) for i in eachindex(hol_sori)]
        ]
    end

    #modify
    function insert(instance, dat_pos, dat_name, hol_value)
        if dat_pos < 1 || dat_pos > length(instance.dat_sori) + 1
            Logging.ERROR("insert position is invalid.")
            result = false
        else
            new_dat_sori = insert!(instance.dat_sori, dat_pos, dat_name)
            new_hol_sori = [insert!(instance.hol_sori[i], dat_pos, hol_value[i]) for i in eachindex(instance.hol_sori)]
            new_sok_dict, new_sok_arr = set_sok(instance, new_dat_sori, new_hol_sori)

            instance.dat_sori, instance.hol_sori, instance._sok_dict, instance._sok_arr = new_dat_sori, new_hol_sori, new_sok_dict, new_sok_arr

            result = true
        end

        return result
    end

    function append(instance, dat_name, hol_value)
        return insert(instance, length(instance.dat_sori) + 1, dat_name, hol_value)
    end

    #simple selection and query
    function select(instance, key)
        key = string(key)
        if !isnull(tryparse(Int, key))
            result = instance._sok_arr[parse(Int, key)]
        else
            result = instance._sok_dict[key]
        end

        return result
    end

    function where(instance, where_string, drop_null)
        new_dat_sori, new_hol_sori = instance.dat_sori, []
        for row_set in instance._sok_arr
            where_result = row_set.where(where_string, drop_null = drop_null)
            if where_result != nothing
                append!(new_hol_sori, [where_result])
            end
        end

        new_instance = Base.deepcopy(instance)
        new_sok_dict, new_sok_arr = set_sok(new_instance, new_dat_sori, new_hol_sori)
        new_instance.dat_sori, new_instance.hol_sori, new_instance._sok_dict, new_instance._sok_arr = new_dat_sori, new_hol_sori, new_sok_dict, new_sok_arr
        new_instance.name = string(new_instance.name, ": ", where_string)

        return new_instance
    end

    function equation(instance, equation_string, inplace)
        equation_to = strip(split(equation_string, "=")[1])
        new_dat_sori, new_hol_sori = in(equation_to, instance.dat_sori) ? instance.dat_sori : append!(instance.dat_sori, equation_to), []
        for row_set in instance._sok_arr
            equation_result = row_set.equation(equation_string)
            if equation_result != nothing
                append!(new_hol_sori, [equation_result])
            end
        end

        if inplace
            new_instance = Base.deepcopy(instance)
        else
            new_instance = instance
        end

        new_sok_dict, new_sok_arr = set_sok(new_instance, new_dat_sori, new_hol_sori)
        new_instance.dat_sori, new_instance.hol_sori, new_instance._sok_dict, new_instance._sok_arr = new_dat_sori, new_hol_sori, new_sok_dict, new_sok_arr
        new_instance.name = string(new_instance.name, ": ", equation_string)

        if inplace
            println(new_instance.to_string())
        else
            return new_instance
        end
    end

    #descriptive statistics
    function describe(instance)
        result = DataStructures.OrderedDict()
        for dat in instance.dat_sori
            result[dat] = instance.select(dat).describe()
        end

        return result
    end
end


text:
export dcs_exhaustive
using Devectorize

function dcs_exhaustive{T}(G::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(G)
  n = size(G[1], 1)

  prev_state = zeros(Bool, n)
  prev_nodes = zeros(Int64, n)
  prev_num_nodes = 0

  binary_state = zeros(Bool, n)

  curr_state = zeros(Bool, n)
  curr_nodes = zeros(Int64, n)
  curr_num_nodes = 0

  opt_idx = 0
  opt_deg = 0.0

  curr_ne = zeros(m, 1)
  curr_deg = zeros(m, 1)

  for j=[1:(2^n-1);]
    int_to_binary_bool!(binary_state, j, n)
    binary_to_gray!(curr_state, binary_state)
    curr_num_nodes = find_set_bits!(curr_nodes, curr_state)

    # Find which node has been added or deleted
    changed_idx = 1
    while (curr_state[changed_idx] == prev_state[changed_idx])
      changed_idx = changed_idx + 1
    end

    # Has the node been added or deleted?
    if curr_state[changed_idx]
      m_sign = 1
    else
      m_sign = -1
    end
    @inbounds begin
      for k=[1:m;]
        delta_deg = 0
        if prev_num_nodes > 0
          for ll=[1:prev_num_nodes;]
            other_node = prev_nodes[ll]
            delta_deg = delta_deg + G[k][other_node, changed_idx] + G[k][changed_idx, other_node]
          end
        end
        curr_ne[k] = curr_ne[k] + m_sign  * delta_deg
        curr_deg[k] = curr_ne[k] / curr_num_nodes
        #      println("G[$k] j: $j: $(curr_ne[k]), $(curr_deg[k])")
      end
    end
    @devec curr_min_deg = minimum(curr_deg)
    if curr_min_deg > opt_deg
      opt_deg = curr_min_deg
      opt_idx = j
    end

    for k=[1:n;]
      prev_state[k] = curr_state[k]
      prev_nodes[k] = curr_nodes[k]
    end
    prev_num_nodes = curr_num_nodes
  end
  int_to_binary_bool!(binary_state, opt_idx, n)
  binary_to_gray!(curr_state, binary_state)
  curr_num_nodes = find_set_bits!(curr_nodes, curr_state)
  opt_S = curr_nodes[1:curr_num_nodes]
  # println("S: $opt_S min_deg: $opt_deg")
  return opt_S , opt_deg
end


text:
mutable struct DetailedOrganism
    id::Int

    genotype::Genotype
    parent::Int

    position::Position
    environment::Environment

    time_birth::Int
    time_death::Int

    num_offsprings::Int
    num_identical_offsprings::Int

    population::Dict{Genotype, Int}

    function DetailedOrganism(organism::Organism, num_offsprings, num_identical_offsprings, population)
        new(get_id(organism), organism.genotype, get_id(organism.parent), organism.position, organism.environment,
            organism.time_birth, organism.time_death, num_offsprings, num_identical_offsprings, population)
    end
end

get_id(organism::DetailedOrganism) = organism.id
get_fitness(organism::DetailedOrganism) = organism.num_identical_offsprings

abstract type PropertyType end
struct GenotypeType <: PropertyType end
struct PositionType <: PropertyType end
struct EnvironmentType <: PropertyType end
struct FitnessType <: PropertyType end
struct NumOffspringsType <: PropertyType end
struct PopulationType <: PropertyType end

get(organism::DetailedOrganism, ::GenotypeType) = organism.genotype
get(organism::DetailedOrganism, ::PositionType) = organism.position
get(organism::DetailedOrganism, ::EnvironmentType) = organism.environment
get(organism::DetailedOrganism, ::FitnessType) = organism.num_identical_offsprings
get(organism::DetailedOrganism, ::NumOffspringsType) = organism.num_offsprings
get(organism::DetailedOrganism, ::PopulationType) = organism.population

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using MPCHelper
using Test


text:
module SimulationLogger

using MacroTools
using BangBang

export @log, @onlylog, @Loggable, @nested_log, @nested_onlylog, __LOG_INDICATOR__

include("macros.jl")

end


text:
function Base.show(io::IO, cms::Array{T, 1}) where T <: ConfusionMatrix
    for i in eachindex(cms)
        cm = cms[i]
        println(io, "$i --> TP: $(cm.tp) FP: $(cm.fp) TN: $(cm.tn) FN $(cm.fn) [P = :$(cm.pos_class), N = :$(cm.neg_class)]")
    end
end

function listfiles(path)
    files = readdir(path)
    res = []
    for f in joinpath.(path, files)
        if isfile(f)
            endswith(f, ".json") && push!(res, f)
        else
            res = vcat(res, listfiles(f))
        end
    end
    return res
end

findResultFile(path, hash) = filter(x -> occursin(hash, x), listfiles(path))[1]

getResultFromFile(path, hash) = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parsefile(findResultFile(path, hash)))

function load(res::Result, iteration::Int; solver = with_optimizer(Gurobi.Optimizer, OutputFlag = 0))
    data_file = res.experiment[:data_file]
    data, labels = load_data(data_file)

    queries = res.al_history[:query_history].values[1:iteration]

    pools = Symbol.(res.experiment[:param][:initial_pools])
    pools[queries] = ifelse.(labels[queries] .== :inlier, :Lin, :Lout)

    kernels = eval.(Meta.parse.(res.experiment[:model][:fitted][:kernel]))
    gamma_strategy = FixedGammaStrategy(kernels)
    C_strategy = FixedCStrategy(res.experiment[:model][:fitted][:model_params][:C])
    init_strategy = SVDD.SimpleSubspaceStrategy(gamma_strategy, C_strategy, gamma_scope=Val(:Subspace))

    subspaces = convert.(Vector{Int}, res.experiment[:model][:param][:subspaces])

    model = SubSVDD(data, subspaces, pools)
    model.adjust_K = res.experiment[:param][:adjust_K]

    v_Lin = res.experiment[:model][:fitted][:model_params][:weight_update_strategy][:v_Lin]
    v_Lout = res.experiment[:model][:fitted][:model_params][:weight_update_strategy][:v_Lout]

    model.v[queries] = ifelse.(labels[queries] .== :inlier, v_Lin, v_Lout)

    SVDD.initialize!(model, init_strategy)
    return data, labels, model
end

#### Subspace Analysis

function n_outlying_subspaces(predictions; sorted=true)
    tmp = map(x -> ifelse.(SVDD.classify.(x) .== :outlier, 1, 0), predictions)
    subspace_outliers = foldl((x,y) -> x .+ y, tmp)
    sort(collect(zip(eachindex(subspace_outliers), subspace_outliers)), by = x -> x[2], rev=true)
end

function outlying_subspaces(observation_id, predictions, subspaces)
    classifications = SVDD.classify(predictions, Val(:Subspace))
    subspaces[map(x-> x[observation_id], classifications).==:outlier]
end

### Publication plotting
COLORS = (inlier = :lightblue, outlier = :red)
SHAPES = (inlier = (inlier = :circle, outlier = :pentagon), outlier = (inlier = :cross, outlier = :+))

@recipe function plot_svdd(m::SVDD.SVDDClassifier, labels::Vector{Symbol}; grid_resolution = 100, axis_overhang = 0.2, db_color=:black)
    grid_range, grid_data = OCALPlots.get_grid(extrema(m.data)..., grid_resolution, axis_overhang)
    grid_scores = reshape(SVDD.predict(m, grid_data), grid_resolution, grid_resolution)
    data_class = SVDD.classify.(SVDD.predict(m, m.data))
    title := "Decision Boundary"
    @series begin
        seriestype := :contourf
        seriescolor --> :greens
        levels := range(0, maximum(grid_scores), length=10)
        grid_range, grid_range, grid_scores
    end

    markeralpha --> 0.8
    markersize --> 5

    for l in [:inlier, :outlier]
        markercolor := COLORS[l]
        for dc in [:inlier, :outlier]
            markershape := SHAPES[dc][l]
            @series begin
                seriestype := :scatter
                label := OCALPlots.get_legend_text(l, dc)
                OCALPlots.split_2d_array(m.data, (labels.==l) .& (data_class .== dc))
            end
       end
    end

    @series begin
        seriestype := :contour
        levels := [0]
        linewidth := 3
        seriescolor := [db_color]
        cbar:= false
        grid_range, grid_range, grid_scores
    end
end

@recipe function plot_svdd(m::SVDD.SubSVDD, labels::Vector{Symbol}, subspace_idx; grid_resolution = 100, axis_overhang = 0.2, db_color=:black)
    data = m.data[m.subspaces[subspace_idx], :]
    grid_range, grid_data = OCALPlots.get_grid(extrema(data)..., grid_resolution, axis_overhang)
    grid_scores = reshape(SVDD.predict(m, grid_data, subspace_idx), grid_resolution, grid_resolution)
    data_class = SVDD.classify.(SVDD.predict(m, data, subspace_idx))

    title --> "Subspace $(m.subspaces[subspace_idx])"
    @series begin
        seriestype := :contourf
        seriescolor --> :greens
        levels := range(0, maximum(grid_scores), length=10)
        grid_range, grid_range, grid_scores
    end

    markeralpha --> 0.8
    markersize --> 5

    for l in [:inlier, :outlier]
        markercolor := COLORS[l]
        for dc in [:inlier, :outlier]
            markershape := SHAPES[dc][l]
            @series begin
                seriestype := :scatter
                label := OCALPlots.get_legend_text(l, dc)
                OCALPlots.split_2d_array(data, (labels.==l) .& (data_class .== dc))
            end
       end
    end

    @series begin
        seriestype := :contour
        levels := [0]
        linewidth := 3
        seriescolor := [db_color]
        cbar:= false
        grid_range, grid_range, grid_scores
    end
end

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
## Environment and packages
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr()
using JLD2, FileIO
using Statistics
using DelimitedFiles
# Set a random seed for reproduceable behaviour
using Random
Random.seed!(5443)

svname = "HudsonBay"
## Data Preprocessing
# The data has been taken from https://jmahaffy.sdsu.edu/courses/f00/math122/labs/labj/q3v1.htm
# Originally published in E. P. Odum (1953), Fundamentals of Ecology, Philadelphia, W. B. Saunders
hudson_bay_data = readdlm("hudson_bay_data.dat", '\t', Float32, '\n')
# Measurements of prey and predator
Xₙ = Matrix(transpose(hudson_bay_data[:, 2:3]))
t = hudson_bay_data[:, 1] .- hudson_bay_data[1, 1]
# Normalize the data; since the data domain is strictly positive
# we just need to divide by the maximum
xscale = maximum(Xₙ, dims =2)
Xₙ .= 1f0 ./ xscale .* Xₙ
# Time from 0 -> n
tspan = (t[1], t[end])

# Plot the data
scatter(t, transpose(Xₙ), xlabel = "t", ylabel = "x(t), y(t)")
plot!(t, transpose(Xₙ), xlabel = "t", ylabel = "x(t), y(t)")

## Direct Identification via SINDy + Collocation
# Create a collocation
dx̂,x̂ = collocate_data(Xₙ,t, GaussianKernel())
# Look at the collocation
plot(t, dx̂')
# Perform sindy

# Create a Basis
@variables u[1:2]

# Generate the basis functions, multivariate polynomials up to deg 5
# and sine
b = [polynomial_basis(u, 5); sin.(u)]
basis = Basis(b, u)
# Create an optimizer for the SINDy problem
opt = STRRidge()#SR3(Float32(1e-2), Float32(1e-2))
# Create the thresholds which should be used in the search process
λ = Float32.(exp10.(-7:0.1:3))
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)
# Test on derivative data
Ψ = SINDy(x̂, dx̂, basis, λ,  opt, g = g, maxiter = 50000, normalize = true, denoise = true)
println(Ψ)
print_equations(Ψ) # Fails
b2 = Basis((u,p,t)->Ψ(u,ones(length(parameters(Ψ))),t),u, linear_independent = true)
Ψ = SINDy(x̂, dx̂, b2, λ,  opt, g = g, maxiter = 50000, normalize = true, denoise = true)
println(Ψ)
print_equations(Ψ) # Fails
parameters(Ψ)
## UDE Approach
# Subsample the data in y -> initial fitting strategy (batching)
# We assume we have only 5 measurements in y, evenly distributed
ty = collect(t[1]:Float32(t[end]/5):t[end])
# Create datasets for the different measurements
XS = zeros(Float32, length(ty)-1, floor(Int64, mean(diff(ty))/mean(diff(t)))+1) # All x data
TS = zeros(Float32, length(ty)-1, floor(Int64, mean(diff(ty))/mean(diff(t)))+1) # Time data
YS = zeros(Float32, length(ty)-1, 2) # Just two measurements in y

for i in 1:length(ty)-1
    idxs = ty[i].<= t .<= ty[i+1]
    XS[i, :] = Xₙ[1, idxs]
    TS[i, :] = t[idxs]
    YS[i, :] = [Xₙ[2, t .== ty[i]]'; Xₙ[2, t .== ty[i+1]]]
end

## Define the network
# Gaussian RBF as activation
rbf(x) = exp.(-(x.^2))

# Define the network 2->5->5->5->2
U = FastChain(
    FastDense(2,5,rbf), FastDense(5,5, rbf), FastDense(5,5, tanh), FastDense(5,2)
)

# Get the initial parameters, first two is linear birth / decay of prey and predator
p = [rand(Float32,2); initial_params(U)]

# Define the hybrid model
function ude_dynamics!(du,u, p, t)
    û = U(u, p[3:end]) # Network prediction
    # We assume a linear birth rate for the prey
    du[1] = p[1]u[1] + û[1]
    # We assume a linear decay rate for the predator
    du[2] = -p[2]*u[2] + û[2]
end

# Define the problem
prob_nn = ODEProblem(ude_dynamics!,Xₙ[:, 1], tspan, p)

## Function to train the network
# Define a predictor
function predict(θ, X = Xₙ[:,1], T = t)
    Array(solve(prob_nn, Vern7(), u0 = X, p=θ,
                tspan = (T[1], T[end]), saveat = T,
                abstol=1e-6, reltol=1e-6,
                sensealg = ForwardDiffSensitivity()
                ))
end

# Multiple shooting like loss
function shooting_loss(θ)
    # Start with a regularization on the network
    l = convert(eltype(θ), 1e-3)*sum(abs2, θ[3:end]) ./ length(θ[3:end])
    for i in 1:size(XS,1)
        X̂ = predict(θ, [XS[i,1], YS[i,1]], TS[i, :])
        # Full prediction in x
        l += sum(abs2, XS[i,:] .- X̂[1,:])
        # Add the boundary condition in y
        l += abs2(YS[i, 2] .- X̂[2, end])
    end

    return l
end

function loss(θ)
    X̂ = predict(θ)
    sum(abs2, Xₙ - X̂) + convert(eltype(θ), 1e-3)*sum(abs2, θ[3:end]) ./ length(θ[3:end])
end

# Container to track the losses
losses = Float32[]

# Callback to show the loss during training
callback(θ,l) = begin
    push!(losses, l)
    if length(losses)%5==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

## Training -> First shooting / batching to get a rough estimate

# First train with ADAM for better convergence -> move the parameters into a
# favourable starting positing for BFGS
res1 = DiffEqFlux.sciml_train(shooting_loss, p, ADAM(0.1f0), cb=callback, maxiters = 100)
println("Training loss after $(length(losses)) iterations: $(losses[end])")
# Train with BFGS to achieve partial fit of the data
res2 = DiffEqFlux.sciml_train(shooting_loss, res1.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 200)
println("Training loss after $(length(losses)) iterations: $(losses[end])")
# Full L2-Loss for full prediction
res3 = DiffEqFlux.sciml_train(loss, res2.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 10000)
println("Final training loss after $(length(losses)) iterations: $(losses[end])")

pl_losses = plot(1:101, losses[1:101], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "ADAM (Shooting)", color = :blue)
plot!(102:302, losses[102:302], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "BFGS (Shooting)", color = :red)
plot!(302:length(losses), losses[302:end], color = :black, label = "BFGS (L2)")
savefig(pl_losses, joinpath(pwd(), "plots", "$(svname)_losses.pdf"))

# Rename the best candidate
p_trained = res3.minimizer

## Analysis of the trained network
# Interpolate the solution
tsample = t[1]:0.5:t[end]
X̂ = predict(p_trained, Xₙ[:,1], tsample)
# Trained on noisy data vs real solution
pl_trajectory = scatter(t, transpose(Xₙ), color = :black, label = ["Measurements" nothing], xlabel = "t", ylabel = "x(t), y(t)")
plot!(tsample, transpose(X̂), color = :red, label = ["UDE Approximation" nothing])
savefig(pl_trajectory, joinpath(pwd(), "plots", "$(svname)_trajectory_reconstruction.pdf"))

# Neural network guess
Ŷ = U(X̂,p_trained[3:end])

pl_reconstruction = scatter(tsample, transpose(Ŷ), xlabel = "t", ylabel ="U(x,y)", color = :red, label = ["UDE Approximation" nothing])
plot!(tsample, transpose(Ŷ), color = :red, lw = 2, style = :dash, label = [nothing nothing])
savefig(pl_reconstruction, joinpath(pwd(), "plots", "$(svname)_missingterm_reconstruction.pdf"))
pl_missing = plot(pl_trajectory, pl_reconstruction, layout = (2,1))
savefig(pl_missing, joinpath(pwd(), "plots", "$(svname)_reconstruction.pdf"))
## Symbolic regression via sparse regression ( SINDy based )

# Create a Basis
@variables u[1:2]

# Generate the basis functions, multivariate polynomials up to deg 5
# and sine
b = [polynomial_basis(u, 5); sin.(u)]
basis = Basis(b, u)

# Create an optimizer for the SINDy problem
opt = STRRidge()
# Create the thresholds which should be used in the search process
λ = Float32.(exp10.(-7:0.1:3))
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)

# Test on uode derivative data
println("SINDy on learned, partial, available data")
Ψ = SINDy(X̂, Ŷ, basis, λ,  opt, g = g, maxiter = 50000, normalize = true, denoise = true)

@info "Found equations:"
print_equations(Ψ)
# Extract the parameter
p̂ = parameters(Ψ)
println("First parameter guess : $(p̂)")

# Define the recovered, hyrid model with the rescaled dynamics
function recovered_dynamics!(du,u, p, t)
    û = Ψ(u, p[3:end]) # Network prediction
    du[1] = p[1]*u[1] + û[1]
    du[2] = -p[2]*u[2] + û[2]
end

p_model = [p_trained[1:2];p̂]
estimation_prob = ODEProblem(recovered_dynamics!, Xₙ[:, 1], tspan, p_model)
# Convert for reuse
sys = modelingtoolkitize(estimation_prob);
dudt = ODEFunction(sys);
estimation_prob = ODEProblem(dudt,Xₙ[:, 1], tspan, p_model)
estimate = solve(estimation_prob, Tsit5(), saveat = t)

##  Fit the found model
function loss_fit(θ)
    X̂ = Array(solve(estimation_prob, Tsit5(), p = θ, saveat = t))
    sum(abs2, X̂ .- Xₙ)
end

# Post-fit the model
res_fit = DiffEqFlux.sciml_train(loss_fit, p_model, BFGS(initial_stepnorm = 0.1f0), maxiters = 1000)
p_fitted = res_fit.minimizer

# Estimate
estimate_rough = solve(estimation_prob, Tsit5(), saveat = 0.1*mean(diff(t)), p = p_model)
estimate = solve(estimation_prob, Tsit5(), saveat = 0.1*mean(diff(t)), p = p_fitted)

# Plot
pl_fitted = plot(t, transpose(Xₙ), style = :dash, lw = 2,color = :black, label = ["Measurements" nothing], xlabel = "t", ylabel = "x(t), y(t)")
plot!(estimate_rough, color = :red, label = ["Recovered" nothing])
plot!(estimate, color = :blue, label = ["Recovered + Fitted" nothing])
savefig(pl_fitted,joinpath(pwd(),"plots","$(svname)recovery_fitting.pdf"))

## Simulation

# Look at long term prediction
t_long = (0.0f0, 50.0f0)
estimate_long = solve(estimation_prob, Tsit5(), saveat = 0.25f0, tspan = t_long,p = p_fitted)
plot(estimate_long.t, transpose(xscale .* estimate_long[:,:]), xlabel = "t", ylabel = "x(t),y(t)")

## Save the results
save(joinpath(pwd(),"results","Hudson_Bay_recovery.jld2"),
    "X", Xₙ, "t" , t, "neural_network" , U, "initial_parameters", p, "trained_parameters" , p_trained, # Training
    "losses", losses, "result", Ψ, "recovered_parameters", p̂, # Recovery
    "model", recovered_dynamics!, "model_parameter", p_model, "fitted_parameter", p_fitted,
    "long_estimate", estimate_long) # Estimation

## Post Processing and Plots

c1 = 3 # RGBA(174/255,192/255,201/255,1) # Maroon
c2 = :orange # RGBA(132/255,159/255,173/255,1) # Red
c3 = :blue # RGBA(255/255,90/255,0,1) # Orange
c4 = :purple # RGBA(153/255,50/255,204/255,1) # Purple

p3 = scatter(t, transpose(Xₙ), color = [c1 c2], label = ["x data" "y data"],
             title = "Recovered Model from Hudson Bay Data",
             titlefont = "Helvetica", legendfont = "Helvetica",
             markersize = 5)

plot!(p3,estimate_long, color = [c3 c4], lw=1, label = ["Estimated x(t)" "Estimated y(t)"])
plot!(p3,[19.99,20.01],[0.0,maximum(Xₙ)*1.25],lw=1,color=:black, label = nothing)
annotate!([(10.0,maximum(Xₙ)*1.25,text("Training \nData",12 , :center, :top, :black, "Helvetica"))])
savefig(p3,joinpath(pwd(),"plots","$(svname)full_plot.pdf"))


text:
using TriMeshes
#using SVG
using LinearAlgebra

debug = false

n = Net()

function profile(path::Vertex)
	[(1,2), (2,1.6), (3,1.75), (3.5, 2), (4,2.2), (5,2), (6,1), (7,0.2), (8, 0.2), (9, 0.8), (9.3, 1), (10, 1.6), (10.3, 2), (10.7, 3), (10.8, 4), (10.9,5), (10.7, 6), (10.5, 7), (10,7.2), (9,7.3), (8, 7.1), (7.8, 7), (7.3, 6), (7,5.8), (6.4, 5), (6,4.75), (5,4), (4,4.5), (3.5, 5), (3.2, 6), (3,6.3), (2.3,7), (2,7.01), (1,6), (0.7, 5), (0.5,4), (0.6, 3)], Vertex(0,0,1)
end

function arc(t::Real)
	Vertex(30 * cos(2pi*t), 0, 30 * sin(2pi*t))
end

function normal(t::Real, tstep::Real, fpath)
	t0 = fpath(t - tstep, tstep)
	t1 = fpath(t, tstep)
	t2 = fpath(t + tstep, tstep)
	v = (t1 - t0) + (t2 - t1)
	n = normalize([v.x, v.y, v.z])
	Vertex(n[1], n[2], n[3])
end

function solid(origin::Vertex, n::Net, slices::Real, fpath, slicer)
	path = Vertex(0,0,0) # to return end point
	tstep = 1/(slices-1)
	ve = 0
	for t in 0.0:tstep:1.0

		if debug println("T: $t") end

		path = origin + fpath(t, tstep)
		pathv = vertex!(n, path)
		slicepts, slicenorm = slicer(path)
		steps = length(slicepts)
		pathnorm = normal(t, tstep, fpath)

		if debug println("pathnorm $pathnorm") end
		if debug println("slicenorm $slicenorm") end

		pathay = angleXZ(pathnorm)
		sliceay = angleXZ(slicenorm)
		if sliceay < 0
			sliceay += 2pi
		end

		pathax = angleYZ(pathnorm)
		sliceax = angleYZ(slicenorm)

		rot = Vertex(0., sliceay-pathay, 0.)

		if debug println("pathax:$(round(Int, rad2deg(pathax))) sliceax:$(round(Int, rad2deg(sliceax))) pathay:$(round(Int, rad2deg(pathay))) sliceay:$(round(Int(rad2deg(sliceay))))") end
		if debug println("rotating x:$(round(Int,rad2deg(rot.x))) y:$(round(Int,rad2deg(rot.y))) z:$(round(Int,rad2deg(rot.z)))") end

		for (x,y) in slicepts
			v = rotate(Vertex(x, y, 0), rot)
			ve = vertex!(n, v + path)
		end

if t < tstep || t >= 1.0
	face!(n, pathv, pathv+steps, pathv+1)
	for v in pathv+2:pathv+steps
		face!(n, pathv, v-1, v)
	end
end

		if t > 0
			face!(n, pathv+1, pathv+steps, pathv-1)
			face!(n, pathv+1, pathv-steps, pathv-1)
			for s in 1:steps-1
				face!(n, pathv+1+s, pathv+s, pathv-1+s-steps)
				face!(n, pathv-1+s-steps, pathv+s-steps, pathv+1+s )
			end
		end
	end
	path, ve
end

function sqr(path::Vertex)
	[(-1,-1), (1,-1), (1,1), (-1,1), ], Vertex(0,0,1)
end

function circ(path::Vertex)
	pts = Vector{Tuple{Real, Real}}()
	r = 1.5
	steps = 8
	for a = 2pi/steps:2pi/steps:2pi
		push!(pts, (r * sin(a), r * cos(a)))
	end
	pts, Vertex(0,0,1)
end

function stockp(t::Float64, tstep::Float64)
	r = 10
	x = y = z = 0.0
	function q1(u)
		r+r*cos(pi - pi*u), r*sin(pi - pi*u)
	end
	function q2(u)
		x1,z1 = q1(1.)
		x1 - 0.5 * r * u, z1 -r * u
	end
	function q3(u)
		x2,z2 = q2(1.)
		x2 + r + r*cos(pi - pi*u), z2 + r*-sin(pi - pi*u)
	end
	function q4(u)
		x3, z3 = q3(1.)
		x3 - 0.5r*u, z3 + r*u
	end

	q4m = 1.00
	q3m = q4m - tstep
	q2m = 0.50
	q1m = q2m - tstep

	if t < 0
		x,z = q1((1/q1m)*abs(t))
		z = -z
	elseif t > 1
		x,z = q4((1/(q4m-q3m))*((2-t)-q3m))
		z = -z
	elseif t < q1m
		x,z = q1((1/q1m)t)
	elseif t < q2m
		x,z = q2((1/(q2m-q1m))*(t-q1m))
	elseif t < q3m
		x,z = q3((1/(q3m-q2m))*(t-q2m))
	else
		x,z = q4((1/(q4m-q3m))*(t-q3m))
	end
	z = -z

	Vertex(x,y,z)
end

#==
function pathsvg(fname, tstep, fpath)
	pts = Vector{Tuple{Real,Real}}()
	st = fpath(1., tstep)
	for stch = 1:1
		for t=-tstep:tstep:1+tstep
			v = fpath(t, tstep)
			push!(pts, (stch*st.x+v.x, stch*st.z+v.z))
		end
	end

	write(fname, Svg([Polyline(pts; style=Style(;strokewidth=3))]), 500, 500; inhtml=true)

end
==#

function sweep(n::Net, steps, multiples=1)
	e = o = Vertex(0,0,0)

	for k in 1:multiples
		o, e = solid(o, n, steps, stockp, circ)
	end
	o, e
end

function knit(n, stitches::Integer, rowpairs::Integer)
	ty = 15
	vs = 1
	for rp = 1:rowpairs
		pe, ve = sweep(n, 16, stitches)
		transformVerts!(n, v->rotate(v, Vertex(deg2rad(40), 0, 0)), vs, ve)
		transformVerts!(n, v->translate(v, Vertex(0,ty,0)), vs, ve)

		vs = ve + 1
		pe, ve = sweep(n, 16, stitches)
		transformVerts!(n, v->rotate(v, Vertex(deg2rad(-40), 0, 0)), vs, ve)
		transformVerts!(n, v->translate(v, Vertex(15,ty,0)), vs, ve)
		vs = ve + 1
		ty += 15
	end

	#STL_ASCII(n, "stk.stl")
end

function go()
	#pathsvg("stk.svg", 1/38, stockp)
	knit(n, 3, 2)
end


text:
abstract type AbstractPoint end

function iselliptic(𝑥::Number,𝑦::Number,𝑎::Number,𝑏::Number)
    𝑦^2 == 𝑥^3 + 𝑎*𝑥 + 𝑏
end

POINTTYPES = Union{Integer,PrimeField}

"""
Represents a point with coordinates (𝑥,𝑦) on an elliptic curve where 𝑦² = 𝑥³ + 𝑎𝑥 + 𝑏
Optional parameter 𝑝 represents finite field 𝐹ₚ and will convert all other parameter to PrimeField
Point(∞,∞,𝑎,𝑏) represents point at infinity
Returns an error if elliptic curve equation isn't satisfied
"""
struct Point{T<:Number,S<:Number} <: AbstractPoint
    𝑥::T
    𝑦::T
    𝑎::S
    𝑏::S
    Point{T,S}(𝑥,𝑦,𝑎,𝑏) where {T<:Number,S<:Number} = new(𝑥,𝑦,𝑎,𝑏)
end

Point(𝑥::Infinity,𝑦::Infinity,𝑎::T,𝑏::T) where {T<:POINTTYPES} = Point{Infinity,T}(𝑥,𝑦,𝑎,𝑏)
Point(𝑥::T,𝑦::T,𝑎::T,𝑏::T) where {T<:POINTTYPES} = !iselliptic(𝑥,𝑦,𝑎,𝑏) ? throw(DomainError("Point is not on curve")) : Point{T,T}(𝑥,𝑦,𝑎,𝑏)
Point(𝑥::Infinity,𝑦::Infinity,𝑎::T,𝑏::T,𝑝::T) where {T<:Integer} = Point(𝑥,𝑦,FieldElement(𝑎,𝑝),FieldElement(𝑏,𝑝))
Point(𝑥::T,𝑦::T,𝑎::T,𝑏::T,𝑝::T) where {T<:Integer} = Point(FieldElement(𝑥,𝑝),FieldElement(𝑦,𝑝),FieldElement(𝑎,𝑝),FieldElement(𝑏,𝑝))

"Formats AbstractPoint as (𝑥, 𝑦) on 𝑦² = 𝑥³ + 𝑎𝑥 + 𝑏 (: 𝐹ₚ)"
function show(io::IO, z::AbstractPoint)
    if typeof(z.𝑥) <: PrimeField
        x, y = z.𝑥.𝑛, z.𝑦.𝑛
    else
        x, y = z.𝑥, z.𝑦
    end

    if typeof(z.𝑎) <: PrimeField
        a, b = z.𝑎.𝑛, z.𝑏.𝑛
        field = string(" : 𝐹", z.𝑎.𝑝)
    else
        a, b = z.𝑎, z.𝑏
        field = ""
    end
    print(io, "(", x, ", ", y, ") on 𝑦² = 𝑥³ + ", a, "𝑥 + ", b, field)
end

"""
Returns the point resulting from the intersection of the curve and the
straight line defined by the points P and Q
"""
function +(𝑃::AbstractPoint,𝑄::AbstractPoint)
    T = typeof(𝑃)
    S = typeof(𝑃.𝑎)
    if 𝑃.𝑎 != 𝑄.𝑎 || 𝑃.𝑏 != 𝑄.𝑏
        throw(DomainError("Points are not on the same curve"))

    # Case 0
    elseif 𝑃.𝑥 == ∞
        return 𝑄
    elseif 𝑄.𝑥 == ∞
        return 𝑃
    elseif 𝑃.𝑥 == 𝑄.𝑥 && 𝑃.𝑦 != 𝑄.𝑦
        # something more elegant should exist to return correct point type
        if T <: Point
            return Point{Infinity,S}(∞, ∞, 𝑃.𝑎, 𝑃.𝑏)
        elseif T <: S256Point
            return S256Point{Infinity}(∞, ∞, 𝑃.𝑎, 𝑃.𝑏)
        end

    # Case 1
    elseif 𝑃.𝑥 != 𝑄.𝑥
        λ = (𝑄.𝑦 - 𝑃.𝑦) ÷ (𝑄.𝑥 - 𝑃.𝑥)
        𝑥 = λ^2 - 𝑃.𝑥 - 𝑄.𝑥
    # Case 2
    else
        λ = (3 * 𝑃.𝑥^2 + 𝑃.𝑎) ÷ (2 * 𝑃.𝑦)
        𝑥 = λ^2 - 2 * 𝑃.𝑥
    end
    𝑦 = λ * (𝑃.𝑥 - 𝑥) - 𝑃.𝑦
    return T(S(𝑥), S(𝑦), 𝑃.𝑎, 𝑃.𝑏)
end

"Scalar multiplication of a Point"
function *(λ::Integer,𝑃::Point)
    𝑅 = Point(∞, ∞, 𝑃.𝑎, 𝑃.𝑏)
    while λ > 0
        𝑅 += 𝑃
        λ -= 1
    end
    return 𝑅
end


text:
@testset "4_ckks_basics" begin
  @testset "EncryptionParameters" begin
    @test_nowarn EncryptionParameters(SchemeType.ckks)
  end
  enc_parms = EncryptionParameters(SchemeType.ckks)

  @testset "polynomial modulus degree" begin
    @test_nowarn set_poly_modulus_degree!(enc_parms, 8192)
    @test poly_modulus_degree(enc_parms) == 8192
  end

  @testset "coefficient modulus" begin
    @test_nowarn coeff_modulus_create(8192, [60, 40, 40, 60])
    @test_nowarn set_coeff_modulus!(enc_parms, coeff_modulus_create(8192, [60, 40, 40, 60]))
  end

  @testset "SEALContext" begin
    @test_nowarn SEALContext(enc_parms)
  end
  context = SEALContext(enc_parms)

  @testset "extract parameters" begin
    @test_nowarn key_context_data(context)
    context_data = key_context_data(context)
    @test_nowarn parms(context_data)
    ec = parms(context_data)
    @test scheme(ec) == SchemeType.ckks
    @test total_coeff_modulus_bit_count(context_data) == 200
    @test_nowarn coeff_modulus(ec)
    bit_counts = [bit_count(modulus) for modulus in coeff_modulus(ec)]
    @test bit_counts == [60, 40, 40, 60]
  end

  @testset "KeyGenerator" begin
    @test_nowarn KeyGenerator(context)
  end
  keygen = KeyGenerator(context)

  @testset "PublicKey" begin
    @test_nowarn PublicKey()
  end
  public_key_ = PublicKey()

  @testset "create_public_key" begin
    @test_nowarn create_public_key!(public_key_, keygen)
  end

  @testset "SecretKey" begin
    @test_nowarn secret_key(keygen)
  end
  secret_key_ = secret_key(keygen)

  @testset "RelinKeys" begin
    @test_nowarn RelinKeys()
  end
  relin_keys_ = RelinKeys()

  @testset "create_relin_keys" begin
    @test_nowarn create_relin_keys!(relin_keys_, keygen)
  end

  @testset "Encryptor" begin
    @test_nowarn Encryptor(context, public_key_)
  end
  encryptor = Encryptor(context, public_key_)

  @testset "Evaluator" begin
    @test_nowarn Evaluator(context)
  end
  evaluator = Evaluator(context)

  @testset "Decryptor" begin
    @test_nowarn Decryptor(context, secret_key_)
  end
  decryptor = Decryptor(context, secret_key_)

  @testset "CKKSEncoder" begin
    @test_nowarn CKKSEncoder(context)
  end
  encoder = CKKSEncoder(context)

  slot_count_ = 4096
  @testset "slot_count" begin
    @test slot_count(encoder) == slot_count_
  end

  @testset "Plaintext" begin
    @test_nowarn Plaintext()
  end
  plain_coeff3 = Plaintext()
  plain_coeff1 = Plaintext()
  plain_coeff0 = Plaintext()
  x_plain = Plaintext()

  input = collect(range(0.0, 1.0, length=slot_count_))
  initial_scale = 2.0^40
  @testset "encode!" begin
    @test_nowarn encode!(plain_coeff3, 3.14159265, initial_scale, encoder)
    @test_nowarn encode!(plain_coeff1, 0.4, initial_scale, encoder)
    @test_nowarn encode!(plain_coeff0, 1.0, initial_scale, encoder)
    @test_nowarn encode!(x_plain, input, initial_scale, encoder)
  end

  @testset "Ciphertext" begin
    @test_nowarn Ciphertext()
  end
  x1_encrypted = Ciphertext()

  @testset "encrypt!" begin
    @test_nowarn encrypt!(x1_encrypted, x_plain, encryptor)
  end

  x3_encrypted = Ciphertext()
  @testset "square!" begin
    @test_nowarn square!(x3_encrypted, x1_encrypted, evaluator)
  end

  @testset "relinearize_inplace!" begin
    @test_nowarn relinearize_inplace!(x3_encrypted, relin_keys_, evaluator)
  end

  @testset "scale (before rescaling)" begin
    @test isapprox(log2(scale(x3_encrypted)), 80)
  end

  @testset "rescale_to_next_inplace!" begin
    @test rescale_to_next_inplace!(x3_encrypted, evaluator) == x3_encrypted
  end

  @testset "scale (after rescaling)" begin
    @test isapprox(log2(scale(x3_encrypted)), 40)
  end

  x1_encrypted_coeff3 = Ciphertext()
  @testset "multiply_plain! and rescale" begin
    @test_nowarn multiply_plain!(x1_encrypted_coeff3, x1_encrypted, plain_coeff3, evaluator)
    @test_nowarn rescale_to_next_inplace!(x1_encrypted_coeff3, evaluator)
  end

  @testset "multiply_inplace!" begin
    @test multiply_inplace!(x3_encrypted, x1_encrypted_coeff3, evaluator) == x3_encrypted
  end

  @testset "relinearize_inplace! and rescale" begin
    @test_nowarn relinearize_inplace!(x3_encrypted, relin_keys_, evaluator)
    @test_nowarn rescale_to_next_inplace!(x3_encrypted, evaluator)
  end

  @testset "multiply_plain_inplace! and rescale" begin
    @test_nowarn multiply_plain_inplace!(x1_encrypted, plain_coeff1, evaluator)
    @test_nowarn rescale_to_next_inplace!(x1_encrypted, evaluator)
  end

  @testset "parms_id" begin
    @test parms_id(x3_encrypted) == UInt64[0x2af33fda5cee1476, 0xe1a78ed1ec9d76b3,
                                           0xed2adee911ba7c4d, 0x94f949f4f9055b1a]
    @test parms_id(x1_encrypted) == UInt64[0xdc264f695a81d156, 0x890465c10f20b410,
                                           0xf03d86f8bc932745, 0x2e67b09e17d4a44c]
    @test parms_id(plain_coeff0) == UInt64[0x12a009457dbf8b0f, 0x6364d5f3d4c92b3c,
                                           0x96a841ccd88e440c, 0x1255677018089458]
  end
  pid1 = parms_id(x3_encrypted)
  pid2 = parms_id(x1_encrypted)
  pid3 = parms_id(plain_coeff0)

  @testset "get_context_data" begin
    @test_nowarn get_context_data(context, pid1)
    @test_nowarn get_context_data(context, pid2)
    @test_nowarn get_context_data(context, pid3)
  end
  cd1 = get_context_data(context, pid1)
  cd2 = get_context_data(context, pid2)
  cd3 = get_context_data(context, pid3)

  @testset "chain_index" begin
    @test chain_index(cd1) == 0
    @test chain_index(cd2) == 1
    @test chain_index(cd3) == 2
  end

  @testset "scale!" begin
    @test_nowarn scale!(x3_encrypted, 2.0^40)
    @test_nowarn scale!(x1_encrypted, 2.0^40)
  end

  last_parms_id = parms_id(x3_encrypted)
  @testset "mod_switch_to_inplace!" begin
    @test mod_switch_to_inplace!(x1_encrypted, last_parms_id, evaluator) == x1_encrypted
    @test mod_switch_to_inplace!(plain_coeff0, last_parms_id, evaluator) == plain_coeff0
  end

  encrypted_result = Ciphertext()
  @testset "add!" begin
    @test_nowarn add!(encrypted_result, x3_encrypted, x1_encrypted, evaluator)
  end

  @testset "add_plain_inplace!" begin
    @test_nowarn add_plain_inplace!(encrypted_result, plain_coeff0, evaluator)
  end

  plain_result = Plaintext()
  @testset "decrypt!" begin
    @test_nowarn decrypt!(plain_result, encrypted_result, decryptor)
  end

  result = similar(input)
  @testset "decode!" begin
    @test_nowarn decode!(result, plain_result, encoder)
  end

  @testset "compare results" begin
    true_result = similar(input)
    for (i, x) in enumerate(input)
      true_result[i] = (3.14159265 * x * x + 0.4) * x + 1
    end
    @test isapprox(result, true_result, atol=1e-4)
  end
end


text:
using PlotlyJS

function linescatter1()
    trace1 = scatter(;x=1:4, y=[10, 15, 13, 17], mode="markers")
    plot([trace1])
end
linescatter1()

testio = open("test.jpg", "w")
PlotlyJS.savefig(testio, linescatter1(); format="jpeg");
close(testio)


text:
using Pkg
using Test

using MPI

# load test packages to trigger precompilation
using DoubleFloats
if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CuArrays
end

# Code coverage command line options; must correspond to src/julia.h
# and src/ui/repl.c
const JL_LOG_NONE = 0
const JL_LOG_USER = 1
const JL_LOG_ALL = 2
const coverage_opts =
    Dict{Int, String}(JL_LOG_NONE => "none",
                      JL_LOG_USER => "user",
                      JL_LOG_ALL => "all")

function runtests()
    nprocs = clamp(Sys.CPU_THREADS, 2, 4)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = dirname(@__FILE__)
    istest(f) = endswith(f, ".jl") && startswith(f, "test_")
    testfiles = sort(filter(istest, readdir(testdir)))

    nfail = 0
    printstyled("Running MPI.jl tests\n"; color=:white)

    for f in testfiles
        coverage_opt = coverage_opts[Base.JLOptions().code_coverage]
        mpiexec() do cmd
            if f == "test_spawn.jl"
                run(`$cmd -n 1 $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
            elseif f == "test_threads.jl"
                withenv("JULIA_NUM_THREAD" => "4") do
                    run(`$cmd -n $nprocs $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
                end
            else
                run(`$cmd -n $nprocs $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
            end
        end
        Base.with_output_color(:green,stdout) do io
            println(io,"\tSUCCESS: $f")
        end
    end
    return nfail
end

exit(runtests())


text:
using fci
using JLD2

@load "_testdata_h4_triplet_ymatrix.jld2"

norbs = 4
nalpha = 3
nbeta = 1

@testset "make xy for grms" begin
    n_unocc_a = (norbs-nalpha)+1
    n_unocc_b = (norbs-nbeta)+1

    #make x matricies
    xalpha = zeros(n_unocc_a, nalpha+1)
    xbeta = zeros(n_unocc_b, nbeta+1)
    #fill first row and columns
    xalpha[:,1] .= 1
    xbeta[:,1] .= 1
    xalpha[1,:] .= 1
    xbeta[1,:] .= 1

    for i in 2:nalpha+1
        for j in 2:n_unocc_a
            xalpha[j, i] = xalpha[j-1, i] + xalpha[j, i-1]
        end
    end

    for i in 2:nbeta+1
        for j in 2:n_unocc_b
            xbeta[j, i] = xbeta[j-1, i] + xbeta[j, i-1]
        end
    end

    #make y matrices
    copya = deepcopy(xalpha)
    copyb = deepcopy(xbeta)
    arraya = zeros(size(xalpha)[2])
    arrayb = zeros(size(xbeta)[2])
    yalpha_test = vcat(transpose(arraya), xalpha[1:size(xalpha)[1]-1, :])
    ybeta_test = vcat(transpose(arrayb), xbeta[1:size(xbeta)[1]-1, :])#=}}}=#
    @test isapprox(ya, yalpha_test, atol=0.05)
    @test isapprox(yb, ybeta_test, atol=0.05)
end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
################################################################################
#
# Run inversion for the Gaussian lens model
#
################################################################################

R = Float32

### Module loading
using LinearAlgebra, JLD2, PyPlot, Dates, Distributed
using Optim, LineSearches
using JUDI

### Load or model synthetic data
base_path = dirname(pathof(JUDI))*"/../examples/twri/"
try
    @load string(base_path*"data/GaussLens_data_acou.jld") model_true dat fsrc
catch e
    @info "Data not found, modeling true data"
    include(base_path*"data/gen_data_gausslens_acou.jl")
    @load string(base_path*"data/GaussLens_data_acou.jld") model_true dat fsrc
end

### Background model
m0 = 1f0./2f0.^2*ones(R, model_true.n)
n = model_true.n
d = model_true.d
o = model_true.o
model0 = Model(n, d, o, m0)

# Pre- and post-conditioning
mask = zeros(Float32, model0.n)
mask[6:end-5, 6:end-5] .= 1
function proj_bounds(m, mmin, mmax)
    m[m[:] .< mmin] .= mmin
    m[m[:] .> mmax] .= mmax
    return m
end
mmin = 1f0/4f0^2
mmax = 1f0/1f0^2
preproc(x) = proj_bounds(x, mmin, mmax)

# [FWI]
function funFWI!(F, G, x)
    i = rand(1:15, 4)
    model0.m .= preproc(x)
    fun, g = fwi_objective(model0, fsrc[i], dat[i])
    F = fun
    isnothing(G) && return fun
    isnothing(g_const) && (global g_const = .1/maximum(abs.(g)))
    G .= g_const .* g[:] .* vec(mask)
    return fun
end

# # [TWRIdual]
ε = zeros(Float32, fsrc.nsrc)
v_bg = sqrt(1/m0[1])
freq_peak = 0.006f0
δ = R(sqrt(2)/2)*v_bg/freq_peak
weight_fun_pars = ("srcfocus", δ)
optwri = TWRIOptions(;grad_corr=true, comp_alpha=true, weight_fun=weight_fun_pars, eps=ε, params=:m)

function funWRI(F, G, x, optwri)
    i = rand(1:15, 4)
    model0.m .= preproc(x)
    fun, g = twri_objective(model0, fsrc[i], dat[i], nothing; optionswri=optwri)
    F = fun
    isnothing(G) && return fun
    isnothing(g_const) && (global g_const = .1/maximum(abs.(g)))
    G .= g_const .* g[:].* vec(mask)
    return fun
end
funWRI!(F, G, x) = funWRI(F, G, x, optwri)

### Optimization
linesearch = HagerZhang(display=1, linesearchmax=10)
# Options
method = LBFGS(;linesearch=linesearch)
niter = 20
optimopt = Optim.Options(g_calls_limit=200, iterations = niter, store_trace = true, show_trace = true, show_every = 1)

# ## Run FWI
g_const = nothing
x0_fwi = vec(m0)
result = optimize(Optim.only_fg!(funFWI!), x0_fwi, method, optimopt)
x0_fwi .= Optim.minimizer(result)
m_inv_fwi = preproc(x0_fwi)
loss_log_fwi = Optim.f_trace(result)
@save "./GaussLens_result_fwi.jld" m_inv_fwi loss_log_fwi

## Run WRI
g_const = nothing
x0_wri = vec(m0)
result = optimize(Optim.only_fg!(funWRI!), x0_wri, method, optimopt)
x0_wri .= Optim.minimizer(result)
m_inv_wri = preproc(x0_wri)
loss_log_wri = Optim.f_trace(result)
@save "./GaussLens_result_wri.jld" m_inv_wri loss_log_wri

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# ---
# title: 999. Available Captures for Rook
# id: problem999
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/available-captures-for-rook/description/>
# hidden: true
# ---
#
# On an 8 x 8 chessboard, there is one white rook.  There also may be empty
# squares, white bishops, and black pawns.  These are given as characters 'R',
# '.', 'B', and 'p' respectively. Uppercase characters represent white pieces,
# and lowercase characters represent black pieces.
#
# The rook moves as in the rules of Chess: it chooses one of four cardinal
# directions (north, east, west, and south), then moves in that direction until
# it chooses to stop, reaches the edge of the board, or captures an opposite
# colored pawn by moving to the same square it occupies.  Also, rooks cannot
# move into the same square as other friendly bishops.
#
# Return the number of pawns the rook can capture in one move.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/02/20/1253_example_1_improved.PNG)
#
#
#
#     Input: [[".",".",".",".",".",".",".","."],[".",".",".","p",".",".",".","."],[".",".",".","R",".",".",".","p"],[".",".",".",".",".",".",".","."],[".",".",".",".",".",".",".","."],[".",".",".","p",".",".",".","."],[".",".",".",".",".",".",".","."],[".",".",".",".",".",".",".","."]]
#     Output: 3
#     Explanation:
#     In this example the rook is able to capture all the pawns.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/02/19/1253_example_2_improved.PNG)
#
#
#
#     Input: [[".",".",".",".",".",".",".","."],[".","p","p","p","p","p",".","."],[".","p","p","B","p","p",".","."],[".","p","B","R","B","p",".","."],[".","p","p","B","p","p",".","."],[".","p","p","p","p","p",".","."],[".",".",".",".",".",".",".","."],[".",".",".",".",".",".",".","."]]
#     Output: 0
#     Explanation:
#     Bishops are blocking the rook to capture any pawn.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2019/02/20/1253_example_3_improved.PNG)
#
#
#
#     Input: [[".",".",".",".",".",".",".","."],[".",".",".","p",".",".",".","."],[".",".",".","p",".",".",".","."],["p","p",".","R",".","p","B","."],[".",".",".",".",".",".",".","."],[".",".",".","B",".",".",".","."],[".",".",".","p",".",".",".","."],[".",".",".",".",".",".",".","."]]
#     Output: 3
#     Explanation:
#     The rook can capture the pawns at positions b5, d6 and f5.
#
#
#
#
# **Note:**
#
#   1. `board.length == board[i].length == 8`
#   2. `board[i][j]` is either `'R'`, `'.'`, `'B'`, or `'p'`
#   3. There is exactly one cell with `board[i][j] == 'R'`
#
#
## @lc code=start
using LeetCode

function num_rook_captures(board::Matrix{Char})::Int
    res = 0
    i, j = findfirst(==('R'), board).I

    for v in (
        @view(board[i, (j + 1):end]),       ## right
        @view(board[i, (j - 1):-1:1]),      ## left
        @view(board[(i - 1):-1:1, j]),      ## up
        @view(board[(i + 1):end, j]),       ## down
    )
        x_ind = findfirst(!=('.'), v)
        if !isnothing(x_ind) && v[x_ind] == 'p'
            res += 1
        end
    end
    return res
end
## @lc code=end


text:
account = getaccount!(client)

@testset "Get User Information" begin
    @test account.email == "sammy@digitalocean.com"
    @test account.dropletlimit == 25
    @test account.floatingiplimit == 5
    @test account.uuid == "b6fr89dbf6d9156cace5f3c78dc9851d957381ef"
    @test account.emailverified
    @test account.status == "active"
    @test account.statusmessage == ""
end;


text:
function calc_phi(x::Array{Float64,2},xi::Array{Float64,2},n::Int64,nn::Int64,c::Float64,vars::Int64)
  phi=zeros(nn,n)
  for i=1:nn
      for j=1:n
          # 2 r^4
          dr=0
          for k=1:vars
            dr = dr + (xi[i,k] - x[j,k])^2
          end
          # r = (dr)
          # phi[i,j] = -(dr^4)/4/c

          # 1 Gaussian
          # phi[i,j] = exp(-sum((xi[i,:] - x[j,:]).^2)/c)
          phi[i,j] = exp(-dr/c)
      end
      # if i in 1:100:nn println(Dates.format(now(), "HH:MM:SS "),i," of ",nn) end
  end

  # xi_all=Array{Float64,3}(undef, nn,1,vars)
  # xi_all[:,1,:]=xi
  # # println(size(xi_all))
  # xi_all=repeat(xi_all,1,n,1)
  # # println(size(xi_all))

  # x_all=Array{Float64,3}(undef, n,1,vars)
  # x_all[:,1,:]=x
  # # println(size(x_all))
  # x_all=repeat(x_all,1,nn,1)
  # # println(size(x_all))
  # x_all = permutedims(x_all, [2, 1, 3])
  # # println(size(x_all))

  # # @time m1=sum((xi_all - x_all).^2,dims=3)
  # # # println(size(m1))
  # phi2=exp.(-sum((xi_all - x_all).^2,dims=3)/c)[:,:,1]
  # # println(size(phi2))
  # # println(mean(abs.(phi-phi2)))

  return phi
end

# phi[i,j] = exp(-sqrt(dr)/c)
# phi[i,j] = exp(-r^2)/2 + (r*pi^(1/2)*erf(r))/2;
# phi[i,j] = 2*(1./(1+exp.(-r)))-1;
# phi[i,j] = exp(-r^2/c[i,j])
# phi[i,j] = exp(-dr/c)
# phi[i,j]=tanh.(r^2/c[i,j])
# phi[i,j] = -(r^3)/3
# phi[i,j] = exp(-r/c)
# phi[i,j] =(c^2*exp(-r/c^2))/2 + (c*sqrt(pi)*erf(sqrt(r)/c)*sqrt(r))/2;
# phi[i,j] = -sum((xi[i,:] - x[j,:]).^2)