text:
# -----------------------------------------------------------------------------

doc"Represents a `Constant Distribution`"
immutable Constant{T} <: DiscreteUnivariateDistribution
  c::T
end
StatsBase.sample(c::Constant) = c.c

immutable HiddenLayerSampler <: DiscreteMultivariateDistribution
  numlayers::AVecI
  numnodes::AVecI
end
function StatsBase.sample(hls::HiddenLayerSampler)
  numlayers = sample(hls.numlayers)
  Int[sample(hls.numnodes) for i in 1:numlayers]
end

immutable VectorSampler{T} <: DiscreteUnivariateDistribution
  vec::AVec{T}
  probs::VecF
end
VectorSampler(vec::AVec) = VectorSampler(vec, fill(1.0 / length(vec), length(vec)))
VectorSampler{T}(vec::AVec{Tuple{T,Float64}}) = VectorSampler(unzip(vec)...)
function StatsBase.sample(vs::VectorSampler)
  r = rand()
  totp = 0.0
  for (i,p) in enumerate(vs.probs)
    totp += p
    if r <= totp
      return vs.vec[i]
    end
  end
  error("probs don't sum to 1?? $r $totp $vs")
end

# -----------------------------------------------------------------------------

StatsBase.sample(d::Distribution) = rand(d)

doc"""
A `ParameterSampler` holds vectors of keyword parameter names and the distributions to sample from
when randomly choosing a parameter set.  Example:

```
  ps = ParameterSampler(Dict(:x => Constant(1), :y => Normal()))
```
"""
type ParameterSampler
  d::Dict
end
StatsBase.sample(ps::ParameterSampler) = [(s,sample(dist)) for (s,dist) in ps.d]

function Base.print(io::IO, ps::ParameterSampler)
  println("ParameterSampler{")
  for k in sort(keys(ps.d))
    println("  $k => $(typeof(ps.d[k]))")
  end
  print("} ")
end
Base.show(io::IO, ps::ParameterSampler) = print(io, ps)

# -----------------------------------------------------------------------------

doc"""
Creates a `VectorTransformer` by sampling (N times) both a transform type and the index that it maps to.

Note that we call `unique` so as not to duplicate data, so we may actually return less than N samples.
"""
function generateTransformer(transformSampler::VectorSampler{DataType}, indices::AVecI, N::Int)
  VectorTransformer(unique(Transformation[sample(transformSampler)(sample(indices)) for i in 1:N]))
end

doc"""
Create numModels models, all using the `buildFn`, which handles creating the model with the constant args/kwargs,
and a random sampling of the parameters from the parameter sampler.
Note we assume all parameters we care about setting are keyword params.
"""
function generateModels(buildFn::Function, ps::ParameterSampler, numModels::Integer)
  OnlineStat[buildFn(; sample(ps)...) for i in 1:numModels]
end

# -----------------------------------------------------------------------------

# TODO: ensemble tools, parameter optimization (grid search, etc)

doc"""
The Ensemble type represents a set of models which have been fit to the same data.
A call to fit! will update all the models with that data.  Then the outputs of
the models (yhatᵢⱼ) are stacked and get passed through to the `agg` aggregator. `agg`
is another OnlineStat which regresses yhats (model estimates) on y (targets).

Note: predict first generates estimates from each model, then uses those to predict y
with the `agg` model aggregator.
"""
type Ensemble{STAT} <: OnlineStat
  models::Vector{OnlineStat}
  aggs::Vector{STAT}  # one aggregation model for each target variable... each agg maps nummodels vars to numtargets vars
end

# the build function should take only keyword arguments to build, which can be populated with a ParameterSampler
function Ensemble{AGG<:OnlineStat}(::Type{AGG}, buildFn::Function, ps::ParameterSampler, numModels::Int, nout::Int, args...; kwargs...)
  Ensemble(generateModels(buildFn, ps, numModels), AGG[AGG(args...; kwargs...) for i in 1:nout])
end

# OnlineStats.statenames(o::Ensemble) = [:nobs]
# OnlineStats.state(o::Ensemble) = Any[nobs(o)]
nobs(o::Ensemble) = 0

StatsBase.predict(models::Vector{OnlineStat}, x::AVecF) = [predict(model, x) for model in models]

# get estimates from each model.  then for each i, grab the estimates of the iᵗʰ target variable
# and use those as the inputs to the iᵗʰ aggregation model
function StatsBase.predict(ensemble::Ensemble, x::AVecF)
  ests = predict(ensemble.models)
  nout = length(ensemble.aggs)
  yhat = zeros(nout)
  for i in 1:nout
    esti = Float64[est[i] for est in ests]
    yhat[i] = predict(ensemble.aggs[i], esti)
  end
end

# similar to predict, but fit! instead
function OnlineStats.fit!(ensemble::Ensemble, x::AVecF, y::AVecF)
  for model in ensemble.models
    fit!(model, x, y)
  end

  ests = predict(ensemble.models)

  yhat = similar(y)
  for i in 1:length(y)
    esti = Float64[est[i] for est in ests]
    fit!(ensemble.aggs[i], esti, y)
  end
end


text:
using ExamplePackage
using Test

@testset "ExamplePackage.jl" begin
    # Write your own tests here.
end

ext: jl


text:
# make QuasiNewton{SR1, approx} callable and then pass that as the
# "update" parameter
# the minimize can take a single SR1 input as -> QuasiNewton{SR1, Approx}
struct SR1{T1} <: QuasiNewton{T1}
   approx::T1
end
function update(H, s, y, scheme::SR1{<:InverseApprox})
   sHy = s - H*y
   d_sHy_y = inv(dot(sHy, y))
   if d_sHy_y > 1e7
      H = H + d_sHy_y*sHy*sHy'
   end
   H
end
function update(B, s, y, scheme::SR1{<:DirectApprox})
   yBs = y - B*s
   d_yBs_s = inv(dot(yBs, s))
   if d_yBs_s > 1e7
      B = B + d_yBs_s*yBs*yBs'
   end
   B
end
function update!(H, s, y, scheme::SR1{<:InverseApprox})
   sHy = s - H*y
   d_sHy_y = inv(dot(sHy, y))
   if d_sHy_y > 1e7
      H .= H .+ d_sHy_y*sHy*sHy'
   end
   H
end
function update!(B, s, y, scheme::SR1{<:DirectApprox})
   yBs = y - B*s
   d_yBs_s = inv(dot(yBs, s))
   if d_yBs_s > 1e7
      B .= B .+ d_yBs_s*yBs*yBs'
   end
   B
end
function update!(A::UniformScaling, s, y, scheme::SR1{<:InverseApprox})
   update(A, s, y, scheme)
end
function update!(A::UniformScaling, s, y, scheme::SR1{<:DirectApprox})
   update(A, s, y, scheme)
end

function find_direction(A, scheme::SR1, ::DirectApprox, ∇f)
   -(A\∇f)
end
function find_direction!(d, B, scheme::SR1, ::DirectApprox, ∇f)
   d .= -(B\∇f)
   d
end


text:
abstract State

immutable StatePair
    current::State
    proposal::State
end

immutable ConcreteState <: State  # move out
    value::Float64
end


text:
#!/usr/bin/env julia

using YAML

# config = JSON.parsefile("config_exp.json")
# config =  YAML.load(open("config_exp.yml"))
config =  YAML.load(open(ARGS[1]))

io = stdout  #  open("config_exp.h","w")
println(io,"# this file has been generated automatically")
for (name,subconfig) in config
    println(io,"# $name")
    for (key,value) in subconfig
        println(io,"export $key=$value")
    end
end


text:
using Flux, CuArrays, Metalhead

struct ResidualBlock
  conv_layers
  norm_layers
  shortcut
end

Flux.treelike(ResidualBlock)

# ResidualBlock Function allows us to define a Residual Block having any number of Convolution and Batch Normalization Layers
function ResidualBlock(filters, kernels::Array{Tuple{Int,Int}}, pads::Array{Tuple{Int,Int}}, strides::Array{Tuple{Int,Int}}, shortcut = identity)
  local conv_layers = []
  local norm_layers = []
  for i in 2:length(filters)
    push!(conv_layers, Conv(kernels[i-1], filters[i-1]=>filters[i], pad = pads[i-1], stride = strides[i-1]))
    push!(norm_layers, BatchNorm(filters[i]))
  end
  ResidualBlock(Tuple(conv_layers),Tuple(norm_layers),shortcut)
end

# Function converts the Array of scalar kernel, pad and stride values to tuples
function ResidualBlock(filters, kernels::Array{Int}, pads::Array{Int}, strides::Array{Int}, shortcut = identity)
  ResidualBlock(filters, [(i,i) for i in kernels], [(i,i) for i in pads], [(i,i) for i in strides], shortcut)
end

function (block::ResidualBlock)(input)
  local value = copy.(input)
  for i in 1:length(block.conv_layers)-1
    value = relu.((block.norm_layers[i])((block.conv_layers[i])(value)))
  end
  relu.(((block.norm_layers[end])((block.conv_layers[end])(value))) + block.shortcut(input))
end

# Function to generate the residual blocks for ResNet18 and ResNet34
function BasicBlock(filters::Int, downsample::Bool = false, res_top::Bool = false)
  if(!downsample || res_top)
    return ResidualBlock([filters for i in 1:3], [3,3], [1,1], [1,1])
  end
  shortcut = Chain(Conv((3,3), filters÷2=>filters, pad = (1,1), stride = (2,2)), BatchNorm(filters))
  ResidualBlock([filters÷2, filters, filters], [3,3], [1,1], [1,2], shortcut)
end

# Function to generate the residual blocks used for ResNet50, ResNet101 and ResNet152
function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
  if(!downsample && !res_top)
    return ResidualBlock([4 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1])
  elseif(downsample && res_top)
    return ResidualBlock([filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1], Chain(Conv((1,1), filters=>4 * filters, pad = (0,0), stride = (1,1)), BatchNorm(4 * filters)))
  else
    shortcut = Chain(Conv((1,1), 2 * filters=>4 * filters, pad = (0,0), stride = (2,2)), BatchNorm(4 * filters))
    return ResidualBlock([2 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,2], shortcut)
  end
end

# Function to build Standard Resnet models as described in the paper "Deep Residual Learning for Image Recognition"
function _ResNet(Block, layers; in_chs = 1, initial_filters::Int = 64, nclasses::Int = 160)

  local top = []
  local residual = []
  local bottom = []

  push!(top, Conv((7,7), in_chs=>initial_filters, pad = (3,3), stride = (2,2)))
  push!(top, x -> maxpool(x, (3,3), pad = (1,1), stride = (2,2)))

  for i in 1:length(layers)
    push!(residual, Block(initial_filters, true, i==1))
    for j in 2:layers[i]
      push!(residual, Block(initial_filters))
    end
    initial_filters *= 2
  end

  push!(bottom, x -> meanpool(x, (7,7)))
  push!(bottom, x -> reshape(x, :, size(x,4)))
  if(Block == Bottleneck)
    push!(bottom, (Dense(2048, nclasses)))
  else
    push!(bottom, (Dense(512, nclasses)))
  end
  push!(bottom, softmax)

  Chain(top..., residual..., bottom...)
end

"""
  get_densenet_model(depth, pretrained::Bool = false)

Creates the Resnet model as per the arguments passed.

Arguments:
1. `depth`: The depth of the Resnet Model. The current allowed
            values are `18`, `34`, `50`, `101` and `152`.
2. `pretrained`: Gets the network with pretrained Imagenet weights.
                 `Metalhead.jl` currently has only the weights for
                 ResNet 50 so sending `depth` anything other than
                 `50` results in an **ERROR**.
"""

