text:
#=
Proof of concept of 2D localization against a scalar field.
In this example we use a sequence of 1D measurements and a known map with a
correlator-like approach to localize.

Case 4: Square helix northward, with overlapping east/west segments
no prior map

2021-03-24,31
=#

using Images, FileIO           # load elevation from PNG
using ImageFiltering
using DSP
# using DataInterpolations       # handles the terrain
using Interpolations
using Distributions
using Statistics
using Caesar, RoME
using Gadfly, Colors
using Cairo
using RoMEPlotting
using Plots
using JSON2
# using ImageView

Gadfly.set_default_plot_size(35cm,25cm)
st=style( major_label_font="CMU Serif",
          minor_label_font="CMU Serif",
          major_label_font_size=20pt,
          minor_label_font_size=20pt)

prjPath = dirname(dirname(pathof(Caesar)))
# this includes Sequences module
include( joinpath(prjPath, "examples","dev","scalar","Sequences.jl") )

# load dem (18x18km span)
img = load(joinpath(prjPath, "examples","dev","scalar","dem.png")) .|> Gray
img = Float64.(img)
# imshow(img)
# interpolant object for querying
x = range(-9e3, 9e3, length = size(img)[1]) # North
y = range(-9e3, 9e3, length = size(img)[2]) # East
dem = Interpolations.LinearInterpolation((x,y), img)

# not to scale on DEM image, but sequence correlation approach does not matter
terrE_ = 1e3*Float64.(@view img[200,:])
terrW_ = 1e3*Float64.(@view img[300,:])
x = range(-100, 100, length = length(terrE_))

global terrE = Interpolations.LinearInterpolation(x, terrE_)
global terrW = Interpolations.LinearInterpolation(x, terrW_)

##

"""
    Point2Point2Northing

IIF will deal with partials,

DevNotes
- FIXME ongoing (2021Q2), see #1206, likely also #1010
"""
struct Point2Point2Northing{B <: IIF.SamplableBelief, T <: Tuple} <: IIF.AbstractRelativeMinimize
  Z::B
  partial::T # which dimension the partial applies to
end

function IIF.getSample(s::CalcFactor{<:Point2Point2Northing}, N::Int=1)
  return ( reshape(rand(s.factor.Z,N),1,N), )
end

(s::CalcFactor{<:Point2Point2Northing})(z, x1, x2) = z .- (x2[1] - x1[1])

##

# matchLeg!(fg, [:x1, :x5], :POSITIVE_X)
# matchLeg!(fg, [:x3, :x7], :NEGATIVE_X)
function matchLeg!( fg::AbstractDFG,
                    legs::AbstractVector{Symbol},
                    direction::Symbol;
                    dofactor::Bool=false,
                    odoPredictedAlign::Real=0,
                    kappa::Real=2)
  #

  dataEntry, dataBytes = getData(fg, legs[1], :elevationSequences)
  myData_a = JSON2.read(IOBuffer(dataBytes), Dict{Symbol, Vector{Float64}})

  dataEntry, dataBytes = getData(fg, legs[2], :elevationSequences)
  myData_b = JSON2.read(IOBuffer(dataBytes), Dict{Symbol, Vector{Float64}})

  xseq_a = myData_a[:x_seq]
  zseq_a = myData_a[:z_seq]

  xseq_b = myData_b[:x_seq]
  zseq_b = myData_b[:z_seq]

  # correlate z_seq0 against z_seq_m4
  s_a = Sequences.MeasurementSequence(xseq_a .- xseq_a[1], zseq_a) # full relative
  s_b = Sequences.MeasurementSequence(xseq_b .- xseq_b[1] .- odoPredictedAlign, zseq_b)
  # FIXME ON FIRE dont have odo ground truth
  # see also: [x,i] = ssdcorr((x,z)_a, (x,z)_b)

  # closure on a unique correlator with s_a and s_b
  _ssdCorr(x) = Sequences.ssd(s_a, Sequences.displace(s_b,x))

  # qr is the range of displacements over which to compute the  SSD
  qr = collect(LinRange(-10.0,10.0,1001))
  intensity = _ssdCorr.(qr)
  # @assert length(intensity) == length(qr)
  # strip nans when no overlap
  mask::BitVector = isnan.(intensity) # .|> x->xor(x, true)
  intensity[mask] .= maximum(intensity[xor.(mask, true)]) #intensity[mask]
  k = abs(maximum(intensity))
  pv = intensity./k  # normalize
  pv = exp.(-pv)
  w = DSP.Windows.tukey(length(intensity), 0.25)
  pv .*= w
  pv .^= kappa
  # hand fix, to emphasize peaks

  # pv ./= sum(pv) # sum to 1, already in AliasingScalarSampler

  # intensity = _mySSDCorr(zseq_a, zseq_b)

  bel = AliasingScalarSampler(qr, pv)
  plr = Point2Point2Northing(bel, (1,))

  f_ = if dofactor
    addFactor!(fg, legs, plr, tags=[:TERRAIN,:MATCH,direction])
  end

  return (qr, pv), f_
end

##

# callback to structure sequences
function newcallback_ex4(fg::AbstractDFG, lastPose::Symbol)
  # temp workaround for terrain interpolators

  @show lastPose
  @show tags = getVariable(fg, lastPose).tags
  @show prevPose = incrSuffix(lastPose, -1)

  # don't do anything extra on first pose :x0
  if lastPose == :x0
    return nothing
  end

  # true location of poses in world
  simPos = getPPE(fg, lastPose, :simulated).suggested
  # generate xseq and zseq using this and previous pose data
  simPosPrev = getPPE(fg, prevPose, :simulated).suggested
  # linear for Point2Point2
  @show simOdo = simPos - simPosPrev

  # add a partial prior in Y if no full prior
  if 0 == length( intersect(ls(fg, lastPose), lsf(fg, PriorPoint2)) )
    pp = PartialPrior(Normal(simPos[2], 0.1),(2,))
    @info "adding partial prior on simulated y position" lastPose simPos[2]
    addFactor!(fg, [lastPose;], pp, graphinit=false, tags=[:EX4_PARTIALPRIOR_Y])
  end

  # assumed 100 measurements per leg
  Nseq = 100

  # avert your eyes: this assumes v is axis-aligned
  # To pay this technical debt: generate a 2D linspace and use it to query a 2D DEM interpolant
  terr, x_seq = if abs(simOdo[2]) < abs(simOdo[1])
    global terrW
    # assume north south
    xs = LinRange(simPos[1], simPosPrev[1], Nseq)
    terrW, xs
  elseif abs(simOdo[1]) < abs(simOdo[2])
    global terrE
    # assume east west
    xs = LinRange(simPos[2], simPosPrev[2], Nseq)
    terrE, xs
  end

  z_seq = terr.(x_seq)
  # add noise to x_seq (odometry) and z_seq (meas noise)
  x_seq_n = 0.01*randn(Nseq)
  z_seq_n = 1.0*randn(Nseq)

  @assert !isapprox( Statistics.median(diff(x_seq)), 0, atol=1e-6) "why is x_seq[1:5]=$(x_seq[1:5])"

  someDict = Dict(:x_seq => x_seq, :z_seq => z_seq, :x_seq_n => x_seq_n, :z_seq_n => z_seq_n)
  addData!( fg, :default_folder_store, lastPose, :elevationSequences,
            Vector{UInt8}(JSON2.write( someDict )),
            mimeType="application/json/octet-stream"  )
  #

  nothing
end

##

# trajectory parameters
NS = 15
skew_x = 2/3

# setup the factor graph
fg = initfg()

# adding a blog store
# getSolverParams(fg).logpath = pwd()
storeDir = joinLogPath(fg,"data")
mkpath(storeDir)
datastore = FolderStore{Vector{UInt8}}(:default_folder_store, storeDir)
addBlobStore!(fg, datastore)

##

# drive the first box, no correlations
RoME.generateCanonicalFG_Boxes2D!(8, dfg=fg, skew_x=skew_x, postpose_cb=newcallback_ex4)

##

# 1st match: north-bound legs on first and second loop
xsq, f_ = matchLeg!(fg, [:x0, :x4], :POSITIVE_X, odoPredictedAlign=0, kappa=3, dofactor=false);
# xsq, f_ = matchLeg!(fg, [:x1, :x5],:NORTH, odoPredictedAlign=0, kappa=3, dofactor=true)
# p=Gadfly.plot(st)
# push!(p,layer(x=xsq[1], y=xsq[2], Geom.line))

# 2nd match: south-bound legs on first and second loop
xsq, f_ = matchLeg!(fg, [:x2, :x6], :NEGATIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true);

# 1st match: north-bound legs on first and second loop
xsq, f_ = matchLeg!(fg, [:x4, :x8],:POSITIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true);
# xsq, f_ = matchLeg!(fg, [:x1, :x5],:POSITIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true)
# p=Gadfly.plot(st)
# push!(p,layer(x=xsq[1], y=xsq[2], Geom.line))

# 2nd match: south-bound legs on first and second loop
xsq, f_ = matchLeg!(fg, [:x6, :x10],:NEGATIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true);

xsq, f_ = matchLeg!(fg, [:x8, :x12],  :POSITIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true);
xsq, f_ = matchLeg!(fg, [:x10, :x14], :NEGATIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true);

# xsq, f_ = matchLeg!(fg, [:x12, :x16],  :POSITIVE_X, odoPredictedAlign=0, kappa=3, dofactor=true)

##

pts = approxConv(fg, :x0x4f1, :x4)

# should peak at +5 to make the relative factor represent the change between x0 dn x4
Gadfly.plot(x=pts[1,:], Geom.histogram)

## TEMPORARY CORRELATION DEV

tree, _, _ = solveTree!(fg)

##

plotKDE(fg, [:x0;:x1;:x2;:x3;:x4], levels=1)
plotKDE(fg, [:x4;:x5;:x6;:x7;:x8], levels=1)

plotKDE(fg, [:x8;:x9;:x10;:x11;:x12], levels=1)
plotKDE(fg, [:x12;:x13;:x14;:x15;:x16], levels=1)

##

# get pose four poses back
# get terrain sequence between (newPose - NW) up to (newPose)
# curr_seq = terrE(curr_x_seq)
# prev_seq = terrE(prev_x_seq)
# # compute SSDcorr wrt sequence 4 poses before
# # add Point2Point2Northing between the two
# intensity = _mySSDCorr(curr_seq, prev_seq)

plr = ScalarFieldSequenceFactorNorthSouth(x_seq_, z_seq, (1,))
addFactor!(fg, [:x0; :x4], plr, tags=[:TERRAIN,:MATCH])
    # plr = ScalarFieldSequenceFactorNorthSouth(x_seq, z_seq, (1,))
    # # # plr = Point2Point2Northing(grid, AliasingScalarSampler(intensity))
    # addFactor!(fg, [prevPose newPose], plr, tags=[:TERRAIN,:MATCH])

##


text:
using TaylorModels

# import Documenter
# Documenter.makedocs(
#     modules = [TaylorModels],
#     format = :html,
#     sitename = "TaylorModels.jl",
#     root = joinpath(dirname(dirname(@__FILE__)), "docs"),
#     pages = Any["Home" => "index.md"],
#     strict = true,
#     linkcheck = true,
#     checkdocs = :exports,
#     authors = "David Sanders"
# )

# using Base.Test

include("TM1.jl")
include("TMN.jl")

ext: jl


text:
module TestDiffMap
    using ManifoldLearning
    using Test

    @testset "DiffMap" begin

	d = 2
	t = 1
	X, L = swiss_roll()
	Y = transform(DiffMap, X; d=d, t=t)

	@test outdim(Y) == d
	@test size(projection(Y), 2) == size(X, 2)
    end
end


text:
using .MiniOmega

function simpletest()
  a = 1 ~ unif
  b = 2 ~ a
  sample(b)
end

function test()
  uniform(a, b) = ω -> unif(ω) * (b - a) + a

  x = 1 ~ uniform(0, 1)
  y = 2 ~ uniform(0, 1)
  d = 7 ~ y
  z = ω -> (x(ω), x(ω), y(ω), d(ω))
  @show sample(z)
  function a(ω)
    x_ = x(ω)
    d = 3 ~ uniform(0, 4)
    e = 4 ~ uniform(0, 4)
    x_ + d(ω) + e(ω)
  end

  @show sample(rt(x, y, z, a))

  # Conditional Independence
  x = 1 ~ uniform(0, 19)
  measure(ω) = x(ω) + uniform(0, 1)(ω)
  m1 = 3 ~ measure <| (x,)
  m2 = 4 ~ measure <| (x,)
  m3 = 5 ~ measure

  @show sample(rt(m1, m2, m3))

  # Linaer regression
  α = 1 ~ uniform(0, 1)
  β = 2 ~ uniform(0, 1)
  f(x, i) = i ~ (ω -> x*α(ω) + β(ω) + uniform(0.0, 1.0)(ω)) <| (α, β)
  y1 = f(0.3, 3)
  y2 = f(0.3, 4)
  @show sample(rt(y1, y2))

  y1b, y2b = pointwise() do
    f2(x, i) = i ~ (x * α + β + uniform(0.0, 1.0)) <| (α, β)
    y1b = f2(0.3, 3)
    y2b = f2(0.3, 4)
    y1b, y2b
  end

  @show sample(rt(y1, y2, y1b, y2b))

  # A bunch of random variables with a shared parents
  pointwise() do
    x = 1 ~ unifQD
    y = x * 100
    function a(ω)
      y(ω) + 5
    end
    function b(ω)
      y(ω) + 10
    end
    c = 2 ~ b <| (y,)
    d = 3 ~ b
    @show sample(rt(y, a, b, c, d))
  end
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
    likelihood(events::Array{<:Number, 1}, background::Float64, kappa::Float64, kernel::Distributions.Distribution, maxT::Number)
Calculate the log likelihood of a Hawkes process for given parameters.

# Arguments

* `events` Vector of events to calculate the likelihood for.
* `background` Background rate.
* `kappa` Kappa parameter.
* `kernel` Function or distribution of the kernel.
* `maxT` Maximum time of the process.

