text:
"""
    DR_TreeNode

Tree node stores information of stage problem

    - `id`: index of node
    - `stage_builder`: adds variables and constraints to model with input (tree::Tree, subtree::SubTree, node::SubTreeNode)
    - `parent`: index of parent node
    - `children`: indices of child nodes
    - `stage`: current stage
    - `ξ`: current scenario
    - `set`: ambiguity set
    - 'cost': current scenario of cost
"""

mutable struct DR_TreeNode <: DD.AbstractTreeNode
    id::Int
    stage_builder::Union{Nothing,Function}
    coupling_variables::Dict{Any,Vector{DD.CouplingVariableRef}}
    parent::Int
    children::Vector{Int}
    stage::Int
    ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}

    set::Union{AbstractAmbiguitySet, Nothing}
    cost::Dict{String, Float64}

    function DR_TreeNode(id::Int, parent::Int, stage::Int, ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}},
            set::Union{AbstractAmbiguitySet, Nothing})
        tn = new()
        tn.id = id
        tn.stage_builder = nothing
        tn.coupling_variables = Dict{Any,Vector{DD.CouplingVariableRef}}()
        tn.parent = parent
        tn.children = Vector{Tuple{Int, Float64}}()
        tn.stage = stage
        tn.ξ = ξ
        tn.set = set
        tn.cost = Dict{Any, Float64}()
        return tn
    end
end

function DR_TreeNode(ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}, set::AbstractAmbiguitySet)
    return DR_TreeNode(1, 0, 1, ξ, set)
end

get_set(node::DR_TreeNode) = node.set
get_cost(node::DR_TreeNode, var_id::Any) = node.cost[var_id]
function set_cost!(node::DR_TreeNode, var_id::String, coeff::Float64)
    node.cost[var_id] = coeff
end

function DD.set_stage_builder!(node::DR_TreeNode, func::Function)
    node.stage_builder = func
    dummynode = DD.SubTreeNode(node, 1.0)

    func(JuMP.Model(), dummynode)

    label = DD.get_id(node)
    if dummynode.obj isa JuMP.AffExpr
        terms = dummynode.obj.terms

        for (ref, var) in dummynode.out
            name = "n$(label)_" * ref
            if haskey(terms, var)
                set_cost!(node, name, terms[var])
            else
                set_cost!(node, name, 0.0)
            end
        end

        for (ref, var) in dummynode.control
            name = "n$(label)_" * ref
            if haskey(terms, var)
                set_cost!(node, name, terms[var])
            end
        end
    else
        for (ref, var) in dummynode.out
            name = "n$(label)_" * ref
            set_cost!(node, name, 0.0)
        end
    end
end

function DD.Tree(ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}, set::AbstractAmbiguitySet)
    return DD.Tree(Dict{Int,DR_TreeNode}(1 => DR_TreeNode(ξ, set)))
end

function DD.add_child!(tree::DD.Tree{DR_TreeNode}, pt::Int, ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}, set::Union{AbstractAmbiguitySet, Nothing})::Int
    #   adds child node to tree.nodes[pt]
    @assert haskey(tree.nodes, pt)                                  # check if pt is valid
    stage = DD.get_stage(tree, pt) + 1                                 # get new stage value
    id = length(tree.nodes) + 1                                     # get node id
    DD.add_node!(tree, DR_TreeNode(id, pt, stage, ξ, set ))            # create node and add to tree
    push!(DD.get_children(tree, pt), id)                               # push id to parent node children
    return id
end

get_set(tree::DD.Tree{DR_TreeNode}, id) = get_set(tree.nodes[id])

"""
    create_Wasserstein_deterministic!

creates deterministic model

# Arguments
    - `tree`: Tree{DR_TreeNode}
"""

function create_Wasserstein_deterministic!(tree::DD.Tree{DR_TreeNode})
    subtree = DD.SubTree(0)
    # add nodes to subtree
    NDict = Dict{Int, Int}()
    for (id, node) in tree.nodes
        subnode = DD.SubTreeNode(node, 0.0)
        DD.add_node!(subtree, subnode)
        if !DD.check_leaf(node)
            NDict[id] = node.set.N
        end
    end

    @variable(subtree.model, lα__[id = keys(NDict)] >= 0)
    @variable(subtree.model, lβ__[id = keys(NDict), s=1:NDict[id]])

    obj = 0
    for (id, subnode) in subtree.nodes
        subnode.treenode.stage_builder(subtree.model, subnode)
        DD.unregister_all!(subtree.model)

        node = subnode.treenode
        if DD.check_root(subnode)
            next_set = get_set(node)
            obj += subnode.obj + next_set.ϵ * lα__[id] + sum(next_set.samples[s].p * lβ__[id,s] for s = 1:next_set.N )
        elseif DD.check_leaf(subnode)
            this_set = get_set(tree, DD.get_parent(node))
            @constraint(subtree.model, [s=1:this_set.N], this_set.norm_func(this_set.samples[s].ξ, node.ξ) * lα__[DD.get_parent(node)] + lβ__[DD.get_parent(node), s] >= subnode.obj)
        else
            this_set = get_set(tree, DD.get_parent(node))
            next_set = get_set(node)
            @constraint(subtree.model, [s=1:this_set.N], this_set.norm_func(this_set.samples[s].ξ, node.ξ) * lα__[DD.get_parent(node)] + lβ__[DD.get_parent(node), s] >=
                                                            subnode.obj + next_set.ϵ * lα__[id] + sum(next_set.samples[s].p * lβ__[id,s] for s = 1:next_set.N) )
        end
    end
    JuMP.set_objective(subtree.model, MOI.MIN_SENSE, obj)

    #
    for (id, subnode) in subtree.nodes
        parent = DD.get_parent(subnode)
        if parent!=0 && haskey(subtree.nodes, parent)
            DD.add_links!(subtree, id, parent)
        end
    end

    subtree.nodelabels = sort!(collect(keys(subtree.nodes)))
    return subtree
end