function get_resnet_model(depth, pretrained::Bool = false)
  if pretrained && depth != 50
    error("Pretrained weights are only available for ResNet 50")
  end
  if depth==18
    _ResNet(BasicBlock, [2, 2, 2, 2]) |> gpu
  elseif depth==34
    _ResNet(BasicBlock, [3, 4, 6, 3]) |> gpu
  elseif depth==50
    if pretrained
      m = ResNet() |> gpu
      Flux.testmode!(m, false)
      m
    else
      _ResNet(Bottleneck, [3, 4, 6, 3]) |> gpu
    end
  elseif depth==101
    _ResNet(Bottleneck, [3, 4, 23, 3]) |> gpu
  elseif depth==152
    _ResNet(Bottleneck, [3, 8, 36, 3]) |> gpu
  end
end


text:
#! /home/sci/dperry/build/julia/julia/julia

using MAT

using RandomFourierFeatures
using Kernel
using Sketch

X = []
Labels = []

name = "gaussian"
recompute = true

d = 0
n = 0
nt = 0

if name == "gaussian"

  n = round(Int64,1e3)
  d = round(Int64,1e2)
	m = 50 # signal dimension

	eta = 10 # signal to noise ratio

  S = randn(n,d)
  D = eye(d)
  for i=1:m
    D[i,i] = max(1.0 - (i-eta*1.0)/m, 0)
  end
  U,R = qr(randn(d,d))

  X = S * D * U

	if eta > 0
		X += randn(size(X)) / eta
	end

	half = round(Int64, ceil(.9*n))
	X[half:end,:] += 1000*randn((n-half+1),d)

	Xtest = X[n/2+1:end,:]
	X = X[1:n/2,:]

	println("X: ", size(X))
	println("Xtest: ", size(Xtest))

end

n = size(X,1)
X_centered = X .- (ones(n,1)*mean(X,1))
nt = size(Xtest, 1)
Xtest_centered = Xtest .- (ones(nt,1)*mean(X,1))