# Notes

* The kernel function must be a proper probability distribution.

# Examples

"""
function likelihood(events::Array{<:Number, 1}, background::Float64, kappa::Float64, kernel::Distributions.Distribution, maxT::Number)
    kernel_f, kernel_cdf = generate_kernel_functions(kernel)
    likelihood(events, background, kappa, kernel_f, kernel_cdf, maxT)
end

function likelihood(events::Array{<:Number, 1}, background::Float64, kappa::Float64, kernel::Function, kernel_cdf::Function, maxT::Number)

    logIntensitySum = sum(log.(intensity(events, events, background, kappa, kernel)))
    backgroundIntegral = background * maxT
    kernelIntegral = sum(kernel_cdf(maxT .- events))

    logIntensitySum - backgroundIntegral - kappa*kernelIntegral
end


text:
################################################################################
# MergingRoadwayOptions
################################################################################

@with_kw mutable struct MergingRoadwayOptions{T} <: RoadwayOptions
    # Options
	"Length of the lanes."
	lane_length::T=6.0

	"Width of the lanes."
	lane_width::T=0.20

	"Position of the incoming lane merging point."
	merging_point::T=1.2

	"Orientation of the incoming lane."
	merging_angle::T=pi/12
end

function build_roadway(opts::MergingRoadwayOptions{T}) where {T}
	ll = opts.lane_length
	lw = opts.lane_width
	mp = opts.merging_point
	θ  = opts.merging_angle
	#     y
	#     ^          x3---------------------------------------x-----x2
	# 	  |							  merging          -      x9       w1
	# 	  |							   point x7   -            x10     w2
	# (0,0) -->x     x4-----------------x-----------------------x---x1
	# 		    merging angle θ  -       \              -       w8  w3
	# 					  -               \      -      w7          w4
	# 		x11    - x5--------------------x------------------------x6
	#       -                       -      x8           w9
	#       w5               -      w6
	#                 -
	#                 x12
	x1  = [ ll,  0.]
	x2  = [ ll,  lw]
	x3  = [ 0.,  lw]
	x4  = [ 0.,  0.]
	x5  = [ 0., -lw]
	x6  = [ ll, -lw]

	x7  = [ mp,  0.]
	x8  = [ mp+lw*tan(θ), -lw]
	x9  = [ mp+lw/tan(θ),  lw]
	x10 = [ mp+lw*tan(θ)+lw/tan(θ), 0.]
	x11 = [ mp-mp*cos(θ), -mp*sin(θ)]
	x12 = [ mp+lw*tan(θ)-mp*cos(θ), -lw-mp*sin(θ)]

	v = [0., 1.]
	w = [-sin(θ), cos(θ)]
	w1 = Wall(x3, x2,  v)
	w2 = Wall(x4, x1, -v)
	w3 = Wall(x4, x1,  v)
	w4 = Wall(x5, x6, -v)

	w5 = Wall(x11, x7,   w)
	w6 = Wall(x12, x8,  -w)
	w7 = Wall(x12, x10, -w)
	@assert x10[1] < x1[1]
	w8 = Wall(x10, x1,  -v)
	@assert x8[1] < x1[1]
	w9 = Wall(x8, x6,   -v)

	# Centerlines
	f_straight_left(s,v=0.)  = VehicleState(s,  lw/2, 0., v)
	f_straight_right(s,v=0.) = VehicleState(s, -lw/2, 0., v)
	f_straight_both(s,v=0.)  = VehicleState(s, -lw/2, 0., v)
	function f_merging_left(s,v=0.)
		if s <= mp
			return VehicleState(mp+3lw/2*tan(θ)-(mp-s)*cos(θ), lw/2-(mp-s)*sin(θ), θ, v)
		else
			return VehicleState(3lw/2*tan(θ)+s,  lw/2, 0., v)
		end
	end
	function f_merging_right(s,v=0.)
		if s <= mp
			return VehicleState(mp+lw/2*tan(θ)-(mp-s)*cos(θ), -lw/2-(mp-s)*sin(θ), θ, v)
		else
			return VehicleState(lw/2*tan(θ)+s, -lw/2, 0., v)
		end
	end
	function f_merging_both(s,v=0.)
		if s <= mp
			return VehicleState(mp+lw/2*tan(θ)-(mp-s)*cos(θ), -lw/2-(mp-s)*sin(θ), θ, v)
		else
			return VehicleState(lw/2*tan(θ)+s, -lw/2, 0., v)
		end
	end

	# Straight Left lane
	straight_left_lane = Lane(1,
					 :straight_left_lane,
					 [w1,w2],
					 Vector{CircularWall}(),
					 StartingArea(
						 VehicleState(0.,  lw/2, 0.0,  0.00),
						 VehicleState(0.,  lw/2, 0.0, -0.05),
						 VehicleState(ll,  lw/2, 0.4,  0.05),
						 ),
					 f_straight_left,
					 )

	# Straight Right lane
	straight_right_lane = Lane(2,
					  :straight_right_lane,
					  [w3,w4],
					  Vector{CircularWall}(),
					  StartingArea(
						  VehicleState(0., -lw/2, 0.0,  0.00),
						  VehicleState(0., -lw/2, 0.0, -0.05),
					      VehicleState(ll, -lw/2, 0.4,  0.05),
					      ),
					  f_straight_right,
					  )

	# Straight Both lane
	straight_both_lane = Lane(3,
					 :straight_both_lane,
					 [w1,w4],
					 Vector{CircularWall}(),
					 StartingArea(
						 VehicleState(0., -lw/2, 0.0,  0.00),
						 VehicleState(0., -lw/2, 0.0, -0.05),
						 VehicleState(ll, -lw/2, 0.4,  0.05),
						 ),
					 f_straight_both,
					 )

 	# Merging Left lane
 	merging_left_lane = Lane(4,
 					 :merging_left_lane,
					 [w1,w5,w7,w8],
 					 Vector{CircularWall}(),
 					 StartingArea(
						 VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ),
						 VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ-0.05),
						 VehicleState(mp+lw/2*tan(θ),           -lw/2,           0.4, θ+0.05),
 						 ),
 					 f_merging_left,
 					 )

 	# Merging Right lane
 	merging_right_lane = Lane(5,
 					  :merging_right_lane,
					  [w3,w5,w6,w9],
 					  Vector{CircularWall}(),
 					  StartingArea(
						  VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ),
						  VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ-0.05),
						  VehicleState(mp+lw/2*tan(θ),           -lw/2,           0.4, θ+0.05),
 					      ),
 					  f_merging_right,
 					  )

 	# Merging Both lane
 	merging_both_lane = Lane(6,
 					 :merging_both_lane,
					 [w1,w5,w6,w9],
 					 Vector{CircularWall}(),
 					 StartingArea(
 						 VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ),
 						 VehicleState(mp+lw/2*tan(θ)-mp*cos(θ), -lw/2-mp*sin(θ), 0.0, θ-0.05),
 						 VehicleState(mp+lw/2*tan(θ),           -lw/2,           0.4, θ+0.05),
 						 ),
 					 f_merging_both,
 					 )

	lanes = [straight_left_lane, straight_right_lane, straight_both_lane,
		merging_left_lane, merging_right_lane, merging_both_lane]
	roadway = Roadway(lanes, opts)
    return roadway
end


text:
import ..Embeddings: AbstractEmbedding

abstract type EmbeddingTriangulationPartition{D, T} <: AbstractTriangulationPartition{D, T} end
abstract type MutableEmbeddingTriangulationPartition{D, T} <: AbstractMutableTriangulationPartition{D, T} end

abstract type EmbeddingTriangulationPartitionFull{D, T} <: AbstractTriangulationPartitionFull{D, T} end
abstract type MutableEmbeddingTriangulationPartitionFull{D, T} <: AbstractMutableTriangulationPartitionFull{D, T} end

struct EmbeddingTriangulation{D, T} <: EmbeddingTriangulationPartition{D, T}
	data::AbstractEmbedding
	simplexindices::DelaunayTriangulation
end

mutable struct MutableEmbeddingTriangulation{D, T} <: MutableEmbeddingTriangulationPartition{D, T}
	data::AbstractEmbedding
	simplexindices::DelaunayTriangulation
end

# Full info
struct EmbeddingTriangulationFull{D, T} <: EmbeddingTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

# Full info
struct MutableEmbeddingTriangulationFull{D, T} <: MutableEmbeddingTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

MT = Union{EmbeddingTriangulationPartition, EmbeddingTriangulationPartitionFull,
    MutableEmbeddingTriangulationPartition, MutableEmbeddingTriangulationPartitionFull}
getpoints(mt::MT) = [mt.data[:, i] for i = 1:maximum(size(mt.data))]

getsimplices(mt::Union{EmbeddingTriangulationPartitionFull, MutableEmbeddingTriangulationPartitionFull}) =
    mt.simplices
getsimplices(mt::Union{EmbeddingTriangulationPartition, MutableEmbeddingTriangulationPartition}) =
    [Simplex(mt.data[:, mt.simplexindices[i]]) for i = 1:length(mt.simplexindices)]

####################
# Pretty printing
####################
function summarise(tp::Union{EmbeddingTriangulation, MutableEmbeddingTriangulation})
    _type = typeof(tp)
    n_simplices = nsimplices(tp)
    n_pts = size(tp.data, 2)
    D = dimension(tp)
    summary = "$_type from $n_simplices simplices and $n_pts points"
end

Base.show(io::IO, tp::Union{EmbeddingTriangulation, MutableEmbeddingTriangulation}) =
    println(io, summarise(tp))

export
EmbeddingTriangulation,
MutableEmbeddingTriangulation,
EmbeddingTriangulationFull,
MutableEmbeddingTriangulationFull


text:
@testset "IndexableGenerator" begin
    @testset "basic" begin
        generator = TimeZones.IndexableGenerator(Char(i) for i in (1:26) .+ 96)

        @test length(generator) == 26
        @test size(generator) == (26,)
        @test axes(generator) == (Base.OneTo(26),)
        @test ndims(generator) == 1

        @test generator[26] == 'z'
        @test lastindex(generator) == 26

        @test collect(generator) == 'a':'z'
    end

    @testset "constructors" begin
        generator = TimeZones.IndexableGenerator(i -> Char(i), (1:26) .+ 96)
        @test collect(generator) == 'a':'z'

        generator = TimeZones.IndexableGenerator(Char, (1:26) .+ 96)
        @test collect(generator) == 'a':'z'
    end
end


text:
tolerance = 10^-6

S = randomGenericSymp(4)

@test size(S) == (8, 8)
@test typeof(S[1, 1]) == BigFloat
@test typeof(S) <: Symp
@test typeof(S) <: AbstractMatrix
@test typeof(Symp(S)) <: Symp
@test typeof(Matrix(S)) <: AbstractMatrix
@test isGeneric(S) == true
@test nModes(S) == 4
@test typeof(transpose(S)) <: Symp
@test typeof(S') <: Symp
@test typeof(-S) <: Symp
@test typeof(+S) <: Symp
@test isapprox( inv(S) * S, Id(4); atol = tolerance )
@test Omega(2) == Symp([0 1 0 0; -1 0 0 0; 0 0 0 1; 0 0 -1 0])
@test dsum(Id(2), Id(2), Id(2)) == Id(6)
@test S^-1 == inv(S)


text:
module MultiFieldTests

using Gridap
using Gridap.FESpaces
using GridapDistributed
using PartitionedArrays
using Test

function main(parts)

  output = mkpath(joinpath(@__DIR__,"output"))

  domain = (0,4,0,4)
  cells = (4,4)
  model = CartesianDiscreteModel(parts,domain,cells)
  Ω = Triangulation(model)

  k = 2
  reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},k)
  reffe_p = ReferenceFE(lagrangian,Float64,k-1,space=:P)

  u((x,y)) = VectorValue((x+y)^2,(x-y)^2)
  p((x,y)) = x+y
  f(x) = - Δ(u,x) + ∇(p,x)
  g(x) = tr(∇(u,x))

  V = TestFESpace(model,reffe_u,dirichlet_tags="boundary")
  Q = TestFESpace(model,reffe_p,constraint=:zeromean)
  U = TrialFESpace(V,u)
  P = TrialFESpace(Q,p)

  VxQ = MultiFieldFESpace([V,Q])
  UxP = MultiFieldFESpace([U,P]) # This generates again the global numbering
  UxP = TrialFESpace([u,p],VxQ) # This reuses the one computed
  @test length(UxP) == 2

  uh, ph = interpolate([u,p],UxP)
  eu = u - uh
  ep = p - ph

  dΩ = Measure(Ω,2*k)
  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9
  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9

  a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - q*(∇⋅u) - (∇⋅v)*p )*dΩ
  l((v,q)) = ∫( v⋅f - q*g )*dΩ

  op = AffineFEOperator(a,l,UxP,VxQ)
  solver = LinearFESolver(BackslashSolver())
  uh, ph = solve(solver,op)

  eu = u - uh
  ep = p - ph

  writevtk(Ω,"Ω",nsubcells=10,cellfields=["uh"=>uh,"ph"=>ph])

  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9
  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9

end

end # module


text:
samplers = [gibbs, langevinprior, langevinposterior, gibbslazy]

sample_pop(f::Function, iters::Int, n::Int=10) = begin
    states = State[State() for i = 1:n]
    results = map(f, states)
    ress = Array(Any, iters*n)
    for i = 1:iters
        state_index = rand(Categorical(n))
        if rand() < 0.5
            (proposed_state, proposed_result) =
                sample(langevinprior, f, states[state_index])
            logp = logpdf_diff(proposed_state, states[state_index])
        else
            (proposed_state, proposed_result) =
                elliptic(f, states[state_index], states[rand(Categorical(n))])
            logp = logpdf_diff(proposed_state, states[state_index])
        end
        (state, res) = sample(samplers[rand(Categorical(current_ps))], f, state, res)
        ress[i] = res
    end
    ress
end

import StatsBase.sample
sample(f::Function, iters::Int; ps=ones(length(samplers))) = begin
    state = State()
    res = f(state)
    ress = Array(Any, iters)
    for i = 1:iters
        current_ps = ps .* map(sampler -> is_sampler_applicable(sampler, state), samplers)
        current_ps = current_ps/sum(current_ps)
        (state, res) = sample(samplers[rand(Categorical(current_ps))], f, state, res)
        ress[i] = res
    end
    ress
end

is_sample{W, S <: Sampler}(::S, ::Sample{W, S}) = true
is_sample(::Sampler, ::TraceElement) = false

is_sampler_applicable_elements(s::Sampler, state::State) =
    map(el -> is_sample(s, el), state.trace)
is_sampler_applicable(s::Sampler, state::State) =
    any(is_sampler_applicable_elements(s, state))

sample(sampler::MH, f::Function, old_state::State, old_res) = begin
    #ngso = n_gibbs_samples(old_state)
    (proposed_state, proposed_res, proposal_logpdf_diff) = sample(sampler, f, old_state)
    #@assert ngso == n_gibbs_samples(old_state)
    log_ap = value(logpdf_diff(sampler, proposed_state, old_state)) + proposal_logpdf_diff
    if rand() < exp(log_ap)
        (proposed_state, proposed_res)
    else
        (old_state, old_res)
    end
end

abstract Count
immutable Lf <: Count
    count::Int
end
immutable Br <: Count
    count::Int
    leaves::Vector{Count}
end
count(sampler::Sampler, si::SI) = begin
    total = 0
    leaves = Array(Count, length(si.trace))
    for i = 1:length(si.trace)
        leaves[i] = count(sampler, si.trace[i])
        total += leaves[i].count
    end
    Br(total, leaves)
end
count{W, S <: Sampler}(::S, ::Sample{W, S}) = Lf(1)
count{S <: Sampler}(s::S, sample::Sample) = Lf(0)

resample_addr(c::Count) = begin
    counter = Int[]
    resample_addr(counter, c)
    counter
end
resample_addr(counter::Vector{Int}, c::Lf) = nothing
resample_addr(counter::Vector{Int}, c::Count) = begin
    ns = map(x -> x.count, c.leaves)
    push!(counter, rand(Categorical(ns/sum(ns))))
    resample_addr(counter, c.leaves[counter[end]])
    nothing
end

import Base.deepcopy
dc(s::State) = State(dc(s.trace))
dc(t::Vector{TraceElement}) = TraceElement[dc(x) for x in t]
dc(i::If) = If(i.cond, dc(i.trace))
#Don't copy when not necessary
dc(s::Sample{Whitened}) = s
dc{S}(s::Sample{UnWhitened, S}) = Sample(unwhitened, S(), s.value, s.logpdf)
dc(s::Sample{Whitened, Conditioned}) = error()
sample(::Gibbs, f::Function, old_state::State) = begin
    old_count = count(gibbs, old_state)
    counter = resample_addr(old_count)
    proposed_state = dc(old_state)
    proposed_state[counter] = resample
    proposed_res = f(proposed_state)
    proposed_count = count(gibbs, proposed_state)
    proposal_diff = log(old_count.count) - log(proposed_count.count) +
                    old_state[counter].logpdf - proposed_state[counter].logpdf
    (proposed_state, proposed_res, proposal_diff)
end
#Gibbs
isgibbs{W}(::Sample{W, Gibbs}) = true
isgibbs(_) = false
isunwhitened(::Sample{UnWhitened}) = true
isunwhitened(u) = false

#GibbsLazy
sample(::GibbsLazy, f::Function, old_state::State) = begin

end

#Langevin prior
sample(::LangevinPrior, f::Function, old_state::State) = begin
    proposed_state = State(TraceElement[langevinprior_(x) for x in old_state.trace])
    (proposed_state, f(proposed_state), 0.)
end
langevinprior_(s::Sample{Whitened, LangevinPrior}) = begin
    new_val = (1-1/tau)*s.value + sqrt(2/tau)*randn(size(s.value)...)
    Sample(whitened, langevinprior, new_val, white_logpdf(new_val))
end
langevinprior_(i::If) =
    If(i.cond, TraceElement[langevinprior_(x) for x in i.trace])
langevinprior_(s::Sample{UnWhitened, LangevinPrior}) = error()
langevinprior_(el::TraceElement) = el

sample(::LangevinPosterior, f::Function, old_state::State) = begin
    #Propagate differentiated variables.
    old_state = diff_rand(old_state)
    f(old_state)
    #Differentiate logpdf.
    logpdf = sum(map(el -> el.logpdf, old_state.trace))
    ReverseDiffOverload.diff(logpdf)
    #Take a gradient step
    proposed_state = diff_rand(grad_step(old_state))#State(TraceElement[diff_rand(grad_step(x)) for x in old_state.trace])
    proposed_res = value(f(proposed_state))
    #Differentiate logpdf.
    logpdf = sum(map(el -> el.logpdf, proposed_state.trace))
    ReverseDiffOverload.diff(logpdf)
    #println()
    pp = - proposal_p(old_state, proposed_state) + proposal_p(proposed_state, old_state)
    (proposed_state, proposed_res, pp)
end

printlnval(x) = (println(x); x)
diff_rand{W}(s::Sample{W, LangevinPosterior}) =
    Sample(W(), langevinposterior, ReverseDiffOverload.Call(value(s)), nothing)
diff_rand{W, S}(s::Sample{W, S}) = s
diff_rand(si::SI) = map(diff_rand, si)
import Base.map
map(f::Function, state::State) = State(TraceElement[f(x) for x in state.trace])
map(f::Function, i::If) = If(i.cond, TraceElement[f(x) for x in i.trace])

const tau = 100.
proposal_mean{W}(s::Sample{W, LangevinPosterior}) = s.value.val + s.value.dval/tau
#proposal_mean(s::Sample{UnWhitened, LangevinPosterior}) = s.value.val + s.value.dval/tau
#proposal_mean(s::Sample{Whitened, LangevinPosterior}) = (1-1/tau)*s.value.val + s.value.dval/tau
grad_step{W}(s::Sample{W, LangevinPosterior}) = begin
    res = proposal_mean(s) + sqrt(2/tau)*randn(size(s.value.val)...)
    Sample(W(), langevinposterior, res, nothing)
end
grad_step(si::SI) = map(grad_step, si)
grad_step(s::Sample{UnWhitened}) = dc(s)
grad_step(s::Sample{Whitened}) = s

proposal_p(state_from::State, state_to::State) = sum(map(proposal_p, state_from.trace, state_to.trace))
proposal_p{W}(sample_from::Sample{W, LangevinPosterior}, sample_to::Sample{W, LangevinPosterior}) =
    logpdf(Normal, value(sample_to), proposal_mean(sample_from), sqrt(2/tau))
proposal_p{W, V}(::Sample{W, LangevinPosterior}, ::Sample{V, LangevinPosterior}) = error()
proposal_p(sample_from::Sample, sample_to::Sample) = 0.

logpdf_diff(sampler::Sampler, proposed_trace::Vector{TraceElement}, old_trace::Vector{TraceElement}) = begin
    total = 0.
    @assert length(proposed_trace) == length(old_trace)
    for i = 1:length(proposed_trace)
        total += logpdf_diff(sampler, proposed_trace[i], old_trace[i])
    end
    total
end
logpdf_diff(sampler::Sampler) = (xy) -> logpdf_diff(sampler, xy[1], xy[2])
logpdf_diff(sampler::Sampler, proposed_si::State, old_si::State) =
    logpdf_diff(sampler, proposed_si.trace, old_si.trace)
logpdf_diff(sampler::Sampler, proposed_si::If, old_si::If) =
    if proposed_si.cond == old_si.cond
        logpdf_diff(sampler, proposed_si.trace, old_si.trace)
    else
        0.
    end
logpdf_diff(sampler::LangevinPrior, proposed_sample::Sample{UnWhitened}, old_sample::Sample{UnWhitened}) =
    proposed_sample.logpdf - old_sample.logpdf
logpdf_diff(sampler::LangevinPrior, proposed_sample::Sample{Whitened}, old_sample::Sample{Whitened}) = 0.
logpdf_diff(sampler::Sampler, proposed_sample::Sample, old_sample::Sample) =
    proposed_sample.logpdf - old_sample.logpdf


text:
include("create_artificial_dataset.jl")

function main()
    output_filepath = "../../data/trajectories/artificial.h5"
    num_scenarios = 10000
    num_col = 20
    pcol = build_parallel_dataset_collector(
        num_scenarios,
        num_col,
        output_filepath
    )
    generate_dataset(pcol)
    analyze_risk_dataset(output_filepath)
end

@time main()


text:
struct TextTree
    children::Vector{Vector{Int}}
    text::Vector{String}
end

function TextTree(D::DESPOT)
    lenb = length(D.children)
    lenba = length(D.ba_children)
    len = lenb + lenba
    children = Vector{Vector{Int}}(len)
    text = Vector{String}(len)
    for b in 1:lenb
        children[b] = D.children[b] .+ lenb
        text[b] = @sprintf("o:%-5s U:%6.2f, l:%6.2f, μ:%6.2f, l₀:%6.2f, |Φ|:%3d",
                           b==1 ? "<root>" : string(D.obs[b]),
                           D.U[b],
                           D.l[b],
                           D.mu[b],
                           D.l_0[b],
                           length(D.scenarios[b]))
    end
    for ba in 1:lenba
        children[ba+lenb] = D.ba_children[ba]
        text[ba+lenb] = @sprintf("a:%-5s l:%6.2f μ:%6.2f, ρ:%6.2f", D.ba_action[ba], ba_l(D, ba), D.ba_mu[ba], D.ba_rho[ba])
    end
    return TextTree(children, text)
end

struct TreeView
    t::TextTree
    root::Int
    depth::Int
end

TreeView(D::DESPOT, b::Int, depth::Int) = TreeView(TextTree(D), b, depth)

Base.show(io::IO, tv::TreeView) = shownode(io, tv.t, tv.root, tv.depth, "", "")

function shownode(io::IO, t::TextTree, n::Int, depth::Int, item_prefix::String, prefix::String)
    print(io, item_prefix)
    print(io, @sprintf("[%-4d]", n))
    print(io, " $(t.text[n])")
    if depth <= 0
        println(io, " ($(length(t.children[n])) children)")
    else
        println(io)
        if !isempty(t.children[n])
            for c in t.children[n][1:end-1]
                shownode(io, t, c, depth-1, prefix*"├──", prefix*"│  ")
            end
            shownode(io, t, t.children[n][end], depth-1, prefix*"└──", prefix*"   ")
        end
    end
end


text:
using IndependentHypothesisWeighting
using Documenter

DocMeta.setdocmeta!(IndependentHypothesisWeighting, :DocTestSetup, :(using IndependentHypothesisWeighting); recursive=true)

makedocs(;
    modules=[IndependentHypothesisWeighting],
    authors="Nikos Ignatiadis <nikos.ignatiadis01@gmail.com> and contributors",
    repo="https://github.com/nignatiadis/IndependentHypothesisWeighting.jl/blob/{commit}{path}#{line}",
    sitename="IndependentHypothesisWeighting.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nignatiadis.github.io/IndependentHypothesisWeighting.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nignatiadis/IndependentHypothesisWeighting.jl",
    devbranch="main",
)


text:
using Distances

using Test
using LinearAlgebra
using Random
using Statistics

@test isempty(detect_ambiguities(Distances))

include("F64.jl")
include("test_dists.jl")


text:
@testset "@col" begin
    t = @Table(A=[1,2,3],B=[2.0,4.0,6.0])

    @test @col(t, A) == [1,2,3]
    @test @col(t, A, B) == t
end

@testset "@select" begin
    t = @Table(A=[1,2,3],B=[2.0,4.0,6.0])

    @test @select(t, A) == @Table(A=[1,2,3])
    @test @select(t, A, C = B) == @Table(A=[1,2,3],C=[2.0,4.0,6.0])
    @test @select(t, C = (A,B) -> A+B) == @Table(C=[3.0,6.0,9.0])
end

@testset "@filter etc" begin
    t = @Table(A=[1,2,3],B=[2.0,4.0,6.0])

    @test @filter(t, A -> A>1) == @Table(A=[2,3], B=[4.0,6.0])
    @test @filter(t, A -> A>1, B -> B < 4.5) == @Table(A=[2],B=[4.0])

    @test (t2 = copy(t); @filter!(t2, A -> A>1, B -> B < 4.5); t2 == @Table(A=[2],B=[4.0]))

    @test @filter_mask(t, A -> A>1, B -> B < 4.5) == [false,true,false]
end


text:
abstract type AbstractObjective end
include("cores/core.jl")
include("maximin_obj.jl")
include("soyster_maximin_obj.jl")
include("de_jong_maximin_obj.jl")
include("robust_maximin_obj.jl")
include("cc_maximin_obj.jl")
include("minimax_obj.jl")
include("custom_obj.jl")
include("no_obj.jl")


text:
module Foo

struct Data
  name::String
  values::Vector{Float64}
end

function doSomething(data::Vector{Data}, io=STDOUT)
  for d in data
    println(io, uppercase(d.name), " ", BLAS.asum(length(d.values), d.values, 1))
  end
end

function bar(data::Vector{Data})
  isempty(data) && return 1
  doSomething(data)
  return 0
end

function baz(data::Vector{Data}, logPath::String)
  isempty(data) && return 1
  open(logPath, "w") do io
    doSomething(data, io)
  end
  return 0
end

end


text:
function ArrayInterface.ismutable(::Type{HybridArray{S,T,N,M,TData}}) where {S,T,N,M,TData}
    return ArrayInterface.ismutable(TData)
end

function ArrayInterface.can_setindex(::Type{HybridArray{S,T,N,M,TData}}) where {S,T,N,M,TData}
    return ArrayInterface.can_setindex(TData)
end

function ArrayInterface.parent_type(::Type{HybridArray{S,T,N,M,TData}}) where {S,T,N,M,TData}
    return TData
end

function ArrayInterface.restructure(x::HybridArray{S}, y) where {S}
    return HybridArray{S}(reshape(convert(Array, y), size(x)...))
end

function ArrayInterface.strides(x::HybridArray)
    return ArrayInterface.strides(parent(x))
end

@generated function ArrayInterface.strides(x::HybridArray{S,T,N,N,Array{T,N}}) where {S,T,N}
    collected_strides = []
    i = 1
    for (argnum, Sarg) in enumerate(S.parameters)
        if i > 0
            push!(collected_strides, ArrayInterface.StaticInt(i))
        else
            push!(collected_strides, :(datastrides[$argnum]))
        end
        if Sarg isa Integer
            i *= Sarg
        else
            i = -1
        end
    end
    return quote
        datastrides = strides(parent(x))
        return tuple($(collected_strides...))
    end
end

@generated function ArrayInterface.size(x::HybridArray{S}) where {S}
    collected_sizes = []
    for (argnum, Sarg) in enumerate(S.parameters)
        if Sarg isa Integer
            push!(collected_sizes, ArrayInterface.StaticInt(Sarg))
        else
            push!(collected_sizes, :(datasize[$argnum]))
        end
    end
    return quote
        datasize = ArrayInterface.size(parent(x))
        return tuple($(collected_sizes...))
    end
end

ArrayInterface.contiguous_axis(::Type{HybridArray{S,T,N,N,TData}}) where {S,T,N,TData} = ArrayInterface.contiguous_axis(TData)
ArrayInterface.contiguous_batch_size(::Type{HybridArray{S,T,N,N,TData}}) where {S,T,N,TData} = ArrayInterface.contiguous_batch_size(TData)
ArrayInterface.stride_rank(::Type{HybridArray{S,T,N,N,TData}}) where {S,T,N,TData} = ArrayInterface.stride_rank(TData)


text:
__precompile__()

module EarthEngineREST

using NPZ
using Retry
using JSON3
using PyCall
using GeoArrays
using EarthEngine
using GeoDataFrames
using DocStringExtensions
import EarthEngine: Initialize
import GeoArrays: geotransform_to_affine

# define a constant pointer to the ee object
const serviceaccount = PyNULL()
const gauth = PyNULL()

# specify different urls and routings based on api version
const baseurl = "https://earthengine.googleapis.com"
const alphaapi = join([baseurl, "v1alpha"],"/")
const betaapi = join([baseurl, "v1beta"], "/")
const latestapi = betaapi

# create a lookup dictionary of api paths so users can use their preference
# add more as different versions become available
const apis = Dict("alpha" => alphaapi, "beta" => betaapi, "latest" => latestapi)

include("eeresttypes.jl")
include("eevalues.jl")
include("eeimages.jl")
include("eetables.jl")

"""
    _sendrequest(session::EESession, endpoint::AbstractString, data::Dict{Symbol,<:Any}; version::AbstractString = "latest", nretries::Integer = 4)

