text:
module SeisNoiseAWS

using Logging
using CSV
using DataFrames
using Dates
using SeisIO
using SeisNoise

# The functions we want to export out of the module
export load_parameter_file
export load_station_file
export build_availability_df
export xcor

# ==============================================================================
# The xcor() parameters
# fs      = df_params.fs[1] # [Hz] resample frequency
# cc_len  = df_params.cc_len[1] # [s]
# freqmin = df_params.freqmin[1] # [Hz] low pass
# cc_step = df_params.cc_step[1] # [s]
# freqmax = df_params.freqmax[1] # [Hz] high pass
# maxlag  = df_params.maxlag[1] # [s]
function xcor(xcor_pairs, df_params, output_dir)

    println("Correlating ", xcor_pairs[1].id," and ", xcor_pairs[2].id)

    # demean, detrend, downsample, time align, taper
    S1 = process_raw( SeisData(xcor_pairs[1]), df_params.fs[1])
    S2 = process_raw( SeisData(xcor_pairs[2]), df_params.fs[1])

    # Split data into windows and try to correlate
    try
        R = RawData.([S1,S2], df_params.cc_len[1], df_params.cc_step[1])

        # Process each window
        detrend!.(R)
        taper!.(R)
        bandpass!.(R, df_params.freqmin[1], df_params.freqmax[1], zerophase=true)

        # Compute the FFT and correlation in frequency domain
        FFT = compute_fft.(R)
        whiten!.(FFT, df_params.freqmin[1], df_params.freqmax[1])
        C = compute_cc(FFT[1], FFT[2], df_params.maxlag[1])

        # demean, detrend, taper, filter, and then save
        clean_up!(C, df_params.freqmin[1], df_params.freqmax[1])
        save_corr(C, output_dir) # Write the complete correlation object

        # abs_max!(C)
        # corrplot(C)
        # stack!( C )

        # deallocate (not sure how much this helps or not)
        R = nothing
        FFT = nothing
        C = nothing

    catch err
        # DomainError: a station has less data than the xcor length
        if isa(err, DomainError)
            @info(err)
            s1_len = length(xcor_pairs[1].x) / xcor_pairs[1].fs
            s2_len = length(xcor_pairs[2].x) / xcor_pairs[2].fs
            @info("Did not correlate: $(xcor_pairs[1].id) len=$s1_len [s], $(xcor_pairs[2].id) len=$s2_len [s]")
        elseif isa(err, ArgumentError) #  - because date folder already exists
            @info(err)
            # Do nothing for now
        else
            @info("Unkown error. May need to look into this!", err)
        end

    finally
        S1 = nothing
        S2 = nothing
    end
end
# ==============================================================================
function load_parameter_file(filename)
    @info("Loading: $filename")
    @assert(isfile(filename),"Parameter file does not exist.")
    df = DataFrame( CSV.File(filename, missingstring="", delim=",", header=1, type=Float64, normalizenames=true) )
    @debug("Done")
    show(df,summary=false); println()
    return df
end
# ==============================================================================
function load_station_file(filename)
    @info("Loading: $filename")
    @assert(isfile(filename),"Station file does not exist.")
    df = DataFrame( CSV.File(filename, missingstring="", delim=",", header=1, normalizenames=true) )
    @debug("Done")
    show(df,summary=false); println()
    return df
end
# ==============================================================================
function build_availability_df(df_st)
    # Compute the station tags
    station_tag = build_station_tags(df_st)
    df_date = fill_days(df_st, station_tag)
    return df_date
end
# ==============================================================================
function build_cor_directory(input_dir)
    # Save correlation files to this directory
    COR_DIR = joinpath(input_dir, "CORR") # local path
    if !isdir(COR_DIR) # make the local directoy if necessary
        mkpath(COR_DIR)
    end
    @debug("Saving correlations to $COR_DIR")
    return COR_DIR
end

# ==============================================================================
# ==============================================================================
# Internal functions to those in this module
# ==============================================================================
# ==============================================================================
function build_station_tags(df_st)

    num_data    = nrow(df_st) # number of data in the IRIS query (channels are separate)
    station_tag = Array{String,1}(undef,num_data) # allocate station_tag

    for ii in 1:num_data
        if ismissing(df_st.Location[ii])
            station_tag[ii] = string(df_st.Network[ii], ".", df_st.Station[ii], ".", ".", df_st.Channel[ii]) # create the station ID
        else
            station_tag[ii] = string(df_st.Network[ii], ".", df_st.Station[ii], ".", df_st.Location[ii], ".", df_st.Channel[ii]) # create the station ID
        end
        @debug("Creating tag: ", station_tag[ii])
    end
    @info("Created tags: $station_tag")
    return station_tag
end
# ==============================================================================
function fill_days(df_st, station_tag)

    # Get the date vector
    date_range = build_date_range(df_st)
    num_days   = length(date_range)
    num_data   = nrow(df_st) # number of data in the IRIS query (channels are separate)

    # Create the station vs. date dataframe
    df_date = DataFrame() # create empy df
    df_date[!, Symbol("station_tag")] = station_tag # column 1 is the station_tag
    # add the other column names and set column values
    for d in date_range
        df_date[!, Symbol(d)] .= false
    end

    st_ends   = Array{Date,1}(df_st.EndDate) # start times of each station_tag
    st_starts = Array{Date,1}(df_st.StartDate) # start times of each station_tag

    # Now populate each row with TRUE/FALSE
    for ii in 1:num_data
        @debug("Adding data days for: ", station_tag[ii])
        idx = ( (st_starts[ii] .<= date_range) + (st_ends[ii] .>= date_range) .== 2)
        df_date[ii, 2:num_days+1] .= idx
    end

    return df_date
end
# ==============================================================================
function build_date_range(df_st)
    # Use oldest to youngest dates for the date df
    date_range = collect( findmin( df_st.StartDate )[1] : Dates.Day(1) : findmax( df_st.EndDate )[1])
    num_days = length(date_range) # number of days of data
    @info("Total number of days of data: $num_days")
    return date_range
end

end


text:
θ_conv_variable_Pr = (
                      Cᵂu★  = 5.18,
                      CᵂwΔ  = 7.37,
                      Cᴷu⁻  = 0.0939,
                      Cᴷc⁻  = 0.835,
                      Cᴷe⁻  = 1.04,
                      Cᴷuʳ  = -0.508,
                      Cᴷcʳ  = +0.256,
                      Cᴷeʳ  = -0.0696,
                      Cᴬc   = 0.544,
                      Cᴬe   = 9.77e-3,
                      Cᴰ    = 0.283,
                      Cᴸᵇ   = 0.0561,
                      Cᴬu   = 0.0,
                      CᴷRiʷ = 0.0,
                      CᴷRiᶜ = 0.0,
                      Cᵟu   = 0.5,
                      Cᵟc   = 0.5,
                      Cᵟe   = 0.5)

θ_variable_Pr = (
                 CᵂwΔ  = 14.1,
                 Cᵂu★  = 2.55,
                 Cᴰ    = 0.354,
                 Cᴸᵇ   = 0.159,
                 Cᴷu⁻  = 0.0963,
                 Cᴷc⁻  = 0.279,
                 Cᴷe⁻  = 1.03,
                 Cᴷuʳ  = -0.525,
                 Cᴷcʳ  = +0.670,
                 Cᴷeʳ  = +1.14,
                 CᴷRiʷ = 0.111,
                 CᴷRiᶜ = 1.98,
                 Cᴬu   = 0.0,
                 Cᴬc   = 0.0,
                 Cᴬe   = 0.0,
                 Cᵟu   = 0.5,
                 Cᵟc   = 0.5,
                 Cᵟe   = 0.5,
                )

θ_constant_Pr = (
                 CᵂwΔ  = 8.66,
                 Cᵂu★  = 2.89,
                 Cᴰ    = 0.489,
                 Cᴸᵇ   = 0.0286,
                 Cᴷu⁻  = 0.0776,
                 Cᴷc⁻  = 0.567,
                 Cᴷe⁻  = 9.00,
                 Cᴷuʳ  = 0.0,
                 Cᴷcʳ  = 0.0,
                 Cᴷeʳ  = 0.0,
                 CᴷRiʷ = 0.0,
                 CᴷRiᶜ = 0.0,
                 Cᴬu   = 0.0,
                 Cᴬc   = 0.0,
                 Cᴬe   = 0.0,
                 Cᵟu   = 0.5,
                 Cᵟc   = 0.5,
                 Cᵟe   = 0.5,
                )

ext: jl

max_stars_repo_name: glwagner/LocalOceanClosureCalibration


text:
function mysum(arr)
    s = 0.0

    for i in arr
        s+=i
    end
    s #最后一行默认返回值
end

println(rand(Bool, 3, 3)) #矩阵

A = rand(1000)
println(mysum(A))

arr = rand(Int,10)
println(mysum(arr))

'''
    is_even(x::Int) -> Bool
判断一个整数`x`是不是偶数

'''
is_even(x::Int)= x % 2 == 0
println(is_even(7))

println(is_even(2.0)) # 强类型

max_stars_repo_name: liushooter/OneDayOneCommit

max_issues_repo_name: liushooter/OneDayOneCommit

max_forks_repo_name: liushooter/OneDayOneCommit


text:
abstract type AbstractKullbackLeiblerProx end

struct FixedMarginal{D<:DiscreteNonParametric} <: AbstractKullbackLeiblerProx
    d::D
end

import Base:length

function length(fm::FixedMarginal)
    length(support(fm.d))
end

# u is prob vector
(fm::FixedMarginal)(u) = probs(fm.d)

function FixedMarginal(v::AbstractVector)
    n = length(v)
    FixedMarginal(DiscreteNonParametric(sort(v), fill(1/n, n)))
end