function comparison(X, Xtest, K_test, center, sigma)
  n = size(X,1)
  d = size(X,2)

	println("starting comparison")
	FDstyles = ["bo-","ro-","go-","mo-","co-"]
	#ells = [100, 250, 500, 1000, 2000]

	ells = [2,5,10,20,30,50]
	epsilons = 1 ./ ells

	firstval = 10
	lastval = 2000
	incrval = 50

	if n >= 1e4
		firstval = 500
		lastval = 1500
		incrval = 100
	elseif n >= 1e3
		firstval = 500
		lastval = 1500
		incrval = 100
	elseif n >= 1e2
		firstval = 1
		lastval = 200
		incrval = 10
	else
		firstval = 1
		lastval = 100
		incrval = 5
	end

	dims = [firstval:incrval:lastval]
	if(ells[end] < firstval)
		dims = [ells; dims]
	end

	println("dims: ", dims)

	counter = length(dims)

	inputdim = size(X,2)
	errs_fro = Float64[]
	errs_l2 = Float64[]
	timings = Float64[]
	timings_test = Float64[]
	timings_fro = Float64[]
	timings_l2 = Float64[]
	errs_rankk_fro = zeros(length(ells),counter)
	errs_rankk_l2 = zeros(length(ells),counter)
	timings_rankk = zeros(length(ells),counter)
	timings_test_rankk = zeros(length(ells),counter)
	timings_rankk_fro = zeros(length(ells),counter)
	timings_rankk_l2 = zeros(length(ells),counter)
	errs_FD_fro = zeros(length(ells),counter)
	errs_FD_l2 = zeros(length(ells),counter)
	timings_FD = zeros(length(ells),counter)
	timings_test_FD = zeros(length(ells),counter)
	timings_FD_fro = zeros(length(ells),counter)
	timings_FD_l2 = zeros(length(ells),counter)
	errs_nystrom_fro = Float64[]
	errs_nystrom_l2 = Float64[]
	timings_nystrom = Float64[]
	timings_test_nystrom = Float64[]
	timings_nystrom_fro = Float64[]
	timings_nystrom_l2 = Float64[]
	errs_nystromwr_fro = Float64[]
	errs_nystromwr_l2 = Float64[]
	timings_nystromwr = Float64[]
	timings_test_nystromwr = Float64[]
	timings_nystromwr_fro = Float64[]
	timings_nystromwr_l2 = Float64[]
	errs_rff_nystrom_fro = Float64[]
	errs_rff_nystrom_l2 = Float64[]
	timings_rff_nystrom = Float64[]
	timings_test_rff_nystrom = Float64[]
	timings_rff_nystrom_fro = Float64[]
	timings_rff_nystrom_l2 = Float64[]

  rankk_toobig = false
	if n >= 1e2
		rankk_toobig = true
	end
	counter = 0
	for dim=dims #firstval:incrval:lastval
		counter += 1

		samples = 10

		err_fro = zeros(samples)
		err_l2 = zeros(samples)
		timing = zeros(samples)
		timing_test = zeros(samples)
		timing_fro = zeros(samples)
		timing_l2 = zeros(samples)
		err_rankk_fro = zeros(samples,length(ells))
		err_rankk_l2 = zeros(samples,length(ells))
		timing_rankk = zeros(samples,length(ells))
		timing_test_rankk = zeros(samples,length(ells))
		timing_rankk_fro = zeros(samples,length(ells))
		timing_rankk_l2 = zeros(samples,length(ells))
		err_FD_fro = zeros(samples,length(ells))
		err_FD_l2 = zeros(samples,length(ells))
		timing_FD = zeros(samples,length(ells))
		timing_test_FD = zeros(samples,length(ells))
		timing_FD_fro = zeros(samples,length(ells))
		timing_FD_l2 = zeros(samples,length(ells))
		err_nystrom_fro = zeros(samples)
		err_nystrom_l2 = zeros(samples)
		timing_nystrom = zeros(samples)
		timing_test_nystrom = zeros(samples)
		timing_nystrom_fro = zeros(samples)
		timing_nystrom_l2 = zeros(samples)
		err_nystromwr_fro = zeros(samples)
		err_nystromwr_l2 = zeros(samples)
		timing_nystromwr = zeros(samples)
		timing_test_nystromwr = zeros(samples)
		timing_nystromwr_fro = zeros(samples)
		timing_nystromwr_l2 = zeros(samples)
		err_rff_nystrom_fro = zeros(samples)
		err_rff_nystrom_l2 = zeros(samples)
		timing_rff_nystrom = zeros(samples)
		timing_test_rff_nystrom = zeros(samples)
		timing_rff_nystrom_fro = zeros(samples)
		timing_rff_nystrom_l2 = zeros(samples)

		for sample=1:samples
			println( "random fourier features" )
			tic()
			rffs = GenerateFunctionsGaussian(Float64, sigma, inputdim, dim)
			generation_time = toc()

			# time how long it takes to compute K using RFF
			tic()
			mean_compute_time = toc()

			tic()
			Xp = RFFProject( X, rffs )
			if center
				Xp = Xp - ones(size(Xp,1),1) * mean(Xp,1)
			end
			# fair comparison:
			Cp = zeros(size(Xp,2),size(Xp,2))
			fakeweights = ones(size(Xp,1))
			for i=1:size(Xp,1)
				Cp += fakeweights[i] * (Xp[i,:]'*Xp[i,:])
			end
			println("compute svd of covariance")
			U,S,V = svd(Cp)

			timing[sample] += toc()
			timing[sample] += generation_time
			timing[sample] += mean_compute_time

			tic()
			Xptest = RFFProject( Xtest, rffs )
			rff_proj_time = toc()

			# full RFF approxmation:
			tic()
			kk = min(ells[end], size(V,2))
			XptestV = Xptest * V
			Kp = XptestV*XptestV'
			timing_test[sample] += toc()
			timing_test[sample] += rff_proj_time

			err_fro[sample] = vecnorm(K_test-Kp) / length(K_test)
			err_l2[sample] = norm(K_test-Kp) / length(K_test[1,:])

			for li=1:length(ells)
				l = ells[li]

				println("RFF rank-",l," - projection time")
				tic()
				k = min(dim,l)
				XpV = Xptest * V[:,1:k]
				Kp = XpV * XpV'
				projection_time = toc()
				timing_test_rankk[sample,li] += projection_time

				# time how long it takes to compute the Fro norm
				println( "random fourier features - Fro norm" )
				tic()
				fro_norm_2 = vecnorm(K_test-Kp) / length(K_test)
				err_rankk_fro[sample,li] = fro_norm_2 # frob norm
				timing_rankk_fro[sample,li] += toc() + projection_time + rff_proj_time

				tic()
				l2_norm = norm(K_test-Kp,2) / length(K_test[1,:])
				err_rankk_l2[sample,li] = l2_norm # l2 norm
				timing_rankk_l2[sample,li] += toc() + projection_time + rff_proj_time

				println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
				println("rank-",l," approximation")
				println("sample ", sample, " / ", samples)
				println("rank(K) = ", rank(K_test))
				println("rank(Khat) = ", rank(Kp) )
				println("RFF dim = ", dim, " fro norm = ", fro_norm_2)
				println("RFF dim = ", dim, " l2 norm = ", l2_norm)
				println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
			end

  		println( "random fourier features with FD" )
			for li = 1:length(ells)
				l = ells[li]
				p = 1
				if name == "forest"
					p = round(Int64,floor(.5*l))
				end
				println("p = ", p)
				tic()
				println("running FD")
				B = zeros(l,dim)
				B = FD(Xp, l , min(dim,round(Int64, ceil(l/2.0))))
				println("svd of B")
				U,S,V = svd(B)

				timing_FD[sample,li] += toc()
				timing_FD[sample,li] += generation_time
				timing_FD[sample,li] += mean_compute_time

				tic()
				kk = min(size(V,2),l)
				XpV = Xptest * V[:,1:kk]
				Kp_FD = XpV * XpV'
				timing_test_FD[sample,li] += toc()
				timing_test_FD[sample,li] += rff_proj_time

  			println( "random fourier features with FD - computing fro " )
				tic()
				fro_norm_2 = vecnorm(K_test-Kp_FD) / length(K_test)
				err_FD_fro[sample,li] = fro_norm_2
				timing_FD_fro[sample,li] += toc()
				tic()
				l2_norm = norm(K_test-Kp_FD,2) / length(K_test[1,:])
				err_FD_l2[sample,li] += l2_norm
				timing_FD_l2[sample,li] += toc()
				println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
				println("sample ", sample, " / ", samples)
				println("RFF dim = ", dim)
				println("rank(K) = ", rank(K_test))
				println("rank(Khat) = ", rank(Kp_FD) )
				println("l = ", l, " fro norm = ", fro_norm_2 )
				println("l = ", l, " l2 norm = ", l2_norm )
				println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
			end

			println( "random without replacement sub matrix (Nystrom method)" )
			# reservoir subsample:
			nystrom_num = min(dim, n)

			UU = zeros(1,1)
			SS = zeros(1,1)
			VV = zeros(1,1)
			X_subset = []
			K_subset = []
			attempts = 0
			while size(UU,1) == 1 && attempts < n
				attempts += 1
				if mod(attempts,1000) == 0
					print(attempts,",")
				end
				# keep sampling until a submatrix that can be factored is found:
				tic()
				subset_ind = collect(1:nystrom_num)
				for j=nystrom_num+1:size(X,1)
					r = round(Int64,ceil(rand()*j))
					if r <= nystrom_num
						# X_subset[r,:] = X[j,:]
						subset_ind[r] = j
					end
				end
				X_subset = X[subset_ind,:]

				K_subset = GaussianKernel(X_subset, sigma)
				try
					UU,SS,VV = svd(K_subset)
				catch(error)
					UU = zeros(1,1)
				end
			end
			K_all = GaussianKernel(Xtest, X_subset, sigma)

      epsilon = 1e-10
			if cond(K_subset) > 1e10
        # TSVD:
				SS[ find( SS .<= epsilon ) ] = 0
        ind = find( SS .> epsilon )
        SS[ ind ] = 1 ./ SS[ ind ]
        SS = sqrt(SS)
			else
				SS = sqrt(SS.^(-1))
			end

      epsilon = 1e-10
			if true || cond(K_subset) > 1e10
				UU,SS,VV = svd(K_subset)
        # TSVD:
				k = min(size(UU,2),dim)
        Utkall = UU[:,1:k]' * K_all'
        SS[ find( SS .< epsilon ) ] = 0
        ind = find( SS .> epsilon )
        SS[ ind ] = 1 ./ SS[ ind ]
        SS = sqrt(SS)
        Utkall = diagm(SS) * Utkall
        K_nystrom_wor = Utkall'*Utkall
        # Tikhonov
        #K_subset = UU * diagm(SS .+ epsilon) * VV'
			  #K_nystrom = K_all * (K_subset \ K_all')
      else
			  K_nystrom_wor = K_all * (K_subset \ K_all')
			end

			timing_nystrom[sample] += toc()

			println( "random sub matrix (Nystrom method) - computing fro" )
			tic()
			fro_norm_2 = vecnorm(K_test-K_nystrom_wor) / length(K_test)
			err_nystrom_fro[sample] = fro_norm_2
			timing_nystrom_fro[sample] += toc()
			tic()
			l2_norm =  norm(K_test-K_nystrom_wor,2) / length(K_test[1,:])
			err_nystrom_l2[sample] =	l2_norm
			timing_nystrom_l2[sample] += toc()
			println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
			println("sample ", sample, " / ", samples)
			println("RFF dim = ", dim)
			println("rank(K) = ", rank(K_test))
			println("rank(Khat) = ", rank(K_nystrom_wor) )
			println( "nystrom sample size ", size(K_subset,1), ", fro norm = ", fro_norm_2)
			println( "nystrom sample size ", size(K_subset,1), ", l2 norm  = ", l2_norm)
			println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

			println( "random sub matrix with replacement (Nystrom method)" )
			# reservoir subsample:
			nystrom_num = min(dim, n)

			# fair comparison: iterate over all data points and compute the kernel space norm
			#                  to use importance sampling (reservoir sampling)
			tic()
			X_subset = zeros(nystrom_num,size(X,2))
			norm_est = 0
			for i=nystrom_num+1:size(X,1)
				kern_norm_sq = GaussianKernel(X[i,:],sigma)[1]
				norm_est += kern_norm_sq
				pc = kern_norm_sq / norm_est
				scale = 1 #1/sqrt(nystrom_num*pc)
				for j=1:nystrom_num
					if rand() <= pc
						X_subset[j,:] = scale * X[i,:]
					end
				end
			end

			# do this the same as before to make sure technique above isn't providing biased sample..
			subset_ind = round(Int64, ceil(rand(nystrom_num)*n))
			X_subset = X[subset_ind,:]

			K_subset = GaussianKernel(X_subset, sigma)

			try
				UU,SS,VV = svd(K_subset)
			catch(error)
				UU = zeros(1,1)
			end

			timing_nystromwr[sample] += toc()

			attempts = 0
			while size(UU,1) == 1 && attempts < n # retry the fast way if we picked a bad subset
				attempts += 1
				if mod(attempts,1000) == 0
					print(attempts,",")
				end
				# keep sampling until a submatrix that can be factored is found:
				subset_ind = round(Int64,ceil(rand(nystrom_num)*n))
				X_subset = X[subset_ind,:]

				K_subset = GaussianKernel(X_subset, sigma)
				try
					UU,SS,VV = svd(K_subset)
				catch(error)
					UU = zeros(1,1)
				end
			end

			tic()
			K_all = GaussianKernel(Xtest, X_subset, sigma)

      epsilon = 1e-10
			if cond(K_subset) > 1e10
        # TSVD:
				SS[ find( SS .<= epsilon ) ] = 0
        ind = find( SS .> epsilon )
        SS[ ind ] = 1 ./ SS[ ind ]
        SS = sqrt(SS)
			else
				SS = sqrt(SS.^(-1))
			end

      epsilon = 1e-10
			if true || cond(K_subset) > 1e10
				UU,SS,VV = svd(K_subset)
        # TSVD:
				k = min(size(UU,2),dim)
        Utkall = UU[:,1:k]' * K_all'
        SS[ find( SS .< epsilon ) ] = 0
        ind = find( SS .> epsilon )
        SS[ ind ] = 1 ./ SS[ ind ]
        SS = sqrt(SS)
        Utkall = diagm(SS) * Utkall
        K_nystrom_wr = Utkall'*Utkall
        # Tikhonov
        #K_subset = UU * diagm(SS .+ epsilon) * VV'
			  #K_nystrom = K_all * (K_subset \ K_all')
      else
			  K_nystrom_wr = K_all * (K_subset \ K_all')
			end

			timing_test_nystromwr[sample] += toc()

			println( "random sub matrix (Nystrom method) - computing fro" )
			tic()
			fro_norm_2 = vecnorm(K_test-K_nystrom_wr) / length(K_test)
			err_nystromwr_fro[sample] = fro_norm_2
			timing_nystromwr_fro[sample] += toc()
			tic()
			l2_norm =  norm(K_test-K_nystrom_wr,2) / length(K_test[1,:])
			err_nystromwr_l2[sample] =	l2_norm
			timing_nystromwr_l2[sample] += toc()
			println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
			println("sample ", sample, " / ", samples)
			println("rank(K) = ", rank(K_test))
			println("rank(Khat) = ", rank(K_nystrom_wr) )
			println("size(K_subset) = ", size(K_subset))
			println( "nystrom sample size ", size(K_subset,1), ", fro norm = ", fro_norm_2)
			println( "nystrom sample size ", size(K_subset,1), ", l2 norm  = ", l2_norm)
			println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

		end

		push!(errs_fro, median(err_fro))
		push!(errs_l2, median(err_l2))
		push!(timings, mean(timing))
		push!(timings_test, mean(timing_test))

		errs_rankk_fro[:,counter] = median(err_rankk_fro,1)
		errs_rankk_l2[:,counter] = median(err_rankk_l2,1)
		timings_rankk[:,counter] = mean(timing_rankk,1)
		timings_test_rankk[:,counter] = mean(timing_test_rankk,1)
		timings_rankk_fro[:,counter] = mean(timing_rankk_fro,1)

		errs_FD_fro[:,counter] = median(err_FD_fro,1)
		errs_FD_l2[:,counter] = median(err_FD_l2,1)
		timings_FD[:,counter] = mean(timing_FD,1)
		timings_test_FD[:,counter] = mean(timing_test_FD,1)
		timings_FD_fro[:,counter] = mean(timing_FD_fro,1)

		push!(errs_nystrom_fro, median(err_nystrom_fro))
		push!(errs_nystrom_l2, median(err_nystrom_l2))
		push!(timings_nystrom, mean(timing_nystrom))
		push!(timings_test_nystrom, mean(timing_test_nystrom))
		push!(timings_nystrom_fro, mean(timing_nystrom_fro))

		push!(errs_nystromwr_fro, median(err_nystromwr_fro))
		push!(errs_nystromwr_l2, median(err_nystromwr_l2))
		push!(timings_nystromwr, mean(timing_nystromwr))
		push!(timings_test_nystromwr, mean(timing_test_nystromwr))
		push!(timings_nystromwr_fro, mean(timing_nystromwr_fro))

		push!(errs_rff_nystrom_fro, median(err_rff_nystrom_fro))
		push!(errs_rff_nystrom_l2, median(err_rff_nystrom_l2))
		push!(timings_rff_nystrom, mean(timing_rff_nystrom))
		push!(timings_test_rff_nystrom, mean(timing_test_rff_nystrom))
		push!(timings_rff_nystrom_fro, mean(timing_rff_nystrom_fro))

		#push!(dims, dim)

		println("dim = ", dim, ", FD erro = ", errs_FD_fro[:,counter])
		println("dim = ", dim, ", Nys erro = ", err_nystrom_fro)
		println("dim = ", dim, ", RFF erro = ", err_fro)

		file = matopen(@sprintf("fd_results_%s_n%d_nt%d_d%d_s%f.mat",name,n,nt,d,sigma),"w")
		write(file, "dims", dims)
		write(file, "ells", ells)
		write(file, "epsilons", epsilons)
		write(file, "errs_fro", errs_fro)
		write(file, "errs_l2", errs_l2)
		write(file, "timings", timings)
		write(file, "timings_test", timings_test)
		write(file, "timings_fro", timings_fro)
		write(file, "timings_l2", timings_l2)
		write(file, "errs_rankk_fro", errs_rankk_fro)
		write(file, "errs_rankk_l2", errs_rankk_l2)
		write(file, "timings_rankk", timings_rankk)
		write(file, "timings_test_rankk", timings_test_rankk)
		write(file, "timings_rankk_fro", timings_rankk_fro)
		write(file, "timings_rankk_l2", timings_rankk_l2)
		write(file, "errs_FD_fro", errs_FD_fro)
		write(file, "errs_FD_l2", errs_FD_l2)
		write(file, "timings_FD", timings_FD)
		write(file, "timings_test_FD", timings_test_FD)
		write(file, "timings_FD_fro", timings_FD_fro)
		write(file, "timings_FD_l2", timings_FD_l2)
		write(file, "errs_nystrom_fro", errs_nystrom_fro)
		write(file, "errs_nystrom_l2", errs_nystrom_l2)
		write(file, "timings_nystrom", timings_nystrom)
		write(file, "timings_test_nystrom", timings_test_nystrom)
		write(file, "timings_nystrom_fro", timings_nystrom_fro)
		write(file, "timings_nystrom_l2", timings_nystrom_l2)
		write(file, "errs_nystromwr_fro", errs_nystromwr_fro)
		write(file, "errs_nystromwr_l2", errs_nystromwr_l2)
		write(file, "timings_nystromwr", timings_nystromwr)
		write(file, "timings_test_nystromwr", timings_test_nystromwr)
		write(file, "timings_nystromwr_fro", timings_nystromwr_fro)
		write(file, "timings_nystromwr_l2", timings_nystromwr_l2)
		write(file, "errs_rff_nystrom_fro", errs_rff_nystrom_fro)
		write(file, "errs_rff_nystrom_l2", errs_rff_nystrom_l2)
		write(file, "timings_rff_nystrom", timings_rff_nystrom)
		write(file, "timings_test_rff_nystrom", timings_test_rff_nystrom)
		write(file, "timings_rff_nystrom_fro", timings_rff_nystrom_fro)
		write(file, "timings_rff_nystrom_l2", timings_rff_nystrom_l2)
		close(file)
	end

	file = matopen(@sprintf("fd_results_%s_n%d_nt%d_d%d_s%f.mat",name,n,nt,d,sigma),"w")
	write(file, "timing_orig", timing_orig)
	write(file, "dims", dims)
	write(file, "ells", ells)
	write(file, "epsilons", epsilons)
	write(file, "errs_fro", errs_fro)
	write(file, "errs_l2", errs_l2)
	write(file, "timings", timings)
	write(file, "timings_test", timings_test)
	write(file, "timings_fro", timings_fro)
	write(file, "timings_l2", timings_l2)
	write(file, "errs_rankk_fro", errs_rankk_fro)
	write(file, "errs_rankk_l2", errs_rankk_l2)
	write(file, "timings_rankk", timings_rankk)
	write(file, "timings_test_rankk", timings_test_rankk)
	write(file, "timings_rankk_fro", timings_rankk_fro)
	write(file, "timings_rankk_l2", timings_rankk_l2)
	write(file, "errs_FD_fro", errs_FD_fro)
	write(file, "errs_FD_l2", errs_FD_l2)
	write(file, "timings_FD", timings_FD)
	write(file, "timings_test_FD", timings_test_FD)
	write(file, "timings_FD_fro", timings_FD_fro)
	write(file, "timings_FD_l2", timings_FD_l2)
	write(file, "errs_nystrom_fro", errs_nystrom_fro)
	write(file, "errs_nystrom_l2", errs_nystrom_l2)
	write(file, "timings_nystrom", timings_nystrom)
	write(file, "timings_test_nystrom", timings_test_nystrom)
	write(file, "timings_nystrom_fro", timings_nystrom_fro)
	write(file, "timings_nystrom_l2", timings_nystrom_l2)
	write(file, "errs_nystromwr_fro", errs_nystromwr_fro)
	write(file, "errs_nystromwr_l2", errs_nystromwr_l2)
	write(file, "timings_nystromwr", timings_nystromwr)
	write(file, "timings_test_nystromwr", timings_test_nystromwr)
	write(file, "timings_nystromwr_fro", timings_nystromwr_fro)
	write(file, "timings_nystromwr_l2", timings_nystromwr_l2)
	write(file, "errs_rff_nystrom_fro", errs_rff_nystrom_fro)
	write(file, "errs_rff_nystrom_l2", errs_rff_nystrom_l2)
	write(file, "timings_rff_nystrom", timings_rff_nystrom)
	write(file, "timings_test_rff_nystrom", timings_test_rff_nystrom)
	write(file, "timings_rff_nystrom_fro", timings_rff_nystrom_fro)
	write(file, "timings_rff_nystrom_l2", timings_rff_nystrom_l2)
	close(file)

	#PyPlot.plt.show()

end

centerit = false
K = []
Ktest = []
timing_orig = 0
if recompute
	tic()

	# compute sigma from data:
	sigma = 0
	println("computing sigma parameter...")
	if n > 1e3
		ss = round(Int64,1e3)
		sampled = X[1:ss,:]
		for i=ss+1:n
			rr = round(Int64,ceil(rand()*n))
			if rr <= ss
				sampled[rr,:] = X[i,:]
			end
		end
		for i=1:ss
			sigma += sum(sqrt(sum((sampled - ones(ss,1)*sampled[i,:]).^2,2)))
		end
		sigma /= ss^2
	else
		for i=1:n
			sigma += sum(sqrt(sum((X - ones(n,1)*X[i,:]).^2,2)))
		end
		sigma /= n^2
	end

	println("computed sigma: ", sigma)

	println("computing kernel matrix...")
	K_test = GaussianKernel(Xtest_centered,sigma)
	nn = size(K,1)
	if centerit
		println("centering kernel matrix...")
		K_test = K_test - (1/nn)*ones(size(K))*K - K*(1/nn)*ones(size(K)) + (1/nn)*ones(size(K))*K*(1/nn)*ones(size(K))
	end
	timing_orig = toc()

	# write out K and timing:
	file = matopen(@sprintf("fd_results_K_%s_%d.mat",name,n),"w")
	write(file,"X",X)
	write(file,"X_centered",X_centered)
	write(file,"Xtest",Xtest)
	write(file,"Xtest_centered",Xtest_centered)
	write(file,"K_test",K_test)
	write(file,"sigma",sigma)
	write(file,"timing_orig",timing_orig)
	close(file)
else
	file = matopen(@sprintf("fd_results_K_%s_%d.mat",name,n))
	K_test = read(file,"K_test")
	X = read(file,"X")
	X_centered = read(file,"X_centered")
	Xtest = read(file,"Xtest")
	Xtest_centered = read(file,"Xtest_centered")
	sigma = read(file,"sigma")
	timing_orig = read(file,"timing_orig")
	println("timing to compute K: ", timing_orig)
	close(file)
end

comparison(X_centered, Xtest_centered, K_test, centerit, sigma)


text:
using MLKernels

SUITE["MLKernels"] = BenchmarkGroup()

mlkernelnames = ["SquaredExponentialKernel"]
kernels=Dict{String,MLKernels.Kernel}()
for k in mlkernelnames
    SUITE["MLKernels"][k] = BenchmarkGroup()
    kernels[k] = eval(Meta.parse("MLKernels."*k*"(alpha)"))
end

for k in mlkernelnames
    SUITE["MLKernels"][k]["k(X,Y)"] = @benchmarkable MLKernels.kernelmatrix($(kernels[k]),$X,$Y)
    # SUITE["MLKernels"][k][kt]["k!(X,Y)"] = @benchmarkable MLKernels.kernelmatrix!(KXY,$(kernels[k][kt]),$X,$Y) setup=(KXY=copy($KXY))
    SUITE["MLKernels"][k]["k(X)"] = @benchmarkable MLKernels.kernelmatrix($(kernels[k]),$X)
    # SUITE["MLKernels"][k][kt]["k!(X)"] = @benchmarkable MLKernels.kernelmatrix!(KX,$(kernels[k][kt]),$X) setup=(KX=copy($KX))
    # SUITE["MLKernels"][k][kt]["kdiag(X)"] = @benchmarkable MLKernels.kerneldiagmatrix($(kernels[k][kt]),$X)
    # SUITE["MLKernels"][k][kt]["kdiag!(X)"] = @benchmarkable MLKernels.kerneldiagmatrix!(kX,$(kernels[k][kt]),$X) setup=(kX=copy($kX))
end
# results = run(SUITE)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module AutoGrad
importall Base  # defining getindex, sin, etc.
export grad, check_grads, @primitive, @zerograd
export Value, getval, recorder, Grad, fixdomain  # these are required for the macros to work
datapath = joinpath(dirname(@__FILE__),"..","data")

include("core.jl")
include("util.jl")
include("interfaces.jl")
include("base/reduce.jl")
include("base/number.jl")
include("base/float.jl")
include("base/broadcast.jl")
include("base/math.jl")
include("base/abstractarray.jl")
include("base/arraymath.jl")
include("linalg/matmul.jl")
include("linalg/dense.jl")
include("linalg/generic.jl")
include("special/bessel.jl")
include("special/erf.jl")
include("special/gamma.jl")
include("special/trig.jl")

end # module


text:
@testset "base interface" begin

############################################################################################

@testset "method caching" begin

# #17057 fallout
@eval @noinline post17057_child(i) = i+1
@eval function post17057_parent(arr::Ptr{Int64})
    i = post17057_child(0)
    unsafe_store!(arr, i, i)
    return nothing
end

# bug 1: emit_invoke performed dynamic call due to NULL child function
#        (hooked module activation because of bug 2 below)
hook_module_activation(ref::Ptr{Void}) = nothing
hooks = Base.CodegenHooks(module_activation=hook_module_activation)
params = Base.CodegenParams(cached=false, runtime=false, hooks=hooks)
Base._dump_function(post17057_parent, Tuple{Ptr{Int64}},
                    #=native=#false, #=wrapper=#false, #=strip=#false,
                    #=dump_module=#true, #=syntax=#:att, #=optimize=#false,
                    params)

# bug 2: default module activation segfaulted on NULL child function if cached=false
params = Base.CodegenParams(cached=false)
Base._dump_function(post17057_parent, Tuple{Ptr{Int64}},
                    #=native=#false, #=wrapper=#false, #=strip=#false,
                    #=dump_module=#true, #=syntax=#:att, #=optimize=#false,
                    params)

end

############################################################################################

end


text:
# function Base.merge!(trace::GenericTrace, db)
#     for pair in db # pairs
#         trace[pair[1]] = pair[2]
#     end
# end

# TODO: Use plotly_trace(; label="WBHP")
# function plotly_trace(
#     label;
#     history_alpha = 1.0, model_alpha = 0.3,
#     #history_alpha = 0.3, model_alpha = 1.0,

#     history_mode = "lines", model_mode = "lines",
#     #history_mode = "lines", model_mode = "markers",
#     #history_mode = "markers", model_mode = "lines",
#     history_line_width = 3, model_line_width = 2,
#     water_color = "30,144,255", oil_color="0,128,0", bhp_color = "0,0,0",
#     db...
#     )

#     trace = nothing
#     if     "WBHP"  == label
#         trace = scatter(
#             mode = model_mode,
#             line = Dict(
#                 :color => "rgba($bhp_color,$model_alpha)",
#                 :width => model_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($bhp_color,$model_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($bhp_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y3",
#         )
#     elseif "WBHPH" == label
#         trace = scatter(
#             mode = history_mode,
#             line = Dict(
#                 :color => "rgba($bhp_color,$history_alpha)",
#                 #:width => history_line_width,
#                 :width => model_line_width, # no need to be thick
#             ),
#             marker = Dict(
#                 :color => "rgba($bhp_color,$history_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($bhp_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y3",
#         )
#     elseif "WOPR"  == label
#         trace = scatter(
#             mode = model_mode,
#             line = Dict(
#                 :shape => "vh",
#                 :color => "rgba($oil_color,$model_alpha)",
#                 :width => model_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($oil_color,$model_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($oil_color)",
#                     :width => 0,
#                     #:width => 1,
#                 ),
#             ),
#         )
#     elseif "WOPRH" == label
#         trace = scatter(
#             #mode = "markers",
#             mode = history_mode,
#             line = Dict(
#                 :shape => "vh",
#                 :color => "rgba($oil_color,$history_alpha)",
#                 :width => history_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($oil_color,$history_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($oil_color)",
#                     :width => 0,
#                     #:width => 1,
#                 ),
#             ),
#         )
#     elseif "WOPT"  == label
#         trace = scatter(
#             mode = model_mode,
#             line = Dict(
#                 :color => "rgba($oil_color,$model_alpha)",
#                 :width => model_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($oil_color,$model_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($oil_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y2",
#          )
#     elseif "WOPTH" == label
#         trace = scatter(
#             mode = history_mode,
#             line = Dict(
#                 :color => "rgba($oil_color,$history_alpha)",
#                 :width => history_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($oil_color,$history_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($oil_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y2",
#         )
#     elseif "WWPR"  == label
#         trace = scatter(
#             mode = model_mode,
#             line = Dict(
#                 :shape => "vh",
#                 :color => "rgba($water_color,$model_alpha)",
#                 :width => model_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($water_color,$model_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($water_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#         )
#     elseif "WWPRH" == label
#         trace = scatter(
#             mode = history_mode,
#             line = Dict(
#                 :color => "rgba($water_color,$history_alpha)",
#                 :width => history_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($water_color,$history_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($water_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#         )
#     elseif "WWPT"  == label
#         trace = scatter(
#             mode = model_mode,
#             line = Dict(
#                 :color => "rgba($water_color,$model_alpha)",
#                 :width => model_line_width,
#              ),
#             marker = Dict(
#                 :color => "rgba($water_color,$model_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($water_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y2",
#         )
#     elseif "WWPTH"  == label
#         trace = scatter(
#             mode = history_mode,
#             line = Dict(
#                 :color => "rgba($water_color,$history_alpha)",
#                 :width => history_line_width,
#             ),
#             marker = Dict(
#                 :color => "rgba($water_color,$history_alpha)",
#                 :size => 4,
#                 :line => Dict(
#                     :color => "rgb($water_color)",
#                     #:width => 1,
#                     :width => 0,
#                 ),
#             ),
#             yaxis = "y2",
#         )
#     else
#         throw(ErrorException("Unknown label: $label"))
#     end

#     for pair in db # pairs
#         trace[pair[1]] = pair[2]
#     end

#     return trace
# end

function plotly_trace(
    ;
    label = "WBHP",
    history_alpha = 1.0, model_alpha = 0.3,
    #history_alpha = 0.3, model_alpha = 1.0,

    history_mode = "lines", model_mode = "lines",
    #history_mode = "lines", model_mode = "markers",
    #history_mode = "markers", model_mode = "lines",
    history_line_width = 3, model_line_width = 2,
    water_color = "30,144,255", oil_color="0,128,0", bhp_color = "0,0,0",
    db...
    )

    trace = nothing
    if     "WBHP"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($bhp_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($bhp_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($bhp_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y3",
        )
    elseif "WBHPH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($bhp_color,$history_alpha)",
                #:width => history_line_width,
                :width => model_line_width, # no need to be thick
            ),
            marker = Dict(
                :color => "rgba($bhp_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($bhp_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y3",
        )
    elseif "WOPR"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($oil_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    :width => 0,
                    #:width => 1,
                ),
            ),
        )
    elseif "WOPRH" == label
        trace = scatter(
            #mode = "markers",
            mode = history_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($oil_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    :width => 0,
                    #:width => 1,
                ),
            ),
        )
        if "markers" == history_mode
            trace[:marker][:color] = "rgba(255,255,255,0)"
            trace[:marker][:size]  = 8
            trace[:marker][:line][:color] = "rgba($oil_color, $history_alpha)"
            trace[:marker][:line][:width] = 1
        elseif "lines" == history_mode
        end
    elseif "WOPT"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y2",
         )
    elseif "WOPTH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y2",
        )
    elseif "WWPR"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($water_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WWPRH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($water_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgba($water_color, $history_alpha)",
                    :width => 0,
                ),
            ),
        )
        if "markers" == history_mode
            trace[:marker][:color] = "rgba(255,255,255,0)"
            trace[:marker][:size]  = 8
            trace[:marker][:line][:color] = "rgba($water_color, $history_alpha)"
            trace[:marker][:line][:width] = 1
        elseif "lines" == history_mode
        end
    elseif "WWPT"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :width => model_line_width,
             ),
            marker = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y2",
        )
    elseif "WWPTH"  == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
            yaxis = "y2",
        )
    else
        throw(ErrorException("Unknown label: $label"))
    end

    for pair in db # pairs
        trace[pair[1]] = pair[2]
    end

    return trace