function DD.create_subtree!(block_id::Int, coupling_variables::Vector{DD.CouplingVariableRef}, nodes::Vector{DR_TreeNode})::DD.SubTree
    subtree = DD.SubTree(block_id)
    # add nodes to subtree
    for node in nodes
        subnode = DD.SubTreeNode(node, 1.0) # dummy weight
        DD.add_node!(subtree, subnode)
    end

    for (id, subnode) in subtree.nodes
        subnode.treenode.stage_builder(subtree.model, subnode)
        DD.unregister_all!(subtree.model)
    end
    JuMP.set_objective_sense(subtree.model, MOI.MIN_SENSE)

    #
    for (id, subnode) in subtree.nodes
        DD.couple_common_variables!(coupling_variables, block_id, subnode)
        parent = DD.get_parent(subnode)
        if parent!=0 && haskey(subtree.nodes, parent)
            DD.add_links!(subtree, id, parent)
        elseif parent!=0 # assuming 1st stage node is 1
            subtree.parent = parent
            subtree.root = id
            DD.couple_incoming_variables!(coupling_variables, block_id, subnode)
        end
        couple_variables_with_cost!(coupling_variables, block_id, subnode)
    end

    subtree.nodelabels = sort!(collect(keys(subtree.nodes)))
    return subtree
end

function couple_variables_with_cost!(coupling_variables::Vector{DD.CouplingVariableRef}, block_id::Int, subnode::DD.SubTreeNode)
    label = DD.get_id(subnode)
    for (ref, var) in subnode.control
        if subnode.obj isa JuMP.AffExpr
            terms = subnode.obj.terms
            if haskey(terms, var)
                DD.couple_variables!(coupling_variables, block_id, label, ref, var)
            end
        end
    end
end

"""
    decomposition_not

outputs the entire tree

# Arguments
    - `tree`: Tree
"""

function DD.decomposition_not(tree::DD.Tree{DR_TreeNode}):: Vector{Vector{DR_TreeNode}}
    nodes = Vector{DR_TreeNode}()
    for (id, node) in tree.nodes
        push!(nodes, node)
    end
    return [nodes]
end

"""
    decomposition_scenario

outputs the scenario decomposition at each leaf nodes

# Arguments
    - `tree`: Tree
"""

function DD.decomposition_scenario(tree::DD.Tree{DR_TreeNode}):: Vector{Vector{DR_TreeNode}}
    node_cluster = Vector{Vector{DR_TreeNode}}()
    for (id, node) in tree.nodes
        if DD.check_leaf(node)
            nodes = Vector{DR_TreeNode}()
            current = node
            while true
                push!(nodes, current)
                pt = DD.get_parent(current)
                current = tree.nodes[pt]
                if DD.check_root(current)
                    push!(nodes, current)
                    break
                end
            end
            push!(node_cluster, nodes)
        end
    end
    return node_cluster
end

"""
    decomposition_temporal

outputs the temporal decomposition at each nodes

# Arguments
    - `tree`: Tree
"""
#need fix

function DD.decomposition_temporal(tree::DD.Tree{DR_TreeNode}):: Vector{Vector{DR_TreeNode}}
    node_cluster = Vector{Vector{TreeNode}}()
    for (id, node) in tree.nodes
        push!(node_cluster,[node])
    end
    return node_cluster
end


text:
module FlameGraphs

using Profile, LeftChildRightSiblingTrees
using Base.StackTraces: StackFrame
using Profile: StackFrameTree
using Colors, FixedPointNumbers, IndirectArrays
using FileIO

# AbstractTree interface for StackFrameTree:
using AbstractTrees
AbstractTrees.children(node::StackFrameTree) = node.down
AbstractTrees.printnode(io::IO, node::StackFrameTree) = print(io, node.frame)

export flamegraph, flamepixels, flametags, FlameColors, StackFrameCategory

include("graph.jl")
include("render.jl")
include("sfcategory.jl")
include("io.jl")

end # module

ext: jl


text:
#abstrcattype.jl

# Metida.jl

abstract type MetidaModel <: StatisticalModel end

abstract type AbstractCovarianceStructure end

abstract type AbstractCovmatMethod end

abstract type AbstractCovarianceType end

abstract type AbstractLMMDataBlocks end

# Data

abstract type AbstractData end

abstract type AbstractDataSet{AbstractData} end

# All have field: id::Dict
abstract type AbstractIdData <: AbstractData end

# All have field:  result::Dict
abstract type AbstractResultData <: AbstractData end

function Base.getindex(a::T, s::Symbol) where T <: AbstractResultData
    return a.result[s]
end

# All have field: result::Dict, data::AbstractIdData
abstract type  AbstractIDResult{AbstractIdData} <: AbstractResultData end

# MetidaFreq.jl

# MetidaNCA.jl
# All have fields: time::Vector, obs::Vector, id::Dict
abstract type AbstractSubject <: AbstractIdData end

# All have field: result::Dict, data::AbstractSubject
#abstract type AbstractSubjectResult{AbstractSubject} <: AbstractResultData end
abstract type AbstractSubjectResult{AbstractSubject} <: AbstractIDResult{AbstractSubject} end

# Descriptive


text:
using Muon, Test

@testset "Index" begin
    include("index.jl")
end
@testset "HDF5 backed matrix" begin
    include("backed_matrix.jl")
end
@testset "aligned mappings" begin
    include("alignedmapping.jl")
end
@testset "AnnData" begin
    include("anndata.jl")
end
@testset "MuData" begin
    include("mudata.jl")
end


text:
using StanSample
using DataFrames, Tables
using Random, Distributions, Test