function sinkhorn(dist_mat, kl_prox1, kl_prox2;
                   ϵ = 1.0,
                   rounds=2,
                   u = fill(1.0, length(kl_prox1)),
                   v = fill(1.0, length(kl_prox2)))

     K = exp.(-dist_mat / ϵ)
     Kv = K*v
     Ktu = K'*u

     for iter=1:rounds
             u = kl_prox1(Kv) ./ Kv
             Ktu = K'*u
             v = kl_prox2(Ktu) ./ Ktu
             Kv = K*v
     end
     (u=u, v=v, plan=u .* K .* v')
end

Base.@kwdef mutable struct FixedMarginalConvolution{ZT <: DiscretizedStandardNormalSamples,
                                            GCAL <: ConvexPriorClass}
    Zs_discr::ZT
    prior_class::GCAL
    model = nothing
end

length(fmc::FixedMarginalConvolution) = nbins(fmc.Zs_discr)

function initialize_dkl_projection(prior_class::ConvexPriorClass,
	                            Zs_discr::DiscretizedStandardNormalSamples)

								mproj = Model(prior_class.solver)

								π = add_prior_variables!(mproj, prior_class)
								f = marginalize(prior_class, Zs_discr, π)

								@variable(mproj, entropy_upper_bd_vec[1:length(f)])
								@variable(mproj, prob_vec_dummy[1:nbins(Zs_discr)])

								for i=1:length(f)
								         @constraint(mproj,
								         [entropy_upper_bd_vec[i]; f[i]; prob_vec_dummy[i]] in MathOptInterface.ExponentialCone())
								end

								@objective(mproj, Min, -sum(entropy_upper_bd_vec ))
								mproj
end

function (fmc::FixedMarginalConvolution)(a)
	fmc(a, fmc.model)
end

# mutating and not clear from type...
function (fmc::FixedMarginalConvolution)(a, ::Nothing)
	mproj = initialize_dkl_projection(fmc.prior_class, fmc.Zs_discr)
	fmc.model = mproj
	fmc(a, fmc.model)
end

function (fmc::FixedMarginalConvolution)(a, mproj::JuMP.Model)
	@show sum(a)
	fix.(mproj[:prob_vec_dummy], a)
	optimize!(mproj)
	prior_vals = JuMP.value.(mproj[:πs])
	@show sum(prior_vals)
	if isapprox(sum(prior_vals) , 1.0; atol=0.001)
		prior_vals = max.(prior_vals, 0.0)
		prior_vals = prior_vals/sum(prior_vals)
	end
	prior_res = fmc.prior_class(prior_vals)
	pdf( marginalize(prior_res, fmc.Zs_discr))
end


text:
int_rules_4_1_10 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.10*(c+d*x)^m*(a+b*sin)^n =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(e') + ~(f') * ~x), ~x) => (-((~c + ~d * ~x) ^ ~m) * cos(~e + ~f * ~x)) / ~f + ((~d * ~m) / ~f) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * cos(~e + ~f * ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && GtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * sin(~(e') + ~(f') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * sin(~e + ~f * ~x)) / (~d * (~m + 1)) - (~f / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * cos(~e + ~f * ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && LtQ(~m, -1)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * Complex(0, ~fz) * ~x) / (~(c') + ~(d') * ~x), ~x) => (im * SinhIntegral((~c * ~f * ~fz) / ~d + ~f * ~fz * ~x)) / ~d <-- FreeQ([~c, ~d, ~e, ~f, ~fz], ~x) && EqQ(~d * ~e - ~c * ~f * ~fz * im, 0)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) / (~(c') + ~(d') * ~x), ~x) => SinIntegral(~e + ~f * ~x) / ~d <-- FreeQ([~c, ~d, ~e, ~f], ~x) && EqQ(~d * ~e - ~c * ~f, 0)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * Complex(0, ~fz) * ~x) / (~(c') + ~(d') * ~x), ~x) => CoshIntegral((-(~c) * ~f * ~fz) / ~d - ~f * ~fz * ~x) / ~d <-- FreeQ([~c, ~d, ~e, ~f, ~fz], ~x) && (EqQ(~d * (~e - pi / 2) - ~c * ~f * ~fz * im, 0) && NegQ((~c * ~f * ~fz) / ~d, 0))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * Complex(0, ~fz) * ~x) / (~(c') + ~(d') * ~x), ~x) => CoshIntegral((~c * ~f * ~fz) / ~d + ~f * ~fz * ~x) / ~d <-- FreeQ([~c, ~d, ~e, ~f, ~fz], ~x) && EqQ(~d * (~e - pi / 2) - ~c * ~f * ~fz * im, 0)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) / (~(c') + ~(d') * ~x), ~x) => CosIntegral((~e - pi / 2) + ~f * ~x) / ~d <-- FreeQ([~c, ~d, ~e, ~f], ~x) && EqQ(~d * (~e - pi / 2) - ~c * ~f, 0)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) / (~(c') + ~(d') * ~x), ~x) => cos((~d * ~e - ~c * ~f) / ~d) * Antiderivative(sin((~c * ~f) / ~d + ~f * ~x) / (~c + ~d * ~x), ~x) + sin((~d * ~e - ~c * ~f) / ~d) * Antiderivative(cos((~c * ~f) / ~d + ~f * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && NeQ(~d * ~e - ~c * ~f, 0)
    @apply_utils Antiderivative(sin(~(e') + pi / 2 + ~(f') * ~x) / sqrt(~(c') + ~(d') * ~x), ~x) => (2 / ~d) * Subst(Antiderivative(cos((~f * (~x) ^ 2) / ~d), ~x), ~x, sqrt(~c + ~d * ~x)) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && (ComplexFreeQ(~f) && EqQ(~d * ~e - ~c * ~f, 0))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) / sqrt(~(c') + ~(d') * ~x), ~x) => (2 / ~d) * Subst(Antiderivative(sin((~f * (~x) ^ 2) / ~d), ~x), ~x, sqrt(~c + ~d * ~x)) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && (ComplexFreeQ(~f) && EqQ(~d * ~e - ~c * ~f, 0))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) / sqrt(~(c') + ~(d') * ~x), ~x) => cos((~d * ~e - ~c * ~f) / ~d) * Antiderivative(sin((~c * ~f) / ~d + ~f * ~x) / sqrt(~c + ~d * ~x), ~x) + sin((~d * ~e - ~c * ~f) / ~d) * Antiderivative(cos((~c * ~f) / ~d + ~f * ~x) / sqrt(~c + ~d * ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && (ComplexFreeQ(~f) && NeQ(~d * ~e - ~c * ~f, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(e') + ~(k') * pi + ~(f') * ~x), ~x) => (im / 2) * Antiderivative((~c + ~d * ~x) ^ ~m * ℯ ^ (-im * ~k * pi) * ℯ ^ (-im * (~e + ~f * ~x)), ~x) - (im / 2) * Antiderivative((~c + ~d * ~x) ^ ~m * ℯ ^ (im * ~k * pi) * ℯ ^ (im * (~e + ~f * ~x)), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~m], ~x) && IntegerQ(2 * ~k)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(e') + ~(f') * ~x), ~x) => (im / 2) * Antiderivative((~c + ~d * ~x) ^ ~m * ℯ ^ (-im * (~e + ~f * ~x)), ~x) - (im / 2) * Antiderivative((~c + ~d * ~x) ^ ~m * ℯ ^ (im * (~e + ~f * ~x)), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~m], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(e') + (~(f') * ~x) / 2) ^ 2, ~x) => (1 / 2) * Antiderivative((~c + ~d * ~x) ^ ~m, ~x) - (1 / 2) * Antiderivative((~c + ~d * ~x) ^ ~m * cos(2 * ~e + ~f * ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~m], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~d * (~b * sin(~e + ~f * ~x)) ^ ~n) / ((~f) ^ 2 * (~n) ^ 2) - (~b * (~c + ~d * ~x) * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~n)) + (((~b) ^ 2 * (~n - 1)) / ~n) * Antiderivative((~c + ~d * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && GtQ(~n, 1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~d * ~m * (~c + ~d * ~x) ^ (~m - 1) * (~b * sin(~e + ~f * ~x)) ^ ~n) / ((~f) ^ 2 * (~n) ^ 2) - (~b * (~c + ~d * ~x) ^ ~m * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~n)) + (((~b) ^ 2 * (~n - 1)) / ~n) * Antiderivative((~c + ~d * ~x) ^ ~m * (~b * sin(~e + ~f * ~x)) ^ (~n - 2), ~x)) - (((~d) ^ 2 * ~m * (~m - 1)) / ((~f) ^ 2 * (~n) ^ 2)) * Antiderivative((~c + ~d * ~x) ^ (~m - 2) * (~b * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (GtQ(~n, 1) && GtQ(~m, 1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * sin(~(e') + ~(f') * ~x) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce((~c + ~d * ~x) ^ ~m, sin(~e + ~f * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~n, 1) && (Not(RationalQ(~m)) || GeQ(~m, -1) && LtQ(~m, 1)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * sin(~(e') + ~(f') * ~x) ^ ~n, ~x) => ((~c + ~d * ~x) ^ (~m + 1) * sin(~e + ~f * ~x) ^ ~n) / (~d * (~m + 1)) - ((~f * ~n) / (~d * (~m + 1))) * Antiderivative(ExpandTrigReduce((~c + ~d * ~x) ^ (~m + 1), cos(~e + ~f * ~x) * sin(~e + ~f * ~x) ^ (~n - 1), ~x), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~n, 1) && (GeQ(~m, -2) && LtQ(~m, -1)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((((~c + ~d * ~x) ^ (~m + 1) * (~b * sin(~e + ~f * ~x)) ^ ~n) / (~d * (~m + 1)) - (~b * ~f * ~n * (~c + ~d * ~x) ^ (~m + 2) * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n - 1)) / ((~d) ^ 2 * (~m + 1) * (~m + 2))) - (((~f) ^ 2 * (~n) ^ 2) / ((~d) ^ 2 * (~m + 1) * (~m + 2))) * Antiderivative((~c + ~d * ~x) ^ (~m + 2) * (~b * sin(~e + ~f * ~x)) ^ ~n, ~x)) + (((~b) ^ 2 * (~f) ^ 2 * ~n * (~n - 1)) / ((~d) ^ 2 * (~m + 1) * (~m + 2))) * Antiderivative((~c + ~d * ~x) ^ (~m + 2) * (~b * sin(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (GtQ(~n, 1) && LtQ(~m, -2))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~c + ~d * ~x) * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) - (~d * (~b * sin(~e + ~f * ~x)) ^ (~n + 2)) / ((~b) ^ 2 * (~f) ^ 2 * (~n + 1) * (~n + 2))) + ((~n + 2) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~c + ~d * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (LtQ(~n, -1) && NeQ(~n, -2))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~c + ~d * ~x) ^ ~m * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) - (~d * ~m * (~c + ~d * ~x) ^ (~m - 1) * (~b * sin(~e + ~f * ~x)) ^ (~n + 2)) / ((~b) ^ 2 * (~f) ^ 2 * (~n + 1) * (~n + 2))) + ((~n + 2) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ ~m * (~b * sin(~e + ~f * ~x)) ^ (~n + 2), ~x) + (((~d) ^ 2 * ~m * (~m - 1)) / ((~b) ^ 2 * (~f) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~c + ~d * ~x) ^ (~m - 2) * (~b * sin(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (LtQ(~n, -1) && (NeQ(~n, -2) && GtQ(~m, 1)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~c + ~d * ~x) ^ ~m, (~a + ~b * sin(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~n, 0) && (EqQ(~n, 1) || (IGtQ(~m, 0) || NeQ((~a) ^ 2 - (~b) ^ 2, 0))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (2 * ~a) ^ ~n * Antiderivative((~c + ~d * ~x) ^ ~m * sin((1 / 2) * (~e + (pi * ~a) / (2 * ~b)) + (~f * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~n) && (GtQ(~n, 0) || IGtQ(~m, 0))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((2 * ~a) ^ IntPart(~n) * (~a + ~b * sin(~e + ~f * ~x)) ^ FracPart(~n)) / sin(~e / 2 + (~a * pi) / (4 * ~b) + (~f * ~x) / 2) ^ (2 * FracPart(~n))) * Antiderivative((~c + ~d * ~x) ^ ~m * sin(~e / 2 + (~a * pi) / (4 * ~b) + (~f * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~n + 1 / 2) && (GtQ(~n, 0) || IGtQ(~m, 0))))
    #= Antiderivative(((~c')+(~d')*(~x))^(~m')*((~a)+(~b')*sin((~e')+(~f')*(~x)))^(~n'),~x) :=  (2*a)^n*Antiderivative((c+d*x)^m*cos(1/2*(e-pi*a/(2*b))+f*x/2)^(2*n),x) <--  FreeQ([a,b,c,d,e,f,m],x) && EqQ(a^2-b^2,0) && IntegerQ(n) &&  (GtQ(n,0) || IGtQ(m,0)) =#
    #= Antiderivative(((~c')+(~d')*(~x))^(~m')*((~a)+(~b')*sin((~e')+(~f')*(~x)))^(~n),~x) :=  (2*a)^IntPart(n)*(a+b*sin(e+f*x))^FracPart(n)/cos(1/2*(e-pi*a/(2*b))+ f*x/2)^(2*FracPart(n))* Antiderivative((c+d*x)^m*cos(1/2*(e-pi*a/(2*b))+f*x/2)^(2*n),x) <--  FreeQ([a,b,c,d,e,f,m],x) && EqQ(a^2-b^2,0) && IntegerQ(n+1/2) &&  (GtQ(n,0) || IGtQ(m,0)) =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * sin(~(e') + ~(k') * pi + ~(f') * Complex(0, ~fz) * ~x)), ~x) => 2 * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (-im * pi * (~k - 1 / 2)) * ℯ ^ (-im * ~e + ~f * ~fz * ~x)) / ((~b + 2 * ~a * ℯ ^ (-im * pi * (~k - 1 / 2)) * ℯ ^ (-im * ~e + ~f * ~fz * ~x)) - ~b * ℯ ^ (-2 * im * ~k * pi) * ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~fz], ~x) && (IntegerQ(2 * ~k) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * sin(~(e') + ~(k') * pi + ~(f') * ~x)), ~x) => 2 * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (im * pi * (~k - 1 / 2)) * ℯ ^ (im * (~e + ~f * ~x))) / ((~b + 2 * ~a * ℯ ^ (im * pi * (~k - 1 / 2)) * ℯ ^ (im * (~e + ~f * ~x))) - ~b * ℯ ^ (2 * im * ~k * pi) * ℯ ^ (2 * im * (~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IntegerQ(2 * ~k) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0)))
    #= Antiderivative(((~c')+(~d')*(~x))^(~m')/((~a)+(~b')*sin((~e')+(~f')*Complex(0,(~fz))*(~x))),(~x)*Symbol) :=  2*im*Antiderivative((c+d*x)^m*ℯ^(-im*e+f*fz*x)/(b+2*im*a*ℯ^(-im*e+f*fz*x)-b*ℯ^(2*(-im* e+f*fz*x))),x) <-- FreeQ([a,b,c,d,e,f,fz],x) && NeQ(a^2-b^2,0) &&  IGtQ(m,0) =#
    #= Antiderivative(((~c')+(~d')*(~x))^(~m')/((~a)+(~b')*sin((~e')+(~f')*(~x))),~x) :=  -2*im*Antiderivative((c+d*x)^m*ℯ^(im*(e+f*x))/(b-2*im*a*ℯ^(im*(e+f*x))-b*ℯ^(2*im*(e+f* x))),x) <-- FreeQ([a,b,c,d,e,f],x) && NeQ(a^2-b^2,0) && IGtQ(m,0) =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * sin(~(e') + ~(f') * Complex(0, ~fz) * ~x)), ~x) => 2 * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (-im * ~e + ~f * ~fz * ~x)) / (-im * ~b + 2 * ~a * ℯ ^ (-im * ~e + ~f * ~fz * ~x) + im * ~b * ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~fz], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => 2 * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (im * (~e + ~f * ~x))) / ((im * ~b + 2 * ~a * ℯ ^ (im * (~e + ~f * ~x))) - im * ~b * ℯ ^ (2 * im * (~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ 2, ~x) => ((~b * (~c + ~d * ~x) ^ ~m * cos(~e + ~f * ~x)) / (~f * ((~a) ^ 2 - (~b) ^ 2) * (~a + ~b * sin(~e + ~f * ~x))) + (~a / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative((~c + ~d * ~x) ^ ~m / (~a + ~b * sin(~e + ~f * ~x)), ~x)) - ((~b * ~d * ~m) / (~f * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative(((~c + ~d * ~x) ^ (~m - 1) * cos(~e + ~f * ~x)) / (~a + ~b * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((-(~b) * (~c + ~d * ~x) ^ ~m * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~n + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (~a / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative((~c + ~d * ~x) ^ ~m * (~a + ~b * sin(~e + ~f * ~x)) ^ (~n + 1), ~x) + ((~b * ~d * ~m) / (~f * (~n + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~n + 1), ~x)) - ((~b * (~n + 2)) / ((~n + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~c + ~d * ~x) ^ ~m * sin(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (ILtQ(~n, -2) && IGtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * (~a + ~b * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * sin(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * sin(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * cos(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * cos(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.10 (c+d x)^m (a+b sin)^n.jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libfyaml_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libfyaml")
JLLWrappers.@generate_main_file("libfyaml", UUID("c8bb72a1-3259-5c06-b5d5-7152be8a520b"))
end  # module libfyaml_jll


text:
########################################################################
function bpca_given_b(X, Σ, b; maxiter=1, S=1)
########################################################################

    #---------------------------------------------
    # Set constants, check dimensions and report
    #---------------------------------------------

    Q = 1

    JITTER = 1e-6

    D, N = size(X)

    @assert(all(size(X) .== size(Σ)))

    @assert(length(b) == D)

    numparam = D*Q + 1

    Σdiag = Diagonal(Σ)

    @printf("BPCA_given_b: %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function logPrior(W, μ)
    #---------------------------------------------

        -1e-10*(sum(W.^2) + sum(μ.^2))

    end

    #---------------------------------------------
    function marginalLogLikelihood(W, μ)
    #---------------------------------------------

        # local ℓ = zero(eltype(W))
        #
        # @inbounds for n in 1:N
        #
        #     local C = W*W' + Matrix(Diagonal(Σ[:,n])) + JITTER*I
        #
        #     ℓ +=  -0.5*logdet(C)-0.5*sum((X[:,n].-μ)'*(C\(X[:,n].-μ))) -0.5*D*log(2π)
        #
        # end

        # return ℓ

        sum(logpdf(MvNormal(μ, W*W' + Σdiag*I + JITTER*I), X))

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local MARK = 0

        local W = reshape(param[MARK+1:MARK+D*Q], D, Q)

        MARK += D*Q

        local s = exp(param[MARK+1])

        MARK += 1

        @assert(MARK == length(param))

        return W, s

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == numparam)

        local W, s = unpack(param)

        return marginalLogLikelihood(W, s*b) + logPrior(W, s*b)

    end

    #---------------------------------------------
    # Initialise with ppca estimates
    #---------------------------------------------

    W, s, _σ = ppca_given_b(X, b)

    #---------------------------------------------
    # Run optimiser
    #---------------------------------------------

    opt = Optim.Options(show_trace = true, show_every = 10, iterations = maxiter, g_tol=1e-6)

    mu, Sigma, neglogev = VI(objective, [[vec(W); log(s)] .+ s*0.001*randn(numparam) for s in 1:30], S = S, optimiser=LBFGS(), optimoptions = opt)

    #---------------------------------------------
    # Return posterior and log evidence
    #---------------------------------------------

    return MvNormal(mu, Sigma), -1*neglogev

end


text:
@testset "Memory" begin

@testset "Pointer-based" begin
    src = 42

    buf1 = Mem.alloc(sizeof(src); coherent=true)

    Mem.set!(buf1, UInt32(57), 1)
    x = Mem.download(UInt32, buf1)
    @test x[1] == UInt32(57)

    GC.@preserve Mem.upload!(buf1, pointer_from_objref(Ref(src)), sizeof(src))

    dst1 = Ref(0)
    GC.@preserve Mem.download!(pointer_from_objref(dst1), buf1, sizeof(src))
    @test src == dst1[]

    buf2 = Mem.alloc(sizeof(src))

    Mem.transfer!(buf2, buf1, sizeof(src))

    dst2 = Ref(0)
    GC.@preserve Mem.download!(pointer_from_objref(dst2), buf2, sizeof(src))
    @test src == dst2[]

    Mem.free(buf2)
    Mem.free(buf1)
end

@testset "Array-based" begin
    src = [42]

    buf1 = Mem.alloc(src)

    Mem.upload!(buf1, src)

    dst1 = similar(src)
    Mem.download!(dst1, buf1)
    @test src == dst1

    buf2 = Mem.upload(src)

    dst2 = similar(src)
    Mem.download!(dst2, buf2)
    @test src == dst2

    Mem.free(buf1)
end

@testset "Type-based" begin
    buf = Mem.alloc(Int)

    # there's no type-based upload, duh
    src = [42]
    Mem.upload!(buf, src)

    dst = Mem.download(eltype(src), buf)
    @test src == dst
end

@testset "Pointer information" begin
    default_agent = get_default_agent()

    N = 1024
    a = rand(N)
    b = Mem.alloc(default_agent, N)

    ptrinfo_host = Mem.pointerinfo(a)
    ptrinfo_hsa = Mem.pointerinfo(b)

    @test ptrinfo_host.type == HSA.EXT_POINTER_TYPE_UNKNOWN
    @test ptrinfo_hsa.type == HSA.EXT_POINTER_TYPE_HSA
    @test_skip ptrinfo_hsa.agentOwner.handle == default_agent.agent.handle

    Mem.free(b)
end

@testset "Page-locked memory (OS allocations)" begin
    a = rand(1024)
    plocked = Mem.lock(a)

    # NOTE - For a single agent, it seems that plocked == pointer(a)
    @test Mem.pointerinfo(pointer(a)).type == HSA.EXT_POINTER_TYPE_LOCKED
    @test Mem.pointerinfo(plocked).type == HSA.EXT_POINTER_TYPE_LOCKED
    @test Mem.pointerinfo(plocked).sizeInBytes == sizeof(a)

    Mem.unlock(a)
    @test Mem.pointerinfo(pointer(a)).type == HSA.EXT_POINTER_TYPE_UNKNOWN
    @test Mem.pointerinfo(plocked).type == HSA.EXT_POINTER_TYPE_UNKNOWN
end

@testset "Memory Region Queries" begin
    @testset "Region API Queries" begin
        for (idx, agent) in enumerate(AMDGPU.get_agents(:gpu))
            regions = AMDGPU.regions(agent)
            regions_global = filter(r->AMDGPU.region_segment(r) == AMDGPU.HSA.REGION_SEGMENT_GLOBAL, regions)
            regions_global_coarse_nohost = filter(r->(AMDGPU.region_global_flags(r) & AMDGPU.HSA.REGION_GLOBAL_FLAG_COARSE_GRAINED > 0) &&
                                                     !AMDGPU.region_host_accessible(r), regions)
            regions_group = filter(r->AMDGPU.region_segment(r) == AMDGPU.HSA.REGION_SEGMENT_GROUP, regions)
            regions_finegrained = filter(r->AMDGPU.region_global_flags(r) & AMDGPU.HSA.REGION_GLOBAL_FLAG_FINE_GRAINED > 0, regions)
            regions_kernarg = filter(r->AMDGPU.region_global_flags(r) & AMDGPU.HSA.REGION_GLOBAL_FLAG_KERNARG > 0, regions)

            @test length(regions_global) > 0
            if idx == 1
                @test length(regions_global_coarse_nohost) >= 1
                @test length(regions_group) == 1
            else
                # BUG: https://github.com/RadeonOpenCompute/ROCR-Runtime/issues/134
                @test length(regions_global_coarse_nohost) == 0
                @test length(regions_group) == 0
            end
            @test length(regions_finegrained) > 0
            @test length(regions_kernarg) > 0

            @test all(r->AMDGPU.region_size(r) > 0, regions)

            @test all(AMDGPU.region_runtime_alloc_allowed, regions_global)
            @test all(r->AMDGPU.region_runtime_alloc_granule(r) > 0, regions_global)
            @test all(r->AMDGPU.region_runtime_alloc_alignment(r) > 0, regions_global)
            @test all(r->AMDGPU.region_alloc_max_size(r) > 0, regions_global)

            @test !any(AMDGPU.region_runtime_alloc_allowed, regions_group)
            @test !any(AMDGPU.region_host_accessible, regions_group)

            @test all(AMDGPU.region_runtime_alloc_allowed, regions_finegrained)
            @test all(AMDGPU.region_host_accessible, regions_finegrained)

            @test all(AMDGPU.region_runtime_alloc_allowed, regions_kernarg)
            @test all(AMDGPU.region_host_accessible, regions_kernarg)

            for region in filter(AMDGPU.region_runtime_alloc_allowed, regions)
                buf = Mem.alloc(agent, region, 8)
                @test buf.ptr != C_NULL
                @test !buf.pool_alloc
                Mem.free(buf)
            end
        end
    end
    @testset "Memory Pool API Queries" begin
        for agent in AMDGPU.get_agents(:gpu)
            pools = AMDGPU.memory_pools(agent)
            pools_global = filter(r->AMDGPU.pool_segment(r) == AMDGPU.HSA.AMD_SEGMENT_GLOBAL, pools)
            pools_group = filter(r->AMDGPU.pool_segment(r) == AMDGPU.HSA.AMD_SEGMENT_GROUP, pools)

            @test length(pools_global) >= 1
            @test length(pools_group) == 1

            @test all(r->AMDGPU.pool_size(r) > 0, pools)
            @test !any(AMDGPU.pool_accessible_by_all, pools)

            @test all(AMDGPU.pool_runtime_alloc_allowed, pools_global)
            @test all(r->AMDGPU.pool_runtime_alloc_granule(r) > 0, pools_global)
            @test all(r->AMDGPU.pool_runtime_alloc_alignment(r) > 0, pools_global)
            @test all(r->AMDGPU.pool_alloc_max_size(r) > 0, pools_global)

            @test !any(AMDGPU.pool_runtime_alloc_allowed, pools_group)

            for pool in filter(AMDGPU.pool_runtime_alloc_allowed, pools)
                buf = Mem.alloc(agent, pool, 8)
                @test buf.ptr != C_NULL
                @test buf.pool_alloc
                Mem.free(buf)
            end
        end
    end
end

@testset "Exceptions" begin
    @test_throws ArgumentError Mem.alloc(Function, 1)   # abstract
    @test_throws ArgumentError Mem.alloc(Array{Int}, 1) # UnionAll
    @test_throws ArgumentError Mem.alloc(Integer, 1)    # abstract
    # TODO: can we test for the third case?
    #       !abstract && leaftype seems to imply UnionAll nowadays...

    # zero-width allocations should be permitted
    null = Mem.alloc(Int, 0)
    Mem.free(null)

    # double-free should throw
    x = Mem.alloc(1)
    Mem.free(x)
    # FIXME: Segfaults... @test_throws HSAError Mem.free(x)
end

@testset "Mutable structs" begin
    @eval mutable struct MutablePtrFree
        foo::Int
        bar::Int
    end
    buf = Mem.alloc(MutablePtrFree)
    Mem.upload!(buf, [MutablePtrFree(0,0)])
    Mem.free(buf)

    @eval mutable struct MutableNonPtrFree
        foo::Int
        bar::String
    end
    @test_throws ArgumentError Mem.alloc(MutableNonPtrFree)
end

end


text:
using HorizonSideRobots
#Основная функция -
function mark_frame_perimetr!(r::Robot)
    #Перемещение робота в Верхний-левый угол -
    num_vert = moves!(r, Nord)
    num_hor = moves!(r, West)

    #Маркировка периметра -
    for side in (Ost, Sud, West,Nord)
        putmarkers!(r, side)
    end

    #Перемещение робота в исходное положение -
    moves!(r, Sud, num_vert)
    moves!(r, Ost, num_hor)

end

function moves!(r::Robot, side:: HorizonSide)
    num_of_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_of_steps+=1
    end
    return num_of_steps
end

function moves!(r::Robot,side::HorizonSide,num_of_steps::Int)
    for _ in 1:num_of_steps
        move!(r,side)
    end
end

function ismarker!(r::Robot, side:: HorizonSide)
    while isborder(r,side)==false
        if ismarker(r,side)==false
            move!(r,side)
            putmarker!(r)
        end
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end


text:
function complain_if_not_supported(obj::AbstractParameterizable, arg0::JString)
    return jcall(obj, "complainIfNotSupported", void, (JString,), arg0)
end

function get_parameters_names(obj::AbstractParameterizable)
    return jcall(obj, "getParametersNames", List, ())
end

function is_supported(obj::AbstractParameterizable, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end


text:
# abstracttensor.jl
#
# Defines AbstractTensor, an abstract tensor to start a type hierarchy for
# representing tensors and extending/uniformizing Julia's built-in
# functionality for working with multilinear objects. A tensor is interpreted
# as a multilinear map whose indices are associated to vector spaces corresponding
# IndexSpace objects.

# Abstract Tensor type
#----------------------
abstract AbstractTensor{S<:IndexSpace,P<:TensorSpace,T,N}
# Any implementation of AbstractTensor should have method definitions for the
# same set of methods which are defined for the dense implementation Tensor
# defined in tensor.jl.

# tensor characteristics
spacetype{S}(::AbstractTensor{S})=S
spacetype{S<:IndexSpace}(::Type{AbstractTensor{S}})=S
spacetype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensor{S,P}})=S
spacetype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensor{S,P,T}})=S
spacetype{S<:IndexSpace,P<:TensorSpace,T,N}(::Type{AbstractTensor{S,P,T,N}})=S
spacetype{TT<:AbstractTensor}(::Type{TT})=spacetype(super(TT))

tensortype{S,P}(::AbstractTensor{S,P})=P
tensortype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensor{S,P}})=P
tensortype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensor{S,P,T}})=P
tensortype{S<:IndexSpace,P<:TensorSpace,T,N}(::Type{AbstractTensor{S,P,T,N}})=P
tensortype{TT<:AbstractTensor}(::Type{TT})=spacetype(super(TT))

Base.eltype{S,P,T}(::AbstractTensor{S,P,T}) = T
Base.eltype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensor{S,P,T}})=T
Base.eltype{S<:IndexSpace,P<:TensorSpace,T,N}(::Type{AbstractTensor{S,P,T,N}})=T
Base.eltype{TT<:AbstractTensor}(::Type{TT}) = eltype(super(TT))