end

function plotly_layout()
    layout = Layout(
        #title = "Pad 6",
        #showlegend = false,
        showlegend = true,
        margin = Dict(
            #:pad => -10,
            #:b => 0,
        ),
        xaxis = Dict(
            :title => "Time",
            :dtick => 86400*1000*30, # unit is ms
            :ticklen => 5,
            :domain => [0.10, 0.90],
            :rangemode => "tozero",
            #:automargin => false,
            #:showspikes => false,
            #:spikemode => "toaxis",
            #:showline => true,
            #:tickformat => "%b %Y",
            #:tickvals = [1, 3, 5, 7, 9, 11],
            #:ticktext = ['One', 'Three', 'Five', 'Seven', 'Nine', 'Eleven']
        ),
        yaxis = Dict(
            :title => "Liquid Flowrate [STB/d]",
            # :title => Dict(
            #     :text => "Liquid Flowrate [STB/d]",
            #     :font => Dict(
            #       :size => 20,
            #     ),
            # ),
            #:range => [0, 35],
            :domain => [0.05, 0.95], # so html does not cut off at the bottom
            :showline => true,
            :ticklen => 5,
            :rangemode => "tozero",

        ),
        yaxis2 = Dict(
            :title => "Liquid Production Volume [STB]",
            # :title => Dict(
            #     :text => "Liquid Production Volume [STB]",
            #     :font => Dict(
            #       :size => 20,
            #     ),
            # ),
            #:range => [0, 2500],
            :rangemode => "tozero",
            :ticklen => 5,
            # titlefont: {color: '#ff7f0e'},
            # tickfont: {color: '#ff7f0e'},
            :anchor => "x", # "free"
            :overlaying => "y",
            :side => "right",
            :position => 1.0,
            :showline => true,
            :zeroline => false,
            :showgrid => false,
        ),
        yaxis3 = Dict(
            :title => "Pressure [psi]",
            # :title => Dict(
            #     :text => "Pressure [psi]",
            #     :font => Dict(
            #       :size => 20,
            #     ),
            # ),
            #:range => [0, 11000],
            :rangemode => "tozero",
            :ticklen => 5,
            :tickformat => ",",
            # titlefont: {color: '#ff7f0e'},
            # tickfont: {color: '#ff7f0e'},
            :anchor => "free", # "free"
            :overlaying => "y",
            :side => "left",
            :position => 0.05,
            :showline => true,
            :zeroline => false,
            :showgrid => false,
        ),
        yaxis4 = Dict(
            # :title => Dict(
            #     :text => "Pressure [psi]",
            #     :font => Dict(
            #       :size => 20,
            #     ),
            # ),
            #:range => [0, 11000],
            :rangemode => "tozero",
            :ticklen => 5,
            :tickformat => ",",
            # titlefont: {color: '#ff7f0e'},
            # tickfont: {color: '#ff7f0e'},
            :anchor => "free", # "free"
            :overlaying => "y",
            :side => "left",
            :position => 1.0,
            :showline => true,
            :zeroline => false,
            :showgrid => false,
        ),
    )

    return layout