Private function for sending requests to EarthEngine REST API. Uses exponential backoff to retry requests if failed with non-fatal error.
"""
function _sendrequest(
    session::EESession,
    endpoint::AbstractString,
    data::Dict{Symbol,<:Any};
    version::AbstractString = "latest",
    nretries::Integer = 4
)

    urlparts = [apis[version], "projects", session.project, endpoint]
    url = join(urlparts, "/")

    @repeat nretries try
        response = session.auth.post(url = url, data = JSON3.write(data))

    catch e
        @delay_retry if e.status_code < 200 && e.status_code >= 500
        end
    end
end

"""
    typetodict(x::T)

Function to convert EarthEngineREST types to dictionary key,value pairs to encode
for API requests
"""
function typetodict(x::T)::Dict{Symbol,Any} where {T}
    fields = fieldnames(T)
    outdict = Dict{Symbol,Any}()
    for fn in fields
        subfield = getfield(x, fn)
        fn = String(fn)
        if length(fieldnames(typeof(subfield))) > 0
            if endswith(String(fn), "of")
                fn = replace(fn, "of" => "")
            end

            outdict[Symbol(fn)] = typetodict(subfield)
        else
            if endswith(fn, "of")
                fn = replace(fn, "of" => "")
            end

            outdict[Symbol(fn)] = subfield
        end
    end
    return outdict
end

"""
    extract_geotransform(x::AffineTransform)