TensorOperations.numind{S,P,T,N}(::AbstractTensor{S,P,T,N})=N
TensorOperations.numind{S,P<:TensorSpace,T,N}(::Type{AbstractTensor{S,P,T,N}})=N
TensorOperations.numind{TT<:AbstractTensor}(::Type{TT})=numind(super(TT))

dim(t::AbstractTensor)=dim(space(t))
sectors(t::AbstractTensor)=sectors(space(t))

order=numind

# check whether a tensor describes a valid state in a Hilbert space encoded as a QuantumSystem object
# by comparing space(t) to the tensor product structure of the Hilbert space
Base.in(t::AbstractTensor,V::VectorSpace)=issubspace(space(t),V)

tensor(t::AbstractTensor)=t
tensor(t::AbstractTensor,P::TensorSpace)= (space(t) == P ? t : throw(SpaceError("tensor not in $P")))

Base.promote_rule{S,P,T1,T2,N}(::Type{AbstractTensor{S,P,T1,N}},::Type{AbstractTensor{S,P,T2,N}})=AbstractTensor{S,P,promote_type(T1,T2),N}
Base.promote_rule{S,P,T1,T2,N1,N2}(::Type{AbstractTensor{S,P,T1,N1}},::Type{AbstractTensor{S,P,T2,N2}})=AbstractTensor{S,P,promote_type(T1,T2)}
Base.promote_rule{S,P,T1,T2}(::Type{AbstractTensor{S,P,T1}},::Type{AbstractTensor{S,P,T2}})=AbstractTensor{S,P,promote_type(T1,T2)}