end

function axis_time(; nargs...)
    db = attr(
        title      = "Time",
        type       = "date",
        tickmode   = "auto",
        nticks     = 12,
        dtick      = "M12",
        ticklen    = 5,
        rangemode  = "tozero",
        automargin = true,
    )

    for pair in nargs # pairs
        db[pair[1]] = pair[2]
    end

    return db
end

function axis_bhp(; nargs...)
    db = attr(
        title = "BHP [psi]",
        ticklen    = 5,
        rangemode  = "tozero",
        automargin = true,
    )

    for pair in nargs # pairs
        db[pair[1]] = pair[2]
    end

    return db
end

function axis_oil_rate(; nargs...)
    db = attr(
        title = "Oil Production Rate, STB/D",
        ticklen    = 5,
        rangemode  = "tozero",
        automargin = true,
    )

    for pair in nargs # pairs
        db[pair[1]] = pair[2]
    end

    return db
end

function axis_water_rate(; nargs...)
    db = attr(
        title = "Water Production Rate, STB/D",
        ticklen    = 5,
        rangemode  = "tozero",
        automargin = true,
    )

    for pair in nargs # pairs
        db[pair[1]] = pair[2]
    end

    return db
end

# Not used
function subplot_trace(; label="WBHP",
    history_alpha = 1.0, model_alpha = 0.3,
    history_mode = "lines", model_mode = "lines",
    #history_mode = "lines", model_mode = "markers",
    #history_mode = "markers", model_mode = "lines",
    history_line_width = 3, model_line_width = 2,
    water_color = "30,144,255", oil_color="0,128,0", bhp_color = "0,0,0",
    db...)
    trace = nothing
    if     "WBHP"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($bhp_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($bhp_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($bhp_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WBHPH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($bhp_color,$history_alpha)",
                #:width => history_line_width,
                :width => model_line_width, # no need to be thick
            ),
            marker = Dict(
                :color => "rgba($bhp_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($bhp_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WOPR"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($oil_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    :width => 0,
                    #:width => 1,
                ),
            ),
        )
    elseif "WOPRH" == label
        trace = scatter(
            #mode = "markers",
            mode = history_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($oil_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    :width => 0,
                    #:width => 1,
                ),
            ),
        )
    elseif "WOPT"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
         )
    elseif "WOPTH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($oil_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($oil_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WWPR"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :shape => "vh",
                :color => "rgba($water_color,$model_alpha)",
                :width => model_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WWPRH" == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WWPT"  == label
        trace = scatter(
            mode = model_mode,
            line = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :width => model_line_width,
             ),
            marker = Dict(
                :color => "rgba($water_color,$model_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    elseif "WWPTH"  == label
        trace = scatter(
            mode = history_mode,
            line = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :width => history_line_width,
            ),
            marker = Dict(
                :color => "rgba($water_color,$history_alpha)",
                :size => 4,
                :line => Dict(
                    :color => "rgb($water_color)",
                    #:width => 1,
                    :width => 0,
                ),
            ),
        )
    else
        throw(ErrorException("Unknown label: $label"))
    end

    for pair in db # pairs
        trace[pair[1]] = pair[2]
    end

    return trace
end

# deprecated

function wbhp_trace(;alpha=0.2, db...)
    trace = scatter(
        mode = "lines",
        line = Dict(
            :color => "rgba(0,0,0,$alpha)",
        ),
        yaxis = "y3",
    )

    merge!(trace, db)

    return trace
end

function wbhph_trace(;alpha=1, db...)
    trace = scatter(
        mode = "markers",
        marker = Dict(
            :color => "rgba(0,0,0,$alpha)",
            :size => 4,
            :line => Dict(
                :color => "rgb(0,0,0)",
                #:width => 1,
                :width => 0,
            ),
        ),
        yaxis = "y3",
    )

    merge!(trace, db)

    return trace
end

"""
    wopr_trace(;db...)
Custom style for well-oil-production-rate

```julia
wopr_trace(name="ADK-1972", x=xvec, y=yvec)
```
"""
function wopr_trace(; alpha=0.2, db...)
    trace = scatter(
        mode = "lines",
        line = Dict(
            :shape => "vh",
            :color => "rgba(0,128,0,$alpha)",
        ),
    )

    merge!(trace, db)

    return trace
end

function woprh_trace(; alpha=1, db...)
    trace = scatter(
        mode = "markers",
        marker = Dict(
            :color => "rgba(0,128,0,$alpha)",
            :size => 4,
            :line => Dict(
                :color => "rgb(0,128,0)",
                :width => 0,
                #:width => 1,
            ),
        ),
    )

    merge!(trace, db)

    return trace
end

function wopt_trace(; alpha=0.2, db...)
    trace = scatter(
       mode = "lines",
       line = Dict(
           :color => "rgba(0,128,0,$alpha)",
        ),
       yaxis = "y2",
    )

    merge!(trace, db)

    return trace
end

function wopth_trace(; alpha=1, db...)
    trace = scatter(
        mode = "markers",
        marker = Dict(
            :color => "rgba(0,128,0,$alpha)",
            :size => 4,
            :line => Dict(
                :color => "rgb(0,128,0)",
                #:width => 1,
                :width => 0,
            ),
        ),
        yaxis = "y2",
    )

    merge!(trace, db)

    return trace
end

function wwpr_trace(; alpha=0.2, db...)
    trace = scatter(
       mode = "lines",
       line = Dict(
           :shape => "vh",
           :color => "rgba(30,144,255,$alpha)",
        ),
    )

    merge!(trace, db)

    return trace
end

function wwprh_trace(; alpha=1, db...)
    trace = scatter(
        mode = "markers",
        marker = Dict(
            :color => "rgba(30,144,255,$alpha)",
            :size => 4,
            :line => Dict(
                :color => "rgb(30,144,255)",
                #:width => 1,
                :width => 0,
            ),
        ),
    )

    merge!(trace, db)

    return trace