Function to extract geotransform vector from AffineTransformation type
"""
function extract_geotransform(x::AffineTransform)
    af = typetodict(x)
    gt = [
        af[:translateX],
        af[:scaleX],
        af[:shearX],
        af[:translateY],
        af[:shearY],
        af[:scaleY],
    ]
end

"""
    extract_geotransform(x::PixelGrid)

Function to extract geotransform vector from PixelGrid type
"""
function extract_geotransform(x::PixelGrid)
    extract_geotransform(x.affineTransform)
end

"""
    extract_bbox(x::AffineTransform, griddims::GridDimensions)

Function to extract bounding box vector from AffineTransform type. Vector will be [W, S, E, N]
"""
function extract_bbox(x::AffineTransform, griddims::GridDimensions)
    geo_t = extract_geotransform(x)
    x_size = griddims.width
    y_size = griddims.height

    xmin = min(geo_t[1], geo_t[1] + x_size * geo_t[2])
    xmax = max(geo_t[1], geo_t[1] + x_size * geo_t[2])
    ymin = min(geo_t[4], geo_t[4] + y_size * geo_t[6])
    ymax = max(geo_t[4], geo_t[4] + y_size * geo_t[6])

    return [xmin, ymin, xmax, ymax]
end

"""
    extract_bbox(x::PixelGrid)

Function to extract bounding box vector from PixelGrid type. Vector will be [W, S, E, N]
"""
function extract_bbox(x::PixelGrid)
    extract_bbox(x.affineTransform,x.dimensions)
end

"""
    extract_affinemap(x::AffineTransform)

Function to extract AffineMap type from AffineTransform type.
"""
function extract_affinemap(x::AffineTransform)
    gt = extract_geotransform(x)
    return geotransform_to_affine(gt)
end

"""
    extract_affinemap(x::PixelGrid)

Function to extract AffineMap type from PixelGrid type.
"""
function extract_affinemap(x::PixelGrid)
    extract_affinemap(x.affineTransform)
end

function extract_lonslats(x::AffineTransform)
    gt = extract_geotransform()

end

"""
    extract_gridcoordinates(x::PixelGrid)

Function to convert PixelGrid to vectors of lat lon coordinates. Returns
a tuple of (lon::Vector{Float64}, lat::Vector{Float64}).
"""
function extract_gridcoordinates(x::PixelGrid)
    bbox = extract_bbox(x)
    gt = extract_geotransform(x)

    lons = collect(bbox[1]:abs(gt[2]):bbox[3])[1:end-1]
    lats = collect(bbox[2]:abs(gt[6]):bbox[4])[2:end]

    return lons, lats
end

"""
    Initialize(session::EESession)

Extends the `Initialize()` function to take an authenticated Earth Engine session
and use those credentials.
"""
function Initialize(session::EESession)
    EarthEngine.Initialize(session.auth.credentials)
end

export Initialize, EESession, PixelGrid, GridDimensions, AffineTransform, computepixels, computetable, computevalue

end # module


text:
"""
    Codecov(; config_file::Union{AbstractString, Nothing}=nothing) -> Codecov

Add `Codecov` to a template's plugins to optionally add a `.codecov.yml` configuration file
to generated repositories, and an appropriate badge to the README. Also updates the
`.gitignore` accordingly.

# Keyword Arguments:
* `config_file::Union{AbstractString, Nothing}=nothing`: Path to a custom `.codecov.yml`.
  If left unset, no file will be generated.
"""
struct Codecov <: GenericPlugin
    gitignore::Vector{String}
    src::Union{String, Nothing}
    dest::String
    badges::Vector{Badge}
    view::Dict{String, Any}

    function Codecov(; config_file::Union{AbstractString, Nothing}=nothing)
        if config_file != nothing
            config_file = if isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError("File $(abspath(config_file)) does not exist"))
            end
        end
        new(
            ["*.jl.cov", "*.jl.*.cov", "*.jl.mem"],
            config_file,
            ".codecov.yml",
            [
                Badge(
                    "Codecov",
                    "https://codecov.io/gh/{{USER}}/{{PKGNAME}}.jl/branch/master/graph/badge.svg",
                    "https://codecov.io/gh/{{USER}}/{{PKGNAME}}.jl",
                ),
            ],
            Dict{String, Any}(),
        )
    end
end
Base.@deprecate_binding CodeCov Codecov

interactive(::Type{Codecov}) = interactive(Codecov; file=nothing)


text:
"""
"""
struct FockSparseState{T,N}
    clusters::Vector{Cluster}
    data::OrderedDict{FockConfig{N}, T }
end

#const ClusteredState2{T,N} = FockSparseState{OrderedDict{ClusterConfig{N}, T}, N}
#const TuckerState2{T,N} = FockSparseState{OrderedDict{TuckerConfig{N}, Array{T}}, N}
#const CompressedState2{T,N} = FockSparseState{OrderedDict{TuckerConfig{N}, Tucker{T}}, N}
#struct ClusteredState2{T,N}
#    clusters::Vector{Cluster}
#    data::FockSparseState{OrderedDict{ClusterConfig{N}, T}}
#end

"""
Glue different Sparse Vector State types together
"""
abstract type AbstractState end

"""
    add_fockconfig!(s::AbstractState, fock::Vector{Tuple{T,T}}) where T<:Integer
"""
function add_fockconfig!(s::AbstractState, fock::Vector{Tuple{T,T}}) where T<:Integer
    add_fockconfig!(s,FockConfig(fock))
end
"""
    setindex!(s::AbstractState, a::OrderedDict, b)
"""
function Base.setindex!(s::AbstractState, a, b)
    s.data[b] = a
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using Glob

"""
    catexamples(output, pattern)

