text:
# # Load some packages

using MinimaxCalibratedEBayes
using RCall
using Plots
using StatsPlots
using LaTeXStrings
using ExponentialFamilies
using MosekTools
using Random

const MCEB = MinimaxCalibratedEBayes
gr() # plotting backend
GR.settextfontprec(232, 3)

# # Prostate data analysis

# ## Preliminaries

# Let us make some analysis choices from the beginning

# ### Empirical Bayes Targets
# Our goal is to estimate the Posterior mean and the local false sign rates:

target_grid = -3:0.2:3
post_means = PosteriorMean.(StandardNormalSample.(target_grid))
lfsrs =  LFSR.(StandardNormalSample.(target_grid));

# ### Empirical Bayes Prior Class
# For the prior class we consider a mixture of Gaussians
# that have standard deviation $0.2$ and with mixing measure
# supported on $[-3,3]$. Note that we also specify the solver
# to be used for solving the convex programming problems;
# here we use Mosek.
gcal = GaussianMixturePriorClass(0.2, -3:0.02:3, Mosek.Optimizer, ["QUIET" =>true])

# ### Discretization of marginal observations
marginal_grid = -4:0.02:4

# ## Load the data

# The dataset is from the following reference:
#   >Dinesh Singh, Phillip G. Febbo, Kenneth Ross, Donald G. Jackson, Judith Manola, Christine Ladd, Pablo Tamayo, Andrew A. Renshaw, Anthony V. D’Amico, Jerome P. Richie, Eric S. Lander, Massimo Loda, Philip W. Kantoff, Todd R. Golub, and William R. Sellers. Gene expression correlates of clinical prostate cancer behavior. Cancer cell, 1(2): 203–209, 2002.

# See the following monograph for further illustration of empirical Bayes methods on this dataset
#   >Bradley Efron. Large-scale inference: Empirical Bayes methods for estimation, testing, and prediction. Cambridge University Press, 2012.
prostz_file = MinimaxCalibratedEBayes.prostate_data_file
R"load($prostz_file)"
@rget prostz;
Zs = StandardNormalSample.(prostz);

# We will implement the preprocessing steps manually here.
# Afterwards, we will use wrapper that automates all of these steps.

# First we split our data into two folds
Random.seed!(200)
Zs_train, Zs_test, _, _ = MCEB._split_data(Zs);

# Let us quickly check the range of values in the folds
extrema(response(Zs_train)), extrema(response(Zs_test))

# Now we use the first fold to estimate an $L_\infty$ neighborhood
# of the marginal density. We also fit a logspline model of the
# prior density to get our pilot estimators.

# ### Neighborhood construction:

infinity_band_options = KDEInfinityBandOptions(a_min=-4.0, a_max=4.0);
fkde_train = fit(infinity_band_options, response.(Zs_train));

# Let us visualize the result by showing a histogram of the datapoints
# the estimated density function and the $L_\infty$ bands.
prostate_marginal_plot = histogram(response.(Zs_train), bins=20, normalize=true,
          label="", alpha=0.4, fillcolor=:lightgray)
plot!(prostate_marginal_plot, fkde_train, fillalpha=0.45, #L"\bar{f}(x)\pm c_m",
      linewidth=1, label=nothing, #L"\bar{f}(x)\pm c_m",
	  legend=:topleft,  bg_legend=:transparent,
	  fg_legend=:transparent)

# ### G-model fitting

# We choose the following class
exp_family = ContinuousExponentialFamilyModel(Uniform(-3.6,3.6),
                         collect(-3.6:0.1:3.6); df=5, scale=true);
exp_family_solver = MCEB.ExponentialFamilyDeconvolutionMLE(cefm = exp_family, c0=0.001);

# Let us fit to the training data
exp_family_fit = fit(exp_family_solver, Zs_train);

# Let us look at the pilot estimates we got
lfsr_logspline_estimate = estimate.(lfsrs , exp_family_fit)
post_mean_logspline_estimate = estimate.(post_means , exp_family_fit);
#src expfamily_fits_plot = plot(plot(target_grid, lfsr_logspline_estimate),
#src                           plot(target_grid, post_mean_logspline_estimate),
#src						   label="")

# ### Implement discretization

# Discretize our test samples to implement minimax calibration methodology
Zs_test_discr = DiscretizedStandardNormalSamples(Zs_test, marginal_grid);
# Set the $L_{\infty}$ band we computed earlier.
Zs_test_discr = set_neighborhood(Zs_test_discr, fkde_train);

# ### Put everything together
# We put all of the above in one struct that will be used in the final step
# of confidence interval construction

# We set `cache_target=true` to speed up computations for the tutorial, cf.
# Appendix C.3.5 of the manuscript for an explanation.

cache_target = true
mceb_setup = MinimaxCalibratorSetup(Zs_train = Zs_train,
                                    Zs_test = Zs_test,
						            prior_class = gcal,
						            fkde_train = fkde_train,
						            Zs_test_discr = Zs_test_discr,
						            pilot_method = exp_family_fit,
									delta_tuner =  MCEB.HalfCIWidth,
									cache_target = cache_target);

# Fit he calibrators for the posterior means.
post_means_fits = fit.(mceb_setup, post_means)
prostate_postmeans_plot = plot(post_means_fits, ylims=(-3.2,3.2),
							fg_legend=:transparent, bg_legend=:transparent,
							legend = (0.35,0.84))
# Now repeat the same for the local false sign rates.
lfsrs_fits = fit.(mceb_setup, lfsrs)
prostate_lfsrs_plot = plot(lfsrs_fits, label=nothing)

# # The impact of Neighborhoods: Moving to opportunity

# We next seek to repeat the analysis in the second dataset. Here we will also
# see how most of the steps manually conducted for the prostate dataset
# can be automated.

# The references for this dataset is the following:
#  >Raj Chetty and Nathaniel Hendren. The impacts of neighborhoods on intergenerational mobility II: County-level estimates. The Quarterly Journal of Economics, 133(3):1163– 1228, 2018.

# ## Loading the dataset
# Let us start by loading the dataset and filtering out missing observations.

using CSV
nbhood_file = MinimaxCalibratedEBayes.chetty_hendren_file
nbhoods_df = CSV.read(nbhood_file, types=Dict(3=>String));
missing_dx = [!ismissing(x) for x in nbhoods_df[!,:p25_coef]]
β_hat = [nbhoods_df[i,:p25_coef] for i in findall(missing_dx)]
se_hat = [nbhoods_df[i,:p25_se] for i in findall(missing_dx)]
nbhood_Zs = StandardNormalSample.(β_hat./se_hat);
prostate_marginal_grid = -4:0.02:4
length(nbhood_Zs)

# ## Setup the minimax calibrator
# The above steps repeat all the setup we conducted for the prostate data
# to compute the `mceb_setup` object.
Random.seed!(34)
mceb_options = MinimaxCalibratorOptions(prior_class = gcal,
                               marginal_grid = prostate_marginal_grid,
							   pilot_options = exp_family_solver,
							   cache_target=cache_target)
nbhood_mceb_setup = fit(mceb_options, nbhood_Zs);

# Again us a sanity check, let us check the range of data in each fold.
extrema(response.(nbhood_mceb_setup.Zs_train)),extrema(response.(nbhood_mceb_setup.Zs_test))

# Also let us plot again the marginal distribution and the estimated $L_{\infty}$ band.
nbhood_marginal_plot = histogram(response.(nbhood_mceb_setup.Zs_train), bins=20, normalize=true,
	            label=nothing, alpha=0.4, fillcolor=:lightgray)
plot!(nbhood_marginal_plot, nbhood_mceb_setup.fkde_train, fillalpha=0.45,
                  label=nothing, linewidth=1)

# Compute the posterior means
nbhood_post_means_fits = fit.(nbhood_mceb_setup, post_means)
nbhood_postmeans_plot = plot(nbhood_post_means_fits,  ylims=(-3.5,3.5), label=nothing)

# Compute the local false sign rates
nbhood_lfsrs_fits = fit.(nbhood_mceb_setup, lfsrs)
nbhood_lfsrs_plot = plot(nbhood_lfsrs_fits,  label=nothing)

# # Construct final plot
using Plots.Measures
dataset_plot = plot(prostate_marginal_plot, prostate_postmeans_plot, prostate_lfsrs_plot,
	 nbhood_marginal_plot, nbhood_postmeans_plot, nbhood_lfsrs_plot,
	 title = ["a) Prostate" "b)" "c)" " d) Neighborhoods" "e)" "f)"],
	 layout= (2,3), size=(750,500))
plot(dataset_plot)

#src savefig(dataset_plot, "dataset_plots.svg")	#src
#src savefig(dataset_plot, "dataset_plots.pdf")	#src


text:
abstract type SUDSStruct end

# STATIONCOMP:  Generic station component information
mutable struct StationComp <: SUDSStruct
  az            ::Int16    # component azimuth, N°E
  inc           ::Int16    # component angle of incidence from vertical
  lat           ::Float64  # latitude, N = +
  lon           ::Float64  # longitude, E = +
  ele           ::Float32  # elevation in meters
  codes         ::Array{UInt8,1}
  # 1       d=dam, n=nuclear plant, v=vault, b=buried, s=surface, etc.
  # 2       annotated comment code
  # 3       type device data recorded on
  # 4       rock type: i=igneous, m=metamorphic, s=sedimentary
  # 5:6     code for type of rock
  # 7       p=permafrost, etc.
  # 8       sensor type: d=displacement, v=velocity, a=acceleration, t=time code
  # 9       data type: see suds_decode.jl
  # 10      data units: d=digital counts, v=mV, n=nanometers (nm/s or nm/s2
  # 11      polarity: n=normal, r=reversed
  # 12      status: d=dead, g=good
  gain          ::Array{Float32,1}  # 1   maximum gain of the amplifier
                                    # 2   abs val at which clipping begins
                                    # 3   conversion to millivolts: mV/count
  a2d           ::Array{Int16,1}    # 1   a2d channel number
                                    # 2   gain of analog-to-digital converter
  t_corr        ::Array{Float32,1}  # 1   clock correction in seconds
                                    # 2   seismological station delay

  StationComp() = new(zero(Int16),
                      zero(Int16),
                      zero(Float64),
                      zero(Float64),
                      zero(Float32),
                      zeros(UInt8, 12),
                      zeros(Float32, 3),
                      zeros(Int16, 2),
                      zeros(Float32, 2)
                      )
end

#=  DESCRIPTRACE:  Descriptive information about a seismic trace.
                   Normally followed by waveform data =#
mutable struct TraceHdr <: SUDSStruct
  net           ::UInt16
  n_ch          ::Int16
  desc          ::UInt8
  ns            ::Int32

  TraceHdr() = new( 0x0000,
                    zero(Int16),
                    0x00,
                    zero(Int32)
                    )
end

# FEATURE:  Observed phase arrival time, amplitude, and period
mutable struct SudsPhase <: SUDSStruct
  pc            ::UInt16
  onset         ::UInt8
  fm            ::UInt8
  snr           ::Int16
  gr            ::Int16
  amp           ::Float32

  SudsPhase() = new(0x0000,
                    0x00,
                    0x00,
                    zero(Int16),
                    zero(Int16),
                    0.0f0)
end

#  ORIGIN: Information about a specific solution for a given event
mutable struct SudsEvtHdr <: SUDSStruct
  evno        ::Int32
  auth        ::Int16
  chars       ::Array{UInt8,1}
  reg         ::Int32
  ot          ::Float64
  lat         ::Float64
  lon         ::Float64
  floats      ::Array{Float32,1}
  model       ::Array{UInt8, 1}
  gap         ::Int16
  d_min       ::Float32
  shorts      ::Array{Int16,1}
  mag         ::Array{Float32,1}

  SudsEvtHdr() = new( zero(Int32),
                      zero(Int16),
                      zeros(UInt8, 6),
                      zero(Int32),
                      zero(Float64),
                      zero(Float64),
                      zero(Float64),
                      zeros(Float32, 4),
                      zeros(UInt8, 6),
                      zero(Int16),
                      zero(Float32),
                      zeros(Int16, 8),
                      zeros(Float32, 3)
                      )

end

# CHANSET
mutable struct ChanSet <: SUDSStruct
  typ     ::Int16             # CHANSET
  n       ::Int16
  sta     ::Array{UInt8,1}
  tu      ::Int32
  td      ::Int32
  inst    ::Int32             # CHANSETENTRY
  stream  ::Int16
  chno    ::Int16

  ChanSet() = new(zero(Int16),
                  zero(Int16),
                  zeros(UInt8, 9),
                  zero(Int32),
                  zero(Int32),
                  zero(Int32),
                  zero(Int16),
                  zero(Int16))
end