end


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

mutable struct IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule <: SwaggerModel
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPort} } # spec name: ports
    to::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPeer} } # spec name: to

    function IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule(;ports=nothing, to=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule, Symbol("to"), to)
        setfield!(o, Symbol("to"), to)
        o
    end
end # type IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule

const _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule = Dict{Symbol,Symbol}(Symbol("ports")=>Symbol("ports"), Symbol("to")=>Symbol("to"))
const _property_types_IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule = Dict{Symbol,String}(Symbol("ports")=>"Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPort}", Symbol("to")=>"Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPeer}")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1NetworkPolicyEgressRule)
    true
end

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


text:
println("Hello World. This is JuliaLang, a language designed for high performance.")

println("It has the speed of C/C++ and the ease of Python. It is revolutionizing scientific coding.")


text:
version = UInt8(131)

# See http://erlang.org/doc/apps/erts/erl_ext_dist.html, tags are in the same
# order.

# Compressed tag is mentioned in introduction
# (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id86906)
compressedtag = UInt8(80)

# ATOM_CACHE_REF (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id87859)
atomcachereftag = UInt8(82)

# SMALL_INTEGER_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id87935)
smallinttag = UInt8(97)

# INTEGER_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id87998)
inttag = UInt8(98)

# FLOAT_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88062)
floattag = UInt8(99)

# ATOM_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88138)
atomtag = UInt8(100)

# REFERENCE_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88237)
reftag = UInt8(101)

# PORT_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88395)
porttag = UInt8(102)

# PID_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88514)
pidtag = UInt8(103)

# SMALL_TUPLE_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88658)
smalltupletag = UInt8(104)

# LARGE_TUPLE_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88750)
largetupletag = UInt8(105)

# MAP_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88840)
maptag = UInt8(116)

# NIL_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id88950)
niltag = UInt8(106)

# STRING_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89001)
stringtag = UInt8(107)

# LIST_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89094)
listtag = UInt8(108)

# BINARY_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89212)
bintag = UInt8(109)

# SMALL_BIG_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89309)
smallbiginttag = UInt8(110)

# LARGE_BIG_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89416)
largebiginttag = UInt8(111)

# NEW_REFERENCE_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89519)
newreftag = UInt8(114)

# SMALL_ATOM_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89692)
smallatomtag = UInt8(115)

# FUN_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id89811)
funtag = UInt8(117)

# NEW_FUN_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90058)
newfuntag = UInt8(112)

# EXPORT_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90435)
exporttag = UInt8(113)

# BIT_BINARY_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90579)
bitbintag = UInt8(77)

# NEW_FLOAT_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90685)
newfloattag = UInt8(70)

# ATOM_UTF8_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90753)
atomutf8tag = UInt8(118)

# SMALL_ATOM_UTF8_EXT (http://erlang.org/doc/apps/erts/erl_ext_dist.html#id90855)
smallatomutf8tag = UInt8(119)


text:
##########################################################################
##########################################################################
## Reproduction code for: Optimal Climate Policy When Damages are Unknown
## Ivan Rudik
## American Economic Journal: Economic Policy
#
# Julia version: 1.1.0
# Operating System: Linux CentOS
#
# Required packages:
# BSON, Roots, DataFrames, CSV, Distributed, BasisMatrices,
# NLopt, Distributions, Dierckx, LinearAlgebra, CompEcon

# 1. All functions in MAIN_SCRIPT.jl must be called in the
# order they appear in the code
# 2. The generate_results() call on line 202 is
# what creates the final results
#########################################################################
#########################################################################

# Load up packages
using Roots, DataFrames, CSV, Distributed, CompEcon
@everywhere using CompEcon
@everywhere using BSON, BasisMatrices, NLopt, Distributions, Dierckx, LinearAlgebra

# Include required .jl files to bring in model functions
@everywhere include("custom_types.jl")
@everywhere include("tools.jl")
@everywhere include("maxbell.jl")
@everywhere include("simulation.jl")
@everywhere include("transitions.jl")
@everywhere include("initialization.jl")
@everywhere include("tax_calc.jl")
@everywhere include("adaptive_grid.jl")
@everywhere include("model.jl")
@everywhere include("model_solvers.jl")
@everywhere include("RESULTS.jl")

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!! SOLVE MODELS !!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

##################################
##################################
# Base models
# Figures 5a, 7, A2, A3
# Tables 2, 3, A4
##################################
##################################

### Solve uncertainty framework
solve_uncert()

### Solve learning framework
solve_learn()

### Solve RC framework
solve_rc(3.9)

### Solve RC+L framework
solve_rcl(3.9)

##################################
##################################
# Base models + alternate damage parameters
# Figures 5b, 5c, 5d
##################################
##################################

simulate_alt_damage_params("results/uncert_results.bson", 0.01, 1.88)
simulate_alt_damage_params("results/uncert_results.bson", .00556486, 3.0)
simulate_alt_damage_params("results/uncert_results.bson", 0.01, 3.0)
simulate_alt_damage_params("results/learning_results_adapt_20.bson", 0.01, 1.88)
simulate_alt_damage_params("results/learning_results_adapt_20.bson", .00556486, 3.0)
simulate_alt_damage_params("results/learning_results_adapt_20.bson", 0.01, 3.0)
simulate_alt_damage_params("results/rc_results_3.9.bson", 0.01, 1.88)
simulate_alt_damage_params("results/rc_results_3.9.bson", .00556486, 3.0)
simulate_alt_damage_params("results/rc_results_3.9.bson", 0.01, 3.0)
simulate_alt_damage_params("results/rcl_results_3.9_adapt_20.bson", 0.01, 1.88)
simulate_alt_damage_params("results/rcl_results_3.9_adapt_20.bson", .00556486, 3.0)
simulate_alt_damage_params("results/rcl_results_3.9_adapt_20.bson", 0.01, 3.0)

##################################
##################################
# Alternate penalty parameter
# Figures 6a and 6b
##################################
##################################

### RC
solve_alt_penalty_parameter_rc()

### RC+L
solve_alt_penalty_parameter_rcl()

##################################
##################################
# Misspecified damage simulations
# Figure 8
##################################
##################################

simulate_misspecified_damages("results/uncert_results.bson", :weitzman)
simulate_misspecified_damages("results/uncert_results.bson", :dietz)
simulate_misspecified_damages("results/uncert_results.bson", :extreme)
simulate_misspecified_damages("results/learning_results_adapt_20.bson", :weitzman)
simulate_misspecified_damages("results/learning_results_adapt_20.bson", :dietz)
simulate_misspecified_damages("results/learning_results_adapt_20.bson", :extreme)
simulate_misspecified_damages("results/rc_results_3.9.bson", :weitzman)
simulate_misspecified_damages("results/rc_results_3.9.bson", :dietz)
simulate_misspecified_damages("results/rc_results_3.9.bson", :extreme)
simulate_misspecified_damages("results/rcl_results_3.9_adapt_20.bson", :weitzman)
simulate_misspecified_damages("results/rcl_results_3.9_adapt_20.bson", :dietz)
simulate_misspecified_damages("results/rcl_results_3.9_adapt_20.bson", :extreme)

##################################
##################################
# Previous iteration state bounds
# Figure A2
##################################
##################################

solve_learn_previous_iter()

##################################
##################################
# Error checks
# Table A4
##################################
##################################

##################################
# Collocation points: Table A4 top panel

### Uncertainty: collocation points
solve_uncert_error_coll()

### Learning: collocation points
solve_learn_error_coll()

### RC: collocation points
solve_rc_error_coll()

### RC+L: collocation points
solve_rcl_error_coll()

##################################
# Quadrature points: Table A4 second panel

### Uncertainty: quadrature points
solve_uncert(quad = 13)

### Learning: quadrature points
solve_learn(load = true, quad = 13)

### RC: quadrature points
solve_rc(3.9, quad = 13)

### RC+L: quadrature points
solve_rcl(3.9, quad = 13)

##################################
# Expand domain: Table A4 third panel

### Uncertainty: expand domain
solve_uncert(expand = 1.1)

### RC: expand domain
solve_rc(3.9, expand = 1.1)

##################################
# Reduce adaptive grid iterations: Table A4 bottom panel

### Learning: adaptive iterations
solve_learn(adapt_iters = 10)

### RC+L: adaptive iterations
solve_rcl(3.9, adapt_iters = 10)

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!! GENERATE RESULTS !!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

generate_results()

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
using DataStructures

struct SpaceIndex
    value::Int
end

struct Spaces
    indices::IntDisjointSets # Union-Find datastructure
    dimensions::Dict{Int, Int}
    polyvars::Dict{Int, Vector{SpaceVariable}}
end
Spaces() = Spaces(IntDisjointSets(0), Dict{Int, Int}(),
                  Dict{Int, Vector{SpaceVariable}}())
function new_space(spaces::Spaces)
    return SpaceIndex(push!(spaces.indices))
end
function new_space(spaces::Spaces,
                   dim::Int)
    space_index = new_space(spaces)
    spaces.dimensions[space_index.value] = dim
    return space_index
end
function new_space(spaces::Spaces,
                   polyvars::Vector{SpaceVariable})
    space_index = new_space(spaces)
    spaces.dimensions[space_index.value] = length(polyvars)
    spaces.polyvars[space_index.value] = polyvars
    return space_index
end
function merge_property(d::Dict, root, key1, key2, name)
    value = nothing
    if key1 in keys(d)
        v1 = d[key1]
        if key2 in keys(d) && v1 != d[key2]
            error("Sets lie on the same spaces but have different $name")
        end
        value = v1
    elseif key2 in keys(d)
        value = d[key2]
    end
    if value !== nothing
        d[root] = value
    end
end
function merge_spaces(spaces::Spaces, a::SpaceIndex, b::SpaceIndex)
    root_a = find_root!(spaces.indices, a.value)
    root_b = find_root!(spaces.indices, b.value)
    root = root_union!(spaces.indices, root_a, root_b)
    # the properties my not be set at a.value when it is set at root_a.value
    merge_property(spaces.dimensions, root, root_a, root_b, "dimension")
    merge_property(spaces.polyvars, root, root_a, root_b,
                   "polynomial variables")
    return SpaceIndex(root)
end
function space_dimension(spaces::Spaces, si::SpaceIndex)
    idx = find_root!(spaces.indices, si.value)
    if !haskey(spaces.dimensions, idx)
        error("Missing dimension information, use Ellipsoid(dimension=...) or PolySet(dimension=...)")
    end
    return spaces.dimensions[idx]
end
function space_polyvars(spaces::Spaces, si::SpaceIndex)
    idx = find_root!(spaces.indices, si.value)
    if !haskey(spaces.polyvars, idx)
        dim = space_dimension(spaces, si)
        DynamicPolynomials.@polyvar x[1:dim]
        spaces.polyvars[idx] = x
    end
    return spaces.polyvars[idx]
end


text:
using BinaryBuilder, Pkg

name = "libtree"
version = v"3.0.1"

sources = [
    ArchiveSource("https://github.com/haampie/libtree/archive/refs/tags/v$(version).tar.gz", "20d3cd66f5c74058de9dd594af8ffd639c795d27ab435c588a3cd43911c1604f")
]

script = raw"""
cd $WORKSPACE/srcdir/libtree-*/
make CFLAGS="-Os -fwhole-program" LDFLAGS="-Wl,-s" "PREFIX=$prefix" install
"""

# Build only on platforms where ELF objects are usually used.
platforms = filter!(p -> Sys.islinux(p) || Sys.isfreebsd(p), supported_platforms(; experimental=true))