Base.convert{S,P,T,N}(::Type{AbstractTensor{S,P,T,N}},t::AbstractTensor{S,P,T,N})=t
Base.convert{S,P,T1,T2,N}(::Type{AbstractTensor{S,P,T1,N}},t::AbstractTensor{S,P,T2,N})=copy!(similar(t,T1),t)
Base.convert{S,P,T}(::Type{AbstractTensor{S,P,T}},t::AbstractTensor{S,P,T})=t
Base.convert{S,P,T1,T2}(::Type{AbstractTensor{S,P,T1}},t::AbstractTensor{S,P,T2})=copy!(similar(t,T1),t)

# Basic algebra
#---------------
Base.copy(t::AbstractTensor)=Base.copy!(similar(t),t)

*(t::AbstractTensor,a::Number)=scale(t,a)
*(a::Number,t::AbstractTensor)=scale(t,a)
/(t::AbstractTensor,a::Number)=scale(t,one(a)/a)
\(a::Number,t::AbstractTensor)=scale(t,one(a)/a)
Base.scale(a::Number,t::AbstractTensor)=scale(t,a)
function Base.scale(t::AbstractTensor,a::Number)
    tnew=similar(t,promote_type(eltype(t),typeof(a)))
    scale!(tnew,t,a)
end
Base.scale!(t::AbstractTensor,a::Number)=scale!(t,t,a)
Base.scale!(a::Number,t::AbstractTensor)=scale!(t,a,t)