mutable struct SUDSBuf
  # struct_tag
  sid::Int16
  nbs::Int32
  nbx::Int32
  nx::Int32
  nz::Int32

  # staident
  hdr       ::Array{UInt8,1}
  id        ::Array{UInt8,1}
  id_str    ::String

  # common variables in many structures
  data_type ::UInt8           # data type code
  sync_code ::UInt8           # sync code
  irig      ::Bool            # whether the time correction is to "IRIG"
  fs        ::Float32         # fs
  rc        ::Float32         # rc
  t_f64     ::Float64         # Begin time, phase time
  t_f32     ::Float32         # Pick time (τ)
  t_i32     ::Int32           # Effective time, time picked (!= pick time)
  t_i16     ::Int16           # Local time correction in minutes
  tc        ::Float64         # Time correction

  # Placeholders
  nx_new    ::Int64
  nx_add    ::Int64

  # too trivial for its own sub-buffer
  comm_i    ::Array{Int16,1}  # Numeric indices associated with comments
  comm_s    ::String

  # sub-buffers
  S05::StationComp
  T::TraceHdr
  P::SudsPhase
  H::SudsEvtHdr
  C::ChanSet

  function SUDSBuf()
    new(
          zero(Int16),        # sid
          zero(Int32),        # nbs
          zero(Int32),        # nbx
          zero(Int32),        # nx
          zero(Int32),        # nz

          zeros(UInt8, 12),   # hdr
          zeros(UInt8, 13),   # id
          "",                 # id_str

          0x00,               # data_type
          0x00,               # sync_code
          false,              # irig
          0.0f0,              # fs
          0.0f0,              # rc
          zero(Float64),      # t_f64
          zero(Float32),      # t_f32
          zero(Int32),        # t_i32
          zero(Int16),        # t_i16
          zero(Float64),      # tc

          KW.nx_new,          # nx_new
          KW.nx_add,          # nx_add

          zeros(Int16, 2),    # comm_i
          "",                 # comm_s

        # structural buffers                Code
        StationComp(),                      # 5
        TraceHdr(),                         # 6-7
        SudsPhase(),                        # 10
        SudsEvtHdr(),                       # 14
        ChanSet(),                          # 32--33

        )
  end
end

const SB = SUDSBuf()

ext: jl


text:
#---------------------------------------------------------------------#
#This function computes the global Mass and Differentiation Matrices.
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function create_rhs(q,qb,M,Dwe,intma,periodicity,gravity,Δ_NL,DFloat)

    #Get lengths of arrays
    (Np,nelem)=size(intma)
    Npoin=size(q,1)

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

    #build differentiation matrix contribution
    fe=zeros(DFloat,Np,2)
    for e=1:nelem
        for i=1:Np
            I=intma[i,e]
            hs=q[I,1]
            hb=qb[I,1]
            h=hs+hb
            u=q[I,2]/h
            fe[i,1]=h*u
            fe[i,2]=(h*u*u + 0.5*(gravity*hs^2))*Δ_NL + gravity*hs*hb
        end
        for i=1:Np
            I=intma[i,e]
            for j=1:Np
                rhs[periodicity[I],:]+=Dwe[i,j]*fe[j,:]
            end
        end
    end

    #build flux matrix contribution (numerical flux)
    flux=zeros(DFloat,2)
    feL=zeros(DFloat,2)
    feR=zeros(DFloat,2)
    for e=1:nelem
        #Left State
        eL=e
        L=intma[Np,eL]
        hsL=q[L,1]
        hbL=qb[L]
        hL=hsL + hbL
        UL=q[L,2]
        uL=UL/hL
        feL[1]=hL*uL
        feL[2]=( hL*uL*uL + 0.5*gravity*hsL^2 )*Δ_NL + gravity*hsL*hbL

        #Right State
        eR=e+1
        if (e == nelem)
            eR=1
        end
        R=intma[1,eR]
        hsR=q[R,1]
        hbR=qb[R]
        hR=hsR + hbR
        UR=q[R,2]
        uR=UR/hR
        feR[1]=hR*uR
        feR[2]=( hR*uR*uR + 0.5*gravity*hsR^2 )*Δ_NL + gravity*hsR*hbR

        #Jump Speed
        λL=( abs(uL) + sqrt(gravity*hL) )*Δ_NL + sqrt(gravity*hbL)*(1.0-Δ_NL);
        λR=( abs(uR) + sqrt(gravity*hR) )*Δ_NL + sqrt(gravity*hbR)*(1.0-Δ_NL);
        λ=max( λL, λR )
        flux[1]=0.5*( feL[1] + feR[1] - λ*(hR - hL) )
        flux[2]=0.5*( feL[2] + feR[2] - λ*(UR - UL) )
        rhs[periodicity[L],:]-=flux[:]
        rhs[periodicity[R],:]+=flux[:]
    end

    #Build RHS vector
    rhs[:,1]=M\rhs[:,1]
    rhs[:,2]=M\rhs[:,2]

    #Return arrays
    return(rhs)
end


text:
function EventState(arg0::ODEEventHandler, arg1::jdouble, arg2::jdouble, arg3::jint, arg4::BracketedUnivariateSolver)
    return EventState((ODEEventHandler, jdouble, jdouble, jint, BracketedUnivariateSolver), arg0, arg1, arg2, arg3, arg4)
end

function do_event(obj::EventState, arg0::ODEStateAndDerivative)
    return jcall(obj, "doEvent", EventState_EventOccurrence, (ODEStateAndDerivative,), arg0)
end

function evaluate_step(obj::EventState, arg0::ODEStateInterpolator)
    return jcall(obj, "evaluateStep", jboolean, (ODEStateInterpolator,), arg0)
end

function get_convergence(obj::EventState)
    return jcall(obj, "getConvergence", jdouble, ())
end

function get_event_handler(obj::EventState)
    return jcall(obj, "getEventHandler", ODEEventHandler, ())
end

function get_event_time(obj::EventState)
    return jcall(obj, "getEventTime", jdouble, ())
end

function get_max_check_interval(obj::EventState)
    return jcall(obj, "getMaxCheckInterval", jdouble, ())
end

function get_max_iteration_count(obj::EventState)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function reinitialize_begin(obj::EventState, arg0::ODEStateInterpolator)
    return jcall(obj, "reinitializeBegin", void, (ODEStateInterpolator,), arg0)
end

function try_advance(obj::EventState, arg0::ODEStateAndDerivative, arg1::ODEStateInterpolator)
    return jcall(obj, "tryAdvance", jboolean, (ODEStateAndDerivative, ODEStateInterpolator), arg0, arg1)
end


text:
mutable struct Document
    texts::Vector{PDText}
    str::String
    char2text::Vector{Int}
end

function Document()
    str = join(map(t -> t.str, texts))
    char2text = Int[]
    for i = 1:length(texts)
        t = texts[i]
        for _ = 1:sizeof(t.str)
            push!(char2text, i)
        end
    end
    Document(texts, str, char2text)
end

function Base.findfirst(doc::Document, query::String, start::Int=1)
    r = findnext(query, doc.str, start)
    if r != nothing
        i = doc.char2text[first(r)]
        j = doc.char2text[last(r)]
    end
    r
end


text:
function gencodeface(filename::String, f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time)

foldername = "app";
opufile = "opu" * filename;
cpufile = "cpu" * filename;
gpufile = "gpu" * filename;

nbc = size(f,2);

stropu = "";
strgpu = "";
for k = 1:nbc
    str1, str2 = gencodebou(filename * string(k), f[:,k], xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time);
    stropu = stropu * str1;
    strgpu = strgpu * str2;
end

tmp = "template <typename T> void " * opufile;
tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uhg, T *nlg, T *tau, T *uinf, T *param, T time, int modelnumber, int ib, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw)\n";
tmp = tmp * "{\n";
for k = 1:nbc
    if k == 1
        tmp = tmp * "\tif (ib == " * string(k) * ")\n";
    else
        tmp = tmp * "\telse if (ib == " * string(k) * ")\n";
    end
    tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw);\n";
end
tmp = tmp * "}\n\n";

tmp = tmp * "template void " * opufile;
tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int);\n";
tmp = tmp * "template void " * opufile;
tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int);\n";
stropu = stropu * tmp;

ioopu = open(foldername * "/" * opufile * ".cpp", "w");
write(ioopu, stropu);
close(ioopu);

tmp = "template <typename T> void " * gpufile;
tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uhg, T *nlg, T *tau, T *uinf, T *param, T time, int modelnumber, int ib, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw)\n";
tmp = tmp * "{\n";
tmp = tmp * "\tint blockDim = 256;\n";
tmp = tmp * "\tint gridDim = (ng + blockDim - 1) / blockDim;\n";
tmp = tmp * "\tgridDim = (gridDim>1024)? 1024 : gridDim;\n";
for k = 1:nbc
    if k == 1
        tmp = tmp * "\tif (ib == " * string(k) * ")\n";
    elseif k<nbc
        tmp = tmp * "\telse if (ib == " * string(k) * ")\n";
    else
        tmp = tmp * "\telse if (ib == " * string(k) * ")\n";
    end
    tmp = tmp * "\t\tkernel"  * gpufile * string(k) * "<<<gridDim, blockDim>>>(f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw);\n";
end
tmp = tmp * "}\n\n";

tmp = tmp * "template void " * gpufile;
tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int);\n";
tmp = tmp * "template void " * gpufile;
tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int);\n";
strgpu = strgpu * tmp;

iogpu = open(foldername * "/" * gpufile * ".cu", "w");
write(iogpu, strgpu);
close(iogpu);

strcpu = replace(stropu, "opu" => "cpu");
strcpu = replace(strcpu, "for (int i = 0; i <ng; i++) {" => "#pragma omp parallel for\n\tfor (int i = 0; i <ng; i++) {");

iocpu = open(foldername * "/" * cpufile * ".cpp", "w");
write(iocpu, strcpu);
close(iocpu);

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
#!/usr/bin/env julia
# See LICENSE.
# Copyright (C) 2019 Akito <akito.kitsune@protonmail.com>

using Dates

let clientConfigDir::String,                ## client-config-dir set in OpenVPN Server config.
    originalIpAddressArray::Array{Int64, 4} ## Starting point in IP address iteration.
end

clientConfigDir = "/var/etc/openvpn/ccd"
# Debugging
#clientConfigDir = "/home/akito/src/julia-serving-hookers/tmp"

if length(ARGS) >= 1 && ARGS[1] == "init"
  ## Creates configuration file.
  ## Backs up old cfg and/or creates new and empty one.
  ##
  ###  Usage
  ## ./ovpn_staticIP_giver.jl init
  println("Create new configuration file?")
  println("If an old configuration file exists it will be backed up.")
  print("Enter YES if you agree: ")
  agreement = readline(stdin)
  if agreement == "YES"
    cT = Dates.format(now(), "yyyymmddHHMMSS")
    filename = "$clientConfigDir/clients.cfg"
    bkfilename = "$clientConfigDir/clients.cfg.$cT"
    try
      oldfile = open(filename, "r")
      newfile = open(bkfilename, "w")
      write(newfile, oldfile)
      close(newfile)
    catch
      file = open(filename, "w")
      write(file, "")
      close(file)
    end
  else
    println("Request denied. Exiting.")
    exit(1)
  end
elseif length(ARGS) >= 1 && ARGS[1] == "emergency"
  ## Pray to Satan for this to never be necessary.
  ##
  ## Creates static IPs for each client manually.
  ## Needed if the config folder is lost and
  ## every config needs to be redone for
  ## assuring the clients their respective IP.
  ##
  ###  Usage
  ## First, create file "clients.cfg" and
  ## add each client line by line, like this:
  ## common-name 10.10.0.13
  ##
  ## Just 1 spelling mistake can make at least
  ## 1 client entry invalid, i.e. this
  ## emergency recovery wouldn't work.
  ##
  ## Once done, run:
  ## ./ovpn_staticIP_giver.jl emergency
  global clientConfigDir
  function emergeClients()
    println("Run emergency mode?!")
    println("If client-config files exist, they will be overwritten!")
    print("Enter YES if you agree: ")
    agreement = readline(stdin)
    if agreement == "YES"
      println("Starting recovery...")
    else
      println("Request denied. Exiting.")
      exit(1)
    end
    clients = open("$clientConfigDir/clients.cfg")
    for line in eachline(clients)
      ## Reads a config called "clients.cfg" which location
      ## is either specified in "ovpn_staticIP_giver.cfg"
      ## or defaults to "/var/etc/openvpn/ccd".
      ##
      ## Iterates through the whole list of client/IP
      ## pairs and generates their respective configs
      ## to re-enable their static IPs.
      lineArray = split(line, " ")
      commonName = lineArray[1]
      ipAddress = lineArray[2]
      clientConfig = open("$clientConfigDir/$commonName", "w")
      write(clientConfig, "ifconfig-push $ipAddress $ipAddress")
      close(clientConfig)
    end
    close(clients)
  end
  try
    cfg = open("ovpn_staticIP_giver.cfg", "r")
    global clientConfigDir = readline(cfg)
    close(cfg)
    emergeClients()
  catch
    emergeClients()
  end