products = [
    ExecutableProduct("libtree", :libtree)
]

dependencies = Dependency[]

# Note: binutils 2.24 has issues with -s and --gc-sections, so use GCC 5 which comes wth a later binutils.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# (@v1.8) pkg> add WebSockets
using WebSockets

function main()
    WebSockets.open("wss://wsaws.okex.com:8443/ws/v5/public") do ws   # 1
        for i = 1:100
            a = time_ns()                   # 2
            write(ws, "ping")               # 3
            data, success = readguarded(ws) # 4
            !success && break               # 5
            b = time_ns()
            println(String(data), " ", (b - a) / 1000000) # 6
            sleep(0.1)
        end
    end
end

main()


text:
mutable struct MCMC{S <: DiseaseStateSequence, M <: ILM}
  event_observations::EventObservations{S, M}
  event_extents::EventExtents{S}
  population::Population
  starting_states::DiseaseStates
  risk_functions::RiskFunctions{S}
  risk_priors::RiskPriors{S}
  substitution_model::Union{Nothing, Type{<:NucleicAcidSubstitutionModel}}
  substitution_model_prior::Union{Nothing, Vector{UnivariateDistribution}}
  transmission_network_prior::Union{Nothing, TNPrior}
  markov_chains::Vector{MarkovChain{S, M}}

  function MCMC(
    obs::EventObservations{S, M},
    ee::EventExtents{S},
    pop::Population,
    states::DiseaseStates,
    rf::RiskFunctions{S},
    rp::RiskPriors{S};
    tnprior::Union{Nothing, TNPrior}=nothing) where {
    S <: DiseaseStateSequence,
    M <: TNILM}
    return new{S, M}(obs, ee, pop, states, rf, rp, nothing, nothing, tnprior, MarkovChain{S, M}[])
  end

  function MCMC(
    obs::EventObservations{S, M},
    ee::EventExtents{S},
    pop::Population,
    states::DiseaseStates,
    rf::RiskFunctions{S},
    rp::RiskPriors{S},
    sm::Type{N},
    smp::Vector{UnivariateDistribution};
    tnprior::Union{Nothing, TNPrior}=nothing) where {
    S <: DiseaseStateSequence,
    M <: PhyloILM,
    N <: NucleicAcidSubstitutionModel}
    return new{S, M}(obs, ee, pop, states, rf, rp, sm, smp, tnprior, MarkovChain{S, M}[])
  end
end

function MCMC(
  obs::EventObservations{S, M},
  ee::EventExtents{S},
  pop::Population,
  rf::RiskFunctions{S},
  rp::RiskPriors{S};
  tnprior::Union{Nothing, TNPrior}=nothing) where {
  S <: DiseaseStateSequence,
  M <: TNILM}
  return MCMC(obs, ee, pop, fill(State_S, individuals(pop)), rf, rp, tnprior=tnprior)
end

function MCMC(
  obs::EventObservations{S, M},
  ee::EventExtents{S},
  pop::Population,
  rf::RiskFunctions{S},
  rp::RiskPriors{S},
  sm::Type{N},
  smp::Vector{UnivariateDistribution};
  tnprior::Union{Nothing, TNPrior}=nothing) where {
  S <: DiseaseStateSequence,
  M <: PhyloILM,
  N <: NucleicAcidSubstitutionModel}
  return MCMC(obs, ee, pop, fill(State_S, individuals(pop)), rf, rp, sm, smp, tnprior=tnprior)
end

function Base.show(io::IO, x::MCMC{S, M}) where {
                   S <: DiseaseStateSequence,
                   M <: ILM}
  return print(io, "$S $M MCMC with $(length(x.markov_chains)) Markov chains")
end

function TNDistribution(x::MCMC; burnin::Int64=0, thin::Int64=1)
  return TNDistribution([TNDistribution(y, burnin=burnin, thin=thin) for y in x.markov_chains])
end


text:
# Obtain the occurrance matrix associated with a UMLS concept in a previously
# obtained pubmed/medline search
# Date: May 6, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using ArgParse
using PubMedMiner
using SQLite
using BioMedQuery.UMLS:Credentials
using JLD
using MySQL

function main(args)

    # initialize the settings (the description is for the help screen)
   s = ArgParseSettings()
   s.description = "This is pubMedOccur"
   s.version = "1.0"

   @add_arg_table s begin
    "--umls_concept"
        help = "UMLS concept for occurrance analysis"
        arg_type = ASCIIString
        required = true
    "--results_dir"
         help = "Path to store the results"
         arg_type = ASCIIString
         required = true
     "mysql"
         help = "Use MySQL backend"
         action = :command
     "sqlite"
          help = "Use SQLite backend"
          action = :command
   end
   @add_arg_table s["mysql"] begin
    "--host"
        help = "Host where you database lives"
        arg_type = ASCIIString
        default = "localhost"
    "--dbname"
        help = "Database name"
        arg_type = ASCIIString
        required = true
    "--username"
        help = "MySQL username"
        arg_type = ASCIIString
        default = "root"
    "--password"
        help = "MySQL password"
        arg_type = ASCIIString
        default = ""
   end

    @add_arg_table s["sqlite"] begin
    "--db_path"
         help = "Path to SQLite database file to store results"
         arg_type = ASCIIString
         required = true
    end

   parsed_args = parse_args(s) # the result is a Dict{String,Any}
   println("-------------------------------------------------------------")
   println(s.description)
   println("Parsed args:")
   for (key,val) in parsed_args
       println("  $key  =>  $(repr(val))")
   end
   println("-------------------------------------------------------------")

   results_dir = parsed_args["results_dir"]

    if !isdir(results_dir)
        mkdir(results_dir)
    end

    occur_path = results_dir*"/occur_sp.jdl"
    labels2ind_path = results_dir*"/labels2ind.jdl"

    umls_concept = parsed_args["umls_concept"]

    db = nothing
    if haskey(parsed_args, "sqlite")
        db_path  = parsed_args["sqlite"]["db_path"]
        db = SQLite.DB(db_path)
    elseif haskey(parsed_args, "mysql")
        host = parsed_args["mysql"]["host"]
        dbname = parsed_args["mysql"]["dbname"]
        username = parsed_args["mysql"]["username"]
        pswd = parsed_args["mysql"]["password"]
        db = MySQL.connect(host, username, pswd, dbname)
    else
        error("Unsupported database backend")
    end

    @time begin
        labels2ind, occur = PubMedMiner.occurance_matrix(db, umls_concept)
    end

    println("-------------------------------------------------------------")
    println("Output Descritor to Index Dictionary")
    println(labels2ind)
    println("-------------------------------------------------------------")

    println("-------------------------------------------------------------")
    println("Output Data Matrix")
    println(occur)
    println("-------------------------------------------------------------")

    # save(occur_path, "occur", occur)
    jldopen(occur_path, "w") do file
        write(file, "occur", occur)
    end
    jldopen(labels2ind_path, "w") do file
        write(file, "labels2ind", labels2ind)
    end

    # file  = jldopen(occur_path, "r")
    # obj2 = read(file, "occur")
    # display(obj2)

    println("-------------------------------------------------------------")
    println("Done computing and saving occurance info to disk")
    println("-------------------------------------------------------------")

end

main(ARGS)


text:
using Printf, Test, Main.MPCC

x0 = ones(6)
m0 = ones(2)
c(x) = [sum(x)]
f1(x) = sum(x[1:6].^2) - sum(x[7:8].^2)
f2(x) = sum(x[7:8].^2)
g(x) = []
ll_lvar = zeros(2)
ll_uvar = ones(2)
ll_lcon, ll_ucon = zeros(0), zeros(0)

bpmpcc = BPMPCCModel(f1,  x0, m0,
                     lvar = fill(-10.0,size(x0)), uvar = fill(10.0,size(x0)),
                     y0 = [0.0], c = c, lcon = [-Inf], ucon = [6.],
                     f2 = f2, g = g, ll_lvar = ll_lvar, ll_uvar = ll_uvar,
                     ll_lcon = ll_lcon, ll_ucon = ll_ucon, no_third_derivative = true)

@test bpmpcc.n == 6
@test bpmpcc.m == 2
@test bpmpcc.llncon == 0

@test bpmpcc.meta.nvar == 12
@test length(bpmpcc.meta.lvar) == 12
@test length(bpmpcc.meta.uvar) == 12
@test bpmpcc.meta.lvar[1:6] == fill(-10.0,size(x0))
@test bpmpcc.meta.uvar[1:6] == fill( 10.0,size(x0))
@test bpmpcc.meta.lvar[7:8] == -Inf * ones(2)
@test bpmpcc.meta.uvar[7:8] ==  Inf * ones(2)

@test bpmpcc.meta.ncon == 3
@test length(bpmpcc.meta.lcon) == 3
@test length(bpmpcc.meta.ucon) == 3
@test length(bpmpcc.meta.y0)   == 3

@test bpmpcc.meta.ncc  == 2 * 2
@test bpmpcc.meta.lccG == zeros(bpmpcc.meta.ncc)
@test bpmpcc.meta.lccG == bpmpcc.meta.lccH

t0 = rand(12)
x0, y0, l1, l2, m1, m2 = Main.MPCC._var_bp(bpmpcc, t0)
@test (x0 == t0[1:6]) && (y0 == t0[7:8]) && (l1 == zeros(0)) && (l2 == zeros(0)) && (m1 == t0[9:10]) && (m2 == t0[11:12])
@test obj(bpmpcc, t0) == f1(t0[1:8])
@test grad(bpmpcc, t0)[1:6]  ==  2*t0[1:6]
@test grad(bpmpcc, t0)[7:8]  == -2*t0[7:8]
@test grad(bpmpcc, t0)[9:12] == zeros(4)

@test length(cons_nl(bpmpcc, t0)) == 3 #d'autres tests?
@test consG(bpmpcc, t0) == vcat(t0[7:8] - ll_lvar, ll_uvar - t0[7:8])
@test consH(bpmpcc, t0) == t0[9:12]

@test size(jac_nl(bpmpcc, t0)) == (3,12)
@test jnlprod(bpmpcc, t0, zeros(12)) == zeros(3)
@test jnltprod(bpmpcc, t0, zeros(3)) == zeros(12)
@test size(jacG(bpmpcc, t0)) == (4,12)
@test jGprod(bpmpcc, t0, zeros(12)) == zeros(4)
@test jGtprod(bpmpcc, t0, zeros(4)) == zeros(12)
@test size(jacH(bpmpcc, t0)) == (4,12)
@test jHprod(bpmpcc, t0, zeros(12)) == zeros(4)
@test jHtprod(bpmpcc, t0, zeros(4)) == zeros(12)

@test hessG(bpmpcc, t0, zeros(4)) == zeros(12,12)
@test hessH(bpmpcc, t0, zeros(4)) == zeros(12,12)
@test hGprod(bpmpcc, t0, zeros(4), t0) == zeros(12)
@test hHprod(bpmpcc, t0, zeros(4), t0) == zeros(12)

@test hess(bpmpcc, t0, zeros(7), obj_weight = 0.0) == zeros(12,12)
@test hess(bpmpcc, t0, obj_weight = 0.0) == zeros(12,12)
@test hprod(bpmpcc, t0, zeros(7), zeros(12), obj_weight = 0.0) == zeros(12)
@test hprod(bpmpcc, t0, zeros(12), obj_weight = 0.0) == zeros(12)