Base.conj(t::AbstractTensor)=Base.conj!(similar(t,conj(space(t))),t)
Base.transpose(t::AbstractTensor)=Base.transpose!(similar(t,space(t).'),t)
Base.ctranspose(t::AbstractTensor)=Base.ctranspose!(similar(t,space(t)'),t)

# Tensor operations
#-------------------
# convenience definition which works for vectors and matrices but also sometimes useful in general case
*{S,P,T1,T2,N1,N2}(t1::AbstractTensor{S,P,T1,N1},t2::AbstractTensor{S,P,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[1:N1-1] ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(1:N1-1,0),'N',t2,vcat(0,numind(t1)-1+(1:N2-1)),'N',0,t3,1:(N1+N2-2)))
Base.At_mul_B{S,P,T1,T2,N1,N2}(t1::AbstractTensor{S,P,T1,N1},t2::AbstractTensor{S,P,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[2:N1].' ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(0,reverse(1:N1-1)),'N',t2,vcat(0,N1-1+(1:N2-1)),'N',0,t3,1:(numind(t1)+numind(t2)-2)))
Base.Ac_mul_B{S,P,T1,T2,N1,N2}(t1::AbstractTensor{S,P,T1,N1},t2::AbstractTensor{S,P,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[2:N1]' ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(0,reverse(1:N1-1)),'C',t2,vcat(0,N1-1+(1:N2-1)),'N',0,t3,1:(numind(t1)+numind(t2)-2)))

⊗{S,P}(t1::AbstractTensor{S,P},t2::AbstractTensor{S,P})=tensorproduct(t1,1:numind(t1),t2,numind(t1)+(1:numind(t2)))
Base.trace{S,P,T}(t::AbstractTensor{S,P,T,2})=scalar(tensortrace(t,[1,1],[]))

# general tensor operations: no error checking, pass to mutating methods
function tensorcopy(A::AbstractTensor,labelsA,outputlabels=labelsA)
    spaceA=space(A)
    spaceC=spaceA[indexin(outputlabels,labelsA)]
    C=similar(A,spaceC)
    tensorcopy!(A,labelsA,C,outputlabels)
    return C
end
function tensoradd{S,P,TA,TB,N}(A::AbstractTensor{S,P,TA,N},labelsA,B::AbstractTensor{S,P,TB,N},labelsB,outputlabels=labelsA)
    spaceA=space(A)
    spaceC=spaceA[indexin(outputlabels,labelsA)]
    T=promote_type(TA,TB)
    C=similar(A,T,spaceC)
    tensorcopy!(A,labelsA,C,outputlabels)
    tensoradd!(1,B,labelsB,1,C,outputlabels)
    return C
end
function tensortrace(A::AbstractTensor,labelsA,outputlabels)
    T=eltype(A)
    spaceA=space(A)
    spaceC=spaceA[indexin(outputlabels,labelsA)]
    C=similar(A,spaceC)
    tensortrace!(1,A,labelsA,0,C,outputlabels)
    return C
end
function tensortrace(A::AbstractTensor,labelsA) # there is no one-line method to compute the default outputlabels
    ulabelsA=unique(labelsA)
    labelsC=similar(labelsA,0)
    sizehint!(labelsC,length(ulabelsA))
    for j=1:length(ulabelsA)
        ind=findfirst(labelsA,ulabelsA[j])
        if findnext(labelsA,ulabelsA[j],ind+1)==0
            push!(labelsC,ulabelsA[j])
        end
    end
    return tensortrace(A,labelsA,labelsC)
end
function tensorcontract{S}(A::AbstractTensor{S},labelsA,conjA,B::AbstractTensor{S},labelsB,conjB,outputlabels=symdiff(labelsA,labelsB);method::Symbol=:BLAS)
    spaceA=space(A)
    spaceB=space(B)
    spaceC=(spaceA ⊗ spaceB)[indexin(outputlabels,vcat(labelsA,labelsB))]
    T=promote_type(eltype(A),eltype(B))
    C=similar(A,T,spaceC)
    tensorcontract!(1,A,labelsA,conjA,B,labelsB,conjB,0,C,outputlabels;method=method)
    return C
end
tensorcontract{S}(A::AbstractTensor{S},labelsA,B::AbstractTensor{S},labelsB,outputlabels=symdiff(labelsA,labelsB);
    method::Symbol=:BLAS)=tensorcontract(A,labelsA,'N',B,labelsB,'N',outputlabels;method=method)

function tensorproduct{S}(A::AbstractTensor{S},labelsA,B::AbstractTensor{S},labelsB,outputlabels=vcat(labelsA,labelsB))
    spaceA=space(A)
    spaceB=space(B)

    spaceC=(spaceA ⊗ spaceB)[indexin(outputlabels,vcat(labelsA,labelsB))]
    T=promote_type(eltype(A),eltype(B))
    C=similar(A,T,spaceC)
    tensorproduct!(1,A,labelsA,B,labelsB,0,C,outputlabels)
    return C
end

Base.eye{S<:ElementarySpace,P<:TensorSpace,T}(t::AbstractTensor{S,P,T},V::S)=eye(T,P,V)

# Factorization
#---------------
Base.svd{S,P,T}(t::AbstractTensor{S,P,T,2})=svd(t,1,2)
leftorth{S,P,T}(t::AbstractTensor{S,P,T,2})=leftorth(t,1,2)
rightorth{S,P,T}(t::AbstractTensor{S,P,T,2})=rightorth(t,1,2)

# general tensor factorizations: permute to correct order and pass to in place methods
Base.svd(t::AbstractTensor,leftind,truncation::TruncationScheme=notrunc())=svd(t,leftind,setdiff(1:numind(t),leftind),truncation)
function Base.svd(t::AbstractTensor,leftind,rightind,truncation::TruncationScheme=notrunc())
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into leftind and rightind.
    N=numind(t)
    p=vcat(leftind,rightind)
    (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
    newt=tensorcopy(t,1:N,p)
    return svd!(newt,length(leftind),truncation)
end

leftorth(t::AbstractTensor,leftind,truncation::TruncationScheme=notrunc())=leftorth(t,leftind,setdiff(1:numind(t),leftind),truncation)
function leftorth(t::AbstractTensor,leftind,rightind=setdiff(1:numind(t),leftind),truncation::TruncationScheme=notrunc())
    # Create orthogonal basis U for left indices, and remainder R for right
    # indices. Decomposition should be unique, such that it always returns the
    # same result for the same input tensor t. QR is fastest but only unique
    # after correcting for phases.
    N=numind(t)
    p=vcat(leftind,rightind)
    (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
    newt=tensorcopy(t,1:N,p)
    return leftorth!(newt,length(leftind),truncation)
end

rightorth(t::AbstractTensor,leftind,truncation::TruncationScheme=notrunc())=rightorth(t,leftind,setdiff(1:numind(t),leftind),truncation)
function rightorth(t::AbstractTensor,leftind,rightind=setdiff(1:numind(t),leftind),truncation::TruncationScheme=notrunc())
    # Create orthogonal basis U for left indices, and remainder R for right
    # indices. Decomposition should be unique, such that it always returns the
    # same result for the same input tensor t. QR is fastest but only unique
    # after correcting for phases.
    N=numind(t)
    p=vcat(leftind,rightind)
    (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
    newt=tensorcopy(t,1:N,p)
    return rightorth!(newt,length(leftind),truncation)
end

svd!(t::AbstractTensor,n::Int) = svd!(t,n,notrunc())
leftorth!(t::AbstractTensor,n::Int) = leftorth!(t,n,notrunc())
rightorth!(t::AbstractTensor,n::Int) = rightorth!(t,n,notrunc())


text:
using BenchmarkTools
using Libtask

macro rep(cnt, exp)
    blk =:(begin end)
    for _ in 1:eval(cnt)
        push!(blk.args, esc(exp))
    end
    blk
end

INTENSITY = 6

indexing(a, x, y) = @rep INTENSITY a[x, y]
setindexing(a, x, y) = @rep INTENSITY a[x, y] = 1
broadcasting(a) = @rep INTENSITY a .+ a

println("= Benchmarks on Arrays =")
A = rand(1000, 1000)
x, y =  abs.(rand(Int, 2) .% 999) .+ 1
print("indexing: ")
@btime indexing($A, $x, $y)
print("set indexing: ")
@btime setindexing($A, $x, $y)
print("broadcast: ")
@btime broadcasting($A)

println("= Benchmarks on TArrays =")
TA = Libtask.localize(deepcopy(A))
print("indexing: ")
@btime indexing($TA, $x, $y)
print("set indexing: ")
@btime setindexing($TA, $x, $y)
print("broadcast: ")
@btime broadcasting($TA)

print("= Benchmarks on task and TArray copying")
function f()
    t = TArray(Int, 1)
    t[1] = 0
    while true
        produce(t[1])
        t[1]
        t[1] = 1 + t[1]
    end
end

@btime begin
    ctask = CTask(f)
    consume(ctask)
    consume(ctask)
    a = copy(ctask)
    consume(a)
    consume(a)
end


text:
"""
## module Potentials

### Summary

This module implements some basic interatomic potentials in pure Julia, as well
as provides building blocks and prototypes for further implementations
The implementation is done in such a way that they can be used either in "raw"
form or withinabstractframeworks.

### Types

### `evaluate`, `evaluate_d`, `evaluate_dd`, `grad`

### The `@D`, `@DD`, `@GRAD` macros

TODO: write documentation

"""
module Potentials

using JuLIP: AbstractAtoms, AbstractNeighbourList, AbstractCalculator,
      JVec, JVecs, mat, vec, JMat, JVecF, SVec, vecs, SMat,
      positions, set_positions!
using StaticArrays: @SMatrix

using NeighbourLists

import JuLIP: energy, forces, cutoff, virial, hessian_pos, hessian, site_energies

export Potential, PairPotential, SitePotential,
     site_energy, site_energy_d, partial_energy, partial_energy_d

"""
`Potential`: genericabstractsupertype for all potential-like things
"""
abstract type Potential <: AbstractCalculator end

"""
`PairPotential`:abstractsupertype for pair potentials

Can also be used as a constructor for analytic pair potentials, e.g.,
```julia
lj = @analytic r -> r^(-12) - 2 * r^(-6)
```
"""
abstract type PairPotential <: Potential end

"""
`SitePotential`:abstractsupertype for generic site potentials
"""
abstract type SitePotential <: Potential end

include("potentials_base.jl")
# * @pot, @D, @DD, @GRAD and related things
# TODO: move this to a separate package???

NeighbourLists.sites(at::AbstractAtoms, rcut::AbstractFloat) =
      sites(neighbourlist(at, rcut))

NeighbourLists.pairs(at::AbstractAtoms, rcut::AbstractFloat) =
      pairs(neighbourlist(at, rcut))

# Implementation of a generic site potential
# ================================================

site_energies(V::SitePotential, at::AbstractAtoms) =
   Float64[ V(r, R) for (_₁, _₂, r, R) in sites(at, cutoff(V)) ]

energy(V::SitePotential, at::AbstractAtoms) = sum_kbn(site_energies(V, at))

evaluate(V::SitePotential, R::AbstractVector{JVecF}) = evaluate(V, norm.(R), R)

evaluate_d(V::SitePotential, R::AbstractVector{JVecF}) = evaluate_d(V, norm.(R), R)

function forces(V::SitePotential, at::AbstractAtoms)
   frc = zerovecs(length(at))
   for (i, j, r, R) in sites(at, cutoff(V))
      dV = @D V(r, R)
      for a = 1:length(j)
         frc[j[a]] -= dV[a]
      end
      frc[i] += sum(dV)
   end
   return frc
end

site_virial(dV, R) =  (
      length(R) > 0 ? (- sum( dVi * Ri' for (dVi, Ri) in zip(dV, R) ))
                    : zero(JMatF) )

virial(V::SitePotential, at::AbstractAtoms) =
      sum(  site_virial((@D V(r, R)), R)
            for (_₁, _₂, r, R) in sites(at, cutoff(V))  )

# TODO: partial_energy and partial_energy_d are not tested properly

function partial_energy(V::SitePotential, at::AbstractAtoms, Idom)
   E = 0.0
   nlist = neighbourlist(at, cutoff(V))
   for i in Idom
      j, r, R = site(nlist, i)
      E += V(r, R)
   end
   return E
end

function partial_energy_d(V::SitePotential, at::AbstractAtoms, Idom)
   F = zeros(JVecF, length(at))
   nlist = neighbourlist(at, cutoff(V))
   for i in Idom
      j, r, R = site(nlist, i)
      dV = @D V(R)
      F[j] += dV
      F[i] -= sum(dV)
   end
   return F
end

partial_energy(V::PairPotential, at::AbstractAtoms, Idom) =
      partial_energy(PairSitePotential(V), at, Idom)
partial_energy_d(V::PairPotential, at::AbstractAtoms, Idom) =
      partial_energy_d(PairSitePotential(V), at, Idom)

site_energy(V::Union{SitePotential, PairPotential}, at::AbstractAtoms, i0::Int) =
      partial_energy(V, at, [i0])
site_energy_d(V::Union{SitePotential, PairPotential}, at::AbstractAtoms, i0::Int) =
      partial_energy_d(V, at, [i0])

include("analyticpotential.jl")
# * AnalyticFunction
# * WrappedPPotential

include("cutoffs.jl")
#   * SWCutoff
#   * ShiftCutoff
#   * SplineCutoff

include("pairpotentials.jl")
# * PairCalculator
# * LennardJonesPotential
# * MorsePotential
# * SimpleExponential

include("adsite.jl")
# * ADPotential : Site potential using ForwardDiff

include("stillingerweber.jl")
# * type StillingerWeber

include("splines.jl")
include("eam.jl")
# EAM, FinnisSinclair

export ZeroSitePotential

@pot type ZeroSitePotential <: SitePotential
end

"a site potential that just returns zero"
ZeroSitePotential

evaluate(p::ZeroSitePotential, r, R) = 0.0
evaluate_d(p::ZeroSitePotential, r, R) = zeros(r)   # TODO: is this a bug?
cutoff(::ZeroSitePotential) = 3.0

"""
`fd_hessian(V, R, h) -> H`

If `length(R) = N` and `length(R[i]) = d` then `H` is an N × N `Matrix{SMatrix}` with
each element a d × d static array.
"""
function fd_hessian{D,T}(V::SitePotential, R::Vector{SVec{D,T}}, h)
   d = length(R[1])
   N = length(R)
   H = zeros( typeof(@SMatrix zeros(d, d)), N, N )
   return fd_hessian!(H, V, R, h)
end

"""
`fd_hessian!(H, V, R, h) -> H`

Fill `H` with the hessian entries; cf `fd_hessian`.
"""
function fd_hessian!{D,T}(H, V::SitePotential, R::Vector{SVec{D,T}}, h)
   N = length(R)
   # convert R into a long vector and H into a big matrix (same part of memory!)
   Rvec = mat(R)[:]
   Hmat = zeros(N*D, N*D)   # reinterpret(T, H, (N*D, N*D))
   # now re-define ∇V as a function of a long vector (rather than a vector of SVecs)
   dV(r) = (evaluate_d(V, r |> vecs) |> mat)[:]
   # compute the hessian as a big matrix
   for i = 1:N*D
      Rvec[i] +=h
      dVp = dV(Rvec)
      Rvec[i] -= 2*h
      dVm = dV(Rvec)
      Hmat[:, i] = (dVp - dVm) / (2 * h)
      Rvec[i] += h
   end
   Hmat = 0.5 * (Hmat + Hmat')
   # convert to a block-matrix
   for i = 1:N, j = 1:N
      Ii = (i-1) * D + (1:D)
      Ij = (j-1) * D + (1:D)
      H[i, j] = SMat{D,D}(Hmat[Ii, Ij])
   end
   return H
end

function fd_hessian(calc::AbstractCalculator, at::AbstractAtoms, h)
   d = 3
   N = length(at)
   H = zeros( typeof(@SMatrix zeros(d, d)), N, N )
   return fd_hessian!(H, calc, at, h)
end

"""
`fd_hessian!{D,T}(H, calc, at, h) -> H`

Fill `H` with the hessian entries; cf `fd_hessian`.
"""
function fd_hessian!(H, calc::AbstractCalculator, at::AbstractAtoms, h)
   D = 3
   N = length(at)
   X = positions(at) |> mat
   x = X[:]
   # convert R into a long vector and H into a big matrix (same part of memory!)
   Hmat = zeros(N*D, N*D)
   # now re-define ∇V as a function of a long vector (rather than a vector of SVecs)
   dE(x_) = (site_energy_d(calc, set_positions!(at, reshape(x_, D, N)), 1) |> mat)[:]
   # compute the hessian as a big matrix
   for i = 1:N*D
      x[i] += h
      dEp = dE(x)
      x[i] -= 2*h
      dEm = dE(x)
      Hmat[:, i] = (dEp - dEm) / (2 * h)
      x[i] += h
   end
   Hmat = 0.5 * (Hmat + Hmat')
   # convert to a block-matrix
   for i = 1:N, j = 1:N
      Ii = (i-1) * D + (1:D)
      Ij = (j-1) * D + (1:D)
      H[i, j] = SMat{D,D}(Hmat[Ii, Ij])
   end
   return H
end

hessian_pos(V::SitePotential, at::AbstractAtoms) =
      _precon_or_hessian_pos(V, at, hess)

# implementation of a generic assembly of a global block-hessian from
# local site-hessians
function _precon_or_hessian_pos(V::SitePotential, at::AbstractAtoms, hfun)
   nlist = neighbourlist(at, cutoff(V))
   I, J, Z = Int[], Int[], JMatF[]
   # a crude size hint
   for C in (I, J, Z); sizehint!(C, 24*npairs(nlist)); end
   for (i, neigs, r, R) in sites(nlist)
      nneigs = length(neigs)
      # [1] the "off-centre" component of the hessian:
      # h[a, b] = ∂_{Ra} ∂_{Rb} V     (this is a nneigs x nneigs block-matrix)
      h = hfun(V, r, R)
      for a = 1:nneigs, b = 1:nneigs
         push!(I, neigs[a])
         push!(J, neigs[b])
         push!(Z, h[a,b])
      end

      # [2] the ∂_{Ri} ∂_{Ra} terms
      # hib = ∂_{Ri} ∂_{Rb} V = - ∑_a ∂_{Ra} ∂_{Rb} V
      # also at the same time we pre-compute the centre-centre term:
      #    hii = ∂_{Ri} ∂_{Ri} V = - ∑_a ∂_{Ri} ∂_{Ra} V
      hii = zero(JMatF)
      for b = 1:nneigs
         hib = -sum( h[a, b] for a = 1:nneigs )
         hii -= hib
         append!(I, (i,         neigs[b] ))
         append!(J, (neigs[b],  i        ))
         append!(Z, (hib,       hib'     ))
      end

      # and finally add the  ∂_{Ri}^2 term, which is precomputed above
      push!(I, i)
      push!(J, i)
      push!(Z, hii)
   end
   return sparse(I, J, Z, length(at), length(at))
end

end


text:
# using Plots
# using PyCall
# using TimerOutputs
# using Statistics
#
# # const to = TimerOutput();
#
# # ht = readHITRAN("/net/fluo/data1/projects/RTM/Mie/hitran/par/02_hit04.par", 2, 1, 6000, 6400)

# using RadiativeTransfer
# include("/home/rjeyaram/RadiativeTransfer/src/RadiativeTransfer.jl")
# using Revise
# using RadiativeTransfer.CrossSection

# # using Test

# # ht = CrossSection.read_hitran("/home/rjeyaram/RadiativeTransfer/test/CO2.data", 2, 1, 6000, 6400)
# ht = CrossSection.read_hitran("/net/fluo/data1/projects/RTM/Mie/hitran/par/02_hit04.par", mol=2, iso=1, ν_min=6000, ν_max=6400)
# modCEF = ErfcHumliErrorFunctionVoigt();
# grid = collect(6000:0.01:6400);

# # @time cs = line_shape(voigt,modCEF,ht,grid,false,1013.0,296,10,0)

# # display(plot(grid, cs, ylims=(0, 3e-17)))

# # @time cs = line_shape(voigt,modCEF,ht,grid,false,1013.15,293,40,0)

# # @time plot(grid, line_shape(doppler,modCEF,ht,grid,false,1013.15,293,40,0), ylims=(0,1e-22))
# # @time plot!(grid, line_shape(lorentz,modCEF,ht,grid,false,1013.15,293,40,0), ylims=(0,1e-22))
# # @time display(plot!(grid, line_shape(voigt,modCEF,ht,grid,false,1013.15,293,40,0), ylims=(0,1e-22)))

# # @time doppler_shift = line_shape(doppler,modCEF,ht,grid,false,1013.15,293,40,0)
# # @time lorentz_shift = line_shape(lorentz,modCEF,ht,grid,false,1013.15,293,40,0)

# # voigt_times = []
# #
# # append!(voigt_times, 1)
# #
# # for i in 1:100
# #     # append!(voigt_times, 1)
# #     time = @elapsed res = line_shape(voigt,modCEF,ht,grid,false,1013.15,293,40,0)
# #     append!(voigt_times, time)
# #     println(time)
# # end
# #
# # println("Mean and stddev")
# # println(mean(voigt_times))
# # println(std(voigt_times))
# #
# # # time = @elapsed res = line_shape(voigt,modCEF,ht,grid,false,1013.15,293,40,0)
# #
# using DelimitedFiles

# CEF = ErfcHumliErrorFunctionVoigt()
# @time cs = absorption_cross_section(Voigt(CEF),ht,grid,false,1013.25,296.0,40,vmr=0)

# using Plots

# maximum(cs)
# plot(grid, cs, ylims=(0, 8e-23))
# plot!(grid, cs, ylims=(0, 8e-23))

# plot

# using Interpolations

# T = 200:10:380;
# p = 1:25:1050;

# cs_matrix = zeros(length(T), length(p), length(grid));

# for iP in eachindex(p)
#     println("iP: ", iP)
#     for iT in eachindex(T)
#         println("iT: ", iT)
#        cs_matrix[iT,iP,:] = line_shape(voigt,modCEF,ht,grid,false,p[iP],T[iT],10.0,0.0);
#     end
#     # println(iP)
# end
# itp = interpolate( cs_matrix, BSpline(Cubic(Line(OnGrid()))))
# sitp = scale(itp,T,p,6000:0.01:6400)

# plot(grid, )

# using Interact

# writedlm( "/home/rjeyaram/RadiativeTransfer/test/voigt_test_julia.csv",  cs, ',')

# using Plots

# b = @manipulate for curr_T = 200:10:380, curr_P = 1:25:1050

#     display(plot(grid, sitp(curr_T, curr_P, grid), ylims=(0, 8e-23)))

# end

# using JLD2
# @save "/home/rjeyaram/RadiativeTransfer/test/savedInterpolation.jld" itp

# function app(req)

#     grid= -10:1:10

#     a = @manipulate for m = 0.1:0.1:2.0#curr_T = 200:10:380, curr_P = 1:25:1050
#         plot(grid, m .* grid)
#     end

#     return a

# end

# webio_serve(page("/", app), port=5055)

# using Mux

# ui = button()
# WebIO.webio_serve(page("/", req -> b), 5055)

# @testset "RadiativeTransfer.jl" begin
#     # Write your tests here.

# end

# using Revise
# using RadiativeTransfer.CrossSection
# using DelimitedFiles
# using Plots

# # pyplot()

# test_ht = CrossSection.read_hitran("helper/CO2.data", ν_min=6000, ν_max=6400)

# temperatures = [100, 175, 250, 325, 400]
# pressures = [250, 500, 750, 1000, 1250]

# grid = collect(6000:0.01:6400);
# CEF = ErfcHumliErrorFunctionVoigt()

# # Threshold -- our value must be within ϵ of the HAPI value
# ϵ = 3.5e-27

# for temp in temperatures
#     println(temp)
#     for pres in pressures
#         jl_cs = absorption_cross_section(Voigt(CEF),test_ht,grid,false,pres,temp,40,vmr=0)
#         py_cs = readdlm("/home/rjeyaram/RadiativeTransfer/test/helper/Voigt_CO2_T" * string(temp) * "_P" * string(pres) * ".csv")
#         Δcs = abs.(jl_cs - py_cs)

#         display(plot(grid, Δcs, ylims=(0, ϵ)))

#         plot(grid, py_cs, ylims=(0, 1.3*maximum(py_cs)))
#         display(plot!(grid, jl_cs, ylims=(0, 1.3*maximum(py_cs))))

#         # println("Max: ", )
#         # println("MAE: ", mean(Δcs))
#         @test maximum(Δcs) < ϵ
#     end

# end

# using Revise
# using Profile
# using RadiativeTransfer
# using RadiativeTransfer.CrossSection

# using BenchmarkTools

# test_ht = read_hitran("/home/rjeyaram/RadiativeTransfer/test/helper/CO2.data", ν_min=6000, ν_max=6400)

# model4 = make_hitran_model(test_ht, Voigt(), architecture = Architectures.CPU())
# model4_GPU = make_hitran_model(test_ht, Voigt(), architecture = Architectures.GPU())

# ν_grid = 6000:0.01:6400
# pressures = 250:250:1250
# temperatures = 100:75:400

# model = make_interpolation_model(test_ht, Voigt(), ν_grid, pressures, temperatures, wing_cutoff = 40, CEF=ErfcHumliErrorFunctionVoigt())

# matrix = Array(model.itp)

# res = absorption_cross_section(model, 6000:0.01:6400, 1000.1, 296.1)

# f(x) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);
# x = rand(4);

# result = DiffResults.HessianResult(x)
# result = ForwardDiff.hessian!(result, f, x);

# const ν_grid = 6000:0.01:6400

# absorption_cross_section(model4, ν_grid, 1000.1, 296.1)
# @btime absorption_cross_section(model4_GPU, )

# @time absorption_cross_section(model4_GPU, ν_grid, 1000.1, 296.1)

# @code_warntype testReturn(3)

# x = [1000.1, 296.1]
# result = DiffResults.HessianResult(x)

# result = ForwardDiff.hessian!(result, CrossSection.acs_shorthand, x);

# x = [0.3, 6.82, 1.3, 0.00001]
# result = DiffResults.JacobianResult(zeros(4568), x)
# ForwardDiff.jacobian!(result, Scattering.phase_shorthand, x);

# using Revise
# using Plots
# using RadiativeTransfer
# using RadiativeTransfer.vSmartMOM
# using CUDA
# using ForwardDiff
# n = 20;
# nSpec = 100;
# RM = rand(n,n,nSpec)
# x = [1.0,2.0]
# function testInv_CPU(x, RM=RM)
#     temp = x[1]* RM .+ x[2]
#     B = similar(temp);
#     vSmartMOM.batch_inv!(B,temp)
#     @show typeof(B)
#     return B
# end
# function testInv_GPU(x, RM=RM)
#     temp = CuArray(x[1]* RM .+ x[2])
#     B = CuArray(similar(temp));
#     vSmartMOM.batch_inv!(B,temp)
#     synchronize()
#     @show typeof(B)
#     return B
# end

# dfdx_CPU = ForwardDiff.jacobian(testInv_CPU, x)
# dfdx_GPU = ForwardDiff.jacobian(testInv_GPU, x)

using RadiativeTransfer, RadiativeTransfer.Absorption, Interpolations, Statistics

function apply_instrument( input_grid::AbstractRange,
                               input_spectra,
                           output_grid::Array{<:Real,1})
    """
down-scales the co-domain of spectral cross-sections grid to the instrument grid
"""

    itp = interpolate(input_spectra, BSpline(Cubic(Line(OnGrid()))))
    sitp = scale(itp, input_grid)
    output_spectra = sitp(output_grid)
       return output_spectra
#    return itp(output_spectral_grid)
end # end of function apply_instrument

ν_min, δν, ν_max = 6000, 0.01, 6400
p, T = 500, 273
molec_num, iso_num = 1, 1
filename = "../H2O_S.data"

hitran_table = CrossSection.read_hitran(artifact("CO2"), mol=molec_num, iso=iso_num, ν_min=ν_min, ν_max=ν_max)
model = make_hitran_model(hitran_table, Voigt(), architecture=GPU());
grid = ν_min:δν:ν_max;
cross_sections = absorption_cross_section(model, grid, p, T)

# define the output grid
out_grid = collect(ν_min:0.1:ν_max)
#in_grid = collect(grid)
in_grid = grid;
out_spectra = apply_instrument(in_grid, cross_sections, out_grid)


text:
using Documenter, TensorBoardLogger

makedocs(
    modules   = [TensorBoardLogger],
    format    = Documenter.HTML(),
    sitename  = "TensorBoardLogger.jl",
    authors   = "Filippo Vicentini",
    pages     = [
            "Home"          => "index.md",
            "Manual"        => Any[
                "Backends"     =>  "custom_behaviour.md",
                "Extending"    =>  "extending_behaviour.md",
                "Explicit Interface" => "explicit_interface.md"
            ],
            "Examples"      => Any[
                "Optim.jl"  => "examples/optim.md"
            ]
    ]
)

deploydocs(
    repo   = "github.com/PhilipVinc/TensorBoardLogger.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
# --- Parsing ---

@doc raw"""
    parse(s::AbstractString)
    parse(io::IO)

Parse a literal string.

```jldoctest
julia> HTM.parse("pineapple: <div class=\\"fruit\\">🍍</div>...")
3-element Vector{Union{String, HTM.Node}}:
 "pineapple: "
 HTM.Node(["div"], [:class => ["fruit"]], Union{String, HTM.Node}["🍍"])
 "..."
```
"""
@inline parse(io::IO) = parseelems(io)
@inline parse(s::AbstractString) = parse(IOBuffer(s))

@doc raw"""
    parseelems(io::IO)

Parse elements.

This function is the entry point for an implementation of a subset of the
[HTML standard](https://html.spec.whatwg.org/multipage/parsing.html#tokenization).

```jldoctest
julia> HTM.parseelems(IOBuffer("pineapple: <div class=\\"fruit\\">🍍</div>..."))
3-element Vector{Union{String, HTM.Node}}:
 "pineapple: "
 HTM.Node(["div"], [:class => ["fruit"]], Union{String, HTM.Node}["🍍"])
 "..."
```
"""
@inline parseelems(io::IO) = parseelems(io -> true, io)
@inline parseelems(predicate, io::IO) = parseelems!(predicate, io, Union{String,HTM.Node}[])
@inline function parseelems!(predicate, io::IO, elems::AbstractVector)
    while !eof(io) && predicate(io)
        skipcomment(io) || pushelem!(elems, parseelem(io))
    end
    return elems
end
@inline pushelem!(elems::AbstractVector, elem) = push!(elems, elem)
@inline pushelem!(elems::AbstractVector, elem::AbstractString) = (isempty(elem) || all(isspace, elem)) ? elems : push!(elems, elem)  # TODO: detect all spaces during reading for performance

@doc raw"""
    parseelem(io::IO)

Parse a single element.

```jldoctest
julia> HTM.parseelem(IOBuffer("pineapple: <div class=\\"fruit\\">🍍</div>..."))
"pineapple: "
```
"""
@inline function parseelem(io::IO)
    startswith(io, '<') && return parsenode(io)
    skipstartswith(io, "\\\$") && return "\$"  # frustrated interp
    return parseinterp(∈("<\$\\"), io)
end

@doc raw"""
    skipcomment(io::IO)

Skip a comment if any.

```jldoctest
julia> io = IOBuffer("<!-- 🍌 -->🍍");

julia> HTM.skipcomment(io)
true

julia> read(io, Char)
'🍍': Unicode U+1F34D (category So: Symbol, other)
```
"""
@inline function skipcomment(io::IO)
    if skipstartswith(io, "<!--")
        while !(eof(io) || skipstartswith(io, "-->"))
            skip(io, 1)
        end
        return true
    end
    return false
end

@doc raw"""
    parsenode(io::IO)

Parse an [`HTM.Node`](@ref) object.

```jldoctest
julia> HTM.parsenode(IOBuffer("<div class=\\"fruit\\">🍍</div>..."))
HTM.Node(["div"], [:class => ["fruit"]], Union{String, HTM.Node}["🍍"])
```
"""
@inline function parsenode(io::IO)
    skipchars(isequal('<'), io)
    tag = parsetag(io)
    attrs = parseattrs(io)
    read(io, Char) === '/' && (skipchars(isequal('>'), io); return Node(tag, attrs))
    endtag = *("</", create_tag(tag), '>')
    children = parseelems(io -> !(startswith(io, endtag) || startswith(io, UENDTAG)), io)
    skipstartswith(io, endtag) || skipstartswith(io, UENDTAG) || error("tag not properly closed")
    return Node(tag, attrs, children)
end

@doc raw"""
    parsetag(io::IO)

Parse a tag.

```jldoctest
julia> HTM.parsetag(IOBuffer("div class=\\"fruit\\">🍍..."))
1-element Vector{String}:
 "div"
```
"""
@inline function parsetag(io::IO)
    🧩 = String[]
    while !(eof(io) || (🍒 = peek(io, Char)) |> isspace || 🍒 ∈ ">/")
        push!(🧩, skipstartswith(io, "\\\$") ? "\$" : parseinterp(isspace ⩔ ∈(">/\$\\"), io))
    end
    return 🧩
end

@doc raw"""
    parseattrs(io::IO)

Parse attributes of a node.

```jldoctest
julia> HTM.parseattrs(IOBuffer("class=\\"fruit\\">🍍..."))
1-element Vector{Pair{Symbol, Vector{String}}}:
 :class => ["fruit"]

julia> HTM.parseattrs(IOBuffer("class=\\"fruit\\" \$(attrs)>🍍..."))
2-element Vector{Pair{Symbol, Vector{String}}}:
     :class => ["fruit"]
 Symbol("") => ["\$(attrs)"]
```
"""
@inline parseattrs(io::IO) = parseattrs!(io, Pair{Symbol,Vector{String}}[])
@inline function parseattrs!(io::IO, attrs::AbstractVector)
    while !eof(io)
        skipchars(isspace, io)
        startswith(io, ('>', '/')) && break
        parseattr!(io, attrs)
    end
    return attrs
end
@inline function parseattr!(io::IO, attrs::AbstractVector)
    startswith(io, '$') && return push!(attrs, Symbol() => [parseinterp(io)])  # spread attributes
    startswith(io, "\\\$") && skip(io, 1)  # no interps in keys: just ignore escaping
    🔑 = parsekey(io)
    let 🍒 = read(io, Char)
        push!(attrs, 🔑 => 🍒 === '=' ? parsevalue(io) : [raw"$(true)"])
        🍒 ∈ ">/" && skip(io, -1)
    end
    return attrs
end

@doc raw"""
    parsekey(io::IO)

Parse an attribute key.

```jldoctest
julia> HTM.parsekey(IOBuffer("class=\\"fruit\\">🍍..."))
:class
```
"""
@inline parsekey(io::IO) = Symbol(readuntil(isspace ⩔ ∈("=>/"), io))

@doc raw"""
    parsevalue(io::IO)

Parse an attribute value.

```jldoctest
julia> HTM.parsevalue(IOBuffer("\\"fruit\\">🍍..."))
1-element Vector{String}:
 "fruit"
```
"""
@inline parsevalue(io::IO) = startswith(io, ('"', '\'')) ? parsequotedvalue(io) : parseunquotedvalue(io)
@inline function parsequotedvalue(io::IO)
    🥝 = read(io, Char)
    🧩 = String[]
    while !(eof(io) || startswith(io, 🥝))
        push!(🧩, skipstartswith(io, "\\\$") ? "\$" : parseinterp(∈((🥝, '$', '\\')), io))
    end
    skipchars(isequal(🥝), io)
    return 🧩
end
@inline function parseunquotedvalue(io::IO)
    startswith(io, "http") && return [parseinterp(isspace ⩔ isequal('>'), io)]
    let f = isspace ⩔ ∈(">/\$\\")
        return skipstartswith(io, "\\\$") ? ["\$", readuntil(f, io)] : [parseinterp(f, io)]
    end
end

@doc raw"""
    parseinterp(io::IO)
    parseinterp(fallback, io::IO)

Parse an interpolation as string, including `$`.

The input must start with `$` if no fallback function is given.
The fallback function is passed to `HTM.readuntil` if the input does not start
with `$`.

```jldoctest
julia> HTM.parseinterp(IOBuffer(raw"$((1, (2, 3)))..."))
"\$((1, (2, 3)))"
```
"""
@inline function parseinterp(io::IO)
    buf = IOBuffer()
    write(buf, read(io, Char))
    (eof(io) || isspace(peek(io, Char))) && return "\$"  # frustrated interp
    # TODO: call Meta.parse here and avoid this while. This will require
    # returning Expr, which will remove some duty from strings. This can lead
    # to better decisions when calling create_element by being able to
    # distinguish between strings and expressions, which is currently
    # impossible.
    if startswith(io, '(')
        n = 1
        write(buf, read(io, Char))
        while n > 0
            🍒 = read(io, Char)
            if 🍒 === '('
                n += 1
            elseif 🍒 === ')'
                n -= 1
            end
            write(buf, 🍒)
        end
    else
        write(buf, readuntil(isspace ⩔ ∈("<>/\"\'=\$\\"), io))
    end
    return String(take!(buf))
end
@inline parseinterp(fallback, io::IO) = startswith(io, '$') ? parseinterp(io) : readuntil(fallback, io)


text:
function _write_array{F<:AbstractFloat}(io::IO, arr::Array{F})
    for v in arr
        @printf(io, "%.16e ", v)
    end
end
function _write_array{I<:Integer}(io::IO, arr::Array{I})
    for v in arr
        @printf(io, "%d ", v)
    end
end

function Base.write(io::IO, gmr::GMR)
    println(io, "Gaussian Mixture Regressor:")
    println(io, "n_targets: ", n_targets(gmr))
    println(io, "n_features: ", n_features(gmr))
    println(io, "n_components: ", n_components(gmr))
    println(io, "vec_A:")
    for (i, arr) in enumerate(gmr.vec_A)
        @printf(io, "\t%d\t", i); _write_array(io, arr); print(io, "\n")
    end
    println(io, "vec_b:")
    for (i, arr) in enumerate(gmr.vec_b)
        @printf(io, "\t%d\t", i); _write_array(io, arr); print(io, "\n")
    end
    println(io, "mixture_Obs:")
    print(io, "\tweights: "); _write_array(io, probs(gmr.mixture_Obs)); print(io, "\n")
    for (i, component) in enumerate(components(gmr.mixture_Obs))
        @printf(io, "\t%d\n", i)
        print(io, "\t\tμ: "); _write_array(io, mean(component)); print(io, "\n")
        print(io, "\t\tΣ: "); _write_array(io, cov(component)); print(io, "\n")
    end
    println(io, "mixture_Act_given_Obs:")
    for (i, component) in enumerate(components(gmr.mixture_Act_given_Obs))
        @printf(io, "\t%d\n", i)
        print(io, "\t\tΣ: "); _write_array(io, cov(component)); print(io, "\n")
    end
end
function Base.read{S<:AbstractString}(io::IO, ::Type{GMR}, lines::Vector{S}=readlines(io), line_index::Int=1)

    function advance!()
        line = strip(lines[line_index])
        line_index += 1
        line
    end
    function advance_until!(str::AbstractString)
        while !contains(lines[line_index], str)
            line_index += 1
        end
        advance!()
    end

    n_targets = parse(Int, split(advance_until!("n_targets"))[2])
    n_features = parse(Int, split(advance_until!("n_features"))[2])
    n_components = parse(Int, split(advance_until!("n_components"))[2])

    vec_A = Array(Matrix{Float64}, n_components)
    vec_b = Array(Vector{Float64}, n_components)
    vec_G = Array(MvNormal, n_components)
    vec_H = Array(MvNormal, n_components)

    advance_until!("vec_A")
    for i in 1 : n_components
        line = advance!()
        tokens = split(line)
        @assert i == parse(Int, tokens[1])
        vec_A[i] = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_targets, n_features), 1:n_targets*n_features)
    end

    advance_until!("vec_b")
    for i in 1 : n_components
        line = advance!()
        tokens = split(line)
        @assert i == parse(Int, tokens[1])
        vec_b[i] = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_targets), 1:n_targets)
    end

    advance_until!("mixture_Obs")
    line = advance!()
    tokens = split(line)
    @assert contains(tokens[1], "weights")
    weights = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_components), 1:n_components)
    components = Array(MvNormal, n_components)
    for i in 1 : n_components
        line = advance!()
        @assert i == parse(Int, line)

        line = advance!()
        tokens = split(line)
        @assert contains(tokens[1], "μ")
        μ = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_features), 1:n_features)

        line = advance!()
        tokens = split(line)
        @assert contains(tokens[1], "Σ")
        Σ = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_features, n_features), 1:n_features*n_features)

        components[i] = MvNormal(μ, Σ)
    end
    mixture_Obs = MixtureModel(components, weights)

    advance_until!("mixture_Act_given_Obs")
    components = Array(MvNormal, n_components)
    for i in 1 : n_components
        line = advance!()
        @assert i == parse(Int, line)

        line = advance!()
        tokens = split(line)
        @assert contains(tokens[1], "Σ")
        Σ = map!(j->parse(Float64, tokens[j+1]), Array(Float64, n_targets, n_targets), 1:n_targets*n_targets)

        components[i] = MvNormal(Array(Float64, n_targets), Σ)
    end
    mixture_Act_given_Obs = MixtureModel(components)

    GMR(vec_A, vec_b, mixture_Obs, mixture_Act_given_Obs)
end

function Base.write{A,F}(io::IO, model::GaussianMixtureRegressionDriver{A,F})
    println(io, "Gaussian Mixture Regressor Driver")
    @printf(io, "rec: %d %.3f\n", record_length(model.rec), model.rec.timestep)
    @printf(io, "context: %.3f %d\n", model.context.Δt, model.context.n_integration_steps)

    gmr = model.gmr
    write(io, model.gmr)

    println(io, "action: ", A)
    print(io, "standardization_means: "); _write_array(io, model.extractor.μ); print(io, "\n")
    print(io, "standardization_stdevs: "); _write_array(io, model.extractor.σ); print(io, "\n")
    print(io, "chosen_indeces: "); _write_array(io, model.extractor.extractor.subset); print(io, "\n")
end
function Base.read(io::IO, ::Type{GaussianMixtureRegressionDriver}, extractor::AbstractFeatureExtractor)
    lines = readlines(io)
    line_index = 1

    function advance!()
        line = strip(lines[line_index])
        line_index += 1
        line
    end
    function advance_until!(str::AbstractString)
        while !contains(lines[line_index], str)
            line_index += 1
        end
        advance!()
    end

    line = advance_until!("rec")
    tokens = split(line)
    rec = SceneRecord(parse(Int, tokens[2]), parse(Float64, tokens[3]))

    line = advance_until!("context")
    tokens = split(line)
    context = IntegratedContinuous(parse(Float64, tokens[2]), parse(Int, tokens[3]))

    gmr = read(io, GMR, lines, line_index)

    line = advance_until!("action")
    if contains(line, "AccelTurnrate")
        A = AccelTurnrate
    elseif contains(line, "AccelDesang")
        A = AccelDesang
    elseif contains(line, "LatLonAccel")
        A = LatLonAccel
    else
        warn("unknown action $line")
        warn("defaulting to AccelTurnrate")
        A = AccelTurnrate
    end

    line = advance_until!("standardization_means")
    tokens = split(line)
    μ = map!(token->parse(Float64, token), Array(Float64, length(tokens)-1), tokens[2:end])

    line = advance_until!("standardization_stdevs")
    tokens = split(line)
    σ = map!(token->parse(Float64, token), Array(Float64, length(tokens)-1), tokens[2:end])

    line = advance_until!("chosen_indeces")
    tokens = split(line)
    chosen_indicators = map!(token->parse(Int, token), Array(Int, length(tokens)-1), tokens[2:end])

    subset_extractor = SubsetExtractor(extractor, chosen_indicators)
    standardizing_extractor = StandardizingExtractor(subset_extractor, μ, σ)

    GaussianMixtureRegressionDriver{A, typeof(standardizing_extractor)}(rec, context,
        gmr, standardizing_extractor,
        Array(Float64, length(standardizing_extractor)),
        Array(Float64, length(A)))
end


text:
using Test
using Random
using Distributions
using SparseArrays

include("../src/LogisticRegression.jl")
lr = LogisticRegression
@testset "basic_functionality" begin
    @test lr.sigmoid(3) == 0.9525741268224334
    @test lr.sigmoid(.5) == 0.6224593312018546
    buncha_vals = [2,3,.5,1]
    res = lr.par_sigmoid(buncha_vals)
    @test res ==  [ 0.8807970779778823,
                    0.9525741268224334,
                    0.6224593312018546,
                    0.7310585786300049 ]
end
@testset "logistic regression" begin
    num_instances = 5000
    # generate dense training data
    x1_v1 = rand(Normal(0, 2), num_instances)
    x1_v2 = rand(Normal(0, .3), num_instances)
    x2_v1 = rand(Normal(4, 1), num_instances)
    x2_v2 = rand(Normal(2, .3), num_instances)

    x1 = [x1_v1 x1_v2]
    x2 = [x2_v1 x2_v2]
    data = [x1 ; x2]
    labels = [zeros(num_instances) ; ones(num_instances)]

    # generate dense test data
    x3_v1 = rand(Normal(0, 2), num_instances)
    x3_v2 = rand(Normal(0, .3), num_instances)
    x4_v1 = rand(Normal(4, 1), num_instances)
    x4_v2 = rand(Normal(2, .3), num_instances)
    x3 = [x3_v1 x3_v2]
    x4 = [x4_v1 x4_v2]
    test_data = [x3 ; x4]
    test_labels = [zeros(num_instances) ; ones(num_instances)]

    model = lr.LogisticRegressor(add_intercept=true)
    lr.fit!(model, data, labels)
    test_preds = lr.predict(model, test_data)
    acc = sum(test_preds .== test_labels)/length(test_labels)
    @test acc > .95 # this should be very easy to do well
    # since the distributions have exactly the same properties for
    # test and train and there's minimal overlap in classes
end
@testset "sparse logistic regression" begin
    num_instances = 5000
    # generate sparse training data
    x1_v1 = sparse(rand(Binomial(4, .2), num_instances))
    x1_v2 = sparse(rand(Binomial(8, .1),
    num_instances))
    x2_v1 = sparse(rand(Binomial(4, .8), num_instances))
    x2_v2 = sparse(rand(Binomial(8, .7), num_instances))

    x1 = [x1_v1 x1_v2]
    x2 = [x2_v1 x2_v2]
    data = [x1 ; x2]
    labels = [zeros(num_instances) ; ones(num_instances)]

    # generate sparse test data
    x3_v1 = sparse(rand(Binomial(4, .2), num_instances))
    x3_v2 = sparse(rand(Binomial(8, .1),
    num_instances))
    x4_v1 = sparse(rand(Binomial(4, .8), num_instances))
    x4_v2 = sparse(rand(Binomial(8, .7), num_instances))
    x3 = [x3_v1 x3_v2]
    x4 = [x4_v1 x4_v2]
    test_data = [x3 ; x4]
    test_labels = [zeros(num_instances) ; ones(num_instances)]

    model = lr.LogisticRegressor(add_intercept=true)
    lr.fit!(model, data, labels)
    test_preds = lr.predict(model, test_data)
    acc = sum(test_preds .== test_labels)/length(test_labels)
    @test acc > .95 # this should be very easy to do well
    # since the distributions have exactly the same properties for
    # test and train and there's minimal overlap in classes
end


text:
Base.names(::DynamicPPL.VarInfo{<:NamedTuple{names}}) where {names} = names
Base.names(::Type{<:DynamicPPL.VarInfo{<:NamedTuple{names}}}) where {names} = names

"""
    bijector(varinfo::DynamicPPL.VarInfo)

Returns a `NamedBijector` which can transform different variants of `varinfo`.

E.g. `ComponentArrays.ComponentArray(varinfo)`, `namedtuple(varinfo)`.
"""
@generated function Bijectors.bijector(varinfo::DynamicPPL.TypedVarInfo; tuplify = false)
    names = Base.names(varinfo)

    expr = Expr(:tuple)
    for n in names
        e = quote
            if tuplify
                $(Bijectors.Stacked)(map($(Bijectors).bijector, tuple(md.$n.dists...)), md.$n.ranges)
            else
                $(Bijectors.Stacked)(map($(Bijectors).bijector, md.$n.dists), md.$n.ranges)
            end
        end
        push!(expr.args, e)
    end

    return quote
        md = varinfo.metadata
        bs = NamedTuple{$names}($expr)
        return $(Bijectors).NamedBijector(bs)
    end
end

#####################################################
### Try to optimize the structure of the bijector ###
#####################################################
iscontiguous(left, right) = maximum(left) + 1 == minimum(right)

Bijectors.up1(b::Bijectors.Stacked) = b

optimize_bijector(b) = b
optimize_bijector(ib::Bijectors.Inverse) = inv(optimize_bijector(ib.orig))

function optimize_bijector(b::Bijectors.NamedBijector{names}) where {names}
    bs = map(b.bs) do b
        optimize_bijector(b)
    end
    return Bijectors.NamedBijector(NamedTuple{names}(bs))
end

# TODO: Maybe also handle the same cases as the evaluation of `TruncatedBijector`,
# e.g. `Bijectors.Log{N}() ∘ Bijectors.Scale{N}(-lb)` if `all(isfinite.(lb))` i.e.
# is lowerbounded.
function optimize_bijector(b::Bijectors.TruncatedBijector{N}) where {N}
    lb, ub = b.lb, b.ub
    return if all(iszero.(lb)) && all(isinf.(ub))
        Bijectors.Log{N}()
    elseif all(isfinite.(lb)) && all(isfinite.(ub))
        Bijectors.Logit{N}(lb, ub)
    else
        b
    end
end

function optimize_bijector(b::Bijectors.Stacked)
    n = length(b.ranges)

    segments = UnitRange[]
    start_idx = 1
    i = start_idx
    while start_idx <= n
        while i < n && iscontiguous(b.ranges[i], b.ranges[i + 1]) && b.bs[i] == b.bs[i + 1]
            i += 1
        end

        push!(segments, start_idx:i)
        i += 1
        start_idx = i
    end

    if b.bs isa Tuple
        segments = tuple(segments...)
    end

    bs = map(segments) do r
        b_r = b.bs[minimum(r)]
        if length(r) > 1 && Bijectors.dimension(b_r) == 0
            optimize_bijector(Bijectors.up1(b_r))
        else
            optimize_bijector(b_r)
        end
    end

    ranges = map(segments) do r
        minimum(b.ranges[minimum(r)]):maximum(b.ranges[maximum(r)])
    end

    @assert length(bs) == length(ranges) "number of bijectors and ranges is different"

    # Nvm; this actually causes issues for certain bijectors which we can't do `up1` to.
    # # If we're left with only a single bijector, don't wrap in `Stacked`.
    # length(bs) == 1 && return optimize_bijector(b.bs[1])

    return Bijectors.Stacked(bs, ranges)
end


text:
abstract type NanosecondTime      <: AbstractTime   end
abstract type NanosecondTimeBase  <: NanosecondTime end  # to define Periods, CompoundPeriods
abstract type NanosecondBasis     <: NanosecondTime end  # a structural trait, inherited

struct TimeDate <: NanosecondBasis
    attime::Time
    ondate::Date

    # ensure other constructors will be give explictly

    function TimeDate(attime::Time, ondate::Date)
        return new(attime, ondate)
    end
    function TimeDate(ondate::Date, attime::Time)
        return new(attime, ondate)
    end
end

@inline at_time(x::TimeDate) = x.attime
@inline on_date(x::TimeDate) = x.ondate

@inline at_time(x::DateTime) = Time(x)
@inline on_date(x::DateTime) = Date(x)

@inline at_time(x::Time) = x
@inline at_time(x::Date) = Time(0)
@inline on_date(x::Date) = x

TimeDate(x::TimeDate) = x

TimeDate(x::DateTime) = TimeDate(at_time(x), on_date(x))
TimeDate(x::Date) = TimeDate(at_time(x), on_date(x))
TimeDate(x::Time) = TimeDate(x, on_date(Dates.now()))

@inline function TimeDate(x::ZonedDateTime)
    datetime = DateTime(x)
    return TimeDate(datetime)
end

@inline Time(x::TimeDate) = at_time(x)
@inline Date(x::TimeDate) = on_date(x)

@inline fastpart(x::Time) = Nanosecond(x) + Microsecond(x)
@inline slowpart(x::Time) = x - (Nanosecond(x) + Microsecond(x))

@inline fastpart(x::Date) = Nanosecond(0)
@inline slowpart(x::Date) = Millisecond(0)

@inline fastpart(x::DateTime) = Nanosecond(0)
@inline slowpart(x::DateTime) = Time(x)

@inline function fastpart(td::TimeDate)
    fast_time = fastpart(at_time(td))
    return isempty(fast_time) ? Nanosecond(0) : fast_time
end
@inline slowpart(x::TimeDate) = slowpart(at_time(x))

timedate(x::TimeDate) = at_time(x), on_date(x)
timedate(x::DateTime) = at_time(x), on_date(x)
timedate(x::Date) = at_time(x), on_date(x)

@inline function timedate(x::ZonedDateTime)
    datetime = DateTime(x)
    return at_time(datetime), on_date(datetime)
end

@inline function DateTime(x::TimeDate)
    tim, dat = at_time(x), on_date(x)
    tim = slowpart(tim)
    return dat + tim
end

function TimeDate(y::Int64, m::Int64=1, d::Int64=1,
                  h::Int64=0, mi::Int64=0, s::Int64=0, ms::Int64=0,
                  us::Int64=0, ns::Int64=0)
  fl, ns = fldmod(ns, NANOSECONDS_PER_MICROSECOND)
  us += fl
  fl, us = fldmod(us, MICROSECONDS_PER_MILLISECOND)
  ms += fl
  fl, ms = fldmod(ms, MILLISECONDS_PER_SECOND)
  s += fl
  fl, s = fldmod(s, SECONDS_PER_MINUTE)
  mi += fl
  fl, mi = fldmod(mi, MINUTES_PER_HOUR)
  h += fl
  fl, h = fldmod(h, HOURS_PER_DAY)
  d += fl
  my, m = fldmod(m, MONTHS_PER_YEAR)
  y += my

  dt = Date(y, m, d)
  tm = Time(h, mi, s, ms, us, ns)
  td = TimeDate(dt, tm)
  return td
end

@inline function ZonedDateTime(tmdt::TimeDate, tz::T) where {T<:AkoTimeZone}
    dtm = DateTime(tmdt)
    return ZonedDateTime(dtm, tz)
end


text:
# Autogenerated wrapper script for Tcl_jll for x86_64-linux-gnu
export libtcl

using Zlib_jll
JLLWrappers.@generate_wrapper_header("Tcl")
JLLWrappers.@declare_library_product(libtcl, "libtcl8.6.so")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libtcl,
        "lib/libtcl8.6.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ModTIFF

using Libtiff_jll

libtiff_wrap_dir = joinpath(@__DIR__, "..", "..", "gen", "libtiff")

using CEnum
include(joinpath(libtiff_wrap_dir, "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t
include(joinpath(libtiff_wrap_dir, "libtiff_fixes.jl"))
include(joinpath(libtiff_wrap_dir, "libtiff_common.jl"))
include(joinpath(libtiff_wrap_dir, "libtiff_api.jl"))

version() = unsafe_string(TIFFGetVersion())

end # module


text:
array_value(e::Identity) = e

array_point(x) = ArrayMPoint(x)
array_point(e::Identity) = e

function inv(M::ArrayManifold, x; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    y = array_point(inv(M.manifold, array_value(x)))
    is_manifold_point(M, y, true; kwargs...)
    return y
end

function inv!(M::ArrayManifold, y, x; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    inv!(M.manifold, array_value(y), array_value(x))
    is_manifold_point(M, y, true; kwargs...)
    return y
end

function identity(M::ArrayManifold, x; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    y = array_point(identity(M.manifold, array_value(x)))
    is_manifold_point(M, y, true; kwargs...)
    return y
end

function identity!(M::ArrayManifold, y, x; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    identity!(M.manifold, array_value(y), array_value(x))
    is_manifold_point(M, y, true; kwargs...)
    return y
end

function compose(M::ArrayManifold, x, y; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    z = array_point(compose(M.manifold, array_value(x), array_value(y)))
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function compose!(M::ArrayManifold, z, x, y; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    compose!(M.manifold, array_value(z), array_value(x), array_value(y))
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function translate(M::ArrayManifold, x, y, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    z = array_point(translate(M.manifold, array_value(x), array_value(y), conv))
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function translate!(M::ArrayManifold, z, x, y, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    translate!(M.manifold, array_value(z), array_value(x), array_value(y), conv)
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function inverse_translate(M::ArrayManifold, x, y, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    z = array_point(inverse_translate(M.manifold, array_value(x), array_value(y), conv))
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function inverse_translate!(M::ArrayManifold, z, x, y, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    inverse_translate!(M.manifold, array_value(z), array_value(x), array_value(y), conv)
    is_manifold_point(M, z, true; kwargs...)
    return z
end

function translate_diff(M::ArrayManifold, x, y, v, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    is_tangent_vector(M, y, v, true; kwargs...)
    vout = ArrayTVector(translate_diff(
        M.manifold,
        array_value(x),
        array_value(y),
        array_value(v),
        conv,
    ))
    xy = translate(M, x, y, conv)
    is_tangent_vector(M, xy, vout, true; kwargs...)
    return vout
end

function translate_diff!(M::ArrayManifold, vout, x, y, v, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    is_tangent_vector(M, y, v, true; kwargs...)
    translate_diff!(
        M.manifold,
        array_value(vout),
        array_value(x),
        array_value(y),
        array_value(v),
        conv,
    )
    xy = translate(M, x, y, conv)
    is_tangent_vector(M, xy, vout, true; kwargs...)
    return vout
end

function inverse_translate_diff(M::ArrayManifold, x, y, v, conv::ActionDirection; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    is_tangent_vector(M, y, v, true; kwargs...)
    vout = ArrayTVector(inverse_translate_diff(
        M.manifold,
        array_value(x),
        array_value(y),
        array_value(v),
        conv,
    ))
    xinvy = inverse_translate(M, x, y, conv)
    is_tangent_vector(M, xinvy, vout, true; kwargs...)
    return vout
end

function inverse_translate_diff!(
    M::ArrayManifold,
    vout,
    x,
    y,
    v,
    conv::ActionDirection;
    kwargs...,
)
    is_manifold_point(M, x, true; kwargs...)
    is_manifold_point(M, y, true; kwargs...)
    is_tangent_vector(M, y, v, true; kwargs...)
    inverse_translate_diff!(
        M.manifold,
        array_value(vout),
        array_value(x),
        array_value(y),
        array_value(v),
        conv,
    )
    xinvy = inverse_translate(M, x, y, conv)
    is_tangent_vector(M, xinvy, vout, true; kwargs...)
    return vout
end


text:
module NewsAPI

import HTTP
using Dates
using Parameters
using JSON3

function __init__()
    if !haskey(ENV, "NEWS_API_KEY")
        @warn "This package isn't usable without ENV[\"NEWS_API_KEY\"] defined."
        @info "Visit https://newsapi.org to get an API key."
        ENV["NEWS_API_KEY"] = "0"
    end
end

const auth_header = HTTP.Header["authorization" => haskey(ENV, "NEWS_API_KEY") ? ENV["NEWS_API_KEY"] : "0"]

abstract type AbstractQuery end

function add_params(query::String, aq::T; verbose=true) where {T <: AbstractQuery}
    for field in fieldnames(T)
        item = getfield(aq, field)
        !isnothing(item) && (query *= "$field=$item&")
    end
    query = rstrip(query, '&')
    query = replace(query, ' ' => "%20")
    query = replace(query, "'" => "%27")
    verbose && @info "Full query:" query
    return query
end

function Base.show(io::IO, q::T) where {T <: AbstractQuery}
    println(io, T)
    for field in fieldnames(T)
        item = getfield(q, field)
        !isnothing(item) && println("  > $field:   $item")
    end
end

#-----------------------------------------------------------------------# Sources
@with_kw_noshow struct Sources <: AbstractQuery
    category    = nothing
    language    = nothing
    country     = nothing
end

function HTTP.get(s::Sources; verbose=true)
    query = add_params("https://newsapi.org/v2/sources?", s)
    JSON3.read(HTTP.get(query, headers=auth_header).body)
end

"""
    sources(; category=nothing, language=nothing, country=nothing, verbose=true)

Reference: https://newsapi.org/docs/endpoints/sources
"""
sources(;verbose=true, kw...) = HTTP.get(Sources(;kw...); verbose=verbose)

#-----------------------------------------------------------------------# TopHeadlines
@with_kw_noshow struct TopHeadlines <: AbstractQuery
    country     = nothing
    category    = nothing
    sources     = nothing
    q           = nothing
    pageSize    = nothing
    page        = nothing
end

function HTTP.get(th::TopHeadlines; verbose=true)
    query = add_params("https://newsapi.org/v2/top-headlines?", th; verbose=verbose)
    return JSON3.read(HTTP.get(query, headers=auth_header).body)
end

"""
    topheadlines(; country, category, sources, q, pageSize, page)

Reference: https://newsapi.org/docs/endpoints/top-headlines
"""
topheadlines(;verbose=true, kw...) = HTTP.get(TopHeadlines(; kw...), verbose=verbose)

#-----------------------------------------------------------------------# Everything
@with_kw_noshow struct Everything <: AbstractQuery
    q               = nothing
    qInTitle        = nothing
    sources         = nothing
    domains         = nothing
    excludeDomains  = nothing
    from            = nothing
    to              = nothing
    language        = nothing
    sortBy          = nothing
    pageSize        = nothing
    page            = nothing
end

function HTTP.get(e::Everything; verbose=true)
    query = add_params("https://newsapi.org/v2/everything?", e; verbose=verbose)
    return JSON3.read(HTTP.get(query, headers=auth_header).body)
end

"""
    everything(q, qInTitle, sources, domains, excludeDomains, from, to, language,
               sortBy, pageSize, page)

Reference: https://newsapi.org/docs/endpoints/everything
"""
everything(;verbose=true, kw...) = HTTP.get(Everything(;kw...); verbose=verbose)

end # module