text:
module ImputationQualityScores
	export dosage_to_likelihoods,likelihoods_to_dosage, rsquared_hat, proper_info,
		bhattacharya_coefficient

	using Distributions

	"The Dosage represents the predicted dosage of the non reference `bb` allele given the data available,
	it will always have a value between 0 and 2.
	The formula is Dosage = Pr(Het|Data) + 2*Pr(Alt|Data)"
	function likelihoods_to_dosage(p_aa, p_ab, p_bb)
		p_bb * 2 + p_ab
	end

	likelihoods_to_dosage{T<:AbstractFloat}(triple::Tuple{T,T,T}) = likelihoods_to_dosage(triple[1],triple[2],triple[3])

	# Note that there is a one-to-infinite correspondence between doages and genotype likelihoods.
	# This is just one of many possible conversions
	function dosage_to_likelihoods(dosage)
		β=Beta(1,2)
		(pdf(β,dosage)/2, pdf(β,abs(dosage-1))/2, pdf(β,abs(dosage-2))/2)
	end

	function minor_allele_frequency(dosages)
		μ=mean(dosages)
		(μ > 1.0) ? 2-μ : μ
	end

	function p_hwe(likelihoods)
	end

	# http://www.pypedia.com/index.php/imputation_r_square_hat
	"Calculates the (estimated) fraction of variance in unobserved 0/1/2 genotype explained by the the individual mean genotypes."
	function rsquared_hat(dosages::Array{Float64,1})
		p = mean(dosages) / 2
		return var(dosages) / 2 * p * (1-p)
	end

	# Adapted from Python code from
	function rsquared_hat(likelihoods::Array{Tuple{Float64,Float64,Float64},1})
		rsquared_hat(likelihoods_to_dosage.(likelihoods))
	end

	# Adapted from http://www.pypedia.com/index.php/imputation_r_square_hat
	function proper_info(likelihoods)
		n = length(likelihoods)
		meanAB = sum([l[2] for l in likelihoods])
		meanBB = sum([l[3] for l in likelihoods])
		sumX2 = sum([x[2] + (4.0 *x[3]) for x in likelihoods])
		sumXbar2 = sum([(x[2] + (2.0 * x[3]))^2.0 for x in likelihoods])
		meanX = (meanAB+2.*meanBB)/n
		rSqHat = (sumXbar2/n-(meanX^2))/(sumX2/n-(meanX^2))
	end

	# check probs
	function check_probs(pr)
		if !isprobvec(pr)
			throw("$pr is not probability distribution")
		end
		pr
	end

	# To simulate probabilities
	probs(pr::Tuple{Float64,Float64,Float64}) = check_probs([pr...])
	probs(pr::Array{Float64,1}) = check_probs(pr)
	probs(pr::Array{Tuple{Float64,Float64,Float64},1}) = normalize(foldl(+, [0,0,0], pr),1)

	# Comparisons of known and imputed genotypes
	"The bhattacharya_coefficient measures overlap between to trinomial probability distributions"
	function bhattacharya_coefficient(known, imputed)
		@assert length(known) == length(imputed)
		sum(sqrt.(probs(known).*probs(imputed)))
	end

	"Hellinger score is a measure of distance between to probability distributions"
	hellinger_score(known,imputed) = 1-sqrt(1-bhattacharya_coefficient(known,imputed))

	"`sen_score` (Scaled Euclidian Norm) score measures distance between genotype-wise expectations of two distributions. Note, can only distinguish between dose differences, and not genotype likelihood differences."
	function sen_score(known,imputed)
		m_observed = 2-(probs(known)[2] + 2*probs(known[1]))
		m_imputed = 2-(probs(imputed)[2] + 2*probs(imputed[1]))
		1 - (((m_observed-m_imputed)^2)/4)
	end

end # module


text:
using Gnuplot
U  = LinRange(-π,π, 100)
V = LinRange(-π,π, 100)
r = 0.5
x = [1 + cos(u) + r * cos(u) * cos(v)  for u in U, v in V]
y = [r * sin(v)                        for u in U, v in V]
z = [sin(u) + r * sin(u) * cos(v)      for u in U, v in V]
bgcp1 = "set object rectangle from screen 0,0 to screen 1,1"
bgcp2 = " behind fillcolor '#212946' fillstyle solid noborder"
bgcolor = bgcp1*bgcp2

@gsp "set pm3d depthorder"
@gsp :- "set style fill solid 1.0"
@gsp :- "set pm3d lighting primary 0.35 specular 0.2"
@gsp :- x y z "w pm3d notit" palette(:Wistia) "set view 108,2"
@gsp :- "unset border" "set xyplane 0" "set view equal xyz"
@gsp :- "unset tics"
@gsp :- "unset colorbox"
@gsp :- bgcolor
save(term="pngcairo size 600,400", output="torus2.png") # hide

ext: jl

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using Test
using SafeTestsets

mkpath("Files")

@safetestset "E1 Tests" begin
    include("e1_test.jl")
end

@safetestset "E2 Tests" begin
    include("e2_test.jl")
end

@safetestset "E3 Tests" begin
    include("e3_test.jl")
end

@safetestset "E4 Tests" begin
    include("e4_test.jl")
end

rm("Files"; force=true, recursive=true)


text:
import SeisIO: sync
export start_end, slide, nearest_start_end, slide_ind, sync
const μs = 1.0e-6
const sμ = 1000000.0
"""

  start_end(C::SeisChannel)

Return start and endtimes of SeisChannel in DateTime format.
"""
function start_end(C::SeisChannel)
    return u2d.(SeisIO.t_win(C.t,C.fs) * 1e-6)
end

"""

  start_end(S::SeisData)

Return start and endtimes of SeisChannel in DateTime format.
"""
function start_end(S::SeisData)
    irr = falses(S.n)
    non_ts = findall(S.fs .== 0)
    irr[non_ts] .= true
    z = zero(Int64)
    start_times = Array{DateTime}(undef,S.n)
    end_times =  Array{DateTime}(undef,S.n)
    for i = 1:S.n
        start_times[i] = u2d(S.t[i][1,2] * μs)
        end_times[i] = u2d(S.t[i][end,2] + (irr[i] ? z : S.t[i][1,2] + sum(S.t[i][2:end,2]) +
        round(Int64, (length(S.x[i])-1)/(μs*S.fs[i]))) * μs)
        # ts data: start time in μs from epoch + sum of time gaps in μs + length of trace in μs
        # non-ts data: time of last sample
    end
    return start_times,end_times
end

"""
    slide(A, cc_len, cc_step, fs, starttime, endtime)

Cut `A` into sliding windows of length `cc_len` points and offset `cc_step` points.

# Arguments
- `A::AbstractArray`: 1D time series.
- `cc_len::Real`: Cross-correlation window length in seconds.
- `cc_step::Real`: Step between cross-correlation windows in seconds.
- `starttime::Float64`: Time of first sample in `A` in unix time.
- `endtime::Float64`: Time of last sample in `A` in unix time.

# Returns
- `out::Array`: Array of sliding windows
- `starts::Array`: Array of start times of each window, in Unix time. E.g to convert
        Unix time to date time, use u2d(starts[1]) = 2018-08-12T00:00:00
"""
function slide(A::AbstractArray, cc_len::Real, cc_step::Real, fs::AbstractFloat,
               starttime::Float64,endtime::Float64)
    N = size(A,1)
    window_samples = Int(cc_len * fs)
    starts = Array(range(starttime,stop=endtime,step=cc_step))
    ends = starts .+ cc_len .- 1. / fs
    ind = findlast(x -> x <= endtime,ends)
    starts = starts[1:ind]
    ends = ends[1:ind]

    # fill array with overlapping windows
    if cc_step == cc_len && N % cc_len == 0
        return Array(reshape(A,window_samples,N ÷ window_samples)),starts
    elseif cc_step == cc_len # disregard data from edge
        return Array(reshape(A[1 : N - N % window_samples], window_samples, N ÷ window_samples)),starts
    else # need overlap between windows
        out = Array{eltype(A),2}(undef, window_samples,length(starts))
        s = convert.(Int,round.((hcat(starts,ends) .- starttime) .* fs .+ 1.))
        @inbounds for ii in eachindex(starts)
            out[:,ii] .= @view(A[s[ii,1]:s[ii,2]])
        end

    end
  return out,starts
end

"""
    slide(C, cc_len, cc_step)

Cut `C` into equal length sliding windows.

# Arguments
- `C::SeisChannel`: SeisChannel.
- `cc_len::Real`: Cross-correlation window length in seconds.
- `cc_step::Real`: Step between cross-correlation windows in seconds.

# Returns
- `out::Array`: Array of sliding windows.
- `starts::Array`: Array of start times of each window, in Unix time. E.g to convert
        Unix time to date time, use u2d(starts[1]) = 2018-08-12T00:00:00
"""
slide(C::SeisChannel, cc_len::Real, cc_step::Real) = slide(C.x,cc_len,cc_step,C.fs,d2u.(start_end(C))...)

"""
    nearest_start_end(C::SeisChannel, cc_len::Float64, cc_step::Float64)

Return best possible start, end times for data in `C` given the `cc_step` and `cc_len`.
"""
function nearest_start_end(C::SeisChannel, cc_len::Float64, cc_step::Float64)
    su,eu = SeisIO.t_win(C.t,C.fs) * μs
    su = round(su,digits=4) # round due to numerical roundoff error
    eu = round(eu,digits=4) # round due to numerical roundoff error
    ideal_start = d2u(DateTime(Date(u2d(su)))) # midnight of same day
    starts = Array(range(ideal_start,stop=eu,step=cc_step))
    ends = starts .+ cc_len .- 1. / C.fs
    return starts[findfirst(x -> x >= su, starts)], ends[findlast(x -> x <= eu,ends)]
end

"""
    nearest_start_end(S::DateTime,E::DateTime, cc_len::Int, cc_step::Int)

Return best possible start, end times for given starttime `S` and endtime `E`.
"""
function nearest_start_end(S::DateTime, E::DateTime, fs::Float64, cc_len::Real, cc_step::Real)
    ideal_start = DateTime(Date(S)) # midnight of same day
    starts = Array(ideal_start:Second(cc_step):E)
    ends = starts .+ Second(cc_len) .- Millisecond(convert(Int,1. / fs * 1e3))
    return d2u(starts[findfirst(x -> x >= S, starts)]), d2u(ends[findlast(x -> x <= E,ends)])
end

function slide_ind(startslice::AbstractFloat,endslice::AbstractFloat,fs::AbstractFloat,t::AbstractArray)
    starttime = t[1,2] * 1e-6
    startind = convert(Int,round((startslice - starttime) * fs,digits=4)) + 1
    endind = convert(Int,round((endslice - starttime) * fs,digits=4)) + 1
    return startind,endind
end

"""
  sync(N,starttime,endtime)

Return all data in `N` that starts no earlier than `starttime` and terminates no later than `endtime`.

# Arguments
- `N::NoiseData`: NoiseData structure (`RawData`, `FFTData`, `CorrData`).
- `starttime::TimeType`: Earliest time to start.
- `endtime::TimeType`: Latest time to start.
"""
function sync(N::NoiseData,starttime::TimeType,endtime::TimeType)
    starttime = DateTime(starttime)
    endtime = DateTime(endtime)
    ind = findall(d2u(starttime) .<= N.t .<= d2u(endtime))
    if length(ind) == 0
        throw(ArgumentError("No data in $(typeof(N)) between $starttime and $endtime."))
    end
    return N[ind]
end


text:
module Board_Test

using Intan, SpikeSorting, DistributedArrays

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

#64 channel amp check

    myamp=RHD2164("PortA1")
    @test myamp == [0,8]
    myamp=RHD2164("PortA2")
    @test myamp == [1,9]
    myamp=RHD2164("PortB1")
    @test myamp == [2,10]
    myamp=RHD2164("PortB2")
    @test myamp == [3,11]
    myamp=RHD2164("PortC1")
    @test myamp == [4,12]
    myamp=RHD2164("PortC2")
    @test myamp == [5,13]
    myamp=RHD2164("PortD1")
    @test myamp == [6,14]
    myamp=RHD2164("PortD2")
    @test myamp == [7,15]

#32 channel amp check

    myamp=RHD2132("PortA1")
    @test myamp == [0]
    myamp=RHD2132("PortA2")
    @test myamp == [1]
    myamp=RHD2132("PortB1")
    @test myamp == [2]
    myamp=RHD2132("PortB2")
    @test myamp == [3]
    myamp=RHD2132("PortC1")
    @test myamp == [4]
    myamp=RHD2132("PortC2")
    @test myamp == [5]
    myamp=RHD2132("PortD1")
    @test myamp == [6]
    myamp=RHD2132("PortD2")
    @test myamp == [7]

myamp=RHD2164("PortA1")
myfpga=FPGA(1,myamp)

    @test myfpga.numDataStreams == 0
    @test myfpga.dataStreamEnabled == zeros(Int64,1,Intan.MAX_NUM_DATA_STREAMS)
    @test myfpga.amps == [0,8]

#=
Single Core Data Structure
=#
myamp=RHD2164("PortA1")
myt=Task_NoTask()
d=Debug(string(dirname(Base.source_path()),"/data/qq.mat"),"qq")

myfpga=FPGA(1,myamp)

(myrhd,ss,myfpgas)=makeRHD([myfpga],debug=d);

    @test myrhd.v == zeros(Int16, Intan.SAMPLES_PER_DATA_BLOCK,64)
    @test length(ss) == 64
    @test typeof(ss) == Array{SpikeSorting.Sorting_1,1}
    @test typeof(myrhd.buf) == Array{SpikeSorting.Spike,2}
    @test size(myrhd.buf,2) == 64
    @test myrhd.nums == zeros(Int64,64)
    @test myrhd.time == zeros(UInt32,Intan.SAMPLES_PER_DATA_BLOCK,1)