@test hessH_structure(bpmpcc) == hessG_structure(bpmpcc)
@test length(hessH_coord(bpmpcc, bpmpcc.meta.x0, zeros(2))) == bpmpcc.meta.nvar * (bpmpcc.meta.nvar + 1) / 2
@test length(hessG_coord(bpmpcc, bpmpcc.meta.x0, zeros(2))) == bpmpcc.meta.nvar * (bpmpcc.meta.nvar + 1) / 2


text:
pushfirst!(LOAD_PATH, joinpath("..", ".."))

using Plots
using Printf
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize

column_ensemble_size = ColumnEnsembleSize(Nz=32, ensemble=(2, 1))

grid = RectilinearGrid(size=column_ensemble_size, z=(-64, 0), topology=(Flat, Flat, Bounded))

closure_1 = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1, background_κz = 1e-5, convective_νz = 1e-3, background_νz = 1e-4)
closure_2 = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1, background_κz = 1e-5, convective_νz = 1e-1, background_νz = 1e-4)

closures = Matrix{typeof(closure_1)}(undef, 2, 1)
closures[1] = closure_1
closures[2] = closure_2

Qᵇ = +1e-8
Qᵘ = +1e-5

b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = :b,
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = FPlane(f=1e-4),
                                    boundary_conditions = (; b=b_bcs, u=u_bcs),
                                    closure = closures)

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

z = znodes(model.tracers.b)

b1 = view(interior(model.tracers.b),    1, 1, :)
u1 = view(interior(model.velocities.u), 1, 1, :)
v1 = view(interior(model.velocities.v), 1, 1, :)

b2 = view(interior(model.tracers.b),    2, 1, :)
u2 = view(interior(model.velocities.u), 2, 1, :)
v2 = view(interior(model.velocities.v), 2, 1, :)

b_plot = plot(b1, z, linewidth = 2, label = "t = 0", xlabel = "Buoyancy", ylabel = "z", legend = :bottomright)

u_plot = plot(u1, z, linewidth = 2, label = "u at t = 0",
              linestyle = :solid, color = :black, xlabel = "Velocity", ylabel = "z", legend = :bottomright)

plot!(u_plot, v1, z, linewidth = 2, label = "v at t = 0", linestyle = :dash, color = :black)

simulation = Simulation(model, Δt = 20.0, stop_time = 48hours)

run!(simulation)

t = prettytime(model.clock.time)

plot!(b_plot, b1, z, linewidth = 2, label = @sprintf("b₁ at t = %s", t))
plot!(b_plot, b2, z, linewidth = 2, label = @sprintf("b₂ at t = %s", t))

plot!(u_plot, u1, z, label = @sprintf("u₁ at t = %s", t), linewidth = 2, linestyle = :dash,  color = :blue)
plot!(u_plot, v1, z, label = @sprintf("v₁ at t = %s", t), linewidth = 2, linestyle = :solid, color = :blue)

plot!(u_plot, u2, z, label = @sprintf("u₂ at t = %s", t), linewidth = 2, linestyle = :dash,  color = :red)
plot!(u_plot, v2, z, label = @sprintf("v₂ at t = %s", t), linewidth = 2, linestyle = :solid, color = :red)

ub_plot = plot(b_plot, u_plot, layout=(1, 2))

display(ub_plot)


text:
using Ipopt, JuMP
import JuMP.MathOptInterface
const MOI = MathOptInterface
using Plots

#include("../experiment/synthetic/synthetic_btg.jl")
include("tdist_mle.jl")

using Ipopt

function optimize_lambda(fun, lambda_range)
    model = Model(Ipopt.Optimizer)
    #lx = hcat([0.01 for i = 1:1], [-5])
    #ux = hcat([100 for i = 1:1], [5])
    lx = [lambda_range[1]]
    ux = [lambda_range[2]]
    register(model, :fun, 1, fun, autodiff = true)
    @variable(model, lx[1] <= au[i=1:1] <= ux[1])
    set_start_value(au[1], 0.0001)
    @NLobjective(model, Min, fun(au[1]))
    JuMP.optimize!(model)
    return value.(au)
end

"""
MLE estimate for single dimensional theta and lambda.
    INPUTS:
        - search bounds
        - 2-variable function fun
    THIS FUNCTION GIVES A WEIRD ERROR.
"""
function optimize_theta_lambda_single_broken(fun, theta_range, lambda_range)
    model = Model(Ipopt.Optimizer)
    lx = [theta_range[1], lambda_range[1]]
    ux = [theta_range[2], lambda_range[2]]
    register(model, :func, 2, fun, autodiff = true)
    @variable(model, lx[i] <= au[i=1:2] <= ux[i])
    set_start_value(au[1], 100.0)
    set_start_value(au[2], 0.001)
    @NLobjective(model, Min, func(au[1], au[2]))
    JuMP.optimize!(model)
    return value.(au)
end

using Optim
function optimize_theta_lambda_single(fun, theta_range, lambda_range)
    x0 = [100.0, 0.01]
    optimize(fun, x0)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
#Demo 1: Transmission line test
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_MAGDB = "Magnitude (dB)"
LBLAX_FREQ = "Frequency (Hz)"
color1 = cons(:a, line = set(color=:red, width=2))
color2 = cons(:a, line = set(color=:blue, width=2))

#==Input data
===============================================================================#
z0 = 50 #Ohms (reference impedance)
ℓ = 5e-3 #Meters
zC = 200 #Ohms (line impedance)
μ0 = 4pi*1e-7 #F/m
ϵ0 = 8.854e-12 #H/m
f = collect(1:.1:100)*1e9

#==Computations
===============================================================================#
ω = 2pi*f
α = 0
β = ω*sqrt(μ0*ϵ0)
γ = α .+ im*β

#Convert γ type to DataF1, function of 1 argument (leverage C-Data toolkit):
γ = DataF1(f, γ)

#ABCD matrix:
A = cosh(γ*ℓ)      ;B = sinh(γ*ℓ)*zC;
C = sinh(γ*ℓ)/zC   ;D = cosh(γ*ℓ);

T = Network(:ABCD, [A B; C D])
S = Network(:S, T, z0=z0)
(s11, s12, s21, s22) = mx2elem(S)

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=2,
	ystrip1 = set(axislabel=LBLAX_MAGDB, striplabel="Reflection Coefficients (𝛤)"),
	ystrip2 = set(axislabel=LBLAX_MAGDB, striplabel="Transmission Coefficients (𝛵)"),
	xaxis = set(label=LBLAX_FREQ),
)
push!(plot,
	cons(:wfrm, dB20(s11), color1, label="s11", strip=1),
	cons(:wfrm, dB20(s22), color2, label="s22", strip=1),
	cons(:wfrm, dB20(s12), color1, label="s12", strip=2),
	cons(:wfrm, dB20(s21), color2, label="s21", strip=2),
)

pcoll = push!(cons(:plotcoll, title="Transmission Line Test"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


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

mutable struct ApplicationGatewayFirewallRuleSetPropertiesFormat <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    ruleSetType::Any # spec type: Union{ Nothing, String } # spec name: ruleSetType
    ruleSetVersion::Any # spec type: Union{ Nothing, String } # spec name: ruleSetVersion
    ruleGroups::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayFirewallRuleGroup} } # spec name: ruleGroups

    function ApplicationGatewayFirewallRuleSetPropertiesFormat(;provisioningState=nothing, ruleSetType=nothing, ruleSetVersion=nothing, ruleGroups=nothing)
        o = new()
        validate_property(ApplicationGatewayFirewallRuleSetPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ApplicationGatewayFirewallRuleSetPropertiesFormat, Symbol("ruleSetType"), ruleSetType)
        setfield!(o, Symbol("ruleSetType"), ruleSetType)
        validate_property(ApplicationGatewayFirewallRuleSetPropertiesFormat, Symbol("ruleSetVersion"), ruleSetVersion)
        setfield!(o, Symbol("ruleSetVersion"), ruleSetVersion)
        validate_property(ApplicationGatewayFirewallRuleSetPropertiesFormat, Symbol("ruleGroups"), ruleGroups)
        setfield!(o, Symbol("ruleGroups"), ruleGroups)
        o
    end
end # type ApplicationGatewayFirewallRuleSetPropertiesFormat

const _property_map_ApplicationGatewayFirewallRuleSetPropertiesFormat = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("ruleSetType")=>Symbol("ruleSetType"), Symbol("ruleSetVersion")=>Symbol("ruleSetVersion"), Symbol("ruleGroups")=>Symbol("ruleGroups"))
const _property_types_ApplicationGatewayFirewallRuleSetPropertiesFormat = Dict{Symbol,String}(Symbol("provisioningState")=>"String", Symbol("ruleSetType")=>"String", Symbol("ruleSetVersion")=>"String", Symbol("ruleGroups")=>"Vector{ApplicationGatewayFirewallRuleGroup}")
Base.propertynames(::Type{ ApplicationGatewayFirewallRuleSetPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayFirewallRuleSetPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayFirewallRuleSetPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayFirewallRuleSetPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayFirewallRuleSetPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayFirewallRuleSetPropertiesFormat[property_name]

function check_required(o::ApplicationGatewayFirewallRuleSetPropertiesFormat)
    (getproperty(o, Symbol("ruleSetType")) === nothing) && (return false)
    (getproperty(o, Symbol("ruleSetVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("ruleGroups")) === nothing) && (return false)
    true
end

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


text:
using LazIO

const c_air, c_water = 2.997*10^8, 2.25*10^8
const n_air_water = c_water/c_air

function correct!(ws::LazIO.LazPoint, uw::LazIO.LazPoint)
    x, y, z = correct((ws.X, ws.Y, ws.Z), (uw.X, uw.Y, uw.Z))
    uw.X = round(Int32, x)
    uw.Y = round(Int32, y)
    uw.Z = round(Int32, z)
    uw
end

"Correct an underwater point according to a watersurface point from the refraction factor by using Snell's Law."
function correct(watersurface::Tuple{Real,Real,Real}, underwater::Tuple{Real,Real,Real})
    x0, y0, z0 = watersurface
    x1, y1, z1 = underwater
    Δx, Δy, Δz = (x1-x0), (y1-y0), (z1-z0)

    if Δz >= 0.
        throw(DomainError(Δz, "The underwater point can not be equal or greater in the z dimension to the watersurface point."))
    end

    xcor, ycor, zcor = correct(Δx, Δy, Δz)

    # Define the new corrected coordinates of the points
    xnew = x0 + xcor
    ynew = y0 + ycor
    znew = z0 - zcor

    return (xnew, ynew, znew)
end

function correct(Δx::Real, Δy::Real, Δz::Real)

    # Calculate the distance of the two points (ρ), the vertical-incidence angle (ϕ) and horizontal angle (θ) between the two points
    ρ, θ, ϕ = cart2pol(Δx, Δy, -Δz)

    # Calculate the underwater angle due to the refraction effect.
    # Using Snell's Law: n(air)*sin(ϕ_air) = n(water)*sin(ϕ_water),
    # where ϕ_air: the angle of incidence (I know that the angle of the scan is 20)
    #      ϕ_water: the angle of refraction
    #      n(air): the refraction idex of medium containing the incident ray, value=1.000293
    #      n(water): the refraction idex of medium containing the transmitted ray, value=1.333
    ϕ_water = asind(n_air_water * sind(ϕ))

    # Distance under the water has been influenced by the different speed of light in the water.
    # In a specific moment, ρ = c_air * t and ρ_water = c_water * t.
    # So, by diving these two equations, the ratio (ρ1/ρ2) = n_air_water
    ρ_water = (ρ * n_air_water)

    # Calculate the correction in the 3d space using the pol2cart function
    pol2cart(ρ_water, θ, ϕ_water)
end