begin
    N = 100
    df = DataFrame(
      :h0 => rand(Normal(10,2 ), N),
      :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
    );
    df[!, :fungus] =
        [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
    df[!, :h1] =
        [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]
    data = Dict(
        :N => nrow(df),
        :h0 => df[:, :h0],
        :h1 => df[:, :h1],
        :fungus => df[:, :fungus],
        :treatment => df[:, :treatment]
    );
end;

stan6_7 = "
data {
  int <lower=1> N;
  vector[N] h0;
  vector[N] h1;
  vector[N] treatment;
  vector[N] fungus;
}
parameters{
  real a;
  real bt;
  real bf;
  real<lower=0> sigma;
}
model {
  vector[N] mu;
  vector[N] p;
  a ~ lognormal(0, 0.2);
  bt ~ normal(0, 0.5);
  bf ~ normal(0, 0.5);
  sigma ~ exponential(1);
  for ( i in 1:N ) {
    p[i] = a + bt*treatment[i] + bf*fungus[i];
    mu[i] = h0[i] * p[i];
  }
  h1 ~ normal(mu, sigma);
}
";

# ╔═╡ 655d0bcb-a4ab-41b4-8525-3e9f066113fd
begin
    m6_7s = SampleModel("m6.7s", stan6_7)
    rc6_7s = stan_sample(m6_7s; data)
end;

if success(rc6_7s)
    nt6_7s = read_samples(m6_7s)
    df6_7s = read_samples(m6_7s, :dataframe)
    a6_7s, cnames = read_samples(m6_7s, :array; return_parameters=true);
end

st = StanSample.convert_a3d(a6_7s, cnames, Val(:table))

# Testing

@test Tables.istable(st) == true

rows = Tables.rows(st)
for row in rows
    rowvals = [Tables.getcolumn(row, col) for col in Tables.columnnames(st)]
end

@test length(Tables.getcolumn(rows, :a)) == 4000
cols = Tables.columns(st)

@test Tables.schema(rows) == Tables.schema(cols)

@test size(DataFrame(st)) == (4000, 4)

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
#=

"""
This module defines the data types related to seismic data:
* `TD` : time domain representation

It also provides methods that apply source and receiver filters onto
seismic data.
"""

using Misfits
using LinearAlgebra
using Statistics
import GeoPhyInv.Interpolation
import GeoPhyInv: AGeom
import GeoPhyInv.Utils
import GeoPhyInv.Coupling
using DSP
using Random

"""
Time domain representation of Seismic Data.

# Fields

* `d::Array{Array{Float64,2},2}` : data
* `fields::Vector{Symbol}` : components recorded at each receiver
* `tgrid` : grid to represent time
* `ageom::AGeom` : acquisition ageometry used to generate the data
"""
mutable struct TD
	tgrid::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
	d::NamedArrays.NamedArray{Array{Float64,2},1,Array{Array{Float64,2},1},Tuple{OrderedCollections.OrderedDict{Symbol,Int64}}}
	nr::Int64
	"adding conditions that are to be false while construction"
	#TD(d, fields, tgrid, ageom) =
#		any([
#		  any([fields[ifield] ∉ [:P, :Vx, :Vz] for ifield in 1:length(fields)]),
#		  length(fields) == 0,
#		  broadcast(size,d) != [(length(tgrid),ageom[iss].nr) for iss=1:length(ageom), ifield=1:length(fields)]
#		  ]) ?
#		error("error in TD construction") : new(d, fields, tgrid, ageom)
end

"""
Initialize TD with zeros
"""
function Base.zero(::Type{TD}, tgrid, r::Recs=Recs(1), fields=[:P])
	nt=length(tgrid); nf=length(fields)
	return TD(tgrid,NamedArray([zeros(nt,r.n) for i in fields], (fields,)), r.n)
end
TD(tgrid,  r::Recs=Recs(1), fields=[:P])=zero(TD, tgrid, s, fields)

Data=Array{NamedD,1}

"""
Method to resample data in time.

# Arguments

* `data::TD` : input data of type `TD`
* `tgrid` : resampling in time according to this time grid

# Return

* data after resampling as `TD`
"""
function interp(data::TD,
		tgrid::StepRangeLen
		)
	nss = length(data.ageom)
	nr = data.ageom.nr
	dataout = TD(
		     [zeros(length(tgrid),data.ageom[iss].nr) for iss=1:nss, ifield=1:length(data.fields)],
	      data.fields,tgrid,data.ageom)
	interp_spray!(data, dataout)
	return dataout
end

function taper!(data::TD, perc=0.0; bperc=perc,eperc=perc)
	nr = data.ageom.nr;	nss = length(data.ageom); nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		Utils.taper!(dd,bperc=bperc,eperc=eperc)
	end
	return data
end

"""
Method to resample data in time.
Can reduce allocations =========

# Arguments

* `data::TD` : input data of type `TD`
* `dataout::TD` : preallocated data of type `TD` that is modified
"""
function interp_spray!(data::TD, dataout::TD, attrib=:interp, Battrib=:B1; pa=nothing)
	# check if datasets are similar
	nss = length(data.ageom)
	nr = data.ageom.nr
	xin=data.tgrid
	xout=dataout.tgrid
	if(pa===nothing)
		pa=Interpolation.Kernel([xin], [xout], :B1)
	end
	for ifield = 1:length(data.fields), iss = 1:nss
		dat=data.d[iss,ifield]
		dato=dataout.d[iss,ifield]
		for ir = 1:nr[iss]
			din=view(dat,:,ir)
			dout=view(dato,:,ir)
	 		Interpolation.interp_spray!(din, dout, pa, attrib)
		end
	end
	return dataout
end

"""
Method used to preallocate `TD` with zeros.

# Arguments

* `fields::Vector{Symbol}` : number of components
* `tgrid` : time domain grid
* `ageom::AGeom` : acquisition ageometry

# Return

* data with zeros as `TD`
"""
function TD_zeros(fields::Vector{Symbol}, tgrid::StepRangeLen, ageom::AGeom)
	return TD([zeros(length(tgrid),ageom[iss].nr) for iss=1:length(ageom), ifield=1:length(fields)],fields,
	   deepcopy(tgrid),deepcopy(ageom))
end
function Base.fill!(data::TD, k::Float64)
	for iss=1:length(data.ageom), ifield=1:length(data.fields)
		fill!(data.d[iss,ifield],k)
	end
end
function TD_zeros(d::TD)
	return TD([zeros(length(d.tgrid),d.ageom[iss].nr) for iss=1:length(d.ageom), ifield=1:length(d.fields)],d.fields,
	   deepcopy(d.tgrid),deepcopy(d.ageom))
end

"Same as `TD_zeros`, except for returning ones"
function TD_ones(fields::Vector{Symbol}, tgrid::StepRangeLen, ageom::AGeom)
	return TD([ones(length(tgrid),ageom[iss].nr) for iss=1:length(ageom), ifield=1:length(fields)],
	   fields,deepcopy(tgrid),deepcopy(ageom))
end

# include rest of the files
for file in ["base", "statistics", "processing", "misfits", "weights"]
	fn=joinpath(@__DIR__, string(file,".jl"))
	include(fn)
	#Revise.track(@__MODULE__,fn)
end

#=
function DDecon(d::TD, wav::AbstractVector{Float64}, ϵ=1e-2)

	dout=TD_zeros(d)
	ntd=length(dout.tgrid)
	wavv=deepcopy(wav);

	paD=Decon.ParamD(ntd=ntd,nts=length(wav), s=wavv)

	paD.ϵ=ϵ

	DDecon!(dout, d, paD)
	return dout
end

function DDecon!(dataout::TD, data::TD, paD)
	nr = data.ageom.nr;	nss = data.ageom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		ddo=dataout.d[iss, ifield]
		for ir = 1:nr[iss]
			for it in 1:nt
				paD.d[it]=dd[it,ir]
			end
			Decon.mod!(paD)
			for it in 1:nt
				ddo[it,ir]=paD.g[it]
			end
		end
	end
	return dataout
end

=#

end # module
	=#


text:
function testinc1()
  sim = LTspiceSimulation("testInc1.asc",executablepath="",tempdir=true)
  @test ~LTspice.does_circuitfilearray_file_match(sim) # because of includes in temp directory
  show(IOBuffer(),sim)
  @test measurementnames(sim) == ("incA","incB","incC","incD","incE","incF")
end
testinc1()


text:
#makes a dataframe with one row per event
#useful for tables 9 and 10 and figure 3
function formcumex(dist::DataFrame;
  refreshcumex::Bool= true,
  cumexpath::String = EVENT_PATH,
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  binextension::String = BIN_EXTENSION,
  cumexname::String = CUMEX_NAME)

  local Fret::Symbol = DIST_RET
  local cumex::DataFrame

  cumexname = "$cumexname-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  #get the appropriate table
  if refreshcumex
    cumex = formcumex(dist,Fret) #generates the table
    outbinstream("$cumexpath\\$cumexname.$binextension", cumex)

  cumex[cumex.exyear.==2015,:] |> CSV.write(
    "output\\cumex-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])-2015.csv")
  println(describe(cumex))
  else #acquire the cached version
    cumex = inbinstream("$cumexpath\\$cumexname.$binextension")
  end

  if occursin("short", string(DIST_TYPE[]))
    cumex |> CSV.write(
      "output\\cumex-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[]).csv")
  end

  return cumex
end

#make the actual table
function formcumex(dist::DataFrame, Fret::Symbol)::DataFrame
  local cumex::DataFrame
  local comp::DataFrame
  local cum::DataFrame
  local ex::DataFrame
  local cumexcols::Vector{Symbol}

  dist.abs = (abs).(dist[!,Fret])
  sdist::SubDataFrame = view(dist,
    dist.match .* #=dist.oneex .* =#dist.dclrexok .* dist.inrange .* dist.primaryyearrange,:)
  #sdist::SubDataFrame = view(dist, dist.match .* dist.dclrexok,:)

  #build the dataframe of yields

  #use this to aggregate over the days where appropriate
  @inline function aggfuncret(v::AbstractVector{<:MFloat64})
    (sum((ismissing).(v)) == length(v)) && return missing

    return prod(skipmissing(1.0 .+ v)) - 1.0
  end

  @inline function aggfuncabs(v::AbstractVector{<:MFloat64})
    (sum((ismissing).(v)) == length(v)) && return missing

    return sum(skipmissing(v))
  end

  #get the cumex returns over the respective windows
  #cumexcols = [:eid, :abs, Fret]
  #cum = aggregateoverdays(aggfunc, sdist, -7:-3, selectcols=cumexcols,
  #  colnames=[:eid, :cumabs, :cumret])
  #ex = aggregateoverdays(aggfunc, sdist, 3:7, selectcols=cumexcols,
  #  colnames=[:eid, :exabs, :exret])
  #cumex = join(cum, ex, on=:eid)

  cumexcols = [:eid, :abs]
  cum = aggregateoverdays(aggfuncabs, sdist, -7:-3, selectcols=cumexcols,
    colnames=[:eid, :cumabs])
  ex = aggregateoverdays(aggfuncabs, sdist, 3:7, selectcols=cumexcols,
    colnames=[:eid, :exabs])
  cumexabs = join(cum, ex, on=:eid)

  cumexcols = [:eid, Fret]
  cum = aggregateoverdays(aggfuncret, sdist, -7:-3, selectcols=cumexcols,
    colnames=[:eid, :cumret])
  ex = aggregateoverdays(aggfuncret, sdist, 3:7, selectcols=cumexcols,
    colnames=[:eid, :exret])
  cumexret = join(cum, ex, on=:eid)
  cumex = join(cumexabs, cumexret, on=:eid)

  #merge in the event-specific fields
  aggnames::Vector{Symbol} = [:permno, :eid, :yield, :exdate, :exyear,
    :primary, :distcd, :days2lastdiv, :ordinary]
  ssdist::SubDataFrame = view(sdist, (!ismissing).(sdist.yield), aggnames)
  agg::DataFrame = aggregate(groupby(ssdist, :eid), v->v[1])

  #fix the names
  @inline nofunction(s::Symbol) = Symbol(replace(string(s), "_function"=>""))
  rename!(agg, (s::Symbol-> s=>nofunction(s)).(names(agg)))
  #println(describe(DataFrame(agg)))
  cumex = join(cumex, agg, on=:eid)

  ssdist = view(sdist, ((r::DataFrameRow)->
    (r.day==0) && (!ismissing(r[Fret]))).(eachrow(sdist)), [:eid, Fret])
  cumex = join(cumex, ssdist, on=:eid)
  rename!(cumex, Fret=>:exdayret)

  #make the remaining time series that we need
  cumex.Dret = cumex.exret .- cumex.cumret
  cumex.Dabs = cumex.exabs .- cumex.cumabs
  cumex.yieldd1d = cumex.yield ./ (1.0 .- cumex.yield)
  finiteormissing!(cumex, :yieldd1d)

  @assert sum(cumex.yieldd1d .> .25 + 10^-8) == 0 #should be true due to prior winsorization

  #use this to create the unit-standardized fields
  function unitize(v::AbstractVector{T}) where T<:Union{Real,Missing}
    local low::Float64
    local high::Float64

    if Missing<:T
      low = minimum(skipmissing(v))
      high = maximum(skipmissing(v))
    else
      low = minimum(v)
      high = maximum(v)
    end
    return (v .- low) ./ (high - low)
  end

  cumex.UDret = unitize(cumex.Dret)
  cumex.UDabs = unitize(cumex.Dabs)
  cumex.Uyield = unitize(cumex.yield)
  cumex.Uyieldd1d = unitize(cumex.yieldd1d)
  cumex.Uexdayret = unitize(cumex.exdayret)

  requiredfields::Vector{Symbol} = [:yield, :Dret, :Dabs]
  cumex = cumex[completecases(view(cumex,:,requiredfields)), :]

  ##now we need the compustat data
  comp = prepcomp(refreshcomp=false)
  cumex = mergecrspcomp!(cumex, comp, joinkind=:left, Fcrspdate = :exdate, allowdatecollisions=true)
  excludedcompfields::Vector{Symbol} = setdiff(names(comp), [:gvkey, :lpermno, :bkliab, :at])
  select!(cumex, Not(excludedcompfields))
  #println(describe(cumex))

  #***now do the various table 10 classifications
  #start with housekeeping
  N::Int = size(cumex,1)
  cumex.Gyield = Vector{MSymbol}(undef, N)
  cumex.Gtaxdiv = Vector{MSymbol}(undef, N)
  cumex.Gfreqdiv = Vector{MSymbol}(undef, N)
  cumex.Gspecialdiv = Vector{MSymbol}(undef, N)
  cumex.Glastdiv = Vector{MSymbol}(undef, N)
  cumex.Gcompustat = Vector{MSymbol}(undef, N)
  cumex.Gbkliab = Vector{MSymbol}(undef, N)
  cumex.Gassets = Vector{MSymbol}(undef, N)

  #The following code handles dividend classifications

  #start with calssifciations based on distributions
  distcds::Vector{Vector{Int}} = (cd::Int->digits(cd)[4:-1:1]).(cumex.distcd) #(Reverse the order)
  #println(distcds)
  cumex.distcd1 = Vector{MInt}(missing, N)
  cumex.distcd2 = Vector{MInt}(missing, N)
  cumex.distcd3 = Vector{MInt}(missing, N)
  cumex.distcd4 = Vector{MInt}(missing, N)
  @mpar for i::Int ∈ 1:N
    (length(distcds[i])≠4) && continue
    cumex[i, :distcd1] = distcds[i][1]
    cumex[i, :distcd2] = distcds[i][2]
    cumex[i, :distcd3] = distcds[i][3]
    cumex[i, :distcd4] = distcds[i][4]
  end

  ###Do the classifications by row
  @mpar for r::DataFrameRow ∈ eachrow(cumex)
    r.Gyield = Gyield(r)
    r.Gtaxdiv = Gtaxdiv(r)
    r.Gfreqdiv = Gfreqdiv(r)
    r.Gspecialdiv = Gspecialdiv(r)
    r.Glastdiv = Glastdiv(r)
    r.Gcompustat = Gcompustat(r)
    r.Gbkliab = Gbkliab(r)
    r.Gassets = Gassets(r)
  end

  #drop some extraneous fields
  select!(cumex, Not([:distcd1, :distcd2, :distcd3, :distcd4]))

  return cumex
end

#checks asset thresholds
@inline function Gassets(r::DataFrameRow)::MSymbol

  (ismissing(r.at)) && (return missing)
  (r.at < 100.0) && (return :assets100)
  #=(100.0 ≤ r.at) && =#(r.at ≤ 3000.0) && (return :assets1003000)
  ( 3000.0 < r.at) && (return :assets3000)

  @assert false #should be MECE
end

#checks bkliab thresholds: bkliab13, bkliab1323, bkliab23
@inline function Gbkliab(r::DataFrameRow)::MSymbol

  (ismissing(r.bkliab)) && (return missing)
  (r.bkliab < 0.33333333) && (return :bkliab13)
  #=(0.33333333 ≤ r.bkliab) && =#(r.bkliab ≤ 0.66666667) && (return :bkliab1323)
  ( 0.66666667 < r.bkliab) && (return :bkliab23)

  @assert false #should be MECE
end

#checks gvkey for as evidence of the compustat merge: noncompustat, :compustat
@inline function Gcompustat(r::DataFrameRow)::MSymbol

  (ismissing(r.gvkey)) && (return :noncompustat)
  (!ismissing(r.gvkey)) && (return :compustat)

  @assert false #should be MECE
end

#classifies the lastdividend date as :lastdiv40, :lastdiv4080, :lastdiv80
@inline function Glastdiv(r::DataFrameRow)::MSymbol

  (ismissing(r.days2lastdiv)) && (return missing)
  (r.days2lastdiv < 40) && (return :lastdiv40)
  #=(40 ≤ r.days2lastdiv) && =#(r.days2lastdiv ≤ 80) && (return :lastdiv4080)
  (80 < r.days2lastdiv) && (return :lastdiv80)

  @assert false #should be MECE
end

#Gets the dividend frequency: :nonspecialdiv, :specialdiv
const SPECIAL_DIV = Dict{Int, MSymbol}(
  0=>missing,
  1=>missing,
  2=>:nonspecialdiv,
  3=>:nonspecialdiv,
  4=>:nonspecialdiv,
  5=>:nonspecialdiv, #for classification purposes
  6=>:nonspecialdiv, #assume yearend is same as annual
  7=>:specialdiv,
  8=>:specialdiv,
  9=>:specialdiv)
@inline function Gspecialdiv(r::DataFrameRow)::MSymbol

  (ismissing(r.distcd3) || ismissing(r.distcd1) || (r.distcd1 ≠ 1)) && (return missing)
  return SPECIAL_DIV[r.distcd3]

  @assert false #should be MECE
end

#Gets the dividend frequency: :monthlydiv, :quarterlydiv, :semiandannualdiv
#codes 0 and 1 set to quarterly in Hartzmark Solomonb
const FREQ_DIV = Dict{Int, MSymbol}(
  0=>missing,
  1=>missing,
  2=>:monthlydiv,
  3=>:quarterlydiv,
  4=>:semiandannualdiv,
  5=>:semiandannualdiv, #for classification purposes
  6=>:semiandannualdiv, #assume yearend is same as annual
  7=>missing,
  8=>missing,
  9=>missing)
@inline function Gfreqdiv(r::DataFrameRow)::MSymbol

  (ismissing(r.distcd3) || ismissing(r.distcd1) || (r.distcd1 ≠ 1)) && (return missing)
  return FREQ_DIV[r.distcd3]

  @assert false #should be MECE
end

#Gets the dividend tax status: :taxdiv, :nontaxdiv
const TAX_STATUS = Dict{Int, MSymbol}(
  0=>missing,
  1=>missing,
  2=>:taxdiv,
  3=>:nontaxdiv,
  4=>:nontaxdiv, #per Ivo pg 27, return of capital dividends are considered non-taxable
  5=>:taxdiv,
  6=>:taxdiv,
  7=>:taxdiv,
  8=>:taxdiv,
  9=>:nontaxdiv)
@inline function Gtaxdiv(r::DataFrameRow)::MSymbol

  (ismissing(r.distcd4) || ismissing(r.distcd1) || (r.distcd1 ≠ 1)) && (return missing)
  return TAX_STATUS[r.distcd4]

  @assert false #should be MECE
end

#Yield grouping: :lowyield, :medyield, :highyield
@inline function Gyield(r::DataFrameRow)::MSymbol

  (ismissing(r.yield) || (!r.primary)) && (return missing)
  (r.yield < 0.0075) && (return :lowyield)
  #=(0.0075≤r.yield) && =#(r.yield≤0.015) && (return :medyield)
  (0.015 < r.yield) && (return :highyield)

  @assert false #should be MECE
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Run principal components analysis (PCA)

export
    projectensemble,
    projectstructure,
    projectpdbfile,
    projectpdbfiles,
    fluctuations,
    findimportantpcs

"""
Runs principal components analysis on an ensemble and returns a `PCA` object.
Assumes the ensemble is aligned sensibly.
"""
function PCA(ensemble::ModelledEnsemble)
    # Calculate the covariance matrix
    strucs = ensemble.strucs
    n_strucs = length(strucs)
    @assert n_strucs > 0 "The ensemble does not contain any structures"
    @assert n_strucs > 1 "The ensemble contains 1 structure - need multiple to do PCA"
    av_coords = centroid(ensemble)
    n_atoms = size(av_coords, 2)
    @assert n_atoms == size(strucs[1].coords, 2) "Number of atoms in the ensemble structures and atoms are different"
    inds_to_use = calphaindices(ensemble.atoms)
    n_to_use = length(inds_to_use)
    @assert n_to_use > 0 "Not calculating PCA on any atoms"
    cov = zeros(3*n_to_use, 3*n_to_use)
    av_to_use = getindex(av_coords, collect(1:3), inds_to_use)
    for struc in strucs
        disps = (getindex(struc.coords, collect(1:3), inds_to_use) - av_to_use)[:]
        cov += disps * transpose(disps)
    end
    cov /= n_strucs
    # Find the eigendecomposition
    eig = eigen(cov)
    # Number of non-zero eigenvalues is limited by number of structures
    n_non_zero = min(3*n_to_use-6, n_strucs-1)
    # Reverse the values in order to put the most important PCs first
    pca = PCA(reverse(eig.values)[1:n_non_zero], eig.vectors[:, end:-1:1][:, 1:n_non_zero], av_to_use)
    # Project each structure onto the PCs to get the principal coordinates
    pca.pcs = projectensemble(ensemble, pca)
    return pca
end

"""
Project the structures in an ensemble onto the eigenvectors of a PCA to get the principal coordinates.
Non-zero eigenvalues only are returned.
Assumes the ensemble is aligned and corresponds exactly to the PCA.
"""
function projectensemble(ensemble::ModelledEnsemble, pca::PCA)
    strucs = ensemble.strucs
    n_strucs = length(strucs)
    @assert n_strucs > 0 "The ensemble does not contain any structures"
    av_coords = pca.av_coords
    n_atoms = size(av_coords, 2)
    inds_to_use = calphaindices(ensemble.atoms)
    @assert n_atoms == length(inds_to_use) "Number of atoms in the ensemble and in the PCA are different"
    n_evecs = size(pca.evecs, 2)
    pcs = zeros(n_evecs, n_strucs)
    for j in 1:n_strucs
        disps = getindex(strucs[j].coords, collect(1:3), inds_to_use) - av_coords
        for i in 1:n_evecs
            pcs[i, j] = dot(pca.evecs[:,i], disps[:])
        end
    end
    return pcs
end

"""
Project a list of atoms onto the eigenvectors of a PCA.
Atom map is the mapping of the input atom indices to the reference atoms.
It is not required if the atoms correspond to the PCA.
`calpha_inds` is only required if atom_map is non-trivial and is a list of C-alpha indices.
0 means the atom is not in the reference atoms.
Assumes the structure is aligned.
"""
function projectstructure(atoms::Array{Atom,1},
                    pca::PCA;
                    atom_map::Array{Int,1}=collect(1:length(atoms)),
                    calpha_inds::Array{Int,1}=Int[])
    @assert length(atoms) > 0 "No atoms in atom list"
    # If not all atoms are present, knowledge of the missing c-alphas is required
    if 0 in atom_map
        @assert length(calpha_inds) > 0
    end
    init_to_use = calphaindices(atoms)
    inds_to_use_in = Int[]
    inds_to_use_ref = Int[]
    evec_inds = Int[]
    j = 0
    for i in 1:length(atom_map)
        if atom_map[i] in init_to_use
            j += 1
            push!(inds_to_use_in, atom_map[i])
            push!(inds_to_use_ref, j)
            push!(evec_inds, 3j-2, 3j-1, 3j)
        elseif atom_map[i] == 0 && i in calpha_inds
            j += 1
        end
    end
    n_evecs = size(pca.evecs, 2)
    pcs = zeros(n_evecs, 1)
    coords_to_use = getindex(atomcoords(atoms), collect(1:3), inds_to_use_in)
    disps = coords_to_use - getindex(pca.av_coords, collect(1:3), inds_to_use_ref)
    for i in 1:n_evecs
        pcs[i] = dot(getindex(pca.evecs, evec_inds, i), disps[:])
    end
    return pcs
end

"""
Projects the structure in a PDB file onto the eigenvectors of a PCA.
The residue numbering between the PDB files must be consistent but missing/extra
residues are allowed.
Reads the PDB file, maps the atoms to the ref atoms, aligns the atoms to the ref
atoms then gets the projections.
Arguments are PDB filepath, reference atoms to get atom map, coords to align to
and the PCA object.
"""
function projectpdbfile(pdb_filepath::AbstractString,
                    atoms_ref::Array{Atom,1},
                    ens_al::Array{Float64},
                    pca::PCA)
    @assert isfile(pdb_filepath) "Not a valid filepath: \"$pdb_filepath\""
    atoms_in = readpdb(pdb_filepath)
    atom_map = atommap(atoms_in, atoms_ref)
    # Input atoms not in the ref atoms are removed to allow alignment
    atoms_in_min = Atom[]
    ens_al_min = Float64[]
    atom_map_new = Int[]
    j = 0
    for i in 1:length(atoms_ref)
        if atom_map[i] != 0
            push!(atoms_in_min, atoms_in[atom_map[i]])
            push!(ens_al_min, ens_al[:,i]...)
            j += 1
            push!(atom_map_new, j)
        else
            push!(atom_map_new, 0)
        end
    end
    alignatoms!(atoms_in_min, reshape(ens_al_min, 3, j))
    pcs = projectstructure(atoms_in_min, pca, atom_map=atom_map_new, calpha_inds=calphaindices(atoms_ref))
    return pcs
end

"Wrapper function for projectpdbfile to project multiple PDBs and write the output."
function projectpdbfiles(out_prefix::AbstractString,
                    pdb_filepaths::Array{<:AbstractString,1},
                    atoms_ref::Array{Atom,1},
                    ens_al::Array{Float64},
                    pca::PCA)
    checkfilepath(out_prefix)
    @assert length(pdb_filepaths) > 0 "No entries in pdb_filepaths"
    pcs_all = Float64[]
    for (i, pdb_filepath) in enumerate(pdb_filepaths)
        pcs = projectpdbfile(pdb_filepath, atoms_ref, ens_al, pca)
        writeprojections("$(out_prefix)$i.tsv", pcs)
        if length(pcs_all) > 0
            pcs_all = hcat(pcs_all, pcs)
        else
            pcs_all = pcs[:]
        end
    end
    return pcs_all
end

"""
Gets mean square fluctuation for each atom in an ensemble.
Assumes the ensemble is aligned sensibly.
Returns a list of mean square fluctuations.
"""
function fluctuations(ensemble::ModelledEnsemble)
    strucs = ensemble.strucs
    n_strucs = length(strucs)
    @assert n_strucs > 0 "The ensemble does not contain any structures"
    @assert n_strucs > 1 "The ensemble contains 1 structure - need multiple to get fluctuations"
    av_coords = centroid(ensemble)
    n_atoms = size(av_coords, 2)
    @assert n_atoms == size(strucs[1].coords, 2) "Number of atoms in the ensemble and in the PCA are different"
    inds_to_use = calphaindices(ensemble.atoms)
    n_to_use = length(inds_to_use)
    @assert n_to_use > 0 "Not calculating fluctuations on any atoms"
    flucs = zeros(1, n_to_use)
    av_to_use = getindex(av_coords, collect(1:3), inds_to_use)
    for struc in strucs
        disps = (getindex(struc.coords, collect(1:3), inds_to_use) - av_to_use)
        flucs += sum(disps .* disps, dims=1)
    end
    flucs /= n_strucs
    return flucs[:]
end

"""
Takes two sets of projections onto the PCs of a PCA and returns the PC indices ordered
such that the first represents the largest difference between the two structures.
"""
function findimportantpcs(pcs_one::Array{Float64}, pcs_two::Array{Float64})
    @assert size(pcs_one) == size(pcs_two) "PC arrays have different sizes"
    @assert size(pcs_one, 2) == 1 "Each set of projections should be of one structure"
    diff = abs.(pcs_one - pcs_two)
    pc_ord = sortperm(diff[:], rev=true)
    return pc_ord
end


text:
using Gumbo: parsehtml
using Cascadia: @sel_str, nodeText
using DataFrames, JDF
using HTTP: get
import HTTP
using TableScraper

export scrape_kifu_depot_table#, download_goratings

function right(str, start)
    String(str[5:end])
end

"""
Scrape table from kifudepot.org
"""
function scrape_kifu_depot_table(; player="", event="", page::Int=1)
    url = "https://kifudepot.net/index.php?page=$page&move=&player=$player&event=$event"
    _scrape_kifu_depot_table(url)
end

function _scrape_kifu_depot_table(url)
    df = DataFrame(scrape_tables(url, identity)[1])
    rename!(df, ["comp", "black", "white", "result", "date"])

    # extract the link to sgf here
    df[!, :kifu_link] = [comp[1].attributes["href"] for comp in df.comp]

    for n in ["comp", "black", "white", "result", "date"]
        df[!, n] = nodeText.(df[!, n])
    end

    tmp = right.(df[!, "result"], 5)

    tmp1 = split.(tmp, "+")

    df[!,:result] .= tmp

    df[!,:who_win] = [String(elem[1]) for elem in tmp1]

    function not_void(elem)
        if length(elem) == 2
            return String(elem[2])
        elseif length(elem) == 1
            return ""
        else
            return "check"
        end
    end

    df[!,:win_by] = [not_void(elem) for elem in tmp1]

    df
end

function extract_sgf(link)
    game = HTTP.get("https://kifudepot.net/"*link)
    eachmatch(sel"#sgf", parsehtml(String(game.body)).root)[1] |> nodeText
end

function extract_komi(sgf::String)
    # from https://stackoverflow.com/questions/66539167/how-to-capture-just-the-next-bracket-using-regular-expression-regex/66539200#66539200
    m = match(r"KM\[(?<komi>[^\]]*)", sgf)
    parse(Float64, m.captures[1])
end

function download_goratings(path)
    df = DataFrame(scrape_tables("https://www.goratings.org/en/")[2])
    rename!(df, [:rank, :player, :sex, :flag, :elo])
    select!(df, [:rank, :player, :elo])
    df[!, :elo] = parse.(Int, df[!, :elo])
    df[!, :rank] = parse.(Int, df[!, :rank])
    JDF.save(path, df)
end


text:
# this is a script to train a basic VAE on the protobuffer pdf data
using GenerativeModels
using GMExtensions
using ConditionalDists
using PdfFraud
using JLD2
using PyPlot
using FileIO
using Flux
using CuArrays
import GenerativeModels: elbo
using Dates
using ValueHistories
using ProgressMeter

"""
	elbo(model, x[; β])

Elbo loss function for convolutional models and 4D inputs.
"""
function elbo(m::GenerativeModels.AbstractVAE, x::AbstractArray{T,4}; β=1) where T
    z = rand(m.encoder, x)
    vecdim = reduce(*, size(x)[1:3])
    llh = mean(-loglikelihood(m.decoder, reshape(x, vecdim, :), z))
    kld = mean(kl_divergence(m.encoder, m.prior, x))
    llh + β*kld
end

"""
	reconstruct(m,x)

Reconstructs a sample `x` through the model `m`.
"""
reconstruct(m,x) = reshape(mean(m.decoder, rand(m.encoder, x)), size(x))
reconstruct(x) = reconstruct(model, x)

"""
	history_progress_cb(prog, model, test_data, history, elbo)

Returns a callback function that stores loss values in the history object and also
prints the current progress.
"""
function history_progress_cb(prog::Progress, model::GenerativeModels.AbstractVAE, test_data::AbstractArray,
							 h::MVHistory, elbox::Function)
	function cb()
		# get losses
		elbo = elbox(test_data) |> cpu
		mse = Flux.mse(test_data, reconstruct(model, test_data)) |> cpu

		# write them to history
		ntuple = GMExtensions.@ntuple elbo mse
		GMExtensions.push_ntuple!(h, ntuple)

		# update the progress bar
		next!(prog, showvalues = [(k, ntuple[k]) for k in keys(ntuple)])
	end
end

# first gather data and their labels
extracted_data = load(PdfFraud.extracted_data)
representations = extracted_data["representations"]
sha256 = extracted_data["sha256"]

# nope ... stack overflow errs
# X = cat(representations..., dims=3)
lens = map(x->size(x,3), representations)
inds = cumsum(lens)
X = zeros(64,48,1,sum(lens))
for (i,representation) in enumerate(representations)
	X[:,:,:,(inds[i]-lens[i]+1):inds[i]] = reshape(representation, 64, 48, 1, :)
end
X = X |> gpu

# for testing purposes, delete after
x = X[:,:,:,1:2]

(h, w, c) = size(X)[1:3]
ldim = 32
kernelsizes = (3, 5)
nchannels = (4, 8)
scalings = (2, 2)
densedims = [256]
vecdim = reduce(*, (h,w,c))
T = Float32

encoder = GMExtensions.conv_encoder((h, w, c), ldim*2, kernelsizes, nchannels, scalings;
	densedims = densedims) |> gpu
enc_dist = CMeanVarGaussian{DiagVar}(encoder)

decoder = GMExtensions.conv_decoder((h, w, c), ldim, reverse(kernelsizes),
	reverse(nchannels), reverse(scalings); densedims = densedims, vec_output = true,
	vec_output_dim = vecdim + 1) |> gpu
dec_dist = CMeanVarGaussian{ScalarVar}(decoder)

model = VAE(ldim, enc_dist, dec_dist)

z = rand(model.encoder, x)
y = mean(model.decoder, z)

β = 1.0
loss(x) = elbo(model, x, β=β)
opt = ADAM(0.0005)
GenerativeModels.update_params!(model, (x,), loss, opt)

_data = [(x,) for _ in 1:500]
h = MVHistory()
prog = Progress(20*length(_data))
cb = history_progress_cb(prog, model, x, h, loss)

Flux.@epochs 20 Flux.train!(loss, params(model), _data, opt, cb = cb)

rx = reconstruct(x)

figure(figsize=(10,5))
subplot(121)
pcolormesh(x[:,:,1,1])
subplot(122)
pcolormesh(rx[:,:,1,1])

figure(figsize=(10,5))
subplot(121)
pcolormesh(x[:,:,1,2])
subplot(122)
pcolormesh(rx[:,:,1,2])

#cmodel = model |> cpu

filename = "$(now()).bson"

save_checkpoint(filename, model, h)

# test loading of the model
filename = filter(x->occursin("bson", x), readdir("."))[1]
model, h = load_checkpoint(filename)
model = model |> gpu
rx = reconstruct(x)


text:
using StochasticDiffEq
using Plots; pyplot()
using LaTeXStrings
f(u,p,t) = -1000.0u*(1-u)*(2-u)
g(u,p,t) = 10.0
prob = SDEProblem(f,g,2.0,(0.0,5.0))

sol = solve(prob,SOSRI(),abstol=1e-2,reltol=1e-2,saveat=0.025,seed=3)
plot(sol,legend=false,yaxis = L"X_t",title = "Pathwise Stiffness Example",
     xtickfont = font(16, "Latin Modern Roman"),titlefont = font(20, "Latin Modern Roman"),
     ytickfont = font(16, "Latin Modern Roman"),guidefont = font(18, "Latin Modern Roman"))

savefig("pathwise_stiffness.png")
savefig("pathwise_stiffness.pdf")