elseif length(ARGS) == 2 && ARGS[1] == "add"
  ## Adds client by common-name and automatically
  ## assigns an available static IP address
  ## permanently to that client by adding an
  ## entry to "clients.cfg".
  ##
  ###  Usage
  ## ./ovpn_staticIP_giver.jl add common-name
  global clientConfigDir
  global originalIpAddressArray = [10, 148, 0, 14]
  commonName = ARGS[2]
  clients = open("$clientConfigDir/clients.cfg", "r+")
  clientConfig = open("$clientConfigDir/$commonName", "w+")
  ipAddressArray = originalIpAddressArray
  newIpAddressArray = originalIpAddressArray
  function assignIP()
    position = 4
    ipArc = 0
    function findNextFreeIP(position, tempIpAddressArray)
      ipArc = tempIpAddressArray[position]
      if position <= -1
        println("Assigning invalid address.")
        exit(1)
      elseif position == 1 && ipArc < 10
        println("Ran out of IP addresses.")
        exit(1)
      end
      if ipArc == ipAddressArray[position] && ipArc < 255
        tempIpAddressArray[position] = ipArc += 1
      elseif ipArc >= 255
        tempIpAddressArray[position] = 1
        position -= 1
        findNextFreeIP(position, tempIpAddressArray)
      end
      return tempIpAddressArray
    end
    function doubleCheckAssignment()
      ## Just before assigning the new static
      ## IP address, it double checks against
      ## every entry, if it is already registered.
      for line in eachline(clients)
        lineArray = split(line, " ")
        cn = lineArray[1]
        ipAddress = lineArray[2]
        tempIpAddressArray = parse.(Int64, split(ipAddress, "."))
        if newIpAddressArray == tempIpAddressArray
          println("Assignment invalid.")
          println("Address already assigned. Exiting.")
          exit(1)
        elseif cn == commonName
          println("Assignment invalid.")
          println("This Common Name is already registered. Exiting.")
          exit(1)
        end
      end
    end
    for line in eachline(clients)
      ## Reads a config called "clients.cfg" which location
      ## is either specified in "ovpn_staticIP_giver.cfg"
      ## or defaults to "/var/etc/openvpn/ccd".
      ##
      ## Iterates through the whole list of client/IP
      ## pairs and generates their respective configs
      ## to re-enable their static IPs.
      global newIpAddressArray
      lineArray = split(line, " ")
      cn = lineArray[1] ## Just to prevent commonName i.e. client duplication.
      ipAddress = lineArray[2] ## This line's IP address read from cofig file.
      ipAddressArray = parse.(Int64, split(ipAddress, ".")) ## Converts read IP address to Array of Ints.
      if ipAddressArray == ""
        ## Ignore empty lines.
        continue
      elseif commonName == cn
        ## Double check to prevent entry duplication.
        println("Client registered already.")
        try
          ## If a file exists and commonName is set in clients.cfg,
          ## then client is definitely assigned a static IP.
          f = open("$clientConfigDir/$commonName", "r")
          close(f)
          println("Client config file present.")
        catch
          ## Both the common-name has to be in the clients.cfg
          ## and the $commonName client config file must be present.
          ## If one is missing, something is not right.
          println("Client registered already, but config file is missing.")
          println("This needs manual fixing.")
          exit(1)
        finally
          println("Client already exists. Exiting.")
          exit(1)
        end
      elseif newIpAddressArray == ipAddressArray
        ## findNextFreeIP continues iterating until an IP
        ## is found that has not been registered before.
        global newIpAddressArray = findNextFreeIP(position, newIpAddressArray)
        continue
      end
    end
    ## Actual assignment of new static IP addresses
    doubleCheckAssignment()
    firstByte = newIpAddressArray[1]
    secondByte = newIpAddressArray[2]
    thirdByte = newIpAddressArray[3]
    fourthByte = newIpAddressArray[4]
    ipAddressArray = "$firstByte.$secondByte.$thirdByte.$fourthByte"
    write(clientConfig, "ifconfig-push $ipAddressArray $ipAddressArray\n")
    println("Assigned $ipAddressArray to $commonName !")
    write(clients, "$commonName $ipAddressArray\n")
    println("Added $commonName to $clientConfigDir/clients.cfg !")
    close(clientConfig)
    close(clients)
  end
  assignIP()
elseif length(ARGS) >= 1 && ARGS[1] == "remove"
  #TODO remove entry from "clients.cfg" and free the static IP address
  println("Not implemented yet.")
elseif length(ARGS) >= 1 && ARGS[1] == "help"
  #TODO print usage help
  println("Not implemented yet.")
elseif length(ARGS) == 0
  #TODO print usage help
  println("No arguments provided.")
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractDensity

Subtypes of `AbstractDensity` must implement the function

* `BAT.eval_logval_unchecked`

For likelihood densities this is typically sufficient, since shape, and
variate bounds will be inferred from the prior.

Densities with a known variate shape may also implement

* `ValueShapes.varshape`

Densities with known variate bounds may also implement

* `BAT.var_bounds`

!!! note

    The function `BAT.var_bounds` is not part of the stable public BAT-API,
    it's name and arguments may change without deprecation.
"""
abstract type AbstractDensity end
export AbstractDensity

"""
    BAT.eval_logval_unchecked(density::AbstractDensity, v::Any)

Compute log of the value of a multivariate density function for the given
variate/parameter-values.

Input:

* `density`: density function
* `v`: argument, i.e. variate / parameter-values

Note: If `eval_logval_unchecked` is called with an argument that is out of bounds,
the behavior is undefined. The result for arguments that are not within
bounds is *implicitly* `-Inf`, but it is the caller's responsibility to handle
these cases.
"""
function eval_logval_unchecked end

"""
    var_bounds(
        density::AbstractDensity
    )::Union{AbstractVarBounds,Missing}

*BAT-internal, not part of stable public API.*

Get the parameter bounds of `density`. See `eval_logval_unchecked` and
`eval_logval` for the implications and handling of bounds.
If bounds are missing, `eval_logval_unchecked` must be prepared to
handle any parameter values.
"""
var_bounds(density::AbstractDensity) = missing

"""
    ValueShapes.totalndof(density::AbstractDensity)::Union{Int,Missing}

Get the number of degrees of freedom of the variates of `density`. May return
`missing`, if the shape of the variates is not fixed.
"""
function ValueShapes.totalndof(density::AbstractDensity)
    shape = varshape(density)
    ismissing(shape) ? missing : ValueShapes.totalndof(shape)
end

"""
    ValueShapes.varshape(
        density::AbstractDensity
    )::Union{ValueShapes.AbstractValueShape,Missing}

    ValueShapes.varshape(
        density::DistLikeDensity
    )::ValueShapes.AbstractValueShape

Get the shapes of the variates of `density`.

For prior densities, the result must not be `missing`, but may be `nothing` if
the prior only supports unshaped variate/parameter vectors.
"""
ValueShapes.varshape(density::AbstractDensity) = missing

"""
    eval_logval(density::AbstractDensity, v::Any, T::Type{<:Real})

*BAT-internal, not part of stable public API.*

Evaluates density log-value via `eval_logval_unchecked`.

Throws an exception on any of these conditions:

* The variate shape of `density` (if known) does not match the shape of `v`.
* The return value of `eval_logval_unchecked` is `NaN`.
* The return value of `eval_logval_unchecked` is an equivalent of positive
  infinity.
"""
function eval_logval end

@inline function eval_logval(density::AbstractDensity, v::Any, T::Type{<:Real})
    v_shaped = fixup_variate(varshape(density), v)
    R = density_logval_type(v_shaped, T)
    v_stripped = stripscalar(v_shaped)
    _generic_eval_logval_impl(density, v_stripped, R)
end

function _generic_eval_logval_impl(density::AbstractDensity, v::Any, T::Type)
    logval = try
        eval_logval_unchecked(density, v)
    catch err
        rethrow(_density_eval_error(density, v, err))
    end

    _check_density_logval(density, v, logval)

    return convert(T, logval)::T
end

function _density_eval_error(density::AbstractDensity, v::Any, err::Any)
    ErrorException("Density evaluation failed at v = $(variate_for_msg(v)) due to exception $err, density has type $(typeof(density))")
end

function _check_density_logval(density::AbstractDensity, v::Any, logval::Real)
    if isnan(logval)
        throw(ErrorException("Log-density must not evaluate to NaN, v = $(variate_for_msg(v)) , density has type $(typeof(density))"))
    end

    if !(logval < float(typeof(logval))(+Inf))
        throw(ErrorException("Log-density must not evaluate to posivite infinity, v = $(variate_for_msg(v)), density has type $(typeof(density))"))
    end

    nothing
end

variate_for_msg(v::Real) = v
# Strip dual numbers to make errors more readable:
variate_for_msg(v::ForwardDiff.Dual) = ForwardDiff.value(v)
variate_for_msg(v::AbstractArray) = variate_for_msg.(v)
variate_for_msg(v::NamedTuple) = map(variate_for_msg, v)

"""
    BAT.density_logval_type(v::Any, T::Type{<:Real})

*BAT-internal, not part of stable public API.*

Determine a suitable return type of log-density functions, given a variate
shape `vs` and a default result type `T`.
"""
function density_logval_type end

@inline function density_logval_type(v::Any, T::Type{<:Real})
    U = float(getnumtype(v))
    promote_type(T, U)
end

default_dlt() = Float32

"""
    BAT.log_zero_density(T::Type{<:Real})

log-density value to assume for regions of implicit zero density, e.g.
outside of variate/parameter bounds/support.

Returns an equivalent of negative infinity.
"""
log_zero_density(T::Type{<:Real}) = float(T)(-Inf)

"""
    BAT.is_log_zero(x::Real, T::Type{<:Real} = typeof(x)}

*BAT-internal, not part of stable public API.*

Check if x is an equivalent of log of zero, resp. negative infinity,
in respect to type `T`.
"""
function is_log_zero(x::Real, T::Type{<:Real} = typeof(x))
    U = typeof(x)

    FT = float(T)
    FU = float(U)

    x_notnan = !isnan(x)
    x_isinf = !isfinite(x)
    x_isneg = x < zero(x)
    x_notgt1 = !(x > log_zero_density(FT))
    x_notgt2 = !(x > log_zero_density(FU))
    x_iseq1 = x ≈ log_zero_density(FT)
    x_iseq2 = x ≈ log_zero_density(FU)

    x_notnan && ((x_isinf && x_isneg) | x_notgt1 | x_notgt2 | x_iseq1 | x_iseq1)
end

# ToDO: Consider deprecation in favor of logdensityof
"""
    logvalof(density::AbstractDensity)::Function

Returns a function that computes the logarithmic value of `density` at a
given point:

```julia
f = logvalof(density)
log_density_at_v = f(v)
```
"""
logvalof(density::AbstractDensity) = logdensityof(density)

"""
    logdensityof(density::AbstractDensity, v)::Real
    logdensityof(density::AbstractDensity)::Function

*Experimental feature, not part of stable public API.*

Computes the logarithmic value of `density` at a given point, resp. returns a
function that does so:

```julia
logy = logdensityof(density, v)
logdensityof(density, v) == logdensityof(density)(v)
```

Note: This function should *not* be specialized for custom density types!
"""
function logdensityof end
export logdensityof

logdensityof(density::AbstractDensity, v::Any) = eval_logval(density, v, default_dlt())

logdensityof(density::AbstractDensity) = LogDensityOf(density)

struct LogDensityOf{D<:AbstractDensity} <: Function
    density::D
end

@inline (lvd::LogDensityOf)(v::Any) = logdensityof(lvd.density, v)

(Base.:-)(lvd::LogDensityOf) = NegLogDensityOf(lvd.density)

ValueShapes.varshape(lvd::LogDensityOf) = varshape(lvd.density)
ValueShapes.unshaped(lvd::LogDensityOf) = LogDensityOf(unshaped(lvd.density))

struct NegLogDensityOf{D<:AbstractDensity} <: Function
    density::D
end

@inline (lvd::NegLogDensityOf)(v::Any) = - logdensityof(lvd.density, v)

(Base.:-)(lvd::NegLogDensityOf) = LogDensityOf(lvd.density)

ValueShapes.varshape(lvd::NegLogDensityOf) = varshape(lvd.density)
ValueShapes.unshaped(lvd::NegLogDensityOf) = NegLogDensityOf(unshaped(lvd.density))

"""
    abstract type DistLikeDensity <: AbstractDensity

A density that implements part of the `Distributions.Distribution` interface.
Such densities are suitable for use as a priors.

Typically, custom priors should be implemented as subtypes of
`Distributions.Distribution`. BAT will automatically wrap them in a subtype of
`DistLikeDensity`.

Subtypes of `DistLikeDensity` are required to support more functionality
than an [`AbstractDensity`](@ref), but less than a
`Distribution{Multivariate,Continuous}`.

A `d::Distribution{Multivariate,Continuous}` can be converted into (wrapped
in) an `DistLikeDensity` via `conv(DistLikeDensity, d)`.

The following functions must be implemented for subtypes:

* `BAT.eval_logval_unchecked`