Generate a module with examples  on a common theme.
"""
function catexamples(output, pattern)
    n,ext = splitext(output)
    if ext == ".jl"
        outputname = n
    else
        outputname = output
        output = output * ".jl"
    end
    List =  Glob.glob(pattern * ".jl")
    Functions = []
    open(output, "w") do o
        @printf o "module %s\n" outputname
        for ex1 in List
            n,ext = splitext(ex1)
            if ext == ".jl" && n != outputname
                println("Example $ex1 in $(pwd())\n")
                open(ex1, "r") do f
                    text = read(f, String)
                    @printf o "\n"
                    @printf o "function %s()\n" n
                    @printf o "%s\n" text
                    @printf o "end # %s\n\n " n
                end
                push!(Functions, n)
            end
        end
        @printf o "function allrun()\n"
        for foo in Functions
            @printf o """println("#####################################################") \n"""
            @printf o """println("# %s ")\n""" foo
            @printf o "    %s()\n" foo
        end
        @printf o "    return true\n"
        @printf o "end # function allrun\n\n"
        @printf o "end # module %s\n" outputname
    end
    println("Wrote $output in $(pwd())\n")
end


text:
function mdl_freq_estimator_ode!(device_states,
        output_ode,
        f0,
        device::DynInverter{C,O,VC,DC,PLL,F}) where {C <: Converter,
                                                   O <: OuterControl,
                                                   VC<: VSControl,
                                                   DC<: DCSource,
                                                   F <: Filter}

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.freq_estimator]
    vod = device_states[external_ix[1]]
    voq = device_states[external_ix[2]]
    δθ_vsm = device_states[external_ix[3]]

    #Obtain inner variables for component
       #vod = device.inner_vars[Vdo_var]
       #voq = device.inner_vars[Vqo_var]
    #δθ_vsm = device.inner_vars[δdqRI_var]

    #Get parameters
      ω_lp = device.freq_estimator.ω_lp
    kp_pll = device.freq_estimator.kp_pll
    ki_pll = device.freq_estimator.ki_pll
        ωb = 2.0*pi*f0
        ωg = 1.0 #TODO:  create getter later

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.freq_estimator]

    #Define internal states for frequency estimator
    internal_states = @view device_states[local_ix]
    vpll_d = internal_states[1]
    vpll_q = internal_states[2]
     ϵ_pll = internal_states[3]
    δθ_pll = internal_states[4]

    #Inputs (control signals)

    #Compute 6 states ODEs (D'Arco EPSR122 Model)
    #Output Voltage LPF (internal state)
    #𝜕vpll_d/𝜕t, D'Arco ESPR122 eqn. 12
    output_ode[local_ix[1]] = ( ω_lp*vod*cos(δθ_pll-δθ_vsm)
                              + ω_lp*voq*sin(δθ_pll-δθ_vsm)
                              - ω_lp*vpll_d)
    #𝜕vpll_q/𝜕t, D'Arco ESPR122 eqn. 12
    output_ode[local_ix[2]] = (- ω_lp*vod*sin(δθ_pll-δθ_vsm)
                               + ω_lp*voq*cos(δθ_pll-δθ_vsm)
                               - ω_lp*vpll_q)
    #PI Integrator (internal state)
    #𝜕dϵ_pll/𝜕t, D'Arco ESPR122 eqn. 13
    output_ode[local_ix[3]] = atan(vpll_q/vpll_d)
    #PLL Frequency Deviation (internal state)
    #𝜕δθ_pll/𝜕t, D'Arco ESPR122 eqn. 15
    output_ode[local_ix[4]] = ( ωb*kp_pll*atan(vpll_q/vpll_d)
                              + ωb*ki_pll*ϵ_pll)

    #Update inner_vars
    #PLL frequency, D'Arco EPSR122 eqn. 16
    device.inner_vars[ω_freq_estimator_var] = ( kp_pll*atan(vpll_q/vpll_d)
                                              + ki_pll*ϵ_pll
                                              + 1.0 )
        #TODO: replace 1.0 w/ grid getter, or not needed if PCC?
end


text:
using Libdl
using PackageCompiler

include("build_IntelOpenMP.jl")
include("build_MKL.jl")
include("deps.jl")

# Read the build_h.jl file
fname = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia", "base", "build_h.jl")
lines = open(fname) do f
    readlines(f, keep = true)
end

# Write path to MKL
lines = map(lines) do l
    if occursin(r"libblas_name", l)
        return "const libblas_name = $(repr(libmkl_rt))\n"
    elseif occursin(r"liblapack_name", l)
        return "const liblapack_name = $(repr(libmkl_rt))\n"
    else
        return l
    end
end

# Write the modified lines to the file
open(fname, "w") do f
    for l in lines
        write(f, l)
    end
end

sysimgpath = PackageCompiler.sysimgbackup_folder("native")
if ispath(sysimgpath)
    rm(sysimgpath, recursive=true)
end

force_native_image!()


text:
#testbug_GB13
include("GB13.jl")
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

rng = Random.seed!(2019)

gbdt = GB13.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                 # Learning rate
  num_iterations = 100,                # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB13.GBLearner(gbdt, :regression)

x, y = GB13.test_data(500; p = 10, sigma = sqrt(10))
x_test, y_test = GB13.test_data(500; p = 10, sigma = sqrt(10.0))
mu = GB13.fried(x)
mu_test = GB13.fried(x_test)

# initialize the hypers
hypers = GB13.Hypers(x, y, 1, 1, 0.8)
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
hypers.α = 0.95
hypers.s = hypers.α .* ones(size(x, 2)) / size(x, 2)
hypers

# initialize opt
opts = GB13.Opts()
opts.num_burn = 400 # 500 in total now

for l in collect(1 : opts.num_burn)

    hypers_copy = deepcopy(hypers)

    # fit the model
    model = GB13.fit!(gbl, x, y, hypers_copy)

    GB13.UpdateS!(hypers_copy) # returns the updated hypers_copy.s
    if (opts.update_s) # true
        hypers.s = hypers_copy.s
    end
    if (opts.update_varcount) # that is, the var count is going to be set to zero again
        hypers.varcount = hypers_copy.varcount
    end

end

hypers

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
# simple iid bootstrap
function bootstrap(data)
    n = size(data,1)
    resampled = similar(data)
    for i = 1:n
        j = rand(1:n)
        resampled[i,:] = data[j,:]
    end
    return resampled
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# ---
# title: 791. Custom Sort String
# id: problem791
# author: Indigo
# date: 2021-01-29
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/custom-sort-string/description/>
# hidden: true
# ---
#
# `S` and `T` are strings composed of lowercase letters. In `S`, no letter
# occurs more than once.
#
# `S` was sorted in some custom order previously. We want to permute the
# characters of `T` so that they match the order that `S` was sorted. More
# specifically, if `x` occurs before `y` in `S`, then `x` should occur before
# `y` in the returned string.
#
# Return any permutation of `T` (as a string) that satisfies this property.
#
#
#
#     **Example :**
#     Input:
#     S = "cba"
#     T = "abcd"
#     Output: "cbad"
#     Explanation:
#     "a", "b", "c" appear in S, so the order of "a", "b", "c" should be "c", "b", and "a".
#     Since "d" does not appear in S, it can be at any position in T. "dcba", "cdba", "cbda" are also valid outputs.
#
#
#
#
# **Note:**
#
#   * `S` has length at most `26`, and no character is repeated in `S`.
#   * `T` has length at most `200`.
#   * `S` and `T` consist of lowercase letters only.
#
#
## @lc code=start
using LeetCode

function custom_sort_string(S::String, T::String)
    letter_rank = fill(0, 128)
    for (idx, ch) in enumerate(codeunits(S))
        letter_rank[ch] = idx
    end
    u = codeunits(T)[:]
    sort!(u, by = x -> letter_rank[x])
    String(u)
end
## @lc code=end


text:
module A

export AType

struct AType
	x::Int
end

end


text:
import Pkg; Pkg.activate("../../")

println("Loading libraries ...")
include("PreProcess.jl")
import CSV
import DataFrames; const DF = DataFrames
println("Finished loading libraries.")

data_dir = "../data"
path_to_data = "$(data_dir)/cb_transformed.csv"

function separatesamples(Y)
  sample_ids = sort(unique(Y.sample_id))
  return Matrix.(DF.select(Y[Y.sample_id .== i, :], DF.Not(:sample_id))
                 for i in sample_ids)
end

# Read transformed CB data into Y as a data frame.
# Replace `missing` with `NaN`.
Y = coalesce.(CSV.read(path_to_data), NaN)
markernames = filter(x -> x != :sample_id, names(Y))
y = separatesamples(Y)
isgoodmarker = PreProcess.preprocess!(y, maxNanOrNegProp=.9,
                                      maxPosProp=.9, rowThresh=-6.0)

includedmarkers = Symbol[]
excludedmarkers = Symbol[]

for j in 1:length(isgoodmarker)
  marker = markernames[j]
  if isgoodmarker[j]
    append!(includedmarkers, [marker])
  else
    append!(excludedmarkers, [marker])
  end
end

# Concatenate y and convert to data frame
Y_reduced = DF.DataFrame(vcat(y...), includedmarkers)

# Convert NaN to missing
Y_reduced = ifelse.(isnan.(Y_reduced), missing, Y_reduced)

# Append a column for sample_id
Y_reduced[!, :sample_id] = vcat([fill(i, size(y[i], 1))
                                 for i in 1:length(y)]...)

# Path to reduced data
path_to_reduced_data = "$(data_dir)/cb_transformed_reduced.csv"

# Write reduced data to disk.
println("Writing reduced-and-transformed CB data to: ")
println(path_to_reduced_data)

CSV.write(path_to_reduced_data, Y_reduced)

println("Done!")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
"""
    trunk_build(trunk,M)

Generate a dictionary of trunk functions evaluated at ``t = 0``

"""
function trunk_build(trunk,M)
    return Dict(i => function (x) return trunk(vcat(0,x))[i] end for i in 1:M)
end

"""
    trunk_ortho_build(utilde,L1,L2,nodes,weights)

Generate a dictionary of orthonormal custom basis functions which use an orthogonal polynomial expansion for evaluation away from the specified quadrature nodes.

"""
function trunk_ortho_build(utilde,L1,L2,nodes,weights)
    trunk_ortho = [];
    W = diagm(0 => weights);
    N = size(utilde,2);
    for i in 1:N
        coeffs = zeros(size(utilde,1));
        coeffs[1] = legendre_norm.(nodes,L1,L2,0)'*W*utilde[:,i];
        for j in 1:(size(utilde,1)-1)
            coeffs[j+1] = legendre_norm.(nodes,L1,L2,j)'*W*utilde[:,i]
        end
        push!(trunk_ortho, (x) -> coeffs'*legendre_norm_collect.(x,L1,L2,(size(utilde,1)-1)));
    end
    return trunk_ortho
end

"""
    build_basis(trunk,L1,L2,M,nodes,weights)

Generate a dictionary of orthonormal custom basis functions based on a specified quadrature rule. An orthogonal polynomial expansion is utilized for the evaluation at points away from the specified quadrature nodes.

"""
function build_basis(trunk,L1,L2,M,nodes,weights)
    trunk_func = trunk_build(trunk,M);

    W = diagm(0 => weights);

    A = zeros(size(nodes,1),M);
    for i in 1:size(nodes,1)
        for k in 1:M
            A[i,k] = trunk_func[k](nodes[i]);
        end
    end

    B = W^(1/2)*A;
    F = svd(B);
    utilde = W^(-1/2)*F.U;

    ortho_trunk_func = trunk_ortho_build(utilde,L1,L2,nodes,weights);

    # Check if the functions are orthonormal to a specified precision at the specified nodes
    orthonormal_check(basis_eval(ortho_trunk_func,nodes),weights;tol=1e-12)

    return ortho_trunk_func, F.S
end

"""
    basis_eval(basis,nodes)

Evaluate the basis functions for a specified spatial grid or at a specified point. Evaluation for a specified grid returns an ``N x M`` matrix with ``N`` respresenting the number of nodes and ``M`` representing the number of custom basis functions. Evalaution for a specified point returns a ``1 x M`` vector.

"""
function basis_eval(basis,nodes::Array)
    apply(f,x) = f(x);
    return permutedims(apply.(basis,permutedims(nodes)));
end
function basis_eval(basis,node)
    apply(f,x) = f(x);
    return permutedims(apply.(basis,node));
end

"""
    expansion_coefficients(basis,fnc,nodes,weights)

Compute the expansion coefficients...

"""
function expansion_coefficients(basis::Array{Any,1},fnc::Function,nodes,weights)
    return basis_eval(basis,nodes)'*diagm(0 => weights)*fnc.(nodes)
end
function expansion_coefficients(basis::Array{Any,1},fnc,nodes,weights)
    return basis_eval(basis,nodes)'*diagm(0 => weights)*fnc
end
function expansion_coefficients(basis::Array{Float64,2},fnc::Function,nodes,weights)
    return basis'*diagm(0 => weights)*fnc.(nodes)
end
function expansion_coefficients(basis::Array{Float64,2},fnc,nodes,weights)
    return basis'*diagm(0 => weights)*fnc
end

"""
    expansion_approximation(basis,coefficients,nodes)

Compute the expansion approximation...

"""
function expansion_approximation(basis::Array{Any,1},coefficients,nodes)
    return basis_eval(basis,nodes)*coefficients
end
function expansion_approximation(basis::Array{Float64,2},coefficients,nodes)
    return basis*coefficients
end

"""
    basis_derivative(basis,nodes)

Compute the derivatives of each basis function using auto differentiation.

"""
function basis_derivative(basis)
  Dbasis = []
  for i in 1:length(basis);
    push!(Dbasis,(x)->ForwardDiff.derivative.(basis[i],x))
  end
  return Dbasis
end

"""
    save_basis(basis,M,pde_function)

"""
function save_basis(basis,M,pde_function)
    @save @sprintf("basis_%i_M_%s.bson",M,pde_function) basis
end

"""
    load_basis(basis,M,pde_function)

"""
function load_basis(M,pde_function)
    @load @sprintf("basis_%i_M_%s.bson",M,pde_function) basis
    return basis
end


text:
using MetaGen
using Gen

#possible_objects = ["person","bicycle","car","motorcycle","airplane"]
possible_objects = [1, 2, 3, 4, 5]
#call it

#include("useful_functions.jl")

dict = @pipe "../Data/data_first_10_detections.json" |> open |> read |> String |> JSON.parse

#try to make objects_observed::Array{Array{Array{Array{Detection2D}}}} of observed objects.
#outer array is for scenes, then frames, the receptive fields, then last is an array of detections

################################################################################
receptive_fields = make_receptive_fields()
objects_observed, camera_trajectories = make_observations(dict, receptive_fields)

#constrain
obs = Gen.choicemap()
for v=1:10
    for f=1:6
        obs[:videos => v => :frame_chain => f => :camera => :camera_location_x] = camera_trajectories[v,f].camera_location.x
    end
end

gt_trace,log_score = Gen.generate(metacog, (possible_objects, 10, 6), obs)
gt_choices = get_choices(gt_trace)

#@show gt_choices
@show log_score

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using Random: AbstractRNG
using Bijectors
import Bijectors: logabsdetjac, bijector
using Distributions: MultivariateDistribution, Continuous
import Distributions: _rand!, _logpdf
using StatsFuns

"""
    Kipping13()