Intan.init_board!(myrhd,myfpgas)

    @test myfpgas[1].numDataStreams == 2
    @test myfpgas[1].dataStreamEnabled == [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    @test myfpgas[1].sampleRate == 30000
    @test myfpgas[1].numWords == Intan.SAMPLES_PER_DATA_BLOCK * (4+2+(2*36) +8+2)
    @test myfpgas[1].numBytesPerBlock == Intan.SAMPLES_PER_DATA_BLOCK * (4+2+(2*36) +8+2) * 2

#=
Multi Core Data Structure
=#
myamp=RHD2132("PortA1")
d=Debug(string(dirname(Base.source_path()),"/data/qq.mat"),"qq")
myt=Task_NoTask()
myfpga1=FPGA(1,myamp)
myfpga2=FPGA(2,myamp)
(myrhd2,ss2,myfpgas)=makeRHD([myfpga1,myfpga2],debug=d,parallel=true);

    @test myrhd2.v == zeros(Int64, Intan.SAMPLES_PER_DATA_BLOCK,64)
    @test length(ss2) == 64
    @test typeof(ss2) == DistributedArrays.DArray{SpikeSorting.Sorting_1,1,Array{SpikeSorting.Sorting_1,1}}
    @test typeof(myrhd2.buf) == SharedArray{SpikeSorting.Spike,2}
    @test size(myrhd2.buf,2) == 64
    @test myrhd2.nums == zeros(Int64,64)

#Intan.init_board!(myrhd2)

#=

    @test myrhd.fpga[1].numDataStreams --> 1
    @test myrhd.fpga[1].dataStreamEnabled --> [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    @test myrhd.fpga[1].sampleRate --> 30000
    @test myrhd.fpga[1].numWords --> 31200
    @test myrhd.fpga[1].numBytesPerBlock --> 62400

=#
#=
Sample Rate Testing
=#

myamp=RHD2164("PortA1")
myt=Task_NoTask()
d=Debug(string(dirname(Base.source_path()),"/data/qq.mat"),"qq")

myfpga=FPGA(1,myamp)

(myrhd,ss,myfpgas)=makeRHD([myfpga],debug=d);

Intan.init_board!(myrhd,myfpgas)

    for i in [1000,1250,1500,2000,2500,3000,3333,4000,5000,6250,8000,10000,12500,15000,20000,25000,30000]
        Intan.setSampleRate(myfpgas[1],i,myrhd.debug.state)
        @test myfpgas[1].sampleRate == i
    end

end


text:
# using GLMakie
import Pkg; Pkg.activate(".")
using Persephone
using LinearAlgebra, TimerOutputs
import Statistics:mean

function main()

    #=========================================================================
        INITIALISE OUTPUT PATHS
    =========================================================================#
    iscluster = false
    if iscluster
        path = "/storage2/unipd/navarro/AnnulusBenchmarks"
    else
        path = "/home/albert/Desktop/output"
    end
    folder = "NoHealingTrue"
    OUT, iplot = setup_output(path, folder)

    #=========================================================================
        RELOAD CHECKPOINT
    =========================================================================#
    load = false
    iload = 131
    i0 = 1

    #=========================================================================
        MAKE GRID
    =========================================================================#
    gr, IDs = init_grid(4; split=2) # 1st argument = N
    PhaseID = 1
    min_inradius = inrectangle(gr)/2
    GlobC = [Point2D{Polar}(gr.θ[i], gr.r[i]) for i in 1:gr.nnod] # → global coordinates

    #=========================================================================
        Boundary conditions
    =========================================================================#
    # T and U boundary conditions
    TBC, UBC = init_BCs(gr, IDs, Ttop = 0.0, Tbot = 1.0, type="free slip" )

    #=========================================================================
        Delete existing statistics file
    =========================================================================#
    ScratchNu, stats_file  = setup_metrics(gr, path, folder, load = load)

    #=========================================================================
        GET DEM STRUCTURE:
    =========================================================================#
    dem_file = joinpath("DEM", "DEM_1e-3_vol30.h5")
    # dem_file = joinpath("newDEM", "Dem_1e-3_vol30.h5")
    Δη, ϕ = 1e-3, 0.3
    D = getDEM(dem_file, Δη, ϕ)

    #=========================================================================
        FIX θ OF ELEMENTS CROSSING π and reshape also r (i.e. periodic boundaries)
    =========================================================================#
    θStokes, rStokes = elementcoordinate(GlobC, @views(gr.e2n[1:3, :]))
    θThermal, rThermal = elementcoordinate(GlobC, gr.e2n_p1)
    fixangles!(θStokes)
    fixangles!(θThermal)
    coordinates = ElementCoordinates(θStokes, rStokes)

    #=========================================================================
        IPs polar coordinates (Some stuff needs to be deprecated)
    =========================================================================#
    θ6, r6 =  elementcoordinate(GlobC, @views(gr.e2n[1:6, :]))
    θ3 = θStokes
    fixangles6!(θ6)
    fixangles!(θ3)
    ipx, ipz = getips(gr.e2n, θ6, r6 )
    ix, iz = polar2cartesian(ipx, ipz)
    transition = 2.22-0.2276
    isotropic_idx = findall(ipz .> transition)
    IntC = [@inbounds(Point2D{Polar}(ipx[i], ipz[i])) for i in CartesianIndices(ipx)] # → ip coordinates

    #=========================================================================
        INITIALISE PARTICLES
    =========================================================================#
    particle_info, particle_weights, particle_fields = particles_generator(
        θThermal, rThermal, IntC, gr.e2n_p1, number_of_particles = 8
    )

    #=========================================================================
        ALLOCATE/INITIALISE FIELDS
    =========================================================================#
    # Allocate velocity and pressure fields
    U, P = init_U_P(gr)
    # Allocate velocity gradient, stress and strain tensors:
    F, _, τ, ε,  = initstress(gr.nel)
    # FSE = [FiniteStrainEllipsoid(1.0, 0.0, 0.0, 1.0, 1.0, 1.0) for _ in 1:gr.nel, _ in 1:6]
    # Allocate viscosity tensor:
    𝓒 = 𝓒init(gr.nel, 7)
    # Allocate nodal velocities
    Ucartesian, Upolar = initvelocity(gr.nnod)
    # Initialise temperature @ nodes
    perturbation = :random
    T = init_temperature(gr, IDs, type = perturbation)
    ΔT = similar(T)
    # Initialise temperature @ particles
    init_particle_temperature!(particle_fields, particle_info, type = perturbation)

    # Annealing rate
    annealing = 1e-3

    # Finite Strain structure
    FSE = FiniteStrain(gr.nel,
        nip = 6,
        ϵ = 0, # a1/a2 at which fabric is destroyed
        annealing_rate = 0, # annealing rate
        r_iso = 0 # depth above which Ω is isotropic
    )

    #= Viscosity type. Options:
        (*) "IsoviscousIsotropic"
        (*) "TemperatureDependantIsotropic"
        (*) "IsoviscousAnisotropic"
        (*) "TemperatureDependantAnisotropic"
    =#
    viscosity_type = :TemperatureDependantAnisotropic

    # Physical parameters for thermal diffusion
    VarT = thermal_parameters(
        κ = 1.0,
        α = 1e-6,
        Cp = 1.0,
        dQdT = 0.0,
    )

    ρ = state_equation(VarT.α, T)
    η = getviscosity(T, viscosity_type, η = 1.81)
    	# η = 1 for isotropic
    	# η = 1.81 for anisotropic with phi = 30%
    	# η = 1/0.6899025321942348 for anisotropic with phi = 20%
    Valη = Val(η)
    g = 1e6
    𝓒 = anisotropic_tensor(FSE.fse, D, Valη)

    #=========================================================================
        SOLVER INVARIANTS (FOR AN IMMUTABLE MESH):
            TODO pre-compute the Jacobians
    =========================================================================#
    # Allocate rotation tensor
    RotationMatrices = rotation_matrix(gr.θ)
    # Allocate sparsity patterns of Stokes block matrices
    KKidx, GGidx, MMidx, = sparsitystokes(gr)
    # Allocate spasity pattern of thermal diffusion stiffness matrix
    CMidx, _ = sparsitythermal(gr.e2n, 6)

    # Stokes immutables
    KS, GS, MS, FS, DoF_U, DoF_P, nn, SF_Stokes, ScratchStokes = stokes_immutables(
        gr, gr.nnod, 2, 3, 6, gr.nel, 7
    )

    SF_Stress = stress_shape_functions()

    # Diffusion_immutables
    KT, MT, FT, DoF_T, valA, SF_Diffusion, ScratchDifussion = diffusion_immutables(
        gr, 2, 3, 6, 7
    )

    # Color elements
    _, color_list = color_mesh(gr.e2n)

    #=========================================================================
        LOAD PREVIOUS MODELS
    =========================================================================#
    if load == true
        i0 = iload + 1
        load_file = string("file_", iload, ".h5")
        toreload = joinpath(pwd(), path, folder, load_file)
        T, F = reloader(toreload)
        FSE, F = getFSE(F, FSE)
        𝓒 = anisotropic_tensor(FSE.fse, D, Valη)
    end

    #=========================================================================
        START SOLVER
    =========================================================================#
    to = TimerOutput()
    Time = 0.0
    T0 = deepcopy(T)

    for iplot in i0:350
        for _ in 1:75
            reset_timer!(to)

            #= Update material properties =#
            state_equation!(ρ, VarT.α, T)
            getviscosity!(η, T)

            #= Stokes solver using preconditioned-CG =#
            Ucartesian, Upolar, U, Ucart, P, to = solveStokes(
                U,
                P,
                gr,
                Ucartesian,
                Upolar,
                g,
                T,
                η,
                𝓒,
                coordinates,
                RotationMatrices,
                PhaseID,
                UBC,
                KKidx,
                GGidx,
                MMidx,
                to,
                solver = :pardiso
            );

            println("min:max Uθ ", extrema(@views U[1:2:end]))
            println("mean speed ", mean(@views @. (√(U[1:2:end]^2 + U[2:2:end]^2))))

            #= Adaptive time-step =#
            Δt = calculate_Δt(Ucartesian, gr.nθ, min_inradius)
            println("Δt = ", Δt)

            #= Stress-Strain postprocessor =#
            @timeit to "F" stress!(F, ε, Ucart, gr.nel, DoF_U, coordinates, SF_Stress, Δt)
            # εII = secondinvariant(ε)

            # isotropic_lithosphere!(F, isotropic_idx)
            # F = healing(F, FSE)
            # @timeit to "FSE" FSE = getFSE(F, FSE)
            # F0, FSE0 = deepcopy(F), deepcopy(FS)
            # @timeit to "FSE" FSE, F = getFSE_healing(F, FSE, ϵ=1e3)
            # @timeit to "FSE" getFSE_annealing!(F, FSE, annealing*(Time+Δt))

            # Finite Strain Ellipsoid calculation
            if FSE.isotropic_domain != 0 # check whether any region of the Earth is always isotropic
                isotropic_lithosphere!(F, isotropic_idx) # force isotropy
            end
            FSE, F = getFSE(F, FSE)

            #= Compute the viscous tensor =#
            @timeit to "viscous tensor" 𝓒 = anisotropic_tensor(FSE.fse, D, Valη)

            #=
                Diffusion solver
            =#
            T, T0, ΔT, to = solveDiffusion_threaded(
                color_list,
                CMidx,
                KT,
                MT,
                FT,
                DoF_T,
                coordinates,
                VarT,
                ScratchDifussion,
                SF_Diffusion,
                valA,
                ρ,
                Δt,
                T,
                T0,
                TBC,
                to,
            )

            @timeit to "Particle advenction" begin
                #=
                Particle advection and mappings
                    Fij : ip   |-> particle
                    T   : node |-> particle
                    Ui  : node |-> particle + advection
                =#
                @sync begin

                    @timeit to "F → particle" particle_fields = F2particle(
                        particle_fields, particle_info, ipx, ipz, F
                    )

                    @timeit to "T → particle" begin
                        interpolate_temperature!(
                            T0,
                            particle_fields,
                            gr,
                            ρ,
                            T,
                            particle_info,
                            particle_weights,
                            VarT,
                            gr.nθ,
                            gr.nr,
                            Δt,
                            ΔT,
                        )
                    end

                end

                @timeit to "advection" particle_info, particle_weights, to = advection_RK2(
                        particle_info,
                        gr,
                        particle_weights,
                        Ucartesian,
                        Δt,
                        θThermal,
                        rThermal,
                        coordinates,
                        IntC,
                        to,
                    )

                println("Min-max particle temperature = ", extrema(particle_fields.T))
            end

            @timeit to "Locate articles" begin
                #= Particle locations =#
                particle_info, particle_weights, found = tsearch_parallel(
                    particle_info, particle_weights, θThermal, rThermal, coordinates, gr, IntC
                )

                lost_particles = length(particle_info) - sum(found)
                check_corruption!(found, particle_fields)
                println("Lost particles: ", lost_particles, " Corrupted particles: ",  length(particle_info) - sum(found) -lost_particles)

                if lost_particles > 0
                    particle_info, particle_weights, particle_fields = purgeparticles(
                        particle_info, particle_weights, particle_fields, found
                    )
                end

            end

            @timeit to "Particle to node/ip" begin
                #=
                    Map back to original immutable locations
                        Fij : particle -> ip
                        T   : particle -> node
                =#
                @timeit to "F → ip" F = F2ip(
                    F, particle_fields, particle_info, particle_weights, gr.nel
                )

                T = T2node(
                    T,
                    particle_fields,
                    particle_info,
                    particle_weights,
                    gr,
                    IDs,
                )

                println("\n Min-max nodal temperature = ", extrema(T))
            end

            @timeit to "Add/reject particles" begin
                particle_info, particle_weights, particle_fields = addreject(
                    T,
                    F,
                    gr,
                    θThermal,
                    rThermal,
                    IntC,
                    particle_info,
                    particle_weights,
                    particle_fields,
                    min_num_particles = 3
                )
            end

            # healing!(F, FSE)

            println("mean T after advection  ", mean(T))

            write_stats(U, T, length(particle_info), gr, Time, ScratchNu, stats_file)

            Time += Δt
            show(to; compact=true)
        end

        #= Save output file =#
        println("\n time = ", Time)
        println("\n Saving output...")
        OUT = IOs(path, folder, "file", iplot)
        savedata(
            OUT,
            Upolar,
            Ucartesian,
            T,
            η,
            𝓒,
            ρ,
            F,
            FSE.fse,
            gr.nθ,
            gr.nr,
            particle_fields,
            particle_info,
            Time,
            Val(η),
        )
        println(" ...done!")

    end
end

main()

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
"""
    abstract type Rotation{N,T} <: StaticMatrix{N,N,T}

An abstract type representing `N`-dimensional rotations. More abstractly, they represent
unitary (orthogonal) `N`×`N` matrices.
"""
abstract type Rotation{N,T} <: StaticMatrix{N,N,T} end

Base.@pure StaticArrays.Size(::Type{Rotation{N}}) where {N} = Size(N,N)
Base.@pure StaticArrays.Size(::Type{Rotation{N,T}}) where {N,T} = Size(N,N)
Base.@pure StaticArrays.Size(::Type{R}) where {R<:Rotation} = Size(supertype(R))
Base.adjoint(r::Rotation) = inv(r)
Base.transpose(r::Rotation{N,T}) where {N,T<:Real} = inv(r)
Base.real(R::Type{<:Rotation}) = R

# Generate zero-matrix with SMatrix
# Note that zeros(Rotation3,dims...) is not Array{<:Rotation} but Array{<:StaticMatrix{3,3}}
Base.zero(::Rotation{N,T}) where {N,T} = @SMatrix zeros(T,N,N)
Base.zero(::Type{Rotation}) = error("The dimension of rotation is not specified.")
Base.zero(::Type{<:Rotation{N}}) where N = @SMatrix zeros(N,N)
Base.zero(::Type{<:Rotation{N,T}}) where {N,T} = @SMatrix zeros(T,N,N)
Base.zeros(::Type{R}) where {R<:Rotation} = zeros(R, ()) # avoid StaticArray constructor
Base.zeros(::Type{R}, dims::Base.DimOrInd...) where {R<:Rotation} = zeros(typeof(zero(R)),dims...)
Base.zeros(::Type{R}, dims::NTuple{N, Integer}) where {R<:Rotation, N} = zeros(typeof(zero(R)),dims)
Base.zeros(::Type{R}, dims::Tuple{}) where {R<:Rotation} = zeros(typeof(zero(R)),dims) # avoid ambiguity

# Generate identity rotation matrix
Base.one(r::Rotation) = one(typeof(r))

# Rotation angles and axes can be obtained by converting to the AngleAxis type
rotation_angle(r::Rotation{3}) = rotation_angle(AngleAxis(r))
rotation_axis(r::Rotation{3}) = rotation_axis(AngleAxis(r))

# `convert` goes through the constructors, similar to e.g. `Number`
Base.convert(::Type{R}, rot::R) where {N,R<:Rotation{N}} = rot
Base.convert(::Type{R}, rot::Rotation{N}) where {N,R<:Rotation{N}} = R(rot)

# Rotation matrices should be orthoginal/unitary. Only the operations we define,
# like multiplication, will stay as Rotations, otherwise users will get an
# SMatrix{3,3} (e.g. rot1 + rot2 -> SMatrix)
Base.@pure StaticArrays.similar_type(::Union{R,Type{R}}) where {R <: Rotation} = SMatrix{size(R)..., eltype(R), prod(size(R))}
Base.@pure StaticArrays.similar_type(::Union{R,Type{R}}, ::Type{T}) where {R <: Rotation, T} = SMatrix{size(R)..., T, prod(size(R))}

@inline function Base.:/(r1::Rotation, r2::Rotation)
    r1 * inv(r2)
end

@inline function Base.:\(r1::Rotation, r2::Rotation)
    inv(r1) * r2
end

@inline function Base.:\(r::Rotation, v::AbstractVector)
    inv(r) * v
end

# This definition is for avoiding anbiguity
@inline function Base.:\(r::Rotation, v::StaticVector)
    inv(r) * v
end

################################################################################
################################################################################
"""
    struct RotMatrix{N,T} <: Rotation{N,T}

A statically-sized, N×N unitary (orthogonal) matrix.

Note: the orthonormality of the input matrix is *not* checked by the constructor.
"""
struct RotMatrix{N,T,L} <: Rotation{N,T} # which is <: AbstractMatrix{T}
    mat::SMatrix{N, N, T, L} # The final parameter to SMatrix is the "length" of the matrix, 3 × 3 = 9
    RotMatrix{N,T,L}(x::AbstractArray) where {N,T,L} = new{N,T,L}(convert(SMatrix{N,N,T,L}, x))
    # fixes #49 ambiguity introduced in StaticArrays 0.6.5
    RotMatrix{N,T,L}(x::StaticArray) where {N,T,L} = new{N,T,L}(convert(SMatrix{N,N,T,L}, x))
end
RotMatrix(x::SMatrix{N,N,T,L}) where {N,T,L} = RotMatrix{N,T,L}(x)

Base.zero(::Type{RotMatrix}) = error("The dimension of rotation is not specified.")

# These functions (plus size) are enough to satisfy the entire StaticArrays interface:
@inline function RotMatrix(t::Tuple)
    n = sqrt(length(t))
    if !isinteger(n)
        throw(DimensionMismatch("The length of input tuple $(t) must be square number."))
    end
    N = Int(n)
    RotMatrix(SMatrix{N,N}(t))
end
@inline (::Type{RotMatrix{N}})(t::Tuple) where N = RotMatrix(SArray{Tuple{N,N}}(t))
@inline RotMatrix{N,T}(t::Tuple) where {N,T} = RotMatrix(SArray{Tuple{N,N},T}(t))
@inline RotMatrix{N,T,L}(t::Tuple) where {N,T,L} = RotMatrix(SArray{Tuple{N,N},T}(t))

# Create aliases RotMatrix2{T} = RotMatrix{2,T,4} and RotMatrix3{T} = RotMatrix{3,T,9}
for N in 2:3
    L = N*N
    RotMatrixN = Symbol(:RotMatrix, N)
    @eval begin
        const $RotMatrixN{T} = RotMatrix{$N, T, $L}
        @inline $RotMatrixN(t::Tuple) = RotMatrix(SArray{Tuple{$N,$N}}(t))
    end
end

Base.@propagate_inbounds Base.getindex(r::RotMatrix, i::Int) = r.mat[i]
@inline Base.Tuple(r::RotMatrix) = Tuple(r.mat)

@inline RotMatrix(θ::Number) = RotMatrix{2}(θ)
@inline function (::Type{<:RotMatrix{2}})(θ::Number)
    s, c = sincos(θ)
    RotMatrix(@SMatrix [c -s; s c])
end
@inline function RotMatrix{2,T}(θ::Number) where T
    s, c = sincos(θ)
    RotMatrix(@SMatrix T[c -s; s c])
end

Base.one(::Type{R}) where {N,R<:RotMatrix{N}} = R(I)

# A rotation is more-or-less defined as being an orthogonal (or unitary) matrix
Base.inv(r::RotMatrix) = RotMatrix(r.mat')

# By default, composition of rotations will go through RotMatrix, unless overridden
@inline Base.:*(r1::Rotation, r2::Rotation) = RotMatrix(r1) * RotMatrix(r2)
@inline Base.:*(r1::RotMatrix, r2::Rotation) = r1 * RotMatrix(r2)
@inline Base.:*(r1::Rotation, r2::RotMatrix) = RotMatrix(r1) * r2
@inline Base.:*(r1::RotMatrix, r2::RotMatrix) = RotMatrix(r1.mat * r2.mat) # TODO check that this doesn't involve extra copying.

# Special case multiplication of 3×3 rotation matrices: speedup using cross product
@inline function Base.:*(r1::RotMatrix{3}, r2::RotMatrix{3})
    ret12 = r1 * r2[:, SVector(1, 2)]
    ret3 = ret12[:, 1] × ret12[:, 2]
    RotMatrix([ret12 ret3])
end

"""
    struct Angle2d{T} <: Rotation{2,T}
        theta::T
    end

A 2×2 rotation matrix parameterized by a 2D rotation by angle.
Only the angle is stored inside the `Angle2d` type, values
of `getindex` etc. are computed on the fly.
"""
struct Angle2d{T} <: Rotation{2,T}
    theta::T
    Angle2d{T}(theta::Number) where T = new{T}(theta)
end

@inline function Angle2d(theta::Number)
    Angle2d{rot_eltype(typeof(theta))}(theta)
end

params(r::Angle2d) = SVector{1}(r.theta)

Angle2d(r::Rotation{2}) = Angle2d(rotation_angle(r))
Angle2d{T}(r::Rotation{2}) where {T} = Angle2d{T}(rotation_angle(r))
@inline (::Type{R})(t::Tuple) where R<:Angle2d = convert(R, RotMatrix2(t))

Base.one(::Type{A}) where {A<: Angle2d} = A(0)

rotation_angle(rot::Angle2d) = rot.theta
function rotation_angle(rot::Rotation{2})
    c = @inbounds rot[1,1]
    s = @inbounds rot[2,1]
    atan(s, c)
end

@inline function Base.:*(r::Angle2d, v::StaticVector)
    if length(v) != 2
        throw(DimensionMismatch("Cannot rotate a vector of length $(length(v))"))
    end
    x,y = v
    s,c = sincos(r.theta)
    T = Base.promote_op(*, typeof(s), eltype(v))
    similar_type(v,T)(c*x - s*y, s*x + c*y)
end

Base.:*(r1::Angle2d, r2::Angle2d) = Angle2d(r1.theta + r2.theta)
Base.:^(r::Angle2d, t::Real) = Angle2d(r.theta*t)
Base.:^(r::Angle2d, t::Integer) = Angle2d(r.theta*t)
Base.inv(r::Angle2d) = Angle2d(-r.theta)

@inline function Base.getindex(r::Angle2d, i::Int)
    if i == 1
        cos(r.theta)
    elseif i == 2
        sin(r.theta)
    elseif i == 3
        -sin(r.theta)
    elseif i == 4
        cos(r.theta)
    else
        throw(BoundsError(r,i))
    end
end

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

_isrotation_eps(T::Type{<:Integer}) = zero(T)
_isrotation_eps(T) = eps(T)

"""
    isrotation(r)
    isrotation(r, tol)

Check whether `r` is a rotation matrix, where `r * r'` is within
`tol` of the identity matrix (using the Frobenius norm). `tol` defaults to
`1000 * eps(float(eltype(r)))` or `zero(T)` for integer `T`.
"""
isrotation

function isrotation(r::StaticMatrix{N,N,T}, tol::Real = 1000 * _isrotation_eps(eltype(T))) where {N,T<:Real}
    m = SMatrix(r)
    d = norm(m*m'-one(SMatrix{N,N}))
    return d ≤ tol && det(m) > 0
end

function isrotation(r::AbstractMatrix{T}, tol::Real = 1000 * _isrotation_eps(eltype(T))) where T<:Real
    if !=(size(r)...)
        return false
    end
    d = norm(r*r'-one(r))
    return d ≤ tol && det(r) > 0
end

function isrotation(r::AbstractMatrix{T}, tol::Real = 1000 * _isrotation_eps(eltype(real(T)))) where T<:Number
    if !isreal(r)
        return false
    end
    return isrotation(real(r), tol)
end

"""
    nearest_rotation(M) -> RotMatrix

Get the nearest special orthonormal matrix from given matrix `M`.
See [Wahba's problem](https://en.wikipedia.org/wiki/Wahba%27s_problem) for more information.
"""
function nearest_rotation(M::StaticMatrix{N,N}) where N
    u, _, v = svd(M)
    s = sign(det(u * v'))
    d = @SVector ones(N-1)
    R = u * Diagonal(push(d,s)) * v'
    return RotMatrix{N}(R)
end

function nearest_rotation(M::AbstractMatrix{T}) where T
    N = size(M,1)
    L = N^2
    M_ = convert(SMatrix{N,N,T,L}, M)
    return nearest_rotation(M_)
end

# A simplification and specialization of the Base.show function for AbstractArray makes
# everything sensible at the REPL.
function Base.show(io::IO, ::MIME"text/plain", X::Rotation)
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end
    summary(io, X)
    if !isa(X, RotMatrix)
        n_fields = length(fieldnames(typeof(X)))
        print(io, "(")
        for i in 1:n_fields
            print(io, getfield(X, i))
            if i < n_fields
                print(io, ", ")
            end
        end
        print(io, ")")
    end
    print(io, ":")
    println(io)
    io = IOContext(io, :typeinfo => eltype(X))
    Base.print_array(io, X)
end


text:
#! /usr/bin/julia

# Rosetta Code, Circles of given radius through two points

immutable Point{T<:FloatingPoint}
    x::T
    y::T
end

immutable Circle{T<:FloatingPoint}
    c::Point{T}
    r::T
end
Circle{T<:FloatingPoint}(a::Point{T}) = Circle(a, zero(T))

using AffineTransforms

function circlepoints{T<:FloatingPoint}(a::Point{T}, b::Point{T}, r::T)
    cp = Circle{T}[]
    r >= 0 || return (cp, "No Solution, Negative Radius")
    if a == b
        if abs(r) < 2eps(zero(T))
            return (push!(cp, Circle(a)), "Point Solution, Zero Radius")
        else
            return (cp, "Infinite Solutions, Indefinite Center")
        end
    end
    ca = Complex(a.x, a.y)
    cb = Complex(b.x, b.y)
    d = (ca + cb)/2
    tfd = tformtranslate([real(d), imag(d)])
    tfr = tformrotate(angle(cb-ca))
    tfm = tfd*tfr
    u = abs(cb-ca)/2
    r-u > -5eps(r) || return(cp, "No Solution, Radius Too Small")
    if r-u < 5eps(r)
        push!(cp, Circle(apply(Point, tfm*[0.0, 0.0]), r))
        return return (cp, "Single Solution, Degenerate Centers")
    end
    v = sqrt(r^2 - u^2)
    for w in [v, -v]
        push!(cp, Circle(apply(Point, tfm*[0.0, w]), r))
    end
    return (cp, "Two Solutions")
end

tp = [Point(0.1234, 0.9876),
      Point(0.0000, 2.0000),
      Point(0.1234, 0.9876),
      Point(0.1234, 0.9876),
      Point(0.1234, 0.9876)]

tq = [Point(0.8765, 0.2345),
      Point(0.0000, 0.0000),
      Point(0.1234, 0.9876),
      Point(0.8765, 0.2345),
      Point(0.1234, 0.9876)]

tr = [2.0, 1.0, 2.0, 0.5, 0.0]

println("Testing circlepoints:")
for i in 1:length(tp)
    p = tp[i]
    q = tq[i]
    r = tr[i]
    (cp, rstatus) = circlepoints(p, q, r)
    println(@sprintf("(%.4f, %.4f), (%.4f, %.4f), %.4f => %s",
                     p.x, p.y, q.x, q.y, r, rstatus))
    for c in cp
        println(@sprintf("    (%.4f, %.4f), %.4f",
                         c.c.x, c.c.y, c.r))
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#=
organizational script for defining the LPBasic module

=#

"""
    LPBasic Module
Julia module implementing the LPBasic line planning problem.

This includes:
- defining structs for containing the relevant problema data
- implementing IO functionality to read relevant data from GTFS data into such structs
- a flexible functionality to construct JuMP models from such structs
- some visualization functionality
"""
module LPBasic

#constains basic struct definitions(Line, ProblemInstance) and constructors
include("definitions.jl")

#constains some utility helper functions
include("utility.jl")

#contains code for validation of instance X solution pairs
include("verifier.jl")

#contains code for reading data from a simplified GTFS feed and loading it into a ProblemInstance
include("DataReader.jl")

#contains code for instance/solution visualization
include("visualization.jl")

#contains code resposible for the optimization models
include("optimization.jl")

end


text:
export
    Query,
    Queryable,
    Marginal,
    Joint,
    ProbValue,
    ProbFunction,
    ProbabilityBounds,
    Expectation,
    Mean,
    Variance,
    answer,
    marginal,
    joint,
    probability,
    probability_bounds,
    expectation,
    mean

"""
    abstract type Query end

General type of query that can be answered after running an algorithm.
"""
abstract type Query end

"""
A query target is either a variable instance or a variable.
Allowing queries to be defined in terms of instances rather than variables makes it possible
to ask queries across multiple instances of a variable at different times.
However, in many cases the current instance of the variable(s) is required and then it is easier
to use variables.
"""
Queryable{O} = Union{VariableInstance{O}, Variable{I,J,O}} where {I,J}

function _resolve(r::Runtime, q::Queryable{O})::VariableInstance{O} where O
    if q isa Instance
        return q
    else
        return current_instance(r, q)
    end
end

"""
    answer(::Query, ::Algorithm, ::Runtime, ::VariableInstance)
    answer(::Query, ::Algorithm, ::Runtime, ::Vector{VariableInstance})
    answer(::Query, ::Algorithm, ::Runtime, ::Queryable)
    answer(::Query, ::Algorithm, ::Runtime, ::Vector{Queryable})

    Answer the query.

    An implementation of an algorithm should implement an `answer` method for any queries
    it can handle. The type hierarchies of `Query` and `Algorithm` will enable
    query answering methods to be used wherever appropriate with the right specialization.
    The implementations of `answer` are differentiated along two dimensions:
    - single or multiple items
    - queryable items in general or specifically instances

    It is expected that an algorithm will implement one of the first two methods for queries it
    can handle. I.e., an algorithm is expected to handle a single instance or a vector of instances.
    If it can handle multiple instances, it should implement a second method and a single instance implementation
    is derived by default using a singleton vector. An algorithm can still override this default
    method if it handles single instances differently from multiple.

    Algorithms will generally not implement the latter two methods, which are provide for convenience.
    Default implementations are provided that delegate to the instance-specific methods.

    Defining a very high-level default implementation that throws an error enables implementations
    to go through sequences of preferences.
"""
function answer(q::Query, a::Algorithm, r::Runtime, i::VariableInstance)
    is = VariableInstance[i]
    answer(q, a, r, is)
end

function answer(q::Query, a::Algorithm, r::Runtime, item::Queryable)
    answer(q, a, r, _resolve(r,item))
end

function answer(q::Query, a::Algorithm, r::Runtime, items::Vector{Queryable})
    insts = VariableInstance[]
    for i in items
        inst::VariableInstance = _resolve(r,i)
        push!(insts, inst)
    end
    answer(q, a, r, insts)
end

answer(::Any, ::Any, ::Any, ::Any) = error("_answer_")

""

struct Marginal <: Query end

"""
    marginal(alg::Algorithm, runtime::Runtime, item::Queryable{O})::Union{Dist{O}, Tuple{Dist{O}, Dist{O}}} where O

Return the marginal distribution over `item`, or lower and upper marginals,
depending on the algorithm.

The returned `Score` assigns a score to each value of `item`.
"""
function marginal(alg::Algorithm, run::Runtime, item::Queryable{O})::Union{Dist{O}, Tuple{Dist{O}, Dist{O}}} where O
    mg = answer(Marginal(), alg, run, item)
    return mg
end

struct Joint <: Query end

"""
    joint(alg::Algorithm, run::Runtime, items::Vector{Queryable})::Union{Score{O}, Tuple{Score{O}, Score{O}}}

Return the joint distribution over `items`, or lower and upper distributions,
depending on the algorithm.

The returned `Score` assigns a score for each Vector of values of the items.
"""
function joint(alg::Algorithm, run::Runtime, items::Vector{Queryable})::Union{Score, Tuple{Score, Score}}
    answer(Joint(), alg, run, items)
end

struct ProbValue{O} <: Query
    value :: O
end

struct ProbFunction <: Query
    fn :: Function
end

"""
    probability(alg::Algorithm, run::Runtime, items::Vector{Queryable}, predicate::Function)::Union{Float64, Tuple{Float64, Float64}}

Return the probability that `items` satisfy `query` or lower and upper probabilities.

`predicate` is a function from Vector{Any} to `Bool`.
"""
function probability(alg::Algorithm, run::Runtime, items::Vector{Queryable}, predicate::Function)::Union{Float64, Tuple{Float64, Float64}}
    answer(ProbFunction(predicate), alg, run, items)
end

function probability(alg::Algorithm, run::Runtime, item::Queryable, predicate::Function)::Union{Float64, Tuple{Float64, Float64}}
    f(vec) = predicate(vec[1])
    insts = VariableInstance[_resolve(run,item)]
    answer(ProbFunction(f), alg, run, insts)
end

"""
    probability(alg::Algorithm, run::Runtime, item::Queryable{O}, value::O)::Union{Float64, Tuple{Float64, Float64}} where O

Return the probability that `item` has `value` or lower and upper probabilities.

The default implementation tries to use the more general probability of a query.
If that fails, it uses the `cpdf` operation on the marginal of `item`.
"""
function probability(alg::Algorithm, run::Runtime, item::Queryable{O}, value::O)::Union{Float64, Tuple{Float64, Float64}} where O
    inst::Instance = _resolve(run, item)
    try
        answer(ProbValue(value), alg, run, inst)
    catch e
        if e == ErrorException("_answer_")
            try
                probability(alg, run, item, x -> x == value)
            catch e
                if e == ErrorException("_answer_")
                    try
                        m = marginal(alg, run, inst)
                        return cpdf(m, (), value)
                    catch e
                        if e == ErrorException("_answer_")
                            error("None of the methods to compute the probability of a value are implemented")
                        else
                            rethrow(e)
                        end
                    end
                else
                    rethrow(e)
                end
            end
        else
            rethrow(e)
        end
    end
end

struct ProbabilityBounds{O} <: Query
    range :: Vector{O}
end

"""
    probability_bounds(alg::Algorithm, run::Runtime, item::Queryable, range::Vector)::Tuple{Vector{Float64}, Vector{Float64}}

    For an algorithm that produces lower and upper bounds, return vectors of lower and upper bounds on probabilities for values in the range.

    The range is important for computing the bounds, because it is assumed that values outside the range have probability zero.
"""
function probability_bounds(alg::Algorithm, run::Runtime, item::Queryable, range::Vector)::Tuple{Vector{Float64}, Vector{Float64}}
    return answer(ProbabilityBounds(range), alg, run, item)
end

struct Expectation <: Query
    fn :: Function
end

struct Mean <: Query end

"""
    expectation(alg::Algorithm, runtime::Runtime, item::Queryable, f::Function)::Float64

Return the expectation of the function `f` over the marginal distribution of `item`.

The default implementation uses the expectation operation on the SFunc representing the
marginal distribution.
"""
function expectation(alg::Algorithm, run::Runtime, item::Queryable, fn::Function)::Float64
    try
        return answer(Expectation(fn), alg, run, item)
    catch e
        if e == ErrorException("_answer_")
            m = marginal(alg, run, item)
            if m isa Tuple
                # marginal produced bounds; support is same for lower and upper bounds
                m = m[1]
            end
            tot = 0.0
            O = output_type(item)
            sup = support(m, (), 1000, O[])
            for x in unique(sup)
                tot += cpdf(m, (), x) * fn(x)
            end
            return tot
            # return Operators.f_expectation(marginal(alg, run, item), (), fn)
        else
            rethrow(e)
        end
    end
end

"""
    mean(alg::Algorithm, runtime::Runtime, item::Queryable)

Return the mean of `item`.
"""
function mean(alg::Algorithm, runtime::Runtime, item::Queryable{O})::Float64 where {O <: Number}
    try
        answer(Mean(), alg, runtime, item)
    catch e
        if e == ErrorException("_answer_")
            return expectation(alg, runtime, item, x -> x)
        else
            rethrow(e)
        end
    end
        return
end

struct Variance <: Query end

"""
    variance(alg::Algorithm, runtime::Runtime, item::Queryable)::Float64

Return the variance of `item`.
"""
function variance(alg::Algorithm, run::Runtime, item::Queryable)::Float64
    #answer(Variance(), alg, run, item)
    try
        return answer(Variance(), alg, run, item)
    catch e
        if e == ErrorException("_answer_")
            exp_x = expectation(alg, run, item, x->x)
            exp_xsq = expectation(alg, run, item, x->x^2)
            variance_x = exp_xsq - exp_x^2
            #return Operators.f_expectation(marginal(alg, run, item), (), fn)
            return variance_x
        else
            rethrow(e)
        end
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module APITests
import ..Pkg # ensure we are using the correct Pkg

using Pkg, Test, REPL
import Pkg.Types.PkgError, Pkg.Resolve.ResolverError
using UUIDs

using ..Utils

@testset "Pkg.activate" begin
    temp_pkg_dir() do project_path
        cd_tempdir() do tmp
            path = pwd()
            Pkg.activate(".")
            mkdir("Foo")
            cd(mkdir("modules")) do
                Pkg.generate("Foo")
            end
            Pkg.develop(Pkg.PackageSpec(path="modules/Foo")) # to avoid issue #542
            Pkg.activate("Foo") # activate path Foo over deps Foo
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            rm("Foo"; force=true, recursive=true)
            Pkg.activate("Foo") # activate path from developed Foo
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("./Foo") # activate empty directory Foo (sidestep the developed Foo)
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("Bar") # activate empty directory Bar
            @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
            Pkg.activate(".")
            Pkg.add("Example") # non-deved deps should not be activated
            Pkg.activate("Example")
            @test Base.active_project() == joinpath(path, "Example", "Project.toml")
            Pkg.activate(".")
            cd(mkdir("tests"))
            Pkg.activate("Foo") # activate developed Foo from another directory
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate() # activate home project
            @test Base.ACTIVE_PROJECT[] === nothing
        end
    end
end

@testset "Pkg.status" begin
    temp_pkg_dir() do project_path
        Pkg.add(PackageSpec(name="Example", version="0.5.1"))
        Pkg.add("Random")
        Pkg.status()
        Pkg.status("Example")
        Pkg.status(["Example", "Random"])
        Pkg.status(PackageSpec("Example"))
        Pkg.status(PackageSpec(uuid = "7876af07-990d-54b4-ab0e-23690620f79a"))
        Pkg.status(PackageSpec.(["Example", "Random"]))
        Pkg.status(; mode=PKGMODE_MANIFEST)
        Pkg.status("Example"; mode=PKGMODE_MANIFEST)
        @test_deprecated Pkg.status(PKGMODE_MANIFEST)
        # issue #1183: Test exist in manifest but not in project
        Pkg.status("Test"; mode=PKGMODE_MANIFEST)
        @test_throws PkgError Pkg.status("Test"; mode=Pkg.Types.PKGMODE_COMBINED)
        @test_throws PkgError Pkg.status("Test"; mode=PKGMODE_PROJECT)
        # diff option
        @test_logs (:warn, r"diff option only available") Pkg.status(diff=true)
        git_init_and_commit(project_path)
        @test_logs () Pkg.status(diff=true)
    end
end

@testset "Pkg.test" begin
    temp_pkg_dir() do tmp
        copy_test_package(tmp, "TestArguments")
        Pkg.activate(joinpath(tmp, "TestArguments"))
        # test the old code path (no test/Project.toml)
        Pkg.test("TestArguments"; test_args=`a b`, julia_args=`--quiet --check-bounds=no`)
        Pkg.test("TestArguments"; test_args=["a", "b"], julia_args=["--quiet", "--check-bounds=no"])
        # test new code path
        touch(joinpath(tmp, "TestArguments", "test", "Project.toml"))
        Pkg.test("TestArguments"; test_args=`a b`, julia_args=`--quiet --check-bounds=no`)
        Pkg.test("TestArguments"; test_args=["a", "b"], julia_args=["--quiet", "--check-bounds=no"])
    end
end

include("FakeTerminals.jl")
import .FakeTerminals.FakeTerminal

@testset "Pkg.precompile_script" begin
    function fake_repl(@nospecialize(f); options::REPL.Options=REPL.Options(confirm_exit=false))
        # Use pipes so we can easily do blocking reads
        # In the future if we want we can add a test that the right object
        # gets displayed by intercepting the display
        input = Pipe()
        output = Pipe()
        err = Pipe()
        Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
        Base.link_pipe!(output, reader_supports_async=true, writer_supports_async=true)
        Base.link_pipe!(err, reader_supports_async=true, writer_supports_async=true)

        repl = REPL.LineEditREPL(FakeTerminal(input.out, output.in, err.in), true)
        repl.options = options

        f(input.in, output.out, repl)
        t = @async begin
            close(input.in)
            close(output.in)
            close(err.in)
        end
        @test read(err.out, String) == ""
        #display(read(output.out, String))
        Base.wait(t)
        nothing
    end

    fake_repl() do stdin_write, stdout_read, repl
        repltask = @async REPL.run_repl(repl)

        for line in split(Pkg.precompile_script, "\n"; keepempty=false)
            sleep(0.1)
            # Consume any extra output
            if bytesavailable(stdout_read) > 0
                copyback = readavailable(stdout_read)
                #@info(copyback)
            end

            # Write the line
            write(stdin_write, line, "\n")

            # Read until some kind of prompt
            readuntil(stdout_read, "\n")
            readuntil(stdout_read, ">")
            #@info(line)
        end

        write(stdin_write, "\x04")
        wait(repltask)
    end
end

end # module APITests


text:
module TestDoctest

import LazyGroupBy
using Documenter: DocMeta, doctest
using Test

DocMeta.setdocmeta!(
    LazyGroupBy,
    :DocTestSetup,
    :(using LazyGroupBy: foldxl, foldxt, foldxd),
)

@testset "doctest" begin
    doctest(LazyGroupBy)
end

end  # module


text: println("Building the Beta Machine Learning Toolkit package (BetaML.jl)...")


text:
using CUDA
using CUDA: cufunc, curand_rng
using CUDA.CUSPARSE: CuSparseMatrix, CuSparseMatrixCSR, CuSparseMatrixCOO
using CUDA.CUSOLVER: CuQR

export cuda_gc, gpu, @gpu!, @cu!

const CuBaseField{B,M,T,A<:CuArray} = BaseField{B,M,T,A}

typealias(::Type{CuArray{T,N}}) where {T,N} = "CuArray{$T,$N}"

# a function version of @cuda which can be referenced before CUDA is
# loaded as long as it exists by run-time (unlike the macro @cuda which must
# exist at compile-time)
function cuda(f, args...; threads=256)
    @cuda threads=threads f(args...)
end

is_gpu_backed(::BaseField{B,M,T,A}) where {B,M,T,A<:CuArray} = true
global_rng_for(::Type{<:CuArray}) = curand_rng()

# handy conversion functions and macros
@doc doc"""

    @gpu! x y

Equivalent to `x = gpu(x)`, `y = gpu(y)`, etc... for any number of
listed variables. See [`gpu`](@ref).
"""
macro gpu!(vars...)
    :(begin; $((:($(esc(var)) = gpu($(esc(var)))) for var in vars)...); nothing; end)
end
@doc doc"""

    gpu(x)

Recursively moves x to GPU, but unlike `CUDA.cu`, doesn't also convert
to Float32. Equivalent to `adapt_structure(CuArray, x)`. Returns nothing.
"""
gpu(x) = adapt_structure(CuArray, x)

@doc doc"""

    @cu! x y

Equivalent to `x = cu(x)`, `y = cu(y)`, etc... for any number of
listed variables. See `CUDA.cu`. Returns nothing.
"""
macro cu!(vars...)
    :(begin; $((:($(esc(var)) = cu($(esc(var)))) for var in vars)...); nothing; end)
end

adapt_structure(::CUDA.Float32Adaptor, proj::ProjLambert) = adapt_structure(CuArray{Float32}, proj)

function Cℓ_to_2D(Cℓ, proj::ProjLambert{T,<:CuArray}) where {T}
    # todo: remove needing to go through cpu here:
    gpu(Complex{T}.(nan2zero.(Cℓ.(cpu(proj.ℓmag)))))
end

### misc
# the generic versions of these trigger scalar indexing of CUDA, so provide
# specialized versions:
pinv(D::Diagonal{T,<:CuBaseField}) where {T} = Diagonal(@. ifelse(isfinite(inv(D.diag)), inv(D.diag), $zero(T)))
inv(D::Diagonal{T,<:CuBaseField}) where {T} = any(Array((D.diag.==0)[:])) ? throw(SingularException(-1)) : Diagonal(inv.(D.diag))
fill!(f::CuBaseField, x) = (fill!(f.arr,x); f)
==(a::CuBaseField, b::CuBaseField) = (==)(promote(a.arr, b.arr)...)
sum(f::CuBaseField; dims=:) = (dims == :) ? sum(f.arr) : (1 in dims) ? error("Sum over invalid dims of CuFlatS0.") : f

# these only work for Reals in CUDA
# with these definitions, they work for Complex as well
CUDA.isfinite(x::Complex) = Base.isfinite(x)
CUDA.sqrt(x::Complex) = CUDA.sqrt(CUDA.abs(x)) * CUDA.exp(im*CUDA.angle(x)/2)
CUDA.culiteral_pow(::typeof(^), x::Complex, ::Val{2}) = x * x
CUDA.pow(x::Complex, p) = x^p

# until https://github.com/JuliaGPU/CUDA.jl/pull/618 (CUDA 2.5)
CUDA.cufunc(::typeof(angle)) = CUDA.angle

# adapting of SparseMatrixCSC ↔ CuSparseMatrixCSR (otherwise dense arrays created)
adapt_structure(::Type{<:CuArray}, L::SparseMatrixCSC)   = CuSparseMatrixCSR(L)
adapt_structure(::Type{<:Array},   L::CuSparseMatrixCSR) = SparseMatrixCSC(L)
adapt_structure(::Type{<:CuArray}, L::CuSparseMatrixCSR) = L
adapt_structure(::Type{<:Array},   L::SparseMatrixCSC)   = L

# CUDA somehow missing this one
# see https://github.com/JuliaGPU/CuArrays.jl/issues/103
# and https://github.com/JuliaGPU/CuArrays.jl/pull/580
ldiv!(qr::CuQR, x::CuVector) = qr.R \ (CuMatrix(qr.Q)' * x)

# some Random API which CUDA doesn't implement yet
Random.randn(rng::CUDA.CURAND.RNG, T::Random.BitFloatType) =
    cpu(randn!(rng, CuVector{T}(undef,1)))[1]

# perhaps minor type-piracy, but this lets us simulate into a CuArray using the
# CPU random number generator
Random.randn!(rng::MersenneTwister, A::CuArray) =
    (A .= adapt(CuArray, randn!(rng, adapt(Array, A))))

# CUDA makes some copies here as a workaround for JuliaGPU/CuArrays.jl#345 &
# NVIDIA/cuFFT#2714055 but it doesn't appear to be needed in the R2C case, and
# in the C2R case we pre-allocate the memory only once (via memoization) as well
# as doing the copy asynchronously
import CUDA.CUFFT: unsafe_execute!
using CUDA.CUFFT: rCuFFTPlan, cufftReal, cufftComplex, CUFFT_R2C, cufftExecR2C, cufftExecC2R, CUFFT_C2R, unsafe_copyto!, CuDefaultStream, pointer

plan_buffer(x) = plan_buffer(eltype(x),size(x))
@memoize plan_buffer(T,dims) = CuArray{T}(undef,dims...)

function unsafe_execute!(plan::rCuFFTPlan{cufftReal,K,false,N},
                            x::CuArray{cufftReal,N}, y::CuArray{cufftComplex,N}
                            ) where {K,N}
    @assert plan.xtype == CUFFT_R2C
    cufftExecR2C(plan, x, y)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftComplex,K,false,N},
                            x::CuArray{cufftComplex,N}, y::CuArray{cufftReal}
                            ) where {K,N}
    @assert plan.xtype == CUFFT_C2R
    cufftExecC2R(plan, unsafe_copyto!(pointer(plan_buffer(x)),pointer(x),length(x),async=true,stream=CuDefaultStream()), y)
end

# monkey-patched version of https://github.com/JuliaGPU/CUDA.jl/pull/436
# until it hits a release
using CUDA.CURAND: curandSetPseudoRandomGeneratorSeed, curandSetGeneratorOffset,
    CURAND_STATUS_ALLOCATION_FAILED, CURAND_STATUS_PREEXISTING_FAILURE, CURAND_STATUS_SUCCESS,
    unsafe_curandGenerateSeeds, throw_api_error, @retry_reclaim, RNG

function Random.seed!(rng::RNG, seed=Base.rand(UInt64), offset=0)
    curandSetPseudoRandomGeneratorSeed(rng, seed)
    curandSetGeneratorOffset(rng, offset)
    res = @retry_reclaim err->isequal(err, CURAND_STATUS_ALLOCATION_FAILED) ||
                              isequal(err, CURAND_STATUS_PREEXISTING_FAILURE) begin
        unsafe_curandGenerateSeeds(rng)
    end
    if res != CURAND_STATUS_SUCCESS
        throw_api_error(res)
    end
    return
end

"""
    cuda_gc()

Gargbage collect and reclaim GPU memory (technically should never be
needed to do this by hand, but sometimes helps with GPU OOM errors)
"""
function cuda_gc()
    GC.gc(true)
    CUDA.reclaim()
end


text:
export simulate, forecast

mutable struct Forecast{T <: AbstractFloat}
    observation_quantiles::Matrix{T}
    observation_forecast::Vector{T}
    observation_scenarios::Matrix{T}
    parameter_quantiles::Array{T, 3}
    parameter_forecast::Matrix{T}
    parameter_scenarios::Array{T, 3}
end

"""
    simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int, kwargs...) where {D, T}

Generate scenarios for the future of a time series by updating the GAS recursion `H` times and taking
a sample of the distribution until it generates `S` scenarios.

By default this method uses the `stationary_initial_params` method to perform the
score driven recursion. If you estimated the model with a different set of `initial_params`
use them here to maintain the coherence of your estimation.
"""
function simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int;
                    initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
    # Filter params estimated on the time series
    params = score_driven_recursion(gas, series; initial_params = initial_params)
    biggest_lag = number_of_lags(gas)
    n_params = num_params(D)
    params_simulation = params[(end - biggest_lag + 1):end, :]
    # Create scenarios matrix
    observation_scenarios = Matrix{T}(undef, H, S)
    parameter_scenarios = Array{T, 3}(undef, H, n_params, S)
    for scenario in 1:S
        # Notice that we know the parameter time_varying parameters for T + 1
        # So the last initial_params is already a part of the future simulation
        # And we must take the
        sim, param = simulate_recursion(gas, H + biggest_lag; initial_params = params_simulation)
        observation_scenarios[:, scenario] = sim[biggest_lag+1:end]
        # The first param is known
        parameter_scenarios[1, :, scenario] = params_simulation[end, :]
        # The last param (param[end, :]) is actually in a time step bigger
        # than H
        parameter_scenarios[2:end, :, scenario] = param[biggest_lag+1:end-1, :]
    end
    return observation_scenarios, parameter_scenarios
end

"""
    forecast(series::Vector{T}, gas::Model{D, T}, H::Int; kwargs...) where {D, T}

Forecast quantiles for future values of a time series by updating the GAS recursion `H` times and
using Monte Carlo method as in Blasques, Francisco, Siem Jan Koopman,
Katarzyna Lasak and Andre Lucas (2016): "In-Sample Confidence Bounds and Out-of-Sample Forecast
Bands for Time-Varying Parameters in Observation Driven Models",
International Journal of Forecasting, 32(3), 875-887.

You can pass the desired quantiles as a `Vector{T}`. The default behavior is to
forecast the median and the `0.025` and `0.975` quantiles.

By default 1000 scenarios are used but you can change this by switching the `S` keyword argument.

By default this method uses the `stationary_initial_params` method to perform the
score driven recursion. If you estimated the model with a different set of `initial_params`
use them here to maintain the coherence of your estimation.
"""
function forecast(series::Vector{T}, gas::Model{D, T}, H::Int;
                    initial_params::Matrix{T} = stationary_initial_params(gas),
                    quantiles::Vector{T} = T.([0.025, 0.5, 0.975]), S::Int = 10_000) where {D, T}

    observation_scenarios, parameter_scenarios = simulate(series, gas, H, S;
                                                    initial_params = initial_params)
    return Forecast(
        get_quantiles(quantiles, observation_scenarios),
        mean(observation_scenarios, dims = 2)[:],
        observation_scenarios,
        get_quantiles(quantiles, parameter_scenarios),
        mean(parameter_scenarios, dims = 3)[:, :, 1],
        parameter_scenarios
        )
end

function get_quantiles(quantiles_probs::Vector{T}, scenarios::Matrix{T}) where T
    @assert all((quantiles_probs .< 1.0) .& (quantiles_probs .> 0.0))
    unique!(sort!(quantiles_probs))
    quantiles = mapslices(x -> quantile(x, quantiles_probs), scenarios; dims = 2)
    return quantiles
end

function get_quantiles(quantiles_probs::Vector{T}, scenarios::Array{T, 3}) where T
    @assert all((quantiles_probs .< 1.0) .& (quantiles_probs .> 0.0))
    unique!(sort!(quantiles_probs))
    quantiles_per_parameter = Array{T, 3}(undef,
                                size(scenarios, 1),
                                size(scenarios, 2),
                                length(quantiles_probs)
                                )
    for j in 1:size(scenarios, 2)
        quantiles_per_parameter[:, j, :] = mapslices(x -> quantile(x, quantiles_probs),
                                                        scenarios[:, j, :]; dims = 2)
    end
    return quantiles_per_parameter
end


text:
import Pkg

Pkg.add("ODBC")
Pkg.add("DataFrames")

using ODBC,DataFrames

home = ENV["HOME"]
ODBC.setunixODBC(;ODBCSYSINI=home, ODBCINSTINI=".odbcinst.ini", ODBCINI=string(home, "/.odbc.ini"))

conn=ODBC.Connection("DSN=DuckDB;")

duckdb_dir = pwd()
DBInterface.execute(conn,"CREATE TABLE test AS SELECT * FROM read_csv_auto('" * duckdb_dir * "/tools/odbc/test/test.csv')")

df1 = DBInterface.execute(conn,"SELECT sum(a) as a, sum(b) as b, sum(c) as c FROM test")|>DataFrame
df2 = DataFrame(a = Float64(3.3), b = Float64(36.0), c = Float64(1.0));

if df1 != df2
    exit(1)
end

println("Julia successful \\o/")

max_stars_repo_name: AldoMyrtaj/duckdb

max_issues_repo_name: AldoMyrtaj/duckdb

max_forks_repo_name: AldoMyrtaj/duckdb


text:
using Base.Test
using Distributions
using HiddenMarkovModel

include("discrete_tests.jl")

include("normal_tests.jl")


text:
"""
    rf = downscale_spaceonly(r,f,weight=1.;fglob=false, fsmooth=false )

Downscale precipitation field `r` using metagaussian spectral field `f`. An optional weights array can be specified. Precipitation can be conserved globally (`fglob`) or using convolution (`fsmooth`).
"""
function downscale_spaceonly(r,f,weight=1.;fglob=false, fsmooth=false )

(nas,nas)=size(r);
(ns,ns)=size(f);

# No negative rain
@compat ii=findall(r[:].<=0);
r[ii].=0;
rs=r;

# Gaussianize
# Alternative: rg=log(rs);
rg=gaussianize(rs);

# Prepare small scale random fields
g=metagauss(f);

# fix the amplitude of f
(pstg,pstgt)=fft3d(g);
pstg=pstg*((ns*ns)^2);
(pstrg,pstrgt)=fft3d(rg);
pstrg=pstrg*((nas*nas)^2);

varg=pstg[div(nas,2)];
varrg=pstrg[div(nas,2)];

c=varrg/varg;
g=g*sqrt(c);

gm=mergespec_spaceonly(rg,g,div(nas,2));

st=Statistics.std(gm);
if st==0.
   st=1.;
end
gm=gm/st;

fm=exp.(gm);

fm=fm.*weight;

# We want the aggregated field to be the same as pa
@compat ii=findall(isnan.(r)); fm[ii].=NaN;

if(fglob)
@compat  ii=findall(.~isnan.(fm));
    ri=interpola(r,ns,1);
    fm=fm*Statistics.mean(ri[ii])/Statistics.mean(fm[ii])
elseif(fsmooth)
@compat    ii=findall(isnan.(fm));
    fmi=interpola(agg(fm,nas,1),ns,1);
    fmi[ii].=NaN;
    fma=smoothconv(fmi,nas);
 # fma=smoothconv(fm,nas);
    ri=interpola(r,ns,1);
    ri[ii].=NaN;
    raa=smoothconv(ri,nas);
    fm=raa./fma.*fm; #pa=p aggregato a L0 e T0;
else
#fm=mergespec_spaceonly(fm,raa,div(nas,2));
    raa=agg(r,nas,1);
    fma=agg(fm,nas,1);
    ca=raa./fma; #pa=p aggregato a L0 e T0;
    cai=interpola(ca,ns,1);
    fm=cai.*fm;
end

return fm
end


text:
include("../../util/src/util_module.jl")

module optitestfun
using ..util
using PyPlot
PyPlot.pygui(true)

function nd_quadratic(x)
    y = x' * x
    return y
end

function quad_with_norm_l1(x, x_sol, lamda = 0.1)
    dx = x - x_sol
    y = dx' * dx + lamda * sum( abs.(x) )
    return y
end

function ackley(x, a = 20, b = 0.2, c = 2*pi )
    d = length(x)
    e1 = -b * sqrt( 1/d * sum(x.^2) )
    e2 = 1/d * sum( cos.(c * x) )
    y = - a * exp(e1) - exp(e2) + a + exp(1.0)
    return y
end

function wheeler_ridge(x, a = 1.5)
	t = (x[1] * x[2] - a)^2 + (x[2] - a)^2;
	y = -exp(-t);
end

function booth(x)
# fcn: R^2 -> R
    y = (x[1] + 2 * x[2] - 7)^2 + (2 * x[1] + x[2] - 5)^2
end

function branin(x, a=1.0, b = 5.1/(4*pi^2), c = 5/pi, r = 6.0, s = 10.0, t = 1/(8*pi))
    y = a * (x[2] - b * x[1]^2 + c * x[1] - r)^2 + s * (1-t) * cos(x[1]) + s
    return y
end

function rosenbrock(x, a = 1, b = 5)
    y = (a - x[1])^2 + b * (x[2] - x[1]^2)^2
    return y
end

function michalewicz(x, dim = 2, m = 10.0)
	y = 0.0
	for dd = dim:-1:1
		xi = x[dd];
		t = sin(xi) * (sin(dd * xi^2 / pi))^(2*m);
		y = y + t;
	end
	y = -y;
end

function flower(x, a = 1.0, b = 1.0, c = 4.0)
	y = a * norm(x) + b * sin( c * atan(x[2], x[1]) );
end

function quadratic(x)
# fcn: R^2 -> R
	y = 0.8 * x[1]^2 + 1.2 * x[2]^2
end

function plot_ackley()

    x_vec = collect(-5:0.1:5)
    y_vec = collect(-6:0.1:6)
    Lx = length(x_vec)
    Ly = length(y_vec)

    (x_mat, y_mat) = util.mesh_grid(x_vec, y_vec)

    f = zeros(Lx, Ly)
    for ii = 1:Lx
        for jj = 1:Ly
            f[ii, jj] = ackley( [x_vec[ii], y_vec[jj]] )
        end
    end

    PyPlot.figure()
    PyPlot.surf(x_mat, y_mat, f)
end

function plot_booth()
    x_vec = collect(-10:0.1:10)
    y_vec = collect(-10:0.1:10)
    Lx = length(x_vec)
    Ly = length(y_vec)

    (x_mat, y_mat) = util.mesh_grid(x_vec, y_vec)

    f = zeros(Lx, Ly)
    for ii = 1:Lx
        for jj = 1:Ly
            f[ii, jj] = booth( [x_vec[ii], y_vec[jj]] )
        end
    end

    PyPlot.figure()
    PyPlot.contour(x_mat, y_mat, f)
end

function plot_branin()
    x_vec = collect(-6:0.1:20)
    y_vec = collect(-5:0.1:22)
    Lx = length(x_vec)
    Ly = length(y_vec)

    (x_mat, y_mat) = util.mesh_grid(x_vec, y_vec)

    f = zeros(Lx, Ly)
    for ii = 1:Lx
        for jj = 1:Ly
            f[ii, jj] = branin( [x_vec[ii], y_vec[jj]] )
        end
    end

    PyPlot.figure()
    PyPlot.contour(x_mat, y_mat, f, [1e-3, 1e-2, 1e-1, 1, 10, 20, 30, 50, 75, 100, 200, 300, 400, 500, 1000]) #, colors=PyPlot.ColorMap

    # PyPlot.figure()
    # PyPlot.plot(f)

end

function plot_rosenbrock(x_domain = [-2, 2], y_domain = [-2, 2])

    x_vec = collect(x_domain[1]:0.1:x_domain[2])
    y_vec = collect(y_domain[1]:0.1:y_domain[2])
    Lx = length(x_vec)
    Ly = length(y_vec)

    (x_mat, y_mat) = util.mesh_grid(x_vec, y_vec)

    f = zeros(Lx, Ly)
    for ii = 1:Lx
        for jj = 1:Ly
            f[ii, jj] = rosenbrock( [x_vec[ii], y_vec[jj]] )
        end
    end

    PyPlot.figure()
    PyPlot.contour(x_mat, y_mat, f, [1e-3, 1e-2, 1e-1, 1, 10, 20, 30, 50, 75, 100, 200, 300, 400, 500, 1000]) #, colors=PyPlot.ColorMap
    # PyPlot.contour(x_mat, y_mat, f) #, colors=PyPlot.ColorMap

end

end # optitestfun

# ================================================================================================ #
# test function visualizations
# ================================================================================================ #
#=

# ------------------------------------------------------------------------------------------------ #
# ackley
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-10.0, 10.0, N_vec) );
y_vec = collect( LinRange(-10.0, 10.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> ackley([x, y]), x_mat, y_mat)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# quadratic
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-10.0, 10.0, N_vec) );
y_vec = collect( LinRange(-10.0, 10.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> quadratic([x, y]), x_mat, y_mat)

min_z = minimum(z_mat)
max_z = maximum(z_mat)
levels = logspace( log10(min_z), log10(max_z), 20 )

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat, levels)

# ------------------------------------------------------------------------------------------------ #
# wheeler_ridge
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-10.0, 25.0, N_vec) );
y_vec = collect( LinRange(-3.0, 6.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> wheeler_ridge([x, y]), x_mat, y_mat)

# min_z = minimum(z_mat)
# max_z = maximum(z_mat)
# levels = logspace( log10(min_z), log10(max_z), 20 )

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# rosenbrock2d
# ------------------------------------------------------------------------------------------------ #
N_vec = 100;
x_vec = collect( LinRange(-2.0, 2.0, N_vec) );
y_vec = collect( LinRange(-2.0, 2.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> rosenbrock2d([x, y], 0.5), x_mat, y_mat)

min_z = minimum(z_mat)
max_z = maximum(z_mat)
levels = logspace( log10(min_z), log10(max_z), 20 )

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat, levels)
# PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# michalewicz
# ------------------------------------------------------------------------------------------------ #
N_vec = 100;
x_vec = collect( LinRange(0.0, 4.0, N_vec) );
y_vec = collect( LinRange(0.0, 4.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> michalewicz([x, y]), x_mat, y_mat)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# flower
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-3.0, 3.0, N_vec) );
y_vec = collect( LinRange(-3.0, 3.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> flower([x, y]), x_mat, y_mat)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# branin
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-6.0, 20.0, N_vec) );
y_vec = collect( LinRange(-6.0, 22.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> branin([x, y]), x_mat, y_mat)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

# ------------------------------------------------------------------------------------------------ #
# booth
# ------------------------------------------------------------------------------------------------ #

N_vec = 100;
x_vec = collect( LinRange(-10.0, 10.0, N_vec) );
y_vec = collect( LinRange(-10.0, 10.0, N_vec) );

(x_mat, y_mat) = meshgrid(x_vec, y_vec);
z_mat = map( (x,y) -> booth([x, y]), x_mat, y_mat)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(x_mat, y_mat, z_mat)

=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#
#  UpdateSystem.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 30/03/2020.
#
#

module UpdateSystem

# Import Julia packages
using LinearAlgebra
using StaticArrays

# Import local program modules

@views function updateSystem!(pos,F,W,t,Δt,D,kT,nParticles)

    pos .+= F[:,1]*(Δt*D/kT) .+ W.*sqrt(2.0*D*Δt)
    fill!(F,@SVector zeros(3))

    return t += Δt

end

export updateSystem!

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
using Requests
import Requests: Response

type GymClient
  remote_base::String
end

function env_list_all(client::GymClient)
  route = "/v1/envs/"
  data = Dict()
  response = get_request(client, route, data)
  all_envs = response["all_envs"]
  return(all_envs)
end

function env_create(client::GymClient, env_id::String)
  route = "/v1/envs/"
  data = Dict("env_id"=>env_id)
  response = post_request(client, route, data)
  instance_id = response["instance_id"]
  return(instance_id)
end

function env_reset(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/reset/"
  data = Dict()
  response = post_request(client, route, data)
  observation = response["observation"]
  return(observation)
end

function env_step(client::GymClient, instance_id::String, action::Number;
                  render::Bool = false)
  route = "/v1/envs/" * instance_id * "/step/"
  data = Dict("instance_id"=>instance_id, "action"=>action, "render"=>render)
  response = post_request(client, route, data)
  observation = response["observation"]
  reward = response["reward"]
  done = response["done"]
  info = response["info"]
  result = Dict("observation"=>observation, "reward"=>reward,
                "done"=>done, "info"=>info)
  return(result)
end

function env_action_space_info(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/action_space/"
  data = Dict()
  response = get_request(client, route, data)
  info = response["info"]
  return(info)
end

function env_action_space_sample(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/action_space/sample"
  data = Dict()
  response = get_request(client, route, data)
  action = response["action"]
  return(action)
end

function env_action_space_contains(client::GymClient, instance_id::String,
                                   action::Number)
  route = "/v1/envs/" * instance_id *  "/action_space/contains/"
  route = route * string(action)
  data = Dict()
  response = get_request(client, route, data)
  member = response["member"]
  return(member)
end

function env_observation_space_info(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/observation_space/"
  data = Dict()
  response = get_request(client, route, data)
  info = response["info"]
  return(info)
end

function env_monitor_start(client::GymClient, instance_id::String,
                           directory::String; force::Bool = false,
                           resume::Bool = false)
  route = "/v1/envs/" * instance_id * "/monitor/start/"
  data = Dict("directory"=>directory, "force"=>force,
              "resume"=>resume)
  response = post_request(client, route, data)
  return(nothing)
end

function env_monitor_close(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/monitor/close/"
  data = Dict()
  response = post_request(client, route, data)
  return(nothing)
end

function env_close(client::GymClient, instance_id::String)
  route = "/v1/envs/" * instance_id * "/close/"
  data = Dict()
  response = post_request(client, route, data)
  return(nothing)
end

function upload(client::GymClient, training_dir::String,
                api_key::String, algorithm_id::String)
  # if (is.null(api_key))
  #   api_key = Sys.getenv("OPENAI_GYM_API_KEY")
  # end
  # route = "/v1/upload/"
  # data = Dict("training_dir"=>training_dir, "algorithm_id"=>algorithm_id,
  #             "api_key"=>api_key)
  # response = post_request(client, route, data)
  return(nothing)
end

function shutdown_server(client::GymClient)
  route = "/v1/shutdown/"
  data = Dict()
  response = post_request(client, route, data)
  return(nothing)
end


text:
module IntegralsModule
export computeIntegralOverlap, computeElectronRepulsionIntegral, computeTensorBlockElectronRepulsionIntegrals, computeIntegralKinetic, computeIntegralNuclearAttraction, computeIntegralThreeCenterOverlap, computeMatrixBlockOverlap, computeMatrixBlockKinetic, computeMatrixBlockNuclearAttraction
using ..BaseModule
using ..BasisFunctionsModule
using ..AtomModule
using ..GeometryModule
using ..DocumentationModule
using ..ShellModule
using ProgressMeter

import Base.normalize!

ProgressMeter.printover(STDOUT," + (GSL........................")
using GSL
ProgressMeter.printover(STDOUT," + IntegralsModule.............")

function GaussianIntegral1D_Valeev(mqn::Int,exponent::Float64)
  m = mqn
  ζ = exponent

  if (mqn%2!=0) # integral over uneven function
    return 0
  else
    return (doublefactorial(m-1)*sqrt(π)) / ((2ζ)^(m/2)*sqrt(ζ))
  end
end
@doc """
I_x = Integrate[x^m Exp[-ζ x^2], {x,-∞,∞}] (acc. to Fundament. of Mol. Integr. Eval. by Fermann, Valeev)
""" GaussianIntegral1D_Valeev
@add_doc GenericCitation("Fundament. of Mol. Integr. Eval. by Fermann, Valeev") GaussianIntegral1D_Valeev

"""
I_x = Integrate[x^m Exp[-ζ x^2], {x,-∞,∞}] (acc. to Mathematica 9)
"""
function GaussianIntegral1D_Mathematica(mqn::Int,exponent::Float64)
  m = mqn
  ζ = exponent

  if (mqn%2!=0) # integral over uneven function
    return 0
  else
    t=(m+1)/2
    return ζ^(-t) * gamma(t)
  end
end

GaussianIntegral1D = GaussianIntegral1D_Mathematica

#function FundamentalIntegral(
#  pgb1::PrimitiveGaussianBasisFunction,
#  pgb2::PrimitiveGaussianBasisFunction,
#  func::Function,
#  pgb3::PrimitiveGaussianBasisFunction,
#  pgb4::PrimitiveGaussianBasisFunction)
#  # (pgb1(x1) pgb2(x1) | func(x1,x2) | pgb3(x2) pgb4(x2) )
#  # pgb1..4 are considered s-Functions
#  # acc. to Molecular Integrals of Gaussian Basis Functions by Gill
#
#  # I = (π G_AB G_CD) / ( 2(ζ η)^(3/2) R^3 ) Integrate[u Sin[u] Exp[-u^2/(4T)] FT[f][u/R],{u,0,∞}]
#end

function GaussProductFundamental(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction)
  # Fundamental means primitives are taken as s-Functions
  # K Exp[-γ r^2_P] = pgb1 * pgb2 (acc. to Fundament. of Mol. Integr. Eval. by Fermann, Valeev)

  α1 = pgb1.exponent
  α2 = pgb2.exponent
  A  = pgb1.center
  B  = pgb2.center

  γ = α1 + α2
  P = (α1*A + α2*B) / γ
  AB = distance(A,B)
  K = exp(-α1*α2*AB^2/γ)

  return (K,P,γ)
end

type PolynomialFactors
  x::Array{Tuple{Float64,Int},1} # Float64*x^Int
  y::Array{Tuple{Float64,Int},1} # Float64*y^Int
  z::Array{Tuple{Float64,Int},1} # Float64*z^Int
end

function GaussProductPolynomialFactor(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction)
  # _Sum_K[f_k r^k]_ K Exp[-γ r^2_P] = pgb1 * pgb2 (acc. to Fundament. of Mol. Integr. Eval. by Fermann, Valeev)

  A  = pgb1.center
  B  = pgb2.center
  (K,P,γ) = GaussProductFundamental(pgb1,pgb2)

  #factors = PolynomialFactors(Array(Tuple{Float64,Int},pgb1.mqn.x+pgb2.mqn.x),Array(Tuple{Float64,Int},pgb1.mqn.y+pgb2.mqn.y),Array(Tuple{Float64,Int},pgb1.mqn.z+pgb2.mqn.z))
  factors = PolynomialFactors(Tuple{Float64,Int}[],Tuple{Float64,Int}[],Tuple{Float64,Int}[])
  for xyz in 1:3
    l1 = getfield(pgb1.mqn,xyz)
    l2 = getfield(pgb2.mqn,xyz)
    sizehint!(getfield(factors,xyz),1+l1+l2)
    for k in 0:l1+l2
      f = 0.
      for q in max(-k,k-2l2):2:min(k,2l1-k)
    i = (k+q) ÷ 2
    j = (k-q) ÷ 2
    f += binomial(l1,i) * binomial(l2,j) * getfield((P-A),xyz)^(l1-i) * getfield((P-B),xyz)^(l2-j)
      end
      #ij = [((k+q) ÷ 2,(k-q) ÷ 2) for q in max(-k,k-2l2):2:min(k,2l1-k)]
      #f = sum([binomial(l1,i) * binomial(l2,j) * (P-A).(xyz)^(l1-i) * (P-B).(xyz)^(l2-j) for(i,j) in ij])
      push!(getfield(factors,xyz),(f,k))
    end
  end

  return factors
end

function computeMatrixBlockOverlap(sh1::Shell,sh2::Shell)
  return [computeIntegralOverlap(cgb1,cgb2) for cgb1 in expandShell(sh1), cgb2 in expandShell(sh2)]
end

function computeMatrixBlockKinetic(sh1::Shell,sh2::Shell)
  return [computeIntegralKinetic(cgb1,cgb2) for cgb1 in expandShell(sh1), cgb2 in expandShell(sh2)]
end

function computeMatrixBlockNuclearAttraction(sh1::Shell,sh2::Shell,atom::Atom)
  return [computeIntegralNuclearAttraction(cgb1,cgb2,atom) for cgb1 in expandShell(sh1), cgb2 in expandShell(sh2)]
end

function computeMatrixBlockNuclearAttraction(sh1::Shell,sh2::Shell,geo::Geometry)
  return mapreduce(atom->computeMatrixBlockNuclearAttraction(sh1,sh2,atom),+,0,geo.atoms)
end

function computeIntegralOverlap(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction)
  # S_ij = Integrate[phi_i(r) phi_j(r),{r el R}] (acc. to Fundament. of Mol. Integr. Eval. by Fermann, Valeev)

  (K,P,γ) = GaussProductFundamental(pgb1,pgb2)
  factors = GaussProductPolynomialFactor(pgb1,pgb2)

  Ix = sum([f*GaussianIntegral1D(i,γ) for (f,i) in factors.x])
  Iy = sum([f*GaussianIntegral1D(i,γ) for (f,i) in factors.y])
  Iz = sum([f*GaussianIntegral1D(i,γ) for (f,i) in factors.z])
  return K*Ix*Iy*Iz::Float64
end

function computeIntegralOverlap(
  cgb1::ContractedGaussianBasisFunction,
  cgb2::ContractedGaussianBasisFunction)

  integral = 0.
  for (coeff1,pgb1) in zip(cgb1.coefficients,cgb1.primitiveBFs)
    for (coeff2,pgb2) in zip(cgb2.coefficients,cgb2.primitiveBFs)
      integral += coeff1*coeff2*computeIntegralOverlap(pgb1,pgb2)
    end
  end
  return integral::Float64
end

function computeIntegralThreeCenterOverlap(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction,
  pgb3::PrimitiveGaussianBasisFunction)
  #
  Ix = 0.
  Iy = 0.
  Iz = 0.
  (K12,center,exponent) = IntegralsModule.GaussProductFundamental(pgb1,pgb2)
  factors12 = IntegralsModule.GaussProductPolynomialFactor(pgb1,pgb2)
  #
  mqn = MQuantumNumber(0,0,0)
  pgb12 = PrimitiveGaussianBasisFunction(center,exponent,mqn)
  (K,P,γ) = IntegralsModule.GaussProductFundamental(pgb12,pgb3)
  #
  for xyz in [:x,:y,:z]
    for (f,i) in getfield(factors12, xyz)
      if (xyz == :x)
        mqn = MQuantumNumber(i,0,0)
        pgb12 = PrimitiveGaussianBasisFunction(center,exponent,mqn)
      elseif (xyz == :y)
        mqn = MQuantumNumber(0,i,0)
        pgb12 = PrimitiveGaussianBasisFunction(center,exponent,mqn)
      elseif (xyz == :z)
        mqn = MQuantumNumber(0,0,i)
        pgb12 = PrimitiveGaussianBasisFunction(center,exponent,mqn)
      end
      factors = IntegralsModule.GaussProductPolynomialFactor(pgb12,pgb3)
      if (xyz == :x)
          Ix += sum([f*g*IntegralsModule.GaussianIntegral1D(j,γ) for (g,j) in factors.x])
      elseif (xyz == :y)
          Iy += sum([f*g*IntegralsModule.GaussianIntegral1D(j,γ) for (g,j) in factors.y])
      elseif (xyz == :z)
          Iz += sum([f*g*IntegralsModule.GaussianIntegral1D(j,γ) for (g,j) in factors.z])
      end
    end
  end
  #
  return K12*K*Ix*Iy*Iz::Float64
end

function computeIntegralThreeCenterOverlap(
  cgb1::ContractedGaussianBasisFunction,
  cgb2::ContractedGaussianBasisFunction,
  cgb3::ContractedGaussianBasisFunction)
  #
  integral = 0.
  for (coeff1,pgb1) in zip(cgb1.coefficients,cgb1.primitiveBFs),
      (coeff2,pgb2) in zip(cgb2.coefficients,cgb2.primitiveBFs),
      (coeff3,pgb3) in zip(cgb3.coefficients,cgb3.primitiveBFs)
    integral += coeff1*coeff2*coeff3*computeIntegralThreeCenterOverlap(pgb1,pgb2,pgb3)
  end
  return integral
end

function incrmqn(mqn::MQuantumNumber,xyz::Symbol,inc::Int)
  mqn_t = [mqn.x,mqn.y,mqn.z]
  xyznum = Dict(:x => 1, :y => 2, :z => 3)
  num = xyznum[xyz]
  mqn_t[num] += inc
  return MQuantumNumber(mqn_t[1],mqn_t[2],mqn_t[3])
end

function computeIntegralKinetic(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction)
  #Ix + Iy + Iz = (pgb1 | 1/2 ∇^2 | pgb2)
  #Ix = 1/2 l1 l2 <-1|-1> + 2 α1 α2 <+1|+1> - α1 l2 <+1|-1> - α2 l1 <-1|+1>
  #(acc. to Fundamentals of Mol. Integr. Eval. by Fermann, Valeev (eq. 4.1 + 4.13))

  integral = 0.

  for xyz in (:x,:y,:z)
    pgb1decr = deepcopy(pgb1)
    pgb1decr.mqn = incrmqn(pgb1.mqn,xyz,-1)
    pgb1incr = deepcopy(pgb1)
    pgb1incr.mqn = incrmqn(pgb1.mqn,xyz,1)

    pgb2decr = deepcopy(pgb2)
    pgb2decr.mqn = incrmqn(pgb2.mqn,xyz,-1)
    pgb2incr = deepcopy(pgb2)
    pgb2incr.mqn = incrmqn(pgb2.mqn,xyz,1)

    l1 = getfield(pgb1.mqn, xyz)
    l2 = getfield(pgb2.mqn, xyz)
    α1 = pgb1.exponent
    α2 = pgb2.exponent

    if(l1*l2!=0) integral += (1/2*l1*l2) * computeIntegralOverlap(pgb1decr,pgb2decr) end
                 integral += (2*α1*α2)   * computeIntegralOverlap(pgb1incr,pgb2incr)
    if(l2   !=0) integral += (-α1*l2)    * computeIntegralOverlap(pgb1incr,pgb2decr) end
    if(l1   !=0) integral += (-l1*α2)    * computeIntegralOverlap(pgb1decr,pgb2incr) end
  end

  return integral::Float64
end

function computeIntegralKinetic(
  cgb1::ContractedGaussianBasisFunction,
  cgb2::ContractedGaussianBasisFunction)

  integral = 0.
  for (coeff1,pgb1) in zip(cgb1.coefficients,cgb1.primitiveBFs),
      (coeff2,pgb2) in zip(cgb2.coefficients,cgb2.primitiveBFs)

    integral += coeff1*coeff2*computeIntegralKinetic(pgb1,pgb2)
  end
  return integral
end

function OverlapFundamental(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction)
  (K,P,γ) = GaussProductFundamental(pgb1,pgb2)
  # Overlap(s-type pgb1, s-type pgb2) (acc. to. Mathematica)
  return K*(π/γ)^(3/2)::Float64
end

function FIntegral(m::Integer,x::Real,APPROX_ZERO::Float64=1e-5)
  # F[m,x] = Integrate[u^2m Exp[-x u^2],{u,0,1}
  #         = 1/2 x^(-0.5-m) ( Gamma[1/2 + m] - Gamma[1/2 + m, x] ) # acc. to Mathematica
  x = max(x,APPROX_ZERO) # for x -> 0
  if (x<=APPROX_ZERO) return 1/(1+2*m) end
  return 1/2 * x^(-1/2-m) * ( GSL.sf_gamma(1/2+m) - GSL.sf_gamma_inc(1/2+m,x) )
end

function ASummand(
  f::Real,
  l1::Integer,
  l2::Integer,
  Ax::Real,
  Bx::Real,
  PCx::Real,
  γ::Real,
  l::Integer,
  r::Integer,
  i::Integer)
  # A[l1,l2,Ax,Bx,Cx,γ] = (-1)^l f_l(l1,l2,PAx,PBx) ((-1)^i l! PCx^(l-2r-2i) (1/(4γ))^(r+i)) / (r! i! (l-2r-2i)!)
  return (-1)^l * f * (-1)^i * factorial(l) * PCx^(l-2r-2i) * (1/(4γ))^(r+i) / (factorial(r) * factorial(i) * factorial(l - 2r - 2i))
end

function computeIntegralNuclearAttraction(
  pgb1::PrimitiveGaussianBasisFunction,
  pgb2::PrimitiveGaussianBasisFunction,
  atom::Atom)
  # V = K * A(l1,l2,Ax,Bx,Cx,γ) * A(m1,m2,Ay,By,Cy,γ) * A(n1,n2,Az,Bz,Cz,γ) * F[l+m+n-2(r+s+t) - (i+j+k),γ PC^2]
  #acc. to. Handbook of Comput. Quant. Chem. by Cook, chap. 7.7.3 final formula
  K,P,γ = GaussProductFundamental(pgb1,pgb2)
  C = atom.position

  result = 0
  for (fx,l) in GaussProductPolynomialFactor(pgb1,pgb2).x,
      r in 0:floor(Int,(l/2)),
      i in 0:floor(Int,((l-2r)/2))

    Ax = ASummand(fx,pgb1.mqn.x,pgb2.mqn.x,pgb1.center.x,pgb2.center.x,(P-C).x,γ,l,r,i)
    for (fy,m) in GaussProductPolynomialFactor(pgb1,pgb2).y,
        s in 0:floor(Int,(m/2))m,
        j in 0:floor(Int,((m-2s)/2))

      Ay = ASummand(fy,pgb1.mqn.y,pgb2.mqn.y,pgb1.center.y,pgb2.center.y,(P-C).y,γ,m,s,j)

      for (fz,n) in GaussProductPolynomialFactor(pgb1,pgb2).z,
          t in 0:floor(Int,(n/2)),
          k in 0:floor(Int,((n-2t)/2))

        Az = ASummand(fz,pgb1.mqn.z,pgb2.mqn.z,pgb1.center.z,pgb2.center.z,(P-C).z,γ,n,t,k)
        result += 2π/γ * K * Ax * Ay * Az * FIntegral(l+m+n-2*(r+s+t)-(i+j+k),γ*distance(P,C)^2)
      end
    end
  end

  return -atom.element.atomicNumber * result
end

function computeIntegralNuclearAttraction(
  cgb1::ContractedGaussianBasisFunction,
  cgb2::ContractedGaussianBasisFunction,
  atom::Atom)
  integral = 0.
  for (coeff1,pgb1) in zip(cgb1.coefficients,cgb1.primitiveBFs),
      (coeff2,pgb2) in zip(cgb2.coefficients,cgb2.primitiveBFs)

      integral += coeff1*coeff2*computeIntegralNuclearAttraction(pgb1,pgb2,atom)
  end
  return integral
end

type θfactors
  x::Array{Tuple{Float64,Int,Int},1} # θ,l,r
  y::Array{Tuple{Float64,Int,Int},1} # θ,l,r
  z::Array{Tuple{Float64,Int,Int},1} # θ,l,r
end

function θFactors(
  μ::PrimitiveGaussianBasisFunction,
  ν::PrimitiveGaussianBasisFunction)

  K,P,γ = IntegralsModule.GaussProductFundamental(μ,ν)

  factors = θfactors(Tuple{Float64,Int,Int}[],Tuple{Float64,Int,Int}[],Tuple{Float64,Int,Int}[])
  for xyz in 1:3
    sizehint!(getfield(factors,xyz),floor(Int,(length(getfield(GaussProductPolynomialFactor(μ,ν),xyz))+1)*5/8))

    for (f,l) in getfield(IntegralsModule.GaussProductPolynomialFactor(μ,ν),xyz),
        r in 0:floor(Int,(l/2))

        θ = f * (factorial(l)*γ^(r-l))/(factorial(r)*factorial(l-2r))
        push!(getfield(factors,xyz),(θ,l,r))
    end
  end
  return factors
end

type Bfactors
  x::Array{Tuple{Float64,Int,Int,Int,Int,Int},1}    # B,l12,r12,i,l34,r34
  y::Array{Tuple{Float64,Int,Int,Int,Int,Int},1}    # B,l12,r12,i,l34,r34
  z::Array{Tuple{Float64,Int,Int,Int,Int,Int},1}    # B,l12,r12,i,l34,r34
end

function BFactors(
  μ::PrimitiveGaussianBasisFunction,
  ν::PrimitiveGaussianBasisFunction,
  λ::PrimitiveGaussianBasisFunction,
  σ::PrimitiveGaussianBasisFunction)

  K1,P,γ1 = IntegralsModule.GaussProductFundamental(μ,ν)
  K2,Q,γ2 = IntegralsModule.GaussProductFundamental(λ,σ)
  δ = 1/(4γ1) + 1/(4γ2)
  #p = (P-Q)    # this might be a typo in the book
  p = (Q-P)

  factors = Bfactors([],[],[])
  for xyz in 1:3,
     (θ12,l12,r12) in getfield(θFactors(μ,ν),xyz),
     (θ34,l34,r34) in getfield(θFactors(λ,σ),xyz)
    #for (i in 0:floor(Integer,((l12-2r12)/2)))    # this might be a typo in the book (otherwise e.g. ERI(s,s,px,px2) != ERI(px,px2,s,s) (where the 2 denotes a second center moved by 0.5 in x direction)
    for i in 0:floor(Int,((l12+l34-2r12-2r34)/2))
      B = (-1)^l34 * θ12 * θ34 * ((-1)^i*(2δ)^(2(r12+r34))*factorial(l12+l34-2r12-2r34)*δ^i*getfield(p,xyz)^(l12+l34-2*(r12+r34+i)))/((4δ)^(l12+l34)*factorial(i)*factorial(l12+l34-2*(r12+r34+i)))
      push!(getfield(factors,xyz),(B,l12,r12,i,l34,r34))
    end
  end
  return factors
end

function computeElectronRepulsionIntegral(
  μ::PrimitiveGaussianBasisFunction,
  ν::PrimitiveGaussianBasisFunction,
  λ::PrimitiveGaussianBasisFunction,
  σ::PrimitiveGaussianBasisFunction)
  # compute (μν|λσ) (Mulliken notation) = Integrate[μ(1) ν(1) 1/Abs(1-2) λ(2) σ(2), d1 d2]
  # in the most straightforward but primitive way (acc. to. Handbook of Comp. Quant. Chem. by Cook eq.7.1)
  A = μ.center
  B = ν.center
  C = λ.center
  D = σ.center
  α1 = μ.exponent
  α2 = ν.exponent
  α3 = λ.exponent
  α4 = σ.exponent
  K1,P,γ1 = IntegralsModule.GaussProductFundamental(μ,ν)
  K2,Q,γ2 = IntegralsModule.GaussProductFundamental(λ,σ)
  δ = 1/(4γ1) + 1/(4γ2)
  p = (P-Q)

  Bfactors = BFactors(μ,ν,λ,σ)
  Ω = 2π^2/(γ1 * γ2) * sqrt(π/(γ1 + γ2)) * exp(-α1*α2*distance(A,B)^2/γ1 - α3*α4*distance(C,D)^2/γ2)

  result = 0.
  for (Bx,l12,r12,i,l34,r34) in Bfactors.x,
      (By,m12,s12,j,m34,s34) in Bfactors.y,
      (Bz,n12,t12,k,n34,t34) in Bfactors.z

    V = (l12+l34+m12+m34+n12+n34) - 2*(r12+r34+s12+s34+t12+t34) - (i+j+k)
    result += Ω * Bx * By * Bz * IntegralsModule.FIntegral(V,distance(P,Q)^2/(4δ))
    #println("Bx*By*Bz*FIntegral[$V] = $Bx * $By * $Bz * $(IntegralsModule.FIntegral(V,distance(P,Q)^2/(4δ)))")
  end
  return result
end

function computeElectronRepulsionIntegral(
  μ::ContractedGaussianBasisFunction,
  ν::ContractedGaussianBasisFunction,
  λ::ContractedGaussianBasisFunction,
  σ::ContractedGaussianBasisFunction)

  integral = 0.
  for (coeff1,pgb1) in zip(μ.coefficients,μ.primitiveBFs),
      (coeff2,pgb2) in zip(ν.coefficients,ν.primitiveBFs),
      (coeff3,pgb3) in zip(λ.coefficients,λ.primitiveBFs),
      (coeff4,pgb4) in zip(σ.coefficients,σ.primitiveBFs)

      integral += coeff1*coeff2*coeff3*coeff4*computeElectronRepulsionIntegral(pgb1,pgb2,pgb3,pgb4)
  end
  return integral
end

function computeTensorBlockElectronRepulsionIntegrals(
  μs::Vector{ContractedGaussianBasisFunction},
  νs::Vector{ContractedGaussianBasisFunction},
  λs::Vector{ContractedGaussianBasisFunction},
  σs::Vector{ContractedGaussianBasisFunction})
  [computeElectronRepulsionIntegral(μ,ν,λ,σ) for μ in μs, ν in νs, λ in λs, σ in σs]
end

function normalize!(cgb::ContractedGaussianBasisFunction)
  N = computeIntegralOverlap(cgb,cgb)
  scale!(cgb.coefficients,1/sqrt(N))
end

end # module


text:
"""
Kim, P.J., Chang, H.J., Song, D.S. and Choi, J.Y., 2007, June.
Fast support vector data description using k-means clustering.
In International Symposium on Neural Networks (pp. 506-514). Springer, Berlin, Heidelberg.

We dub the method K-Means subsamples for SVDD (KMSVDD)
Parameters:
on benchmark data k ∈ [5, 10, 20, 30, 40, 50]
init strategy and solver for SVDD solution
"""
struct KMSVDD <: Sampler
    k::Int
    init_strategy::SVDD.InitializationStrategy
    solver::OptimizerFactory
    function KMSVDD(k::Int, init_strategy::SVDD.InitializationStrategy, solver::OptimizerFactory)
        k > 0 || throw(ArgumentError("Invalid number of clusters $(k)."))
        return new(k, init_strategy, solver)
    end
end

function sample(sampler::KMSVDD, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    clustering = kmeans(data, sampler.k)
    model = SVDD.VanillaSVDD(data)
    SVDD.initialize!(model, sampler.init_strategy)
    SVDD.set_adjust_K!(model, true)

    sample_mask = falses(size(data, 2))
    for i in 1:sampler.k
        cluster_indices = findall(clustering.assignments .== i)
        svs = train_and_find_support_vectors(model, data[:, cluster_indices], sampler.solver)
        sample_mask[cluster_indices[svs]] .= true
    end
    sv_mask = findall(sample_mask)
    svs = train_and_find_support_vectors(model, data[:, sv_mask], sampler.solver)
    sample_mask = falses(size(data, 2))
    sample_mask[sv_mask[svs]] .= true
    return sample_mask
end