* `ValueShapes.varshape`

* `Distributions.sampler`

* `Statistics.cov`

* `BAT.var_bounds`

!!! note

    The function `BAT.var_bounds` is not part of the stable public BAT-API
    and subject to change without deprecation.
"""
abstract type DistLikeDensity <: AbstractDensity end
export DistLikeDensity

"""
    var_bounds(density::DistLikeDensity)::AbstractVarBounds

*BAT-internal, not part of stable public API.*

Get the parameter bounds of `density`. Must not be `missing`.
"""
function var_bounds end

"""
    ValueShapes.totalndof(density::DistLikeDensity)::Int

Get the number of degrees of freedom of the variates of `density`. Must not be
`missing`, a `DistLikeDensity` must have a fixed variate shape.
"""
ValueShapes.totalndof(density::DistLikeDensity) = totalndof(var_bounds(density))

"""
    BAT.AnyDensityLike = Union{...}

Union of all types that BAT will accept as a probability density, resp. that
`convert(AbstractDensity, d)` supports:

* [`AbstractDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyDensityLike = Union{
    AbstractDensity,
    Distributions.ContinuousDistribution
}
export AnyDensityLike

"""
    BAT.AnySampleable = Union{...}

Union of all types that BAT can sample from:

* [`AbstractDensity`](@ref)
* [`DensitySampleVector`](@ref)
* `Distributions.Distribution`
"""
const AnySampleable = Union{
    AbstractDensity,
    DensitySampleVector,
    Distributions.Distribution
}
export AnySampleable

"""
    BAT.AnyIIDSampleable = Union{...}

Union of all distribution/density-like types that BAT can draw i.i.d.
(independent and identically distributed) samples from:

* [`DistLikeDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyIIDSampleable = Union{
    DistLikeDensity,
    Distributions.Distribution
}
export AnyIIDSampleable


text:
rand_in(min::Integer, max::Integer) = begin
    rand(min:max)
end

rand_in(min::AbstractFloat, max::AbstractFloat) = begin
    rand() * (max - min) + min
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _geo2eci{T<:AbstractFloat}(lat::T, long::T, alt::T, jd::T)
    Re    = planets["earth"].eqradius*1e-3
    lat   = deg2rad(lat)
    long  = deg2rad(long)
    gst   = ct2lst(0.0, jd)
    sid_angle = gst*pi/12.0 # Sidereal angle.
    theta = long + sid_angle # Azimuth
    altRe = alt + Re
    r     = altRe*cos(lat)
    return r*cos(theta), r*sin(theta), altRe*sin(lat)
end

"""
    geo2eci(latitude, longitude, altitude, jd) -> x, y, z

### Purpose ###

Convert geographic spherical coordinates to Earth-centered inertial coordinates.

### Explanation ###

Converts from geographic spherical coordinates (latitude, longitude, altitude)
to ECI (Earth-Centered Inertial) (x, y, z) rectangular coordinates.  Julian days
is also needed.

Geographic coordinates assume the Earth is a perfect sphere, with radius equal
to its equatorial radius.  ECI coordinates are in km from Earth center at epoch
TOD (True of Date).

### Arguments ###

* `latitude`: geographic latitude, in degrees.
* `longitude`: geographic longitude, in degrees.
* `altitude`: geographic altitude, in kilometers.
* `jd`: Julian days.

The three coordinates can be passed as a 3-tuple `(latitude, longitude,
altitude)`.  In addition, `latitude`, `longitude`, `altitude`, and `jd` can be
given as arrays of the same length.

### Output ###

The 3-tuple of ECI (x, y, z) coordinates, in kilometers.  The TOD epoch is the
supplied `jd` time.

If geographical coordinates are given as arrays, a 3-tuple of arrays of the same
length is returned.

### Example ###

Obtain the ECI coordinates of the intersection of the equator and Greenwich's
meridian on 2015-06-30T14:03:12.857

``` julia
geo2eci(0, 0, 0, jdcnv("2015-06-30T14:03:12.857"))
# => (-4024.8671780315185,4947.835465127513,0.0)
```

### Notes ###

`eci2geo` converts Earth-centered inertial coordinates to geographic spherical
coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
geo2eci(lat::Real, long::Real, alt::Real, jd::Real) =
    _geo2eci(promote(float(lat), float(long), float(alt), float(jd))...)

geo2eci(lla::Tuple{Real, Real, Real}, jd::Real) =
    geo2eci(lla..., jd)

function geo2eci{LA<:Real, LO<:Real, AL<:Real, JD<:Real}(lat::AbstractArray{LA},
                                                         long::AbstractArray{LO},
                                                         alt::AbstractArray{AL},
                                                         jd::AbstractArray{JD})
    @assert length(lat) == length(long) == length(alt) == length(jd)
    typela = typeof(float(one(LA)))
    x = similar(lat, typela)
    y = similar(lat, typela)
    z = similar(lat, typela)
    for i in eachindex(lat)
        x[i], y[i], z[i] =
            geo2eci(lat[i], long[i], alt[i], jd[i])
    end
    return x, y, z
end


text:
struct MDP
    γ  # discount factor
    𝒮  # state space
    𝒜  # action space
    T  # transition function
    R  # reward function
    TR # sample transition and reward
end

MDP(γ, 𝒮, 𝒜, T, R) = MDP(γ, 𝒮, 𝒜, T, R, nothing)

function MDP(T::Array{Float64, 3}, R::Array{Float64, 2}, γ::Float64)
    MDP(γ, 1:size(R,1), 1:size(R,2), (s,a,s′)->T[s,a,s′], (s,a)->R[s,a], nothing)
end

function MDP(T::Array{Float64, 3}, R, γ::Float64)
    MDP(γ, 1:size(T,1), 1:size(T,2), (s,a,s′)->T[s,a,s′], R, nothing)
end

struct MDPInitialStateDistribution{MDP}
    mdp::MDP
end
Base.rand(S::MDPInitialStateDistribution) = generate_start_state(S.mdp)

function get_mdp_type(mdp; γ::Float64=discount(mdp))
    return MDP(
        γ,
        ordered_states(mdp),
        ordered_actions(mdp),
        (s,a, s′=nothing) -> begin
            S′ = transition(mdp, s, a)
            if s′ == nothing
                return S′
            end
            return pdf(S′, s′)
        end,
        (s,a) -> reward(mdp, s, a),
        (s, a)->begin
            s′ = rand(transition(mdp,s,a))
            r = reward(mdp, s, a)
            return (s′, r)
        end
    )
end


text:
#jl Use `Literate.notebook("lm-lstm.jl", ".", execute=false)` to convert to notebook.

# # Character Language Modeling using LSTM on Penn treebank

# ### Importing packages
using Base.Iterators, IterTools, Knet, Printf, LinearAlgebra, StatsBase, Random
#-

# ### Charset
# This will hold our set of characters that the model will be able to process.
struct Charset
    c2i::Dict{Any,Int}
    i2c::Vector{Any}
    eow::Int
    mask::Int
end

function Charset(charset::String; eow="", mask="-")
    i2c = [ eow; mask; [ c for c in charset ]  ]
    c2i = Dict( c => i for (i, c) in enumerate(i2c))
    return Charset(c2i, i2c, c2i[eow], c2i[mask])
end
#-

# ### TextReader
# Here we will read our input files and split into characters line-by-line.
struct TextReader
    file::String
    charset::Charset
end

function Base.iterate(r::TextReader, s=nothing)
    s === nothing && (s = open(r.file))
    eof(s) && return close(s)
    return [ get(r.charset.c2i, c, r.charset.eow) for c in readline(s)], s
end

Base.IteratorSize(::Type{TextReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{TextReader}) = Base.HasEltype()
Base.eltype(::Type{TextReader}) = Vector{Int}
#-

# ### DataIterator
# In this iterator we will handle iterating over the data, and prepare our models input and output by converting lines of similar lengths into mini training batches.
struct DataIterator
    src::TextReader
    batchsize::Int
    maxlength::Int
    batchmajor::Bool
    bucketwidth::Int
    buckets::Vector
    batchmaker::Function
end

function DataIterator(src::TextReader; batchmaker = arraybatch, batchsize = 128, maxlength = typemax(Int),
                batchmajor = false, bucketwidth = 2, numbuckets = min(128, maxlength ÷ bucketwidth))
    ## buckets[i] is an array of sentence pairs with similar length
    buckets = [ [] for i in 1:numbuckets ]
    DataIterator(src, batchsize, maxlength, batchmajor, bucketwidth, buckets, batchmaker)
end

Base.IteratorSize(::Type{DataIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{DataIterator}) = Base.HasEltype()
Base.eltype(::Type{DataIterator}) = NTuple{2}

function Base.iterate(d::DataIterator, state=nothing)
    ## When file is finished but buckets are partially full
    if state == 0
        for i in 1:length(d.buckets)
            if length(d.buckets[i]) > 0
                batch = d.batchmaker(d, d.buckets[i])
                d.buckets[i] = []
                return batch, state
            end
        end
        ## terminate iteration
        return nothing
    end

    while true
        src_next = iterate(d.src, state)

        if src_next === nothing
            state = 0
            return iterate(d, state)
        end

        (src_word, src_state) = src_next
        state = src_state
        src_length = length(src_word)

        (src_length > d.maxlength) && continue
        (src_length < 2) && continue

        i = Int(ceil(src_length / d.bucketwidth))
        i > length(d.buckets) && (i = length(d.buckets))

        push!(d.buckets[i], src_word)
        if length(d.buckets[i]) == d.batchsize
            batch = d.batchmaker(d, d.buckets[i])
            d.buckets[i] = []
            return batch, state
        end
    end
end

function arraybatch(d::DataIterator, bucket)
    src_eow = d.src.charset.eow
    src_lengths = map(x -> length(x), bucket)
    max_length = max(src_lengths...)
    x = zeros(Int64, length(bucket), max_length + 2)

    for (i, v) in enumerate(bucket)
        to_be_added = fill(src_eow, max_length - length(v) + 1)
        x[i,:] = [src_eow; v; to_be_added]
    end

    ## default d.batchmajor is false
    d.batchmajor && (x = x')

    ## the output in lang. model is same as input sequence but shifted one step
    return (x[:, 1:end-1], x[:, 2:end])
end
#-

# ### Embedding and Linear layers
mmul(w,x) = (w == 1 ? x : w == 0 ? 0 : reshape(w * reshape(x,size(x,1),:), (:, size(x)[2:end]...)))

struct Embed; w; end
Embed(charsetsize::Int, embedsize::Int) = Embed(param(embedsize, charsetsize))
(l::Embed)(x) = l.w[:, x]

struct Linear; w; b; end
Linear(inputsize::Int, outputsize::Int) = Linear(param(outputsize, inputsize), param0(outputsize))
(l::Linear)(x) = mmul(l.w,x) .+ l.b
#-

# ### NLL Loss masking function
# This function masks the training output array by zeros, this way we don't propagate loss from paddings
function mask(a, pad)
    a = copy(a)
    for i in 1:size(a, 1)
        j = size(a,2)
        while a[i, j] == pad && j > 1
            if a[i, j - 1] == pad
                a[i, j] = 0
            end
            j -= 1
        end
    end
    return a
end
#-

# ### RNN Language model
# Our model consists of four layers. Size of their outputs are as the following, where T is sequence length, B is batchsize, H is the hidden size of RNN, and V is vocabulary size, E is the embedding size:
#
# * **(T, B)** - Input
# * **(E, T, B)** - Embedding
# * **(H, T, B)** - RNN
# * **(V, T, B)** - Projection
struct LModel
    srcembed::Embed
    rnn::RNN
    projection::Linear
    dropout::Real
    srccharset::Charset
end

function LModel(hidden::Int, srcembsz::Int, srccharset::Charset; layers=1, dropout=0)

    srcembed = Embed(length(srccharset.i2c), srcembsz)
    rnn = RNN(srcembsz, hidden; bidirectional=false, numLayers=layers, dropout=dropout)
    projection = Linear(hidden, length(srccharset.i2c))

    LModel(srcembed, rnn, projection, dropout, srccharset)
end

function (s::LModel)(src, tgt; average=true)
    s.rnn.h, s.rnn.c = 0, 0
    srcembed = s.srcembed(src)
    rnn_out = s.rnn(srcembed)
    dims = size(rnn_out)
    output = s.projection(dropout(reshape(rnn_out, dims[1], dims[2] * dims[3]), s.dropout))
    scores = reshape(output, size(output, 1), dims[2], dims[3])
    nll(scores, mask(tgt, s.srccharset.eow); dims=1, average=average)
end

function generate(s::LModel; start="", maxlength=30)
    s.rnn.h, s.rnn.c = 0, 0
    chars = fill(s.srccharset.eow, 1)
    start = [ c for c in start ]
    starting_index = 1
    for i in 1:length(start)
        push!(chars, s.srccharset.c2i[start[i]])
        charembed = s.srcembed(chars[i:i])
        rnn_out = s.rnn(charembed)
        starting_index += 1
    end

    for i in starting_index:maxlength
        charembed = s.srcembed(chars[i:i])
        rnn_out = s.rnn(charembed)
        output = s.projection(dropout(rnn_out, s.dropout))
        push!(chars, s.srccharset.c2i[ sample(s.srccharset.i2c, Weights(Array(softmax(reshape(output, length(s.srccharset.i2c)))))) ] )

        if chars[end] == s.srccharset.eow
            break
        end
    end

    join([ s.srccharset.i2c[i] for i in chars ], "")
end
#-

# ### Evaluation metrics
#
# Here we define a `loss(model, data)` which returns a `(Σloss, Nloss)` pair if `average=false` and
# a `Σloss/Nloss` average if `average=true` for a whole dataset.
#
# `report_lm(loss)` calculates character perplexity and bit-per-character metrics.
function loss(model, data; average=true)
    mean([model(x,y) for (x,y) in data])
end

report_lm(loss) = (loss=loss, ppl=exp.(loss), bpc=loss ./ log(2))
#-

# ### Training procedure
# Train our model using Adam optimizer and saving the best performing model on dev set.
function train!(model, steps, trn, dev, tst...)
    bestmodel, bestloss = deepcopy(model), loss(model, dev)
    progress!(adam(model, trn), steps=steps) do y
        devloss = loss(model, dev)
        tstloss = map(d->loss(model,d), tst)
        if devloss < bestloss
            bestloss, bestmodel = devloss, deepcopy(model)
        end
        println(stderr)
        (trn=report_lm(tstloss), dev=report_lm(devloss))
    end
    return bestmodel
end
#-

## seed random generators for reproducability
Random.seed!(123);

## Define char set, batchsize and max sequence length
char_set = " #\$&'*-./0123456789<>N\\abcdefghijklmnopqrstuvwxyz"
datadir = "../data/ptb"
BATCHSIZE, MAXLENGTH = 16, 256

@info "Reading data"
charset = Charset(char_set)
train_reader = TextReader("$datadir/train.txt", charset)
dev_reader = TextReader("$datadir/valid.txt", charset)
test_reader = TextReader("$datadir/test.txt", charset)

dtrn = DataIterator(train_reader, batchsize=BATCHSIZE, maxlength=MAXLENGTH)
ddev = DataIterator(dev_reader, batchsize=BATCHSIZE, maxlength=MAXLENGTH)
dtst = DataIterator(test_reader, batchsize=BATCHSIZE, maxlength=MAXLENGTH)
#-

@info "Initializing Language Model"
epochs = 10
ctrn = collect(dtrn)
trn = collect(flatten(shuffle!(ctrn) for i in 1:epochs))
trnmini = ctrn[1:20]
dev = collect(ddev)
model = LModel(256, 256, charset; layers=2, dropout=0.2)
#-

@info "Starting training ..."
model = train!(model, length(ctrn), trn, dev, trnmini)
#-

@info "Finished training, Starting evaluation ..."
trnloss = loss(model, dtrn);
println("Training set scores:       ", report_lm(trnloss))
devloss = loss(model, ddev);
println("Development set scores:    ", report_lm(devloss))
tstloss = loss(model, dtst);
println("Test set scores:           ", report_lm(tstloss))
#-

# ### Sampling sentences from the model
@info "Sample sentences from the model"
print("enter prompt or leave empty for no prompt, CTRL+C to exit")

while true
    print("prompt:")
    prompt = lowercase(readline())
    println(generate(model; start=prompt, maxlength=1024))
end
#-


text:
using LinearAlgebra
using Measurements
using Plots
using Formatting
using LsqFit
using Printf

## initialisation

zscore(x::Measurement) = Measurements.value(x) / Measurements.uncertainty(x)

include("constants.jl")
include("types.jl")
include("data.jl")
include("residuals.jl")
include("plotting.jl")

## load data

# import list of methyl names - numbering here will be used to select methyls for further analysis
methyls = loadmethyls("data/methyls.txt")
nmethyls = length(methyls)
for (i, name) in zip(1:nmethyls, methyls)
    println("$i\t$name")
end

# import Hahn echo relaxation rates and specify associated field strengths
R2data_raw = loadR2("data/R2data.txt", [600, 700, 800, 950])

# import S2tc and CSA measurements
S2tc = loadS2tc("data/S2tc.txt")
ΔσH = load1Hcsa("data/etaH.txt", S2tc, 950)
ΔσC = load13Ccsa("data/13Ccsa.txt")

# subtract CSA contributions from measured R2 rates
R2data = subtractCSA(R2data_raw, S2tc, ΔσH, ΔσC)

# import CPMG data - file format is based on parsed ChemEx output
# specify experiment type, field strength, relaxation time (0.04 s) and a minimum error (0.3 s-1)
CPMGdata1 = loadCPMG("data/950-MQ-CPMG.exp", methyls, :MQ, 950, 0.04, min_err=0.3)
CPMGdata2 = loadCPMG("data/800-MQ-CPMG.exp", methyls, :MQ, 800, 0.04, min_err=0.3)
CPMGdata3 = loadCPMG("data/800-1H-CPMG.exp", methyls, :SQ1H, 800, 0.04, min_err=0.3)

## select a group of residues for analysis

# # FIRST GROUP
# active_methyls = [4, 6, 13, 34, 35, 3, 7, 33]
# spinpar0 = [
#     # ΔδH, ΔδC, R20(ZQ, DQ, QQ⁻, QQ⁺), R20(CPMG1, CPMG2, CPMG3)
#     -0.015  0.5  5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     0.015  0.4  5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     0.03  0.4  5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     -0.03  0.4  5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     -0.02  0.3  5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     .01 .1 5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     .01 .1 5.0  5.0  50.0  50.0  10.0  10.0  10.0
#     .01 .1 5.0  5.0  50.0  50.0  10.0  10.0  10.0
# ]
# mp0 = modelpars_2state(900, .03)  # initial estimate for kex, pB

# SECOND GROUP
active_methyls = [8, 20, 21, 19, 22]
spinpar0 = [
    # ΔδH, ΔδC, R20(ZQ, DQ, QQ⁻, QQ⁺), R20(CPMG1, CPMG2, CPMG3)
    -0.0754704  1.15793   29.4364   12.1437   156.053   30.1853   6.66643  10.1832   12.3749
    0.0518553  0.572343   6.43592   6.26245   26.4896  52.4839   7.10874   6.46275  12.9216
   -0.1         0.247   7.6476   10.7173    40.3527  40.0117  10.4753    9.84693  16.6869
   0.0518553  0.572343   6.43592   6.26245   26.4896  52.4839   7.10874   6.46275  12.9216
   -0.0855321  0.446828   7.6476   10.7173    95.3527  73.0117  10.4753    9.84693  16.6869
    ]
mp0 = modelpars_2state(6900, .16)  # initial estimate for kex, pB

# apply selection to imported data
methyls = methyls[active_methyls]
R2data = R2data[active_methyls]
CPMGdata1 = CPMGdata1[active_methyls]
CPMGdata2 = CPMGdata2[active_methyls]
CPMGdata3 = CPMGdata3[active_methyls]
nmethyls = length(active_methyls)

## prepare functions for fitting

# get list of field strengths for HE measurements
R2_B0s = R2data[1].B0

# pack model parameters and spin parameters into a single vector
function pack(mp, spinpar)
    vec([mp.kex; mp.pB; vec(spinpar)])
end

# unpack parameter vector into model and spin parameters
function unpack(p)
    mp = modelpars_2state(p[1:2]...)
    spinpars = reshape(p[3:end], nmethyls, :)
    return mp, spinpars
end

# simulate collection of experiments - this may need to be edited, e.g. if using more/less CPMG experiments
function simulate_spins(mp, spinpars)
    nmethyls, ~ = size(spinpars)
    simR2 = []
    simCPMG1 = []
    simCPMG2 = []
    simCPMG3 = []
    for i=1:nmethyls
        ΔδH = spinpars[i, 1]
        ΔδC = spinpars[i, 2]
        R2_ZQ = spinpars[i, 3]
        R2_DQ = spinpars[i, 4]
        R2_QQ⁻ = spinpars[i, 5]
        R2_QQ⁺ = spinpars[i, 6]
        R2cpmg1 = spinpars[i, 7]
        R2cpmg2 = spinpars[i, 8]
        R2cpmg3 = spinpars[i, 9]
        push!(simR2, sim_R2s(mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, R2_B0s))
        push!(simCPMG1, sim_CPMG(CPMGdata1[i], mp, ΔδH, ΔδC, R2cpmg1))
        push!(simCPMG2, sim_CPMG(CPMGdata2[i], mp, ΔδH, ΔδC, R2cpmg2))
        push!(simCPMG3, sim_CPMG(CPMGdata3[i], mp, ΔδH, ΔδC, R2cpmg3))
    end
    return simR2, simCPMG1, simCPMG2, simCPMG3
end

# calculate residuals - this may need to be edited, e.g. if using more/less CPMG experiments
function calcresiduals(mp, spinpars)
    nmethyls, ~ = size(spinpars)

    resid = Vector{Float64}()

    for i=1:nmethyls
        ΔδH = spinpars[i, 1]
        ΔδC = spinpars[i, 2]
        R2_ZQ = spinpars[i, 3]
        R2_DQ = spinpars[i, 4]
        R2_QQ⁻ = spinpars[i, 5]
        R2_QQ⁺ = spinpars[i, 6]
        R2cpmg1 = spinpars[i, 7]
        R2cpmg2 = spinpars[i, 8]
        R2cpmg3 = spinpars[i, 9]
        append!(resid, residuals_R2(R2data[i], mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, R2data[i].B0))
        append!(resid, residuals_CPMG(CPMGdata1[i], mp, ΔδH, ΔδC, R2cpmg1))
        append!(resid, residuals_CPMG(CPMGdata2[i], mp, ΔδH, ΔδC, R2cpmg2))
        append!(resid, residuals_CPMG(CPMGdata3[i], mp, ΔδH, ΔδC, R2cpmg3))
    end

    return resid
end

# calculate χ² = sum of squares of residuals, weighted by uncertainty
χ2(mp, spinpar) = sum(calcresiduals(mp, spinpar).^2)

## least-squares fitting

# run the fit
fit = LsqFit.lmfit(p -> calcresiduals(unpack(p)...), pack(mp0, spinpar0), Float64[], show_trace=true)

# extract best-fit parameters
mpfit, spinparfit = unpack(fit.param)

# calculate parameter uncertainties
mpfit_err, spinparfit_err = unpack(fit.param .± stderror(fit))

# calculate best-fit χ² value
χ2(mpfit, spinparfit)

# calculate best-fit predictions and plot
simR2, simCPMG1, simCPMG2, simCPMG3 = simulate_spins(mpfit, spinparfit)
for i=1:nmethyls
    plot(plotdata(R2data[i], simR2[i]),
        plotdata([CPMGdata1[i], CPMGdata2[i], CPMGdata3[i]],
            [simCPMG1[i], simCPMG2[i], simCPMG3[i]]),
        layout=(1,2),title=methyls[i],size=(800,300)) |> display
end

println()
println("kex = $(mpfit_err.kex) s⁻¹")
println("pB = $(mpfit_err.pB)")
spinparfit_err

## plot a single residue fit (used for main text figures)
font = Plots.font("Helvetica", 6)
plotstyle = Dict(
    :guidefont=>font,
    :xtickfont=>font,
    :ytickfont=>font,
    :legendfont=>font,
    :framestyle=>:box,
    :grid=>false,
    :lw=>0.5)

i = 1  # select residue number
plot(plotdata(R2data[i], simR2[i]),
    plotdata([CPMGdata1[i], CPMGdata2[i], CPMGdata3[i]], [simCPMG1[i], simCPMG2[i], simCPMG3[i]]),
    layout=(1,2),
    legend=:topright,
    size=(500,200);
    plotstyle...)

## combine plots (used for SI figures)
font = Plots.font("Helvetica", 6)
plotstyle = Dict(
    :guidefont=>font,
    :xtickfont=>font,
    :ytickfont=>font,
    :legendfont=>font,
    :framestyle=>:box,
    :grid=>false,
    :lw=>0.5)

plts = []
# order = [2, 7, 3, 8, 4, 6, 1, 5]  # first group
order = [1, 4, 2, 3, 5]  # second group
for j=1:nmethyls
    i = order[j]
    push!(plts, plotdata(R2data[i], simR2[i]))
    push!(plts, plotdata([CPMGdata1[i], CPMGdata2[i], CPMGdata3[i]],
            [simCPMG1[i], simCPMG2[i], simCPMG3[i]]))
end
for i=1:6 # pad second plots
    push!(plts, plot())
end
plot(plts..., layout=(4,4), legend=nothing, size=(800,600), ms=3; plotstyle...)

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
function encode(s)
    count = 1
    ans = []
    for idx in 2:length(s)
        if s[idx] == s[idx - 1]
            count += 1
        else
            if count > 1
                push!(ans, count)
            end
            push!(ans, s[idx - 1])
            count = 1
        end
        if idx == length(s)
            if count > 1
                push!(ans, count)
            end
            push!(ans, s[idx])
        end
    end
    join(ans)
end

function decode(s)

end

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
using Random
const frequencies = Dict("c"=>16.3516,"c#"=>17.32391, "d"=>18.35405, "d#"=>19.44544, "e"=>20.60172,
"f"=>21.82676, "f#"=>23.12465, "g"=>24.49971, "g#" => 25.95654, "a"=>27.5, "a#"=>29.13524,
"b"=>30.86771
)
#check input string to see if it matches convention
function checkInput(input::AbstractString)
    if (length(input) >3 || length(input) < 2 ) return false end
    return !isnothing(match(r"([a-g]#?[\d])"i, input))
end
#returns the correct frequency for the given note  and octave
function getFreq(input::AbstractString)
    if !checkInput(input) return false end
    if length(input)==2
        noteName = lowercase(string(input[1]))
        octave = parse(Int, input[2])
    else
        noteName = lowercase(input[1:2])
        octave = parse(Int, input[3])
    end
    return frequencies[noteName] * 2^octave
end
function getFreq(input::Real)
    return input
end
#returns a vector where fs is sampling rate, from 0 to seconds in increments of 1/fs
function getTimeVector(seconds,fs)
    return 0.0 : inv(fs) : seconds
end
#returns a sin wave at give frequency for givin time
function getSine(seconds, hz, fs = 41000.0)
    t = 0.0 : inv(fs) : seconds
    sin.(2π  *  hz  *  t )
end
#returns a sin wave at given note name and time
function getSine(seconds, freq::AbstractString, fs = 41000.0)
    hz = getFreq(freq)
    t = 0.0 : inv(fs) : seconds
    sin.(2π  *  hz  *  t )
end
#=
get a sin wave with vibrato code adapted from:
https://www.audiolabs-erlangen.de/resources/MIR/FMP/C1/C1S3_Timbre.html
#Meinard Müller and Frank Zalkow: FMP Notebooks: Educational Material for Teaching and Learning Fundamentals of Music Processing. Proceedings of the International Conference on Music Information Retrieval (ISMIR), Delft, The Netherlands, 2019.

=#
function getSine(seconds, hz, fs, vib_amp=1, vib_rate=5)
    t = 0.0 : inv(fs) : seconds
    freq_vib = hz .+ vib_amp .* sin.(t * 2π * vib_rate)
    phase_vib = zeros(length(t))
    for i in 2:length(t)
        phase_vib[i] = phase_vib[i-1] + 2π* freq_vib[i-1] / fs
    end
    sin.(phase_vib)
end
function getSine(seconds, freq::AbstractString, fs, vib_amp=1, vib_rate=5)
    hz = getFreq(freq)
    t = 0.0 : inv(fs) : seconds
    freq_vib = hz .+ vib_amp .* sin.(t * 2π * vib_rate)
    phase_vib = zeros(length(t))
    for i in 2:length(t)
        phase_vib[i] = phase_vib[i-1] + 2π* freq_vib[i-1] / fs
    end
    sin.(phase_vib)
end

#=returns a saw wave at give frequency for givin time
#note to self:change to the series formation of saw wave
to remove aliasing
=#
function getSaw(seconds, hz, fs = 41000.0)
    timeVector = getTimeVector(seconds, fs)
    2 .* mod.(timeVector * hz, 1).-1
end
#returns a sin wave at given note name and time
function getSaw(seconds, freq::AbstractString, fs = 41000.0)
    hz = getFreq(freq)
    timeVector = getTimeVector(seconds, fs)
    2 .* mod.(timeVector * hz, 1).-1
end
function lfo(signal::AbstractVector{T},hz,fun::Function, i::Tuple = (1,40),fs=41000)where
    {T<:Real}
    inc = Float64(0)
    newSignal = T[]
    f(x) = (i[2]-i[1])*sin(2π * x * hz)/2 + (i[2]+i[1])/2
    for i in signal
        push!(newSignal, fun(i,f(inc)))
        inc+=inv(fs)
    end
    return newSignal
end

#returns a vector containing random noise for seconds
function getRand(seconds,fs= 41000)
    rng = MersenneTwister(1234)
    count = length(0.0 : inv(fs) : seconds)
    return rand!(rng, zeros(count)).*2 .-1
end

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
export Linear

"""
**Linear function**

    Linear(c)

Returns the function
```math
f(x) = \\langle c, x \\rangle.
```
"""
struct Linear{R <: Real, A <: AbstractArray{R}} <: ProximableFunction
    c::A
end

is_separable(f::Linear) = true
is_convex(f::Linear) = true
is_smooth(f::Linear) = true

function (f::Linear{R})(x::AbstractArray{R}) where R
    return dot(f.c, x)
end

function prox!(y::AbstractArray{R}, f::Linear{R}, x::AbstractArray{R}, gamma::Union{R, AbstractArray{R}}=1.0) where R
    y .= x .- gamma.*(f.c)
    fy = dot(f.c, y)
    return fy
end

function gradient!(y::AbstractArray{R}, f::Linear{R}, x::AbstractArray{R}) where R
    y .= f.c
    return dot(f.c, x)
end

fun_name(f::Linear) = "Linear function"
fun_expr(f::Linear) = "x ↦ c'x"

function prox_naive(f::Linear, x, gamma)
    y = x - gamma.*(f.c)
    fy = dot(f.c, y)
    return y, fy
end


text:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Test

rate = (u,p,t) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (u,p,t) -> 0.01u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
integrator = init(jump_prob,FunctionMap())
sol = solve(jump_prob,FunctionMap())

jump_prob

using Plots; plotly(); plot(sol)

nums = Int[]
@time for i in 1:1000
  sol = solve(jump_prob,FunctionMap())
  push!(nums,sol[end][1])
end
mean(nums)

using ProfileView
@profile for i in 1:1000; solve(jump_prob,FunctionMap()); end
Profile.clear()
@profile for i in 1:1000; solve(jump_prob,FunctionMap()); end
ProfileView.view()


text:
tests = [
    "streams"
]

for t in tests
    tfile = string(t, ".jl")
    println("  * $tfile ...")
    include(tfile)
end


text:
module AdvancedMH

# Import the relevant libraries.
using AbstractMCMC
using Random
using Requires
using Distributions

# Import specific functions and types to use or overload.
import AbstractMCMC: step!, AbstractSampler, AbstractTransition,
    transition_type, bundle_samples

# Exports
export MetropolisHastings, DensityModel, sample, psample, RWMH, StaticMH, Proposal, Static, RandomWalk

# Abstract type for MH-style samplers.
abstract type Metropolis <: AbstractSampler end
abstract type ProposalStyle end

struct RandomWalk <: ProposalStyle end
struct Static <: ProposalStyle end

# Define a model type. Stores the log density function and the data to
# evaluate the log density on.
"""
    DensityModel{F<:Function} <: AbstractModel

`DensityModel` wraps around a self-contained log-liklihood function `logdensity`.

Example:

```julia
l
DensityModel
```
"""
struct DensityModel{F<:Function} <: AbstractModel
    logdensity :: F
end

# Create a very basic Transition type, only stores the
# parameter draws and the log probability of the draw.
struct Transition{T<:Union{Vector, Real, NamedTuple}, L<:Real} <: AbstractTransition
    params :: T
    lp :: L
end

# Store the new draw and its log density.
Transition(model::M, params::T) where {M<:DensityModel, T} = Transition(params, logdensity(model, params))

# Tell the interface what transition type we would like to use.
transition_type(model::DensityModel, spl::Metropolis) = typeof(Transition(spl.init_params, logdensity(model, spl.init_params)))

# Calculate the density of the model given some parameterization.
logdensity(model::DensityModel, params) = model.logdensity(params)
logdensity(model::DensityModel, t::Transition) = t.lp

# A basic chains constructor that works with the Transition struct we defined.
function bundle_samples(
    rng::AbstractRNG,
    model::DensityModel,
    s::Metropolis,
    N::Integer,
    ts::Type{Any};
    param_names=missing,
    kwargs...
) where {ModelType<:AbstractModel, T<:AbstractTransition}
    return ts
end

function bundle_samples(
    rng::AbstractRNG,
    model::DensityModel,
    s::Metropolis,
    N::Integer,
    ts::Vector{T},
    chain_type::Type{Vector{NamedTuple}};
    param_names=missing,
    kwargs...
) where {ModelType<:AbstractModel, T<:AbstractTransition}
    # Check if we received any parameter names.
    if ismissing(param_names)
        param_names = ["param_$i" for i in 1:length(keys(ts[1].params))]
    else
        # Deepcopy to be thread safe.
        param_names = deepcopy(param_names)
    end

    push!(param_names, "lp")

    # Turn all the transitions into a vector-of-NamedTuple.
    ks = tuple(Symbol.(param_names)...)
    nts = [NamedTuple{ks}(tuple(t.params..., t.lp)) for t in ts]

    return nts
end

function __init__()
    @require MCMCChains="c7f686f2-ff18-58e9-bc7b-31028e88f75d" include("mcmcchains-connect.jl")
    @require StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("structarray-connect.jl")
end

# Include inference methods.
include("proposal.jl")
include("mh-core.jl")

end # module AdvancedMH


text:
type SSLConfig
    data::Ptr{Void}
    rng
    chain::CRT
    dbg
    cert
    key
    alpn_protos

    function SSLConfig()
        conf = new()
        conf.data = Libc.malloc(1000)  # 360
        ccall((:mbedtls_ssl_config_init, MBED_TLS), Void, (Ptr{Void},), conf.data)
        finalizer(conf, conf->begin
            ccall((:mbedtls_ssl_config_free, MBED_TLS), Void, (Ptr{Void},), conf.data)
            Libc.free(conf.data)
        end
        )
        conf
    end
end

type SSLContext <: IO
    data::Ptr{Void}
    config::SSLConfig
    isopen::Bool
    bio

    function SSLContext()
        ctx = new()
        ctx.data = Libc.malloc(1000)  # 488
        ccall((:mbedtls_ssl_init, MBED_TLS), Void, (Ptr{Void},), ctx.data)
        finalizer(ctx, ctx->begin
            ccall((:mbedtls_ssl_free, MBED_TLS), Void, (Ptr{Void},), ctx.data)
            Libc.free(ctx.data)
        end
        )
        ctx
    end
end

function config_defaults!(config::SSLConfig; endpoint=MBEDTLS_SSL_IS_CLIENT,
    transport=MBEDTLS_SSL_TRANSPORT_STREAM, preset=MBEDTLS_SSL_PRESET_DEFAULT)
    @err_check ccall((:mbedtls_ssl_config_defaults, MBED_TLS), Cint,
        (Ptr{Void}, Cint, Cint, Cint),
        config.data, endpoint, transport, preset)
end

function authmode!(config::SSLConfig, auth)
    ccall((:mbedtls_ssl_conf_authmode, MBED_TLS), Void,
        (Ptr{Void}, Cint),
        config.data, auth)
end

function rng!(config::SSLConfig, f_rng::Ptr{Void}, ctx)
    ccall((:mbedtls_ssl_conf_rng, MBED_TLS), Void,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}),
        config.data, f_rng, ctx)
end

function rng!(config::SSLConfig, rng::AbstractRNG)
    config.rng = rng
    rng!(config, c_rng, pointer_from_objref(rng))
end

function ca_chain!(config::SSLConfig, chain=crt_parse_file(TRUSTED_CERT_FILE))
    config.chain = chain
    ccall((:mbedtls_ssl_conf_ca_chain, MBED_TLS), Void,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}),
        config.data, chain.data, C_NULL)
end

function own_cert!(config::SSLConfig, cert::CRT, key::PKContext)
    config.cert = cert
    config.key = key
    @err_check ccall((:mbedtls_ssl_conf_own_cert, MBED_TLS), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}),
        config.data, cert.data, key.data)
end

function setup!(ctx::SSLContext, conf::SSLConfig)
    ctx.config = conf
    @err_check ccall((:mbedtls_ssl_setup, MBED_TLS), Cint,
        (Ptr{Void}, Ptr{Void}),
        ctx.data, conf.data)
end

function set_bio!(ssl_ctx::SSLContext, ctx, f_send::Ptr{Void}, f_recv::Ptr{Void})
    ccall((:mbedtls_ssl_set_bio, MBED_TLS), Void,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
        ssl_ctx.data, ctx, f_send, f_recv, C_NULL)
end

function f_send(c_ctx, c_msg, sz)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    jl_msg = unsafe_wrap(Array, c_msg, sz, false)
    return Cint(write(jl_ctx, jl_msg))
end

function f_recv(c_ctx, c_msg, sz)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    jl_msg = unsafe_wrap(Array, c_msg, sz, false)
    n = readbytes!(jl_ctx, jl_msg, sz)
    return Cint(n)
end

function set_bio!{T<:IO}(ssl_ctx::SSLContext, jl_ctx::T)
    ssl_ctx.bio = jl_ctx
    set_bio!(ssl_ctx, pointer_from_objref(jl_ctx), c_send, c_recv)
    nothing
end

function dbg!(conf::SSLConfig, f::Ptr{Void}, p)
    ccall((:mbedtls_ssl_conf_dbg, MBED_TLS), Void,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}),
        conf.data, f, p)
end

function f_dbg(c_ctx, level, filename, number, msg)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    jl_ctx(level, String(filename), number, String(msg))
    nothing
end

function dbg!(conf::SSLConfig, f)
    conf.dbg = f
    dbg!(conf, c_dbg, pointer_from_objref(f))
    nothing
end

function handshake(ctx::SSLContext)
    @err_check ccall((:mbedtls_ssl_handshake, MBED_TLS), Cint,
        (Ptr{Void},), ctx.data)
    ctx.isopen = true
    nothing
end

function set_alpn!(conf::SSLConfig, protos)
    conf.alpn_protos = protos
    @err_check ccall((:mbedtls_ssl_conf_alpn_protocols, MBED_TLS), Cint,
                     (Ptr{Void}, Ptr{Ptr{Cchar}}), conf.data, protos)
    nothing
end

function alpn_proto(ctx::SSLContext)
    rv = ccall((:mbedtls_ssl_get_alpn_protocol, MBED_TLS), Ptr{Cchar},
               (Ptr{Void},), ctx.data)
    unsafe_string(rv)
end

if Base.VERSION < v"0.5.0-dev+2301"
    import Base: read, write
    const unsafe_read = read
    const unsafe_write = write
    @noinline Base.write(ctx::SSLContext, msg::Base.RefValue{UInt8}) = write(ctx, Base.unsafe_convert(Ptr{UInt8}, msg), UInt(sizeof(UInt8)))
    # the @eval macros are also necessary for this Compat layer to get the symbol name correct for the method
else
    import Base: unsafe_read, unsafe_write
end

@eval function $(Symbol(unsafe_write))(ctx::SSLContext, msg::Ptr{UInt8}, N::UInt)
    nw = 0
    while nw < N
        ret = ccall((:mbedtls_ssl_write, MBED_TLS), Cint,
                    (Ptr{Void}, Ptr{Void}, Csize_t),
                    ctx.data, msg, N - nw)
        ret < 0 && mbed_err(ret)
        nw += ret
        msg += ret
    end
    return Int(nw)
end

Base.write(ctx::SSLContext, msg::UInt8) = write(ctx, Ref(msg))

@eval function $(Symbol(unsafe_read))(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt; err=true)
    nread::UInt = 0
    while nread < nbytes
        n = ccall((:mbedtls_ssl_read, MBED_TLS), Cint,
                  (Ptr{Void}, Ptr{Void}, Csize_t),
                  ctx.data, buf + nread, nbytes - nread)
        if n == MBEDTLS_ERR_SSL_PEER_CLOSE_NOTIFY || n == 0
            ctx.isopen = false
            err ? throw(EOFError()) : return nread
        end
        if n != MBEDTLS_ERR_SSL_WANT_READ
            n < 0 && mbed_err(n)
            nread += n
        end
    end
end

Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes=length(buf)) = readbytes!(ctx, buf, UInt(nbytes))
function Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes::UInt)
    nr = unsafe_read(ctx, pointer(buf), nbytes; err=false)
    if nr !== nothing
        resize!(buf, nr::UInt)
    else
        nr = nbytes
    end
    return Int(nr::UInt)
end

function Base.readavailable(ctx::SSLContext)
    # For unknown reasons, nb_available on SSLContext erroneously returns 0
    # until `read` is called on the context. As a temporary hack, we read one
    # byte from the SSL context to cause nb_available to be accurate.
    # TODO: figure out and fix the root cause of this
    b = IOBuffer()
    write(b, read(ctx, 1))
    write(b, read(ctx, nb_available(ctx)))
    return takebuf_array(b)
end

function Base.eof(ctx::SSLContext)
    # Not quite semantically correct, since nb_available might still be zero
    # when this returns false (ie, the underlying socket has bytes available but not
    # a complete record)
    nb_available(ctx)>0 && return false
    return eof(ctx.bio)
end

function Base.close(ctx::SSLContext)
    if isopen(ctx.bio)
        try  # This is ugly, but a harmless broken pipe exception will be thrown if the peer closes the connection without responding
            ccall((:mbedtls_ssl_close_notify, MBED_TLS), Cint, (Ptr{Void},), ctx.data)
        catch
        end
        close(ctx.bio)
    end
    ctx.isopen = false
    nothing
end

Base.isopen(ctx::SSLContext) = ctx.isopen && isopen(ctx.bio)

function get_peer_cert(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_peer_cert, MBED_TLS), Ptr{Void}, (Ptr{Void},), ctx.data)
    return CRT(data)
end

function get_version(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_version, MBED_TLS), Ptr{UInt8}, (Ptr{Void},), ctx.data)
    return String(data)
end

function get_ciphersuite(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_ciphersuite, MBED_TLS), Ptr{UInt8}, (Ptr{Void},), ctx.data)
    return String(data)
end

function Base.nb_available(ctx::SSLContext)
    n = ccall((:mbedtls_ssl_get_bytes_avail, MBED_TLS), Csize_t, (Ptr{Void},), ctx.data)
    return Int(n)
end

function hostname!(ctx::SSLContext, hostname)
    @err_check ccall((:mbedtls_ssl_set_hostname, MBED_TLS), Cint,
      (Ptr{Void}, Cstring), ctx.data, hostname)
end

function __sslinit__()
    global const c_send = cfunction(f_send, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t))
    global const c_recv = cfunction(f_recv, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t))
    global const c_dbg = cfunction(f_dbg, Void,
        (Ptr{Void}, Cint, Ptr{UInt8}, Cint, Ptr{UInt8}))
    global const TRUSTED_CERT_FILE = joinpath(dirname(@__FILE__), "../deps/cacert.pem")
end


text:
abstract type RegularizerUnsupervised <: Regularizer end

##############################################################################
# todo: replace with function of xhat,x

#########################################
# Trivial Matrix Regularizer
#########################################
struct NoReg <: RegularizerUnsupervised
end

eval(R::NoReg, X) = 0
prox(R::NoReg, X) = X

#########################################
# Nonnegative Matrix Regularizer
#########################################
struct PosReg <: RegularizerUnsupervised
end

eval(R::PosReg, X) = any(X .< 0) ? Inf : 0

function prox(R::PosReg, X)
    Q = copy(X)
    Q[Q.<0] = 0
    return Q
end


text:
export to_array

import Base: show

"""
    to_array()

Creates a `to_array` operator, which reduces all values into a single array and returns this result when the source completes.

# Producing

Stream of type `<: Subscribable{Vector{L}}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> to_array(), logger())
;

# output

[LogActor] Data: [1, 2, 3]
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
to_array() = ToArrayOperator()

struct ToArrayOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{Vector{L}}, operator::ToArrayOperator, source) where L
    return proxy(Vector{L}, source, ToArrayProxy())
end

operator_right(operator::ToArrayOperator, ::Type{L}) where L = Vector{L}

struct ToArrayProxy <: ActorProxy end

actor_proxy!(::Type{Vector{L}}, proxy::ToArrayProxy, actor::A) where { L, A } = ToArrayActor{L, A}(actor)

struct ToArrayActor{L, A} <: Actor{L}
    values :: Vector{L}
    actor  :: A

    ToArrayActor{L, A}(actor::A) where { L, A } = new(Vector{L}(), actor)
end

on_next!(actor::ToArrayActor, data::L) where L = push!(actor.values, data)
on_error!(actor::ToArrayActor, err)            = error!(actor.actor, err)

function on_complete!(actor::ToArrayActor)
    next!(actor.actor, actor.values)
    complete!(actor.actor)
end

Base.show(io::IO, ::ToArrayOperator)         = print(io, "ToArrayOperator()")
Base.show(io::IO, ::ToArrayProxy)            = print(io, "ToArrayProxy()")
Base.show(io::IO, ::ToArrayActor{L}) where L = print(io, "ToArrayActor($L)")


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

using LinearAlgebra

# resultant number of clusters may differ from the cohort parametrization in the likelihood
# this feels more salient than a relabeling algorithm using anchor examples... what if the anchors are weird "outliers"?
function count_coincidences!(coincident::AbstractMatrix{Int}, # can be sparse for huge voter bases?
    voters::Vector{VoterRealization{N,M,T}})::Nothing where {N,M,T}
  n_voters = length(voters)
  for first_index in 1:n_voters
    first_cohort = voters[first_index].membership # don't thrash the cache by indexing this in the inner loop
    for second_index in 1:n_voters
      second_cohort = voters[second_index].membership
      if first_cohort == second_cohort
        coincident[second_index, first_index] += 1
      end
    end
  end
end

# our predictive posterior (distribution of voters marginalized over Gibbs samples)
# no longer has `n_cohorts` distinct components; in fact, it comprises arbitrarily
# many humps integrated to a smooth surface
# so, why relabel? also, I disagree philosophically with most methods that use ordinary
# points as exemplars or anchors for clusters

using Clustering

# try MCL, which normalizes automagically, and then maybe Affinity Propagation.
# or diffusion maps..
# if we have Likert answers, can we score that against a Markov transition matrix?
# maybe remove diagonal and then normalize rows to sum to unit, and do a knn-style
# classifier (not real-valued prediction?)


text:
"""
    reduce_axes(a, dims)

Returns the appropriate axes for a measure that reduces dimensions along the
dimensions `dims`.

## Example
```jldoctest
julia> reduce_axes((Index{:a}(1:4), Index{:b}(1:4)), 2)
(Index{a}(1:4 => Base.OneTo(4)), Index{b}(1:1 => Base.OneTo(1)))

julia> reduce_axes((Index{:a}(1:4), Index{:b}(1:4)), :a)
(Index{a}(1:1 => Base.OneTo(1)), Index{b}(1:4 => Base.OneTo(4)))
```
"""
reduce_axes(x::AbstractArray, dims) = reduce_axes(axes(x), dims)
reduce_axes(x::Tuple, dims) = _reduce_axes(x, to_dims(x, dims))
reduce_axes(x::Tuple, dims::Colon) = ()
_reduce_axes(x::Tuple{Vararg{Any,D}}, dims::Int) where {D} = _reduce_axes(x, (dims,))
function _reduce_axes(x::Tuple{Vararg{Any,D}}, dims::Tuple{Vararg{Int}}) where {D}
    Tuple(map(i -> ifelse(in(i, dims), reduce_axis(x[i]), x[i]), 1:D))
end

"""
    reduce_axis(a)

Reduces axis `a` to single value. Allows custom index types to have custom
behavior throughout reduction methods (e.g., sum, prod, etc.)

See also: [`reduce_axes`](@ref)

## Example
```jldoctest
julia> reduce_axis(Index{:a}(1:4))
Index{a}(1:1 => Base.OneTo(1))

julia> reduce_axis(1:4)
1:1
```
"""
function reduce_axis(x::AbstractIndex)
    if isempty(x)
        error("Cannot reduce empty index.")
    else
        return unsafe_reindex(x, 1:1)
    end
end
reduce_axis(x::OneTo{T}) where {T} = OneTo(one(T))
reduce_axis(x::OneToSRange{T}) where {T} = OneToSRange(one(T))
reduce_axis(x::OneToMRange{T}) where {T} = OneToMRange(one(T))
reduce_axis(x::UnitRange{T}) where {T} = UnitRange{T}(one(T), one(T))
reduce_axis(x::UnitSRange{T}) where {T} = UnitSRange{T}(one(T), one(T))
reduce_axis(x::UnitMRange{T}) where {T} = UnitMRange{T}(one(T), one(T))


text:
using VLMinorityGameKit

# what is my memory?
memory = 3

# build the design -
design_array = design(memory)


text:
@testset "564.Find-the-Closest-Palindrome.jl" begin
    @test nearest_palindromic("1231") == "1221"
    @test nearest_palindromic("99") == "101"
    @test nearest_palindromic("100") == "99"
    @test nearest_palindromic("900") == "898"
    @test nearest_palindromic("101") == "99"
    @test nearest_palindromic("1") == "0"
    @test nearest_palindromic("10") == "9"
    @test nearest_palindromic("123") == "121"
end


text:
struct StableNonLinearLESModel{T} <: AbstractLESModel
    c::T
    Δ²::T
    tau::SymTrTenField{T,3,2,false}
    reduction::Vector{T}
end

function StableNonLinearLESModel(c::T,Δ::Real,dim::NTuple{3,Integer}) where {T<:Real}
    data = SymTrTenField{T}(dim,(LX,LY,LZ))
    #fill!(data,0)
    reduction = zeros(THR ? Threads.nthreads() : 1)
    return StableNonLinearLESModel{T}(c,Δ^2, data,reduction)
end

is_stable_nl_model(a::Type{T}) where {T<:StableNonLinearLESModel} = true
@inline @par is_stable_nl_model(s::Type{T}) where {T<:@par(AbstractSimulation)} = is_stable_nl_model(LESModelType)
@inline is_stable_nl_model(s::T) where {T<:AbstractSimulation} = is_stable_nl_model(T)
is_stable_nl_model(a) = false

msg(a::StableNonLinearLESModel) = "\nLES model: Stable Non-linear Model\nConstant: $(a.c)\nFilter Width: $(sqrt(a.Δ²))\n"

@inline function stable_nl_eddy_viscosity(S::SymTen,W::AntiSymTen,c::Real,Δ2::Real)
    a = S:S
    nu = ifelse(a == 0.,zero(eltype(a)),
         @fastmath max(0.0,c*Δ2*((square(W):S) - (square(S):S))/(2a)))
    mnu = ifelse(a == 0.,zero(eltype(a)),
         @fastmath -min(0.0,c*Δ2*((square(W):S) - (square(S):S))/(2a)))
    return nu,mnu
end


text:
using RiemannHilbert, ApproxFun, ApproxFunBase, LinearAlgebra, SingularIntegralEquations, Test
import ApproxFunBase: pieces

@testset "NLS" begin
    x,t = 0.0,0.0
    r = k -> 0.5exp(-k^2)
    r̄ = k -> conj(r(conj(k)))
    θ = k -> 4t*k^2+2x*k
    Gf = k -> [1-abs2(r(k))                      -conj(r(k))*exp(-im*θ(k));
                    r(k)*exp(im*θ(k))      1  ]

    g = Fun(Gf, -5.5..5.5)
    @test first(g) ≈ last(g) ≈ [1 0; 0 1]
    φ = transpose(rhsolve(transpose(g), 500))
    istieltjes(φ)(0.1)
    @test φ(1+im) ≈ [0.9847854664616504-0.014938097892888065im -0.0778564768587141-0.05302708037085426im;
                0.07824233623082648+0.05395699475793745im 1.0120649351483406+0.006873030868410332im] # mathematica

    x,t = 1.0,1.0

    g = Fun(Gf, -5.5..5.5)
    @test first(g) ≈ last(g) ≈ [1 0; 0 1]
    φ = transpose(rhsolve(transpose(g), 1000))
    @test φ(1+im) ≈ [0.9952708020858667-0.00828828357987084im -0.036773740429971044-0.01579754685573759im;
                0.009328821936648003+0.04514458247615043im 1.0050689684367242+ 0.006553783105674419im]

    @test norm(φ.coefficients[end-50:end]) ≤ 1E-14
    @test φ(0.1⁺) ≈ φ(0.1⁻)g(0.1)

    z₀ = -x/(4t)
    Γ = Segment(z₀,z₀+3exp(3im*π/4)) ∪ Segment(z₀,-5.5) ∪ Segment(z₀,z₀+3exp(-3im*π/4)) ∪ Segment(z₀,z₀+3exp(-im*π/4)) ∪ Segment(z₀,z₀+3exp(im*π/4))

    U = z -> [1 -r̄(z)/(1-r(z)r̄(z))*exp(-im*θ(z)); 0 1]
    D = z -> [(1-r(z)r̄(z)) 0; 0 1/(1-r(z)r̄(z))]
    L = z -> [1 0; r(z)/(1-r(z)r̄(z))*exp(im*θ(z)) 1]
    M = z -> [1 -r̄(z)*exp(-im*θ(z)); 0 1]
    P = z -> [1 0; r(z)*exp(im*θ(z)) 1]
    G = Fun(function(z)
            z ∈ component(Γ,1) && return inv(U(z))
            z ∈ component(Γ,2) && return inv(D(z))
            z ∈ component(Γ,3) && return inv(L(z))
            z ∈ component(Γ,4) && return M(z)
            z ∈ component(Γ,5) && return P(z)
            error("Not in contour")
        end, Γ)

    @test prod(first.(pieces(G))) ≈ [1 0; 0 1]
    @test norm(last.(pieces(G)) .- Ref(I)) ≤ 1E-14

    Φ = transpose(rhsolve(transpose(G), 3000))
    @test norm(Φ.coefficients[end-50:end]) ≤ 10E-14

    @test Φ(2im) ≈ φ(2im)

    z = Fun(ℂ)
    2im*(z*Φ[1,2])(Inf)≈ 2im*(z*φ[1,2])(Inf)

    # make evaluation faster
    φ2 = Matrix(φ)
    φ̃ = z -> [φ2[1,1](z) φ2[1,2](z); φ2[2,1](z) φ2[2,2](z)]

    @time V = Fun(function(z)
            z ∈ component(Γ,1) && return φ̃(z)*(inv(U(z))-I)
            z ∈ component(Γ,2) && return φ̃(z-10eps()im)L(z)-φ̃(z+10eps()im)inv(U(z))
            z ∈ component(Γ,3) && return φ̃(z)*(I-L(z))
            z ∈ component(Γ,4) && return φ̃(z)*(M(z)-I)
            z ∈ component(Γ,5) && return φ̃(z)*(I-inv(P(z)))
            error("Not in contour")
        end, Γ, 1000)

    @test norm(V.coefficients[end-20:end]) ≤ 1E-15
    @test P(z₀+3exp(im*π/4)) ≈ [1 0; 0 1]
    Φ = I+cauchy(V)
    Φ̃ = function(z)
            3π/4  < angle(z-z₀) < π     && return φ(z)inv(U(z))
            -π    < angle(z-z₀) < -3π/4 && return φ(z)L(z)
            -3π/4 < angle(z-z₀) < -π/4  && return φ(z)
            -π/4  < angle(z-z₀) <  0    && return φ(z)M(z)
            0    < angle(z-z₀) <  π/4  && return φ(z)inv(P(z))
            π/4  < angle(z-z₀) < 3π/4  && return φ(z)
            error("Not defined")
        end

    @test Φ̃(0.1+10eps()im) ≈ Φ̃(0.1-10eps()im)

    z = z₀+0.1exp(im*π/4); @test Φ(z-10eps())-Φ(z+10eps()) ≈ Φ̃(z-10eps())-Φ̃(z+10eps()) ≈ φ̃(z)*(I-inv(P(z)))
    z = z₀+0.1exp(3im*π/4); @test Φ(z-10eps())-Φ(z+10eps()) ≈ Φ̃(z-10eps())-Φ̃(z+10eps()) ≈ φ̃(z)*(inv(U(z))-I)
    z = z₀-0.1+0im; @test Φ(z-10im*eps())-Φ(z+10im*eps()) ≈ Φ̃(z-10im*eps())-Φ̃(z+10im*eps()) ≈ φ̃(z-10im*eps())*L(z)-φ̃(z+10im*eps())*inv(U(z))
    z = z₀+0.1exp(-3im*π/4); @test Φ(z+10eps())-Φ(z-10eps()) ≈ Φ̃(z+10eps())-Φ̃(z-10eps()) ≈ φ̃(z)*(I-L(z))
    z = z₀+0.1exp(-im*π/4); @test Φ(z+10eps())-Φ(z-10eps()) ≈ Φ̃(z+10eps())-Φ̃(z-10eps()) ≈ φ̃(z)*(M(z)-I)
    z = z₀+0.1+0im; @test Φ(z-10im*eps()) ≈ Φ(z+10im*eps()) ≈ Φ̃(z-10im*eps()) ≈ Φ̃(z+10im*eps())
end


text:
using PyPlot, Suppressor, Printf, ProgressMeter
using LinearAlgebra, Distributions, Statistics, StatsFuns, Random
using Distances
using Optim
using Interpolations
using PDMats
using ApproximateVI
using ForwardDiff

# Utilities

include("UTIL/gputil.jl")
include("UTIL/smallutil.jl")
include("UTIL/astroutil.jl")
lsstwaves = setlsstwavelengths()
@info(@sprintf("setting LSST wavelengths to \n"))
display(lsstwaves())

# Synthetic data generation

include("SyntheticDataset/simulatedata.jl")

# GP approach

include("GP/complete_lower_bound.jl")
include("GP/exact_expectations.jl")
include("GP/complete_lower_bound_grad.jl")
include("GP/exact_expectations_grad.jl")
include("GP/gpmodel.jl")
include("GP/gpmodelVIexact.jl")

# Line intersection

include("LineIntersection/noisyintersection.jl")
include("LineIntersection/linemarginal.jl")

# Probabilistic principal component algorithms

include("PPCA/bpca.jl")
include("PPCA/ppca.jl")


text:
using PyCall
@pyimport matplotlib.pyplot as plt

@doc """
Initial conditions
""" ->
function from_prior(num::Int64)
	@assert num > 1
	x = rand(num)
	y = x.^2 + 0.1*randn(num)
	return hcat(x, y)
end

@doc """
Measure upper corner count of a single point
""" ->
function ucc(x::Float64, y::Float64, X::Vector{Float64}, Y::Vector{Float64})
	@assert length(X) == length(Y)
	count = 0
	for(i in 1:length(X))
		if((X[i] > x) && (Y[i] > y))
			count += 1
		end
	end
	return count
end

@doc """
Measure upper corner counts of a bunch of points
""" ->
function uccs(X::Vector{Float64}, Y::Vector{Float64})
	@assert length(X) == length(Y)
	counts = Array(Int64, (length(X), ))
	for(i in 1:length(X))
		counts[i] = ucc(X[i], Y[i], X, Y)
	end
	return counts
end

# Generate some points
params = from_prior(1000)

# Assign x and y
xy = -log(rand(size(params)))

# Plot
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(params[:,1], params[:,2], "b.")

plt.subplot(1, 2, 2)
plt.plot(xy[:,1], xy[:,2], "b.")
plt.show()


text:
using Optim

export solve

include("policyiteration.jl") # Optimize over control in each timestep
include("policytimestep.jl")  # Choose from discrete set of controls in each timestep

function solveconstant{T<:Real}(model::HJBOneDim{T}, K::Int, N::Int, M::Int)
    # TODO: better names
    # K = number of points in space domain
    # N+1 = number of points in time domain
    # M   = number of control values
    x = linspace(model.xmin, model.xmax, K)
    avals = linspace(model.amin, model.amax, M)
    Δx = (model.xmax-model.xmin)/(K-1) # TODO: use diff(x) to accomodate non-uniform grid
    Δτ = model.T/N # TODO: introduce non-uniform timesteps?

    vinit = zeros(x)
    for i = 1:length(x)
        vinit[i] = model.g(x[i])
    end

    v, pol = timeloopconstant(model, N, Δτ, vinit, avals, x, Δx)
    return v, pol
end

function solveiteration(model::HJBOneDim, K::Int, N::Int)
    # K = number of points in space domain
    # N+1 = number of points in time domain
    x = linspace(model.xmin, model.xmax, K)
    Δx = (model.xmax-model.xmin)/(K-1) # TODO: use diff(x) to accomodate non-uniform grid
    Δτ = model.T/N # TODO: introduce non-uniform timesteps?

    vinit = zeros(x)
    for i = 1:length(x)
        vinit[i] = model.g(x[i])
    end

    v, pol = timeloopiteration(model, N, Δτ, vinit, x, Δx)
    return v, pol
end

# TODO: Create a better interface than this,
# maybe with a discretisation object that contains the x,a,t-arrays etc?
function solve(model::HJBOneDim, K::Int, N::Int)
    solveiteration(model, K, N)
end

function solve(model::HJBOneDim, K::Int, N::Int, M::Int)
    solveconstant(model, K, N, M)
end

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
using Distributions
using Distances
using Test
using StatsBase
using ApproxBayes
using Plots
using Random

#useful to check what the backend is
println(backend())
tests = ["sampling", "parameter", "parameter_parallel", "bayesfactor", "util"]

println("Running tests ...")

@time for t in tests
    fn = "test_$t.jl"
    println("* $fn ...")
    include(fn)
end


text:
#=
 TODO try to get rid of these two kernels
=#

function _mul_ld_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    x       ::CuDeviceVector{T}
) where T
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    ξ = Ξ[i]
    numelem = length(elements)
    val12 = T(0)
    val3  = T(0)
    for j in 1:numelem
        ld = laplacepot_double(ξ, elements[j])

        val12 = CUDA.fma(ld, x[j], val12)
        val3  = CUDA.fma(ld, x[j + 2numelem], val3)
    end
    dst[i] = -val12
    dst[i + numelem] = val12
    dst[i + 2numelem] = -val3
    nothing
end

function _mul_yd_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    x       ::CuDeviceVector{T},
    pref    ::T,
    yuk     ::T
) where T
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    ξ = Ξ[i]
    numelem = length(elements)
    val1 = T(0)
    val2 = T(0)
    for j in 1:numelem
        yd = regularyukawapot_double(ξ, elements[j], yuk)

        val1 = CUDA.fma(yd, x[j], val1)
        val2 = CUDA.fma(yd, x[j + 2numelem], val2)
    end
    dst[i] = pref * val2 - val1
    nothing
end