A non-informative prior for two-parameter limb-darkening coefficients using *triangular sampling* ([Kipping 2013](https://ui.adsabs.harvard.edu/abs/2013MNRAS.435.2152K/)).

# Examples

```jldoctest
julia> using Random; rng = Random.seed!(10);

julia> rand(rng, Kipping13())
2-element Vector{Float64}:
 0.24716310305467298
 0.08836997249298882

julia> rand(rng, Kipping13(), 5)
2×5 Matrix{Float64}:
 0.0664907  0.124817   1.00732    0.806902  0.74165
 0.520411   0.222718  -0.389412  -0.314755  0.0768429
```

# References

> [Kipping (2013)](https://ui.adsabs.harvard.edu/abs/2013MNRAS.435.2152K/)
>
>   "Efficient, uninformative sampling of limb darkening coefficients for two-parameter laws"
"""
struct Kipping13 <: MultivariateDistribution{Continuous} end

Base.length(::Kipping13) = 2

function _rand!(rng::AbstractRNG, ::Kipping13, x::AbstractVector{T}) where T
    q1, q2 = rand(rng, T, 2)
    sqrtq1 = sqrt(q1)
    twoq2 = 2 * q2
    x[begin] = sqrtq1 * twoq2
    x[end] = sqrtq1 * (1 - twoq2)
    return x
end

_logpdf(::Kipping13, x::AbstractArray{T}) where {T} = zero(T)

struct Kipping13Transform <: Bijector{1} end

function (::Kipping13Transform)(x::AbstractVector)
    usum = sum(x)
    q = [usum^2, 0.5 * first(x) / usum]
    return @. log(q) - log(1 - q)
end

function (::Inverse{<:Kipping13Transform})(y::AbstractVector)
    tmp = map(logistic, y)
    sqrtq1 = sqrt(first(tmp))
    twoq2 = 2 * last(tmp)
    tmp[begin] = sqrtq1 * twoq2
    tmp[end] = sqrtq1 * (1 - twoq2)
    return tmp
end

logabsdetjac(::Kipping13Transform, y::Number) = -2 * softplus(-y) - y
logabsdetjac(k::Kipping13Transform, y) = sum((yi) -> logabsdetjac(k, yi), y)
bijector(::Kipping13) = Kipping13Transform()


text:
#=

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

=#
using QuantEcon: LAE, lae_est
using Distributions

s = 0.2
delta = 0.1
a_sigma = 0.4  # A = exp(B) where B ~ N(0, a_sigma)
alpha = 0.4  # We set f(k) = k**alpha
psi_0 = Beta(5.0, 5.0)  # Initial distribution
phi = LogNormal(0.0, a_sigma)

function p(x, y)
    #=
    Stochastic kernel for the growth model with Cobb-Douglas production.
    Both x and y must be strictly positive.
    =#
    d = s * x.^alpha

    # scipy silently evaluates the pdf of the lognormal dist at a negative
    # value as zero. It should be undefined and Julia recognizes this.
    pdf_arg = clamp.((y .- (1-delta) .* x) ./ d, eps(), Inf)
    return pdf(phi, pdf_arg) ./ d
end

n = 10000  # Number of observations at each date t
T = 30  # Compute density of k_t at 1,...,T+1

# Generate matrix s.t. t-th column is n observations of k_t
k = Array{Float64}(n, T)
A = rand!(phi, Array{Float64}(n, T))

# Draw first column from initial distribution
k[:, 1] = rand(psi_0, n) ./ 2  # divide by 2 to match scale=0.5 in py version
for t=1:T-1
    k[:, t+1] = s*A[:, t] .* k[:, t].^alpha + (1-delta) .* k[:, t]
end

# Generate T instances of LAE using this data, one for each date t
laes = [LAE(p, k[:, t]) for t=T:-1:1]

# Plot
ygrid = linspace(0.01, 4.0, 200)
laes_plot = []
colors = []
for i = 1:T
    psi = laes[i]
    push!(laes_plot, lae_est(psi, ygrid))
    push!(colors,  RGBA(0, 0, 0, 1 - (i - 1)/T))
end
plot(ygrid, laes_plot, color=reshape(colors, 1, length(colors)), lw=2, xlabel="capital", legend=:none)
t=LaTeXString("Density of \$k_1\$ (lighter) to \$k_T\$ (darker) for \$T=$T\$")
plot!(title=t)


text:
"""
Utility data structures and functions
"""
module MCUtility
include("stopwatch.jl")
export StopWatch, check

include("color.jl")
export black, red, green, yellow, blue, magenta, cyan, white

export progressBar
"""
    progressBar(step, total)

Return string of progressBar (step/total*100%)
"""
function progressBar(step, total)
    barWidth = 70
    percent = round(step / total * 100.0, digits=2)
    str = "["
    pos = barWidth * percent / 100.0
    for i = 1:barWidth
        if i <= pos
            str *= "█"
        else
            str *= " "
        end
    end
    str *= "] $step/$total=$percent%"
    return str
end

end


text:
module DualNumbers

using SpecialFunctions
import NaNMath
import Calculus

include("dual.jl")

export
    Dual,
    Dual128,
    Dual64,
    Dual32,
    DualComplex256,
    DualComplex128,
    DualComplex64,
    dual,
    epsilon,
    realpart,
    dualpart,
    isdual,
    dual_show,
    conjdual,
    absdual,
    abs2dual,
    ɛ,
    imɛ

end # module


text:
using DrWatson
@quickactivate "EPExperiments"

using EPExperiments
using Random
using Logging
using LoggingExtras
using JLD2
using FileIO
using CSV
using DataFrames

const LOG_DIR_NAME = "radon_nuts_baseline"
const LOG_FILENAME = "out.log"

@with_kw struct RadonNUTSBaseline
    experiment_name::String = "test"
    num_samples::Int = 100
    num_chains::Int = 4
    nuts_acceptance_ratio::Float64 = 0.65
    seed::Int = 1234
    data_file::String = "radon.csv"
    radon_model::Expectation = radon_unbounded
    cost_fn_name::Symbol = :radon_hierarchical
end

function radon_nuts_baseline(experiment_config)
    @unpack seed, num_samples, nuts_acceptance_ratio, num_chains = experiment_config
    @unpack radon_model, cost_fn_name, data_file = experiment_config

    result_folder = mkpath(joinpath(
        datadir(LOG_DIR_NAME),
        savename(experiment_config)
    ))

    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    Random.seed!(seed)

    radon_df = DataFrame(CSV.File(datadir("exp_raw", data_file)))
    log_radon, county_idx, floor, num_counties = preprocess_radon_df(
        radon_df, data_file
    )

    radon_exp = radon_model(
        log_radon,
        county_idx,
        floor,
        num_counties;
        cost_fn=RADON_COST_FNS[cost_fn_name]
    )
    radon_turing = radon_exp.gamma2

    @time begin
    nuts_samples = sample(
        radon_turing,
        Turing.NUTS(nuts_acceptance_ratio),
        MCMCThreads(),
        num_samples,
        num_chains;
        discard_adapt=false
    )
    end

    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "nuts_samples" => nuts_samples,
        "experiment_config" => experiment_config
    ))

    with_logger(logger) do
        @info "MCMC samples:" nuts_samples
    end
end

radon_nuts_baseline(RadonNUTSBaseline(
    experiment_name="collect_baseline_samples",
    num_samples=10000,
    num_chains=16,
    radon_model=radon_hierarchical,
    seed=200
))


text:
@testset "Spaces" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 2, "a1")
  OK = maximal_order(K)
  p = 3*OK
  Kt, t = K["t"]

  infp = infinite_places(K)

  E, b = NumberField(t^2 + 3)
  s = involution(E)

  F = GF(3)

  Hecke.change_base_ring(::Hecke.NfRel, ::Hecke.gfp_mat) = error("asd")
  @test_throws ErrorException hermitian_space(E, F[1 2; 2 1])

  Hecke.change_base_ring(::Hecke.NfRel, x::Hecke.gfp_mat) = x
  @test_throws ErrorException hermitian_space(E, F[1 2; 2 1])

  V = @inferred hermitian_space(E, 3)
  @test fixed_field(V) == base_field(E)
  @test sprint(show, V) isa String
  @test !isquadratic(V)
  @test ispositive_definite(V)
  @test gram_matrix(V) == identity_matrix(E,3)
  @test V !== hermitian_space(E, 3, cached = false)
  v = matrix(E,1,3,[1,1,1])
  @test inner_product(V, v, v) == matrix(E,1,1,[3])

  @test_throws ArgumentError hermitian_space(E, E[1 b; b 1])
  @test_throws ArgumentError hermitian_space(E, FlintQQ[1 0; 1 1])

  V = @inferred hermitian_space(E, FlintQQ[1 2; 2 1])
  @test V === hermitian_space(E, FlintQQ[1 2; 2 1])
  @test V !== hermitian_space(E, FlintQQ[1 2; 2 1], cached = false)
  @test ishermitian(V)
  @test !isdefinite(V)
  @test involution(V) == s
  @test det(V) == -discriminant(V)

  U = hermitian_space(E, E[a 1; 1 1])
  V = hermitian_space(E, E[0 1; 1 0])
  W = hermitian_space(E, E[0 1 0; 1 0 0; 0 0 1])
  @test !Hecke.islocally_hyperbolic(U,p)
  @test !Hecke.isisotropic(U,p)
  @test Hecke.islocally_hyperbolic(V,p)
  @test Hecke.isisotropic(V,p)
  @test all(v -> Hecke.isisotropic(V, v), infp)
  @test !Hecke.islocally_hyperbolic(W,p)
  @test Hecke.isisotropic(W,p)
  @test_throws AssertionError Hecke.islocally_hyperbolic(V, 2*OK)

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2-3
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(1)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(a), 0, 0, 0, (2 * a), 0, 0, 0, (a + 5)])
  V = hermitian_space(E, D)
  DD = matrix(E, 2, 2, [a + 3, 0, 0, (a + 2)])
  W = hermitian_space(E, DD)
  @test @inferred Hecke.isrepresented_by(W, V)
  @test @inferred Hecke.isrepresented_by(V, V)
  DD = E(a) * identity_matrix(E, 5)
  W = hermitian_space(E, DD)
  @test @inferred Hecke.isrepresented_by(V, W)

  # they are not isometric at some dyadic prime ideal:
  DD = identity_matrix(E, 4)
  W = hermitian_space(E, DD)
  @test !(@inferred Hecke.isrepresented_by(W, V))

  K, a = rationals_as_number_field()
  OK = maximal_order(K)
  Kt, t = K["t"]

  E, b = NumberField(t^2+17)

  p = 2*OK
  q = 17*OK

  V = hermitian_space(E, E[102 b; -b 0])
  H = hermitian_space(E, E[0 1; 1 0])
  W = hermitian_space(E, E[1 1 2; 1 2 3; 2 3 1])

  @test_throws ErrorException hasse_invariant(V, p)
  @test_throws ErrorException witt_invariant(W, q)

  for r in [p,q]
    @test Hecke.islocally_represented_by(V,H,r)
  end
  @test Hecke.isrepresented_by(V,H)
  @test !Hecke.islocally_represented_by(W,V,p)
  @test Hecke.islocally_represented_by(V, W, p)
  @test !Hecke.isrepresented_by(W,V)

  @test !isisometric(V, W, ZZ(2))
  @test all(p -> isisometric(V, H, ZZ(p)), PrimesSet(1, 20))

  Qx, x = QQ["x"]
  K, a = number_field(x^3+3, "a")
  Kt, t = K["t"]
  E, b = number_field(t^2-2, "b")

  U = hermitian_space(E, E[1 1;1 1])
  V1 = hermitian_space(E, E[102 b; -b 0])
  V2 = hermitian_space(E, E[1 0;0 1])
  H = hermitian_space(E, E[0 1; 1 0])
  W = hermitian_space(E, E[1 1 2; 1 2 3; 2 3 1])

  infp = infinite_places(K)

  @test_throws ArgumentError isisometric(U, V)  # U is not regular
  @test isisometric(V1, V1, infp[1])
  @test !isisometric(V1, W, infp[2]) # infp[2] is complex but V and W don't have same rank
  @test count(v -> isisometric(V1, H, v), infp) == 1 # V1 and H are not everywhere locally isometric
  @test !isisometric(V1, H)
  @test isisometric(V1, V2)

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1d36ffb4-47cf-479c-add4-c56ef1b1a860
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 200b3ca5-0648-489d-b10d-62a851cfdb4a
begin
	using Plots
	plotly()
end

# ╔═╡ 2042164a-c8f4-48ba-82e7-330205b1d942
md"
# 找一個點，取代一個群
"

# ╔═╡ aecba716-d27b-11eb-10e0-f529b2367d05
data = vcat(randn(400, 2).-0.5, 2randn(201, 2).+1);

# ╔═╡ 51d14132-f210-4982-87b0-9626bd185ae4
begin
	scatter(data[:, 1], data[:, 2], label="data")

	mean = [sum(data[:, 1])/601, sum(data[:, 2])/601]
	scatter!([mean[1]], [mean[2]], label="mean", markersize=10)

	medoid_i = argmin(
		[sum(sum(abs.(data[j, :]-data[i, :])) for j in 1:601) for i in 601]
	)
	scatter!([data[1, medoid_i]], [data[2, medoid_i]], label="medoid", markersize=10)

	median = [sort(data[:, 1])[301], sort(data[:, 2])[301]]
	scatter!([median[1]], [median[2]], label="median", markersize=10)

	mode = [-0.5, -0.5]
	scatter!([mode[1]], [mode[2]], label="mode", markersize=10)
end

# ╔═╡ fc7216ea-86c2-4822-8ceb-a2f9d332dd48
md"
## 中心 (medoid)：

定義 'swapping cust' = $J = argmin(\sum_{j=1}^m \sum_{j=1}^n d(x_i, c_j))$

中心所在之位置滿足與其他點之距離和最小

即在座標空間中，所有點繞著他轉
"

# ╔═╡ 28fdfd71-6505-4c48-b6df-506eb6385572
md"
## 中位 (median)：

中位數落在資料空間之正中間
"

# ╔═╡ a1ef3f18-d27a-416d-97af-a7efa5299d94
md"
## 眾數 (mode)：

眾數落在座標空間中資料密度最大的地方
"

# ╔═╡ bd8d7cae-107a-4181-a2d5-5ad3815b4995
md"---"

# ╔═╡ d5e47835-0f8a-4dbd-9ac0-904868a0d818
md"
## K-Medoid

1. 在有限步數內收斂
2. O(n²)
3. 初始值不能亂給
4. 因初始值的不確定性而不保證收斂到全域最小值
5. 對於類別型變數無法定義 1-norm
"

# ╔═╡ c0417983-bb49-4797-aa62-29be1db7d4d8
md"
## K-Mode

1. 可處理類別型變數 (K-Means 也不行)
2. 使用 cosine similarity 來構造距離
"

# ╔═╡ f7dd3999-9e32-498a-bd3f-35c2526e452f
md"
## K-Prototypes

1. 可處理類別型變數以及數執行變數
2. 使用歐式距離以及 cosine similarity 的線性組合來構造距離
"

# ╔═╡ Cell order:
# ╟─2042164a-c8f4-48ba-82e7-330205b1d942
# ╟─1d36ffb4-47cf-479c-add4-c56ef1b1a860
# ╠═200b3ca5-0648-489d-b10d-62a851cfdb4a
# ╠═aecba716-d27b-11eb-10e0-f529b2367d05
# ╟─51d14132-f210-4982-87b0-9626bd185ae4
# ╟─fc7216ea-86c2-4822-8ceb-a2f9d332dd48
# ╟─28fdfd71-6505-4c48-b6df-506eb6385572
# ╟─a1ef3f18-d27a-416d-97af-a7efa5299d94
# ╟─bd8d7cae-107a-4181-a2d5-5ad3815b4995
# ╟─d5e47835-0f8a-4dbd-9ac0-904868a0d818
# ╟─c0417983-bb49-4797-aa62-29be1db7d4d8
# ╟─f7dd3999-9e32-498a-bd3f-35c2526e452f


text:
# module ShaleDrillingLikelihood_CounterfactualObjects_Test
using Revise

using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters
using Distributed
using Test

using ShaleDrillingLikelihood: drill, _x, _D, LeaseCounterfactual,
    ObservationGroup, InitialDrilling, ichars, zchars, uniti,
    end_ex0, _data, _i, simulate_lease!, Pprime,
    state_if_never_drilled, tstart, jtstart, first_state,
    update_sparse_state_transition!, vw_revenue,
    split_thetas, SimulationTmp, SimulationPrimitives,
    simulationPrimitives_information, dataFramesForSimulations,
    SharedSimulations, update_sim_dataframes_from_simdata!,
    value_function,check_theta, doSimulations, actionspace, statespace

num_i = 50
M = 10
nz = 13
nψ = 13
Tstop = 150

datafull, thetafull = last(MakeTestData(; num_i=num_i, nz=nz, nψ=nψ))

ddata = drill(datafull)
wpold = statespace(_model(ddata))

newrwrd = DrillReward(
    DrillingRevenue(Unconstrained(), NoTrend(), GathProcess(), NoLearn(), WithRoyalty()),
    DrillingCost_constant(),
    ExtensionCost_Constant()
)
newwp = PerpetualProblem(wpold)

@testset "Revision of statespace" begin
    newd = DataDrill(ddata, newrwrd, newwp)
    for (xo, xn) in zip(_x(ddata), _x(newd))
        Do = _D(statespace(_model(ddata)), xo)
        Dn = _D(statespace(_model(newd)), xn)
        @test Do == Dn
    end
end

ddm_novf = DDM_NoVF(_model(ddata))
data_dso = DataDrillStartOnly(ddata)
data_for_xfer = DataDrill(data_dso, ddm_novf)

@test ShaleDrillingLikelihood.num_i(datafull) ==
    ShaleDrillingLikelihood.num_i(data_for_xfer)

newdatafull = DataSetofSets(datafull, data_for_xfer)
theta = thetafull
sim = SimulationDraws(datafull, M)
update!(sim, ThetaRho())
sharesim = SharedSimulations(data_for_xfer)
simprim = SimulationPrimitives(newdatafull, sim, Tstop, theta, sharesim)
simtmp = SimulationTmp(simprim)

@assert _nparm(newrwrd) == _nparm(drill(newdatafull))

@test value_function(_model(drill(newdatafull))) == nothing
@test SimulationTmp(simprim) === simtmp

simlist = [
    (newrwrd, PerpetualProblem(statespace(data_for_xfer)), theta),
]

df_t, df_Tstop = dataFramesForSimulations(simlist, data_for_xfer, Tstop)

update_sim_dataframes_from_simdata!(df_t, df_Tstop, sharesim, 1)

# make these objects
u = ObservationGroup(ddata, 1)
r = ObservationGroup(u, InitialDrilling())
l = ObservationGroup(r, 1)
lc = LeaseCounterfactual(l)

# using ShaleDrillingLikelihood: actionspace
#
# for (obs, zt) in lc
#     @show obs.ichars, obs.z, obs.y, obs.x, zt, actionspace(statespace(obs.model), obs.x)
# end

@testset "last counterfactual is an undrilled state" begin
    for unit in ddata
        for regime in unit
            for lease in regime
                lc = LeaseCounterfactual(lease)
                @test length(lc) >= 0
                obslast = 0
                for (t,(obs,zt)) in enumerate(lc)
                    if t == length(lc)
                        wp = statespace(_model(obs))
                        @test _x(obs) == end_ex0(wp)+1
                    end
                end
            end
        end
    end
end

# to make sure state if not drilled is OK
@testset "iteration through state if never drilled in lease cntrfact" begin
    @test state_if_never_drilled(lc, 0) == 1

    @test tstart(lc) == firstindex(lc)
    @test length(lc) == length(jtstart(lc):length(zchars(ddata)))
    @test first_state(lc) == first(_x(l))
    @test _x(first(first(lc))) == _x(first(l))

    for (t, (obs,zt)) in enumerate(lc)
        if t == 1
            @test _x(obs) == first_state(lc)
        else
            _x(obs) > first_state(lc)
        end
    end
end # testset

@testset "trying out sparse_state_transition" begin
    thetasvw = split_thetas(newdatafull, theta)
    thet = first(thetasvw)
    for (obs,zt) in lc
        simm = getindex(view(sim, 1), rand(1:M))
        update_sparse_state_transition!(simtmp, obs, simm, thet)
    end
end

@testset "simulate 1 lease 1 time" begin
    simi = view(sim, uniti(l))
    simim = simi[1]

    fill!(sharesim, 0)
    simulate_lease!(simprim, l, simim, 1.0)

    @test abs(sum(sharesim.d0[:,1])) > 0
    @test abs(sum(sharesim.profit[:,1])) > 0
    @test abs(sum(sharesim.surplus[:,1])) > 0
    @test abs(sum(sharesim.revenue[:,1])) > 0
    @test abs(sum(sharesim.D_at_T[:,1])) > 0
end

using ShaleDrillingLikelihood: simulate_unit!

@testset "simulate unit" begin
    simprim2 = SimulationPrimitives(newdatafull, sim,
        newrwrd, wpold, Tstop, theta, sharesim)
        update!(sim, ThetaRho())
    fill!(sharesim, 0)
    do_r = false
    for i in 1:ShaleDrillingLikelihood.num_i(newdatafull)
        simulate_unit!(simprim2, i, do_r)
    end

    @test abs(sum(sharesim.d0        )) > 0
    @test abs(sum(sharesim.d1        )) > 0
    @test abs(sum(sharesim.d0psi     )) > 0
    @test abs(sum(sharesim.d1psi     )) > 0
    @test abs(sum(sharesim.d0eur     )) > 0
    @test abs(sum(sharesim.d1eur     )) > 0
    @test abs(sum(sharesim.d0eursq   )) > 0
    @test abs(sum(sharesim.d1eursq   )) > 0
    @test abs(sum(sharesim.d0eurcub  )) > 0
    @test abs(sum(sharesim.d1eurcub  )) > 0
    @test abs(sum(sharesim.epsdeq1   )) > 0
    @test abs(sum(sharesim.epsdgt1   )) > 0
    @test abs(sum(sharesim.Prdeq1    )) > 0
    @test abs(sum(sharesim.Prdgt1    )) > 0
    @test abs(sum(sharesim.Eeps      )) > 0
    @test abs(sum(sharesim.profit    )) > 0
    @test abs(sum(sharesim.surplus   )) > 0
    @test abs(sum(sharesim.revenue   )) > 0
    @test abs(sum(sharesim.drillcost )) > 0
    @test abs(sum(sharesim.extension )) > 0
    @test abs(sum(sharesim.D_at_T    )) > 0
end

@testset "using global objects" begin

    ddata = drill(datafull)
    wp = statespace(_model(ddata))
    rwrd = DrillReward(
        DrillingRevenue(Unconstrained(), NoTrend(), GathProcess(), NoLearn(), WithRoyalty()),
        DrillingCost_constant(),
        ExtensionCost_Constant()
    )

    ddm_novf = DDM_NoVF(_model(drill(datafull)))
    datadrill_dso = DataDrillStartOnly(drill(datafull))
    datadrill_bare= DataDrill(data_dso, ddm_novf)

    sharesim2 = SharedSimulations(datadrill_bare)
    fill!(sharesim2, 0)

    data_bare = DataSetofSets(datafull, datadrill_bare)
    check_theta(data_bare, theta)
    simMat = SimulationDraws(data_bare, M)
    set_g_BaseDataSetofSets(data_bare)
    set_g_SimulationDrawsMatrix(simMat)
    set_g_SharedSimulations(sharesim2)
    simprim = set_g_SimulationPrimitives(rwrd, wp, Tstop, theta)

    @test simprim.sharedsim === sharesim2 #FIXME

    fill!(sharesim2, 0)
    for i in 1:ShaleDrillingLikelihood.num_i(newdatafull)
        simulate_unit!(i, true)
    end

    @test all(isfinite.(sharesim2.d0))

    g_sharesim = get_g_SharedSimulations()

    @test abs(sum(g_sharesim.d0        )) > 0
    @test abs(sum(g_sharesim.d1        )) > 0
    @test abs(sum(g_sharesim.d0psi     )) > 0
    @test abs(sum(g_sharesim.d1psi     )) > 0
    @test abs(sum(g_sharesim.d0eur     )) > 0
    @test abs(sum(g_sharesim.d1eur     )) > 0
    @test abs(sum(g_sharesim.d0eursq   )) > 0
    @test abs(sum(g_sharesim.d1eursq   )) > 0
    @test abs(sum(g_sharesim.d0eurcub  )) > 0
    @test abs(sum(g_sharesim.d1eurcub  )) > 0
    @test abs(sum(g_sharesim.epsdeq1   )) > 0
    @test abs(sum(g_sharesim.epsdgt1   )) > 0
    @test abs(sum(g_sharesim.Prdeq1    )) > 0
    @test abs(sum(g_sharesim.Prdgt1    )) > 0
    @test abs(sum(g_sharesim.Eeps      )) > 0
    @test abs(sum(g_sharesim.profit    )) > 0
    @test abs(sum(g_sharesim.surplus   )) > 0
    @test abs(sum(g_sharesim.revenue   )) > 0
    @test abs(sum(g_sharesim.drillcost )) > 0
    @test abs(sum(g_sharesim.extension )) > 0
    @test abs(sum(g_sharesim.D_at_T    )) > 0

end

@testset "multiple simulations" begin
    rwrd1 = reward(_model(ddata))
    wp1 = statespace(_model(ddata))
    rwrd2 = DrillReward(
        DrillingRevenue(Unconstrained(), NoTrend(), GathProcess(), NoLearn(), NoRoyalty()),
        DrillingCost_constant(),
        ExtensionCost_Constant()
    )
    wp2 = PerpetualProblem(wp1)

    theta_d = thetafull[1:_nparm(ddata)]

    simlist = [
        (rwrd1, wp1, thetafull),
        (rwrd2, wp1, thetafull),
        (rwrd1, wp2, thetafull),
    ]

    # pids = start_up_workers(ENV)
    # @everywhere using ShaleDrillingLikelihood
    M = 50
    Tstop = 60
    doSimulations(datafull, simlist, Tstop, M)
end

# end # module


text:
module timeevolution_master

export master, master_nh, master_h

using ...bases
using ...operators
using ...states
using ...ode_dopri

"""
Evaluate master equation for diagonal jump operators.
"""
function dmaster_h(rho::DenseOperator, H::Operator,
                    Gamma::Vector{Complex128}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(complex(0,-1.), H, rho, complex(0.), drho)
    operators.gemm!(complex(0,1.), rho, H, complex(1.), drho)
    for i=1:length(J)
        operators.gemm!(Gamma[i], J[i], rho, complex(0.), tmp)
        operators.gemm!(complex(1.), tmp, Jdagger[i], complex(1.), drho)

        operators.gemm!(complex(-0.5), Jdagger[i], tmp, complex(1.), drho)

        operators.gemm!(Gamma[i], rho, Jdagger[i], complex(0.), tmp)
        operators.gemm!(complex(-0.5), tmp, J[i], complex(1.), drho)
    end
    return drho
end

"""
Evaluate master equation for non-diagonal jump operators.
"""
function dmaster_h(rho::DenseOperator, H::Operator,
                    Gamma::Matrix{Complex128}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(complex(0,-1.), H, rho, complex(0.), drho)
    operators.gemm!(complex(0,1.), rho, H, complex(1.), drho)
    for j=1:length(J), i=1:length(J)
        operators.gemm!(Gamma[i,j], J[i], rho, complex(0.), tmp)
        operators.gemm!(complex(1.), tmp, Jdagger[j], complex(1.), drho)

        operators.gemm!(complex(-0.5), Jdagger[j], tmp, complex(1.), drho)

        operators.gemm!(Gamma[i,j], rho, Jdagger[j], complex(0.), tmp)
        operators.gemm!(complex(-0.5), tmp, J[i], complex(1.), drho)
    end
    return drho
end

"""
Evaluate master equation for non-hermitian Hamiltonian and diagonal jump operators.
"""
function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
                    Gamma::Vector{Complex128}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(complex(0,-1.), Hnh, rho, complex(0.), drho)
    operators.gemm!(complex(0,1.), rho, Hnh_dagger, complex(1.), drho)
    for i=1:length(J)
        operators.gemm!(Gamma[i], J[i], rho, complex(0.), tmp)
        operators.gemm!(complex(1.), tmp, Jdagger[i], complex(1.), drho)
    end
    return drho
end

"""
Evaluate master equation for non-hermitian Hamiltonian and non-diagonal jump operators.
"""
function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
                    Gamma::Matrix{Complex128}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(complex(0,-1.), Hnh, rho, complex(0.), drho)
    operators.gemm!(complex(0,1.), rho, Hnh_dagger, complex(1.), drho)
    for j=1:length(J), i=1:length(J)
        operators.gemm!(Gamma[i,j], J[i], rho, complex(0.), tmp)
        operators.gemm!(complex(1.), tmp, Jdagger[j], complex(1.), drho)
    end
    return drho
end

"""
Integrate master equation with the given derivative function.

Arguments
---------

dmaster
    A function f(t, rho, drho) that calculates the time-derivative of rho at
    time t and stores the result in drho.
tspan
    Vector specifying the points of time for which output should be displayed.
rho0
    Initial density operator.

Keyword arguments
-----------------

fout (optional)
    If given this function fout(t, rho) is called every time an output should
    be displayed. To limit copying to a minimum the given density operator rho
    is further used and therefore must not be changed.

kwargs
    Further arguments are passed on to the ode solver.
"""
function integrate_master(dmaster::Function, tspan, rho0::DenseOperator; fout=nothing, kwargs...)
    nl = prod(rho0.basis_l.shape)
    nr = prod(rho0.basis_r.shape)
    N = nl*nr
    as_operator(x::Vector{Complex128}) = DenseOperator(rho0.basis_l, rho0.basis_r, reshape(x, nl, nr))
    as_vector(rho::DenseOperator) = reshape(rho.data, N)
    f = (x->x)
    if fout==nothing
        tout = Float64[]
        xout = DenseOperator[]
        function fout_(t, rho::DenseOperator)
            push!(tout, t)
            push!(xout, deepcopy(rho))
        end
        f = fout_
    else
        f = fout
    end
    f_(t, x::Vector{Complex128}) = f(t, as_operator(x))
    dmaster_(t, x::Vector{Complex128}, dx::Vector{Complex128}) = dmaster(t, as_operator(x), as_operator(dx))
    ode(dmaster_, float(tspan), as_vector(rho0), f_; kwargs...)
    return fout==nothing ? (tout, xout) : nothing
end

function _check_input(rho0::DenseOperator, H::Operator, J::Vector, Jdagger::Vector, Gamma::Union{Vector{Float64}, Matrix{Float64}})
    operators.check_samebases(rho0, H)
    for j=J
        @assert typeof(j) <: Operator
        operators.check_samebases(rho0, j)
    end
    for j=Jdagger
        @assert typeof(j) <: Operator
        operators.check_samebases(rho0, j)
    end
    @assert length(J)==length(Jdagger)
    if typeof(Gamma) == Matrix{Float64}
        @assert size(Gamma, 1) == size(Gamma, 2) == length(J)
    elseif typeof(Gamma) == Vector{Float64}
        @assert length(Gamma) == length(J)
    else
        error()
    end
end

"""
Integrate the master equation with dmaster_h as derivative function.

For further information look at :func:`master(,::DenseOperator,::Operator,)`
"""
function master_h(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                Gamma::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
                Jdagger::Vector=map(dagger, J),
                fout::Union{Function,Void}=nothing,
                tmp::DenseOperator=deepcopy(rho0),
                kwargs...)
    _check_input(rho0, H, J, Jdagger, Gamma)
    Gamma = complex(Gamma)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, Gamma, J, Jdagger, drho, tmp)
    return integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

master_h(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master_h(tspan, tensor(psi0, dagger(psi0)), H, J; kwargs...)

"""
Integrate the master equation with dmaster_nh as derivative function.

In this case the given Hamiltonian is assumed to be the non-hermitian version.
For further information look at :func:`master(,::DenseOperator,::Operator,)`
"""
function master_nh(tspan, rho0::DenseOperator, Hnh::Operator, J::Vector;
                Gamma::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
                Hnhdagger::Operator=dagger(Hnh),
                Jdagger::Vector=map(dagger, J),
                fout::Union{Function,Void}=nothing,
                tmp::DenseOperator=deepcopy(rho0),
                kwargs...)
    _check_input(rho0, Hnh, J, Jdagger, Gamma)
    Gamma = complex(Gamma)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, Gamma, J, Jdagger, drho, tmp)
    return integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

master_nh(tspan, psi0::Ket, Hnh::Operator, J::Vector; kwargs...) = master_nh(tspan, tensor(psi0, dagger(psi0)), Hnh, J; kwargs...)

"""
Integrate the master equation with dmaster_nh as derivative function.

There are two implementations for integrating the master equation:

* ``master_h``: Usual formulation of the master equation.
* ``master_nh``: Variant with non-hermitian Hamiltonian.

The ``master`` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

Arguments
---------

tspan
    Vector specifying the points of time for which output should be displayed.
rho0
    Initial density operator (must be a dense operator). Can also be a
    state vector which is automatically converted into a density operator.
H
    DenseOperator specifying the Hamiltonian.
J
    Vector containing all jump operators.

Keyword Arguments
-----------------

Gamma (optional)
    Vector or matrix specifying the coefficients for the jump operators.
Jdagger (optional)
    Vector containing the hermitian conjugates of the jump operators. If they
    are not given they are calculated automatically.
fout (optional)
    If given this function fout(t, rho) is called every time an output should
    be displayed. To limit copying to a minimum the given density operator rho
    is further used and therefore must not be changed.
kwargs
    Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                Gamma::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
                Jdagger::Vector=map(dagger, J),
                fout::Union{Function,Void}=nothing,
                tmp::DenseOperator=deepcopy(rho0),
                kwargs...)
    _check_input(rho0, H, J, Jdagger, Gamma)
    Gamma = complex(Gamma)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, Gamma, J, Jdagger, drho, tmp)
    return integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

function master(tspan, rho0::DenseOperator, H::DenseOperator, J::Vector{DenseOperator};
                Gamma::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)),
                Jdagger::Vector{DenseOperator}=map(dagger, J),
                fout::Union{Function,Void}=nothing,
                tmp::DenseOperator=deepcopy(rho0),
                kwargs...)
    _check_input(rho0, H, J, Jdagger, Gamma)
    Hnh = deepcopy(H)
    if typeof(Gamma) == Matrix{Float64}
        for i=1:length(J), j=1:length(J)
            Hnh -= 0.5im*Gamma[i,j]*Jdagger[i]*J[j]
        end
    elseif typeof(Gamma) == Vector{Float64}
        for i=1:length(J)
            Hnh -= 0.5im*Gamma[i]*Jdagger[i]*J[i]
        end
    else
        error()
    end
    Hnhdagger = dagger(Hnh)
    Gamma = complex(Gamma)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, Gamma, J, Jdagger, drho, tmp)
    return integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

master(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master(tspan, tensor(psi0, dagger(psi0)), H, J; kwargs...)

end #module


text:
"""
treeknit

# Arguments

- `nwk1`: Newick file for first tree
- `nwk2`: Newick file for second tree

# Options

- `-o, --outdir <arg>`: directory for results; Example `-o=treeknit_results`
- `-g, --gamma <arg>`: value of γ; Example `-g=2`
- `--seq-lengths <arg>`: length of the sequences. Example: `--seq-length "1500 2000"`
- `--n-sa-it <arg>`: number of SA iterations per temperature and per leaf; default 0.2

# Flags

- `--naive`: Naive inference (overrides `-g`).
- `--no-likelihood`: Do not use branch length likelihood test to sort between different MCCs
- `--no-resolve`: Do not attempt to resolve trees before inferring MCCs.
- `-v, --verbose`: verbosity
"""
@main function treeknit(
	nwk1::AbstractString, nwk2::AbstractString;
	# options
	outdir::AbstractString = "treeknit_results",
	gamma::Float64 = 2.,
	seq_lengths::AbstractString = "1 1",
	n_sa_it::Float64 = 0.2,
	# flags
	naive::Bool = false,
	no_likelihood::Bool = false,
	no_resolve::Bool = false,
	verbose::Bool = false,
)

	println("Treeknit: ")
	println("Input trees: $nwk1 \t $nwk2")
	println("Results directory: $outdir")
	println("γ: $gamma")

	# Setting up directories
	mkpath(outdir)

	# Setting loggers
	loggers = []
	## File log.txt
	io = open(outdir*"/log.txt", "w+")
	file_logger = FormatLogger(io) do io, args
		println(io, args.file, ":", args.line, " [", args.level, "] ", args.message)
	end;
	push!(loggers, file_logger)
	## If verbose, stderr
	if verbose
		push!(loggers, ConsoleLogger())
	else
		push!(loggers, MinLevelLogger(ConsoleLogger(), Logging.Warn))
	end

	global_logger(TeeLogger(loggers...))

	# Reading trees
	@info "Input Newick files: $nwk1 \t $nwk2. Reading trees..."
	t1 = read_tree(nwk1)
	t2 = read_tree(nwk2)
	if !TreeTools.share_labels(t1, t2)
		error("Trees must share leaves")
	end

	# Setting up OptArgs
	sl = map(i->parse(Int, i), split(seq_lengths, " "))
	oa = OptArgs(;
		γ = gamma,
		likelihood_sort = !no_likelihood,
		resolve = !no_resolve,
		Md = 1 / n_sa_it,
		seq_lengths = sl,
		verbose=true,
	)

	#
	@info "Parameters: $oa"

	@info "Inferring MCCs...\n"
	out = @timed computeMCCs(t1, t2, oa; naive)
	MCCs = out[1]
	@info "Found $(length(MCCs)) MCCs (runtime $(out[2]))\n"

	verbose && println()

	@info "Resolving trees based on found MCCs..."
	rS = resolve!(t1, t2, MCCs)
	@info "Resolved $(length(rS[1])) splits in $(nwk1) and $(length(rS[1])) splits in $(nwk2)\n"

	verbose && println()

	@info "Building ARG from trees and MCCs..."
	arg, rlm, lm1, lm2 = SRG.arg_from_trees(t1, t2, MCCs)
	@info "Found $(length(arg.hybrids)) reassortments in the ARG.\n"

	verbose && println()

	# Write output
	@info "Writing results in $(outdir)"
	write_mccs(outdir * "/" * "MCCs.dat", MCCs)
	for (nwk, t) in zip((nwk1,nwk2),(t1,t2))
		fn = basename(nwk)
		name, ext = splitext(fn)
		write_newick(outdir * "/" * name * ".resolved" * ext, t)
	end
	write(outdir * "/" * "arg.nwk", arg)
	write_rlm(outdir * "/" * "nodes.dat", rlm)

	close(io)

	println()
end

function write_rlm(filename, rlm)
	open(filename, "w") do io
		for (n, (k,v)) in enumerate(rlm)
			write(io, k)
			write(io, ",")
			for (i,vv) in enumerate(v)
				if isnothing(vv)
					write(io, " ")
				else
					write(io, vv)
				end
				if i != length(v)
					write(io, ",")
				end
			end
			if n != length(rlm)
				write(io, "\n")
			end
		end
	end
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

fbase = fdistribute(x->sin(x), y->sin(y), (x,y)->0, 2, 2, 2)
dbase = distribute(x->sin(x), y->sin(y), (x,y)->0, 2, 2, 2)

@test fbase[[1,1]].value ≈ dbase[[1,1]].value
@test fbase[[1,2]].value ≈ dbase[[1,2]].value
@test fbase[[2,1]].value ≈ dbase[[2,1]].value
@test fbase[[2,2]].value ≈ dbase[[2,2]].value


text:
using REPL
using Test
using TerminalPager

function _create_pagerd(str::AbstractString)
    lines = split(str, '\n')
    matches = NTuple{4, Int}[]
    term = REPL.Terminals.TTYTerminal("", stdin, stdout, stderr)
    iobuf = IOBuffer()
    buf = IOContext(iobuf, :color => get(stdout, :color, true))
    pagerd = TerminalPager.Pager(term = term,
                                 buf = buf,
                                 display_size = displaysize(term.out_stream),
                                 start_row = 1,
                                 start_col = 1,
                                 lines = lines,
                                 num_lines = length(lines))

    return pagerd
end

include("./internals/key_bindings.jl")
include("./internals/key_processing.jl")
include("./internals/view.jl")


text:
using Distributions, PlotRecipes, Statistics

# Order of variables in function does not matter (since cor(y,x) = cor(x,y))
function correlation_ttest(y,x; h0=0.0)
    n = length(y)
    v = n - 2
    X = [ones(n) x]
    b = inv(X'X)*X'y
    res = y .- X*b
    s2 = sum(res.^2)/v
    vb = s2.*inv(transpose(X)*X)
    seb = sqrt(vb[2,2])
#    that = abs(b[2]/seb)
    ts = (seb*std(x)/std(y)).*rand(TDist(v),1000000) .+ (b[2]*std(x)/std(y))
    that = abs(mean(ts) - h0)/std(ts)
    todds = (1.0 + (that^2)/v)^(0.5*(v+1))
    pval = 2.0*minimum([cdf(TDist(v),that) (1.0 - cdf(TDist(v),that))])
	ps = [0.005; 0.025; 0.5; 0.95; 0.975]
	q = quantile(ts,[0.005,0.025,0.5,0.95,0.975])
	qs = [ps q]
	plt = plot(ts,st=:density, label="posterior", fill=(0,0.4,:blue),alpha=0.4, title="Posterior for Correlation Coefficient")
	vline!([cor(y,x) h0],linecolor = [:green :black],label=["mean" "H0"])
	return todds, pval, qs, plt
end

# Example of use (vary coefficient on x to increase/reduce correlation)
#n = 50
#x = randn(n)
#y = 1.0 .+ 1.0.*x .+ randn(50)
#cor(x,y)
# cor(x,y) = bhat*sd(x)/sd(y)

#results = correlation_ttest(x,y, h0=0.5)
#results[1]
#results[2]
#results[3]
#results[4]
# savefig("trash.png")


text:
@userplot Spectrogram

@recipe function f(h::Spectrogram;
    yscale=:identity,
    n=nothing,
    noverlap=0.9,
    seriescolor=:inferno,
    dbscale=false
)
    if  (typeof(h.args[1]) <: SeisData)
        @warn("Only plotting spectrogram of first SeisChannel in SeisData")
        C = h.args[1][1]
    elseif (typeof(h.args[1]) <: SeisChannel)
        C = h.args[1]
    else
        error("First argument must be SeisChannel or SeisData.")
    end

    if size(C.t,2) != 2
        error("Cannot plot SeisChannel with gaps. Please use ungap! first.")
    end

    # set up the subplots

    # framestyle := [:none :axes]
    grid := false
    layout := (2,1)

    # get spectrogram
    if isnothing(n)
        n = nextpow(2,C.fs ^ 2 / 50)
    else
        n = nextpow(2,n)
    end
    newoverlap = floor(Int, n * noverlap)
    spec = DSP.spectrogram(C.x,n,newoverlap,fs=C.fs)
    t = u2d.(SeisIO.t_expand(C.t,C.fs) * SeisIO.μs)

    if dbscale
        power = 10 .* log10.(spec.power)
    else
        power = sqrt.(spec.power)
    end

    # top plot time series
    @series begin
        subplot := 1
        label --> C.id
        legendfontsize --> 8
        ytickfontsize --> 10
        xtickfontsize --> 8
        seriescolor --> :black
        yscale := :identity
        t,C.x
    end

    # main heatmap
    @series begin
        legend := false
        seriestype := :heatmap
        yscale --> yscale
        ytickfontsize --> 10
        xtickfontsize --> 8
        subplot := 2
        xguide --> "Time [s]"
        yguide --> "Frequency [Hz]"
        xlims --> (t[1],t[end])
        seriescolor --> cmap
        t[1] .+ Millisecond.(spec.time .* 1000), spec.freq[2:end], power[2:end,:]
    end
end


text:
module DataSet

function two_squared(n)
    x = rand(0:100,n)*0.01
    y = rand(0:100,n)*0.01
    _t_data = [] #正解データ（上1か下0か）

    for i in 1:n
        if x[i]^2 >= y[i]
            #下にある
            append!(_t_data,0)
        else
            append!(_t_data,1)
        end
    end

    _data = hcat(x,y)

    return _data,_t_data
end

function noise_sin(n)
    data,t_data = normal_sin(n,interval=0.05)
    t_data += rand(size(t_data,1))*0.2

    return data,t_data
end

function normal_sin(n;interval=0.1)
    data = Array(-pi*2*n : interval : pi*2*n)
    t_data = sin.(data)

    return data,t_data
end

end