text:
type WannierParameters
	N1::Int
	N2::Int
	N3::Int
	t1::Array{Float64,1}
	t2::Array{Float64,1}
	t3::Array{Float64,1}
	M::Array{Complex128,4} # nband x nband x Ntot x nntot
	A::Array{Complex128,5} # nband x nwannier x N1 x N2 x N3
	eigs::Array{Float64,2} # Ntot x nband
	neighbors::Array{Int,2} # Ntot x nntot
        kpt_order::Array{Int,1} # permutation such that the order of the kpoints in the file has perm[3] as fastest index
	ijk_to_K :: Array{Int,3}
	K_to_ijk :: Array{Int,2}
	nntot::Int
	nband::Int
        nwannier :: Int
	filename::String
	map::Bool
	logMethod::Bool
        recp_unit_cell :: Array{Float64,2} #3 x 3
        wb :: Float64 #must be inputed manually for now, and be the same for all b vectors
        displ_vecs :: Array{Float64,3} # integer array of displacements between Brillouin zones, 3 x Ntot x nntot
end

# read mmn, win and optionally amn
function read_system(filename, read_amn=true,read_eig=true)
    println("Reading $filename.win")
    win = open("$filename.win")
    nwannier = -1
    nband = -1
    N1,N2,N3 = -1,-1,-1
    unit_cell = zeros(3,3)
    fill!(unit_cell,NaN)
    Ns = (-1,-1,-1)
    kpt_order = [1,2,3]
    while !eof(win)
        line = readline(win)
	line = lowercase(line) # handle case insensitive win files (relic of Fortran)
        line = replace(line, "="," ")
        line = replace(line, ":"," ")
        if(startswith(line,"!"))
            continue
        elseif contains(line,"mp_grid")
            N1,N2,N3 = Base.map(x -> parse(Int,x), split(line)[2:4])
            Ns = [N1,N2,N3]
        elseif contains(line,"num_bands")
            nband = parse(Int, split(line)[2])
        elseif contains(line,"num_wann")
            nwannier = parse(Int, split(line)[2])
        elseif contains(line,"begin unit_cell_cart")
            readline(win) # skip "bohr"
            for i = 1:3
                unit_cell[:,i] = Base.map(x -> parse(Float64,x), split(readline(win))) #weird convention...
            end
        elseif contains(line,"begin kpoints")
            # try to guess the kpoint format
            read_arr() = Base.map(x -> parse(Float64,x), split(readline(win)))
            @assert read_arr() == [0.,0.,0.] #first line
            #second line: detect fastest changing index -> kpt_order[3]
            arr = read_arr()
            nonzero_index = find(x-> x != 0., arr)
            @assert length(nonzero_index) == 1
            @assert isapprox(sum(arr), 1/Ns[nonzero_index[1]],rtol=1e-2)
            kpt_order[3] = nonzero_index[1]

            #read rest of the first part to make sure it's correct
            for i=2:Ns[nonzero_index[1]]-1
                arr = read_arr()
                nonzero_index = find(x-> x != 0., arr)
                @assert length(nonzero_index) == 1
                @assert isapprox(sum(arr), i/Ns[nonzero_index[1]],rtol=1e-2)
            end

            if count(x -> x != 1, Ns) >= 2
                # get the second fastest changing index
                arr = read_arr()
                nonzero_index = find(x-> x != 0., arr)
                @assert length(nonzero_index) == 1
                @assert isapprox(sum(arr), 1/Ns[nonzero_index[1]],rtol=1e-2)
                kpt_order[2] = nonzero_index[1]
            else
                # special 1D case
                kpt_order[2] = kpt_order[3] == 1 ? 2 : 1
            end

            # get last fastest changing index
            kpt_order[1] = 6 - sum(kpt_order[2:3])
        end
    end
    @assert N1 != -1
    @assert N2 != -1
    @assert N3 != -1
    @assert nwannier != -1
    if nband == -1
        nband = nwannier
    end
    @assert all(x->!isnan(x), unit_cell)

    # TODO implement that logic (a point has two neighbors that are the same)
    @assert N1 != 2
    @assert N2 != 2
    @assert N3 != 2

    recp_unit_cell = 2pi * inv(unit_cell)'

    println("Reading $filename.mmn")
    mmn = open("$filename.mmn")
    readline(mmn) # skip header
    line = readline(mmn)
    nband2,Ntot,nntot = Base.map(x -> parse(Int64,x),split(line))
    @assert nband2 == nband
    @assert Ntot == N1*N2*N3

    M = zeros(Complex128,nband,nband,Ntot,nntot) # overlap matrix, (K) representation
    neighbors = zeros(Int64,Ntot,nntot) # for each point, list of neighbors, (K) representation
    displ_vecs = zeros(Int64,3,Ntot,nntot)
    while !eof(mmn)
        for nneighbor = 1:nntot
            line = readline(mmn)
            arr = split(line)
            K = parse(Int64, arr[1])
            Kpb = parse(Int64, arr[2])
            displ_vecs[:,K,nneighbor] = [parse(Int64, arr[3]),parse(Int64, arr[4]),parse(Int64, arr[5])]
            neighbors[K,nneighbor,:] = Kpb
            for mn = 0:nband^2-1
                m,n = mod(mn,nband)+1, div(mn,nband)+1
                line = readline(mmn)
                arr = split(line)
                ol = parse(Float64, arr[1]) + im*parse(Float64, arr[2])
                M[m,n,K,nneighbor] = ol
		@assert !isnan(ol)
            end
        end
    end

    t1 = collect(0:N1-1)/N1
    t2 = collect(0:N2-1)/N2
    t3 = collect(0:N3-1)/N3

    # We switch between a big index K=1:N^3 and three indices i,j,k = 1:N using these arrays
    K_to_ijk = zeros(Int64,Ntot,3)
    ijk_to_K = zeros(Int64,N1,N2,N3)
    Ns = [N1,N2,N3]
    for i=0:N1-1
        for j=0:N2-1
            for k=0:N3-1
                ijk = [i,j,k]
                big_offset = Ns[kpt_order[1]]*Ns[kpt_order[2]]
                med_offset = Ns[kpt_order[2]]
                K = ijk[kpt_order[1]]*big_offset + ijk[kpt_order[2]]*med_offset + ijk[kpt_order[3]]+1
                ijk_to_K[i+1,j+1,k+1] = K
                K_to_ijk[K,:] = [i+1 j+1 k+1]
            end
        end
    end

    #try to guess wb from (B1) of MV
    sumb = zeros(3,3)
    @assert ijk_to_K[1,1,1] == 1
    for ib = 1:nntot
        neighbor = neighbors[1,ib]
        i_n,j_n,k_n = K_to_ijk[neighbor,:]
	b = recp_unit_cell*([t1[i_n];t2[j_n];t3[k_n]] .+ displ_vecs[:,1,ib] .- [t1[1];t2[1];t3[1]])
        sumb += b*b'
    end
    @assert sumb / norm(sumb) ≈ diagm([Int(N1>1), Int(N2>1), Int(N3>1)])
    wb = 1/norm(sumb)
    println("Finite difference condition satisfied, wb=$wb")

    A = zeros(Complex128,nband,nwannier,N1,N2,N3)
    if read_amn
        println("Reading $filename.amn")
        amn = open("$filename.amn")
        readline(amn)
        arr = split(readline(amn))
        @assert parse(Int64, arr[1]) == nband
        @assert parse(Int64, arr[2]) == N1*N2*N3
        @assert parse(Int64, arr[3]) == nwannier
        lines = readlines(amn)
        for line in lines
            arr = split(line)
            m = parse(Int64, arr[1])
            n = parse(Int64, arr[2])
            kpt = parse(Int64, arr[3])
            Aijkmn = parse(Float64, arr[4]) + im*parse(Float64, arr[5])
            A[m,n,K_to_ijk[kpt,1], K_to_ijk[kpt,2], K_to_ijk[kpt,3]] = Aijkmn
        end

        # normalization
        for i=1:N1,j=1:N2,k=1:N3
            A[:,:,i,j,k] = normalize_matrix(A[:,:,i,j,k])
        end
    else
        for n = 1:nwannier
            A[n,n,:,:,:] = 1
        end
    end

    eigs = zeros(Ntot,nband)
    if read_eig
        println("Reading $filename.eig")
        eigfile = open("$filename.eig")
        while(!eof(eigfile))
            arr = split(readline(eigfile))
            iband = parse(Int,arr[1])
            iK = parse(Int,arr[2])
            val = parse(Float64,arr[3])
            eigs[iK,iband] = val
        end
    end

    map = true
    logMethod = true

    println("$nband bands, $nwannier Wannier, $N1 x $N2 x $N3 grid, $nntot neighbors")

    return WannierParameters(N1,N2,N3,t1,t2,t3,M,A,eigs,neighbors,kpt_order,ijk_to_K,K_to_ijk,nntot,nband,nwannier,filename,map,logMethod,recp_unit_cell,wb,displ_vecs)
end

# Computes overlap between two neighboring K points
function overlap(K1,K2,p)
    #fix useful if N3 = 1 (e.g. for 2D models)
    if(K1 == K2)
    	return eye(Complex128,p.nband)
    end
    for nneighbor=1:p.nntot
        if p.neighbors[K1,nneighbor] == K2
            return view(p.M, :, :, K1, nneighbor)
        end
    end
    error("No neighbors found, K1 = $(K1), K2 = $(K2)")
    return eye(Complex128,p.nband)
end
# Computes the overlap between K1 and K2, rotated by A1 and A2 respectively
function overlap(v1,v2,A1,A2,p)
    return (A1')*overlap(v1, v2, p)*A2
end
# Computes the overlap between K1 and K2, using the rotation contained in the array A
function overlap_A(v1,v2,p)
    i1,j1,k1 = v1
    i2,j2,k2 = v2
    K1 = p.ijk_to_K[i1,j1,k1]
    K2 = p.ijk_to_K[i2,j2,k2]
    A1 = p.A[:,:,i1,j1,k1]
    A2 = p.A[:,:,i2,j2,k2]
    return overlap(K1,K2,A1,A2,p)
end

# Power of a unitary (or at least, normal) matrix A
function powm(A,p)
    d,V = eig(A)
    V = normalize_matrix(V)
    accuracy = vecnorm(V*Diagonal(d)*V'-A)
    @assert(accuracy < 1e-10)
    return V*Diagonal(d.^p)*V'
end

# Normalize a matrix A to be unitary. If X is a matrix with orthogonal columns and A a non-singular matrix, then Löwdin-orthogonalizing X*A is equivalent to computing X*normalize_matrix(A)
function normalize_matrix(A)
    U,S,V = svd(A)
    accuracy = vecnorm(U*Diagonal(S)*V'-A)
    @assert(accuracy < 1e-10)
    return U*V'
end

function normalize_matrix_chol(A)
    ovl = A'A
    return A/chol(ovl)
end

# normalize and freeze a block of a matrix
# A = [Uf* Ur*]
# A*A = I, Uf* Uf = I
# UfUf* + UrUr* = I
# From this, obtain Uf*Ur = 0
# Strategy: first orthogonalize Uf, then project Uf out of Ur, then orthogonalize the range of Ur
function normalize_and_freeze(A,frozen,not_frozen)
    # orthogonalize Uf
    Uf = A[frozen,:]'
    U,S,V = svd(Uf)
    Uf = U*V'
    # Uf = normalize_matrix_chol(Uf)

    # project Uf out of Ur
    Ur = A[not_frozen,:]'
    Ur -= Uf*Uf'*Ur

    # # alternative method, maybe more stable but slower
    # ovl = Ur'Ur
    # S, U = eig(Hermitian(ovl))
    # S = real(S)
    # @assert !any(x -> 1e-11 <= x <= 1e-9, S)
    # @assert count(x -> x > 1e-10, S) == size(A,2) - nfrozen
    # Sm12 = map(x-> x < 1e-10 ? 0. : 1/sqrt(x), S)
    # Ur = Ur*(U*diagm(Sm12)*U')

    # renormalize the range of Ur
    U,S,V = svd(Ur)
    eps = 1e-10
    @assert !any(x -> 1e-11 <= x <= 1e-9, S)
    @assert count(x -> x > 1e-10, S) == size(A,2) - count(frozen)
    S[S .> eps] = 1
    S[S .< eps] = 0
    Ur = U*Diagonal(S)*V'

    A[not_frozen,:] = Ur'

    B = vcat(Uf',Ur')
    B[frozen,:] = Uf'
    B[not_frozen,:] = Ur'
    @assert B'B ≈ eye(size(A,2)) rtol=1e-12
    @assert B[frozen,:]*B[frozen,:]' ≈ eye(count(frozen)) rtol=1e-12
    @assert vecnorm(Uf'*Ur)<1e-10
    return B
end

# Propagate A0, defined at the first kpt, to the given list of kpts.
# Those must be neighbors, and only the first kpoint is assumed to have been rotated
function propagate(A0, kpts,p)
    N = length(kpts)
    n,m = size(A0)
    @assert(n==m,"Non square matrix given as argument in propagate")
    As = zeros(Complex128,n,n,N)
    As[:,:,1] = A0
    for i0=2:N
        As[:,:,i0] = normalize_matrix(overlap(kpts[i0],kpts[i0-1],p)*As[:,:,i0-1])
    end
    return As
end

# Plot obstructions
function plot_surface_obstructions(p, suffix="")
    if p.N3 != 1
        phases = zeros(p.N1,p.N2,nwannier)
        for i=1:p.N1,j=1:p.N2
            Obs = normalize_matrix(overlap_A([i,j,p.N3],[i,j,1],p)) #rotation at image point
            phases[i,j,:] = sort(imag(log(eig(Obs)[1])))
        end
        figure()
        xx = [p.t1[i] for i=1:p.N1,j=1:p.N2]
        yy = [p.t2[j] for i=1:p.N1,j=1:p.N2]
        for n=1:nwannier
            plot_surface(xx,yy,phases[:,:,n],rstride=1,cstride=1)
        end
        savefig("wannierize$filename.$suffix.pdf")
        close()
    end
end

function write_amn(p,A,filename)
    ## Output amn file
    out = open("$filename.amn","w")
    write(out, "Created by wannierize.jl ", string(now()), "\n")
    write(out, "$(p.nband) $(p.N1*p.N2*p.N3) $(p.nwannier)\n")
    for K=1:p.N1*p.N2*p.N3
        for n=1:p.nwannier
            for m = 1:p.nband
                coeff = A[m,n,p.K_to_ijk[K,1], p.K_to_ijk[K,2], p.K_to_ijk[K,3]]
                write(out, "$m $n $K $(real(coeff)) $(imag(coeff))\n")
            end
        end
    end
    close(out)
end

max_stars_repo_name: VVitale/wannier

max_issues_repo_name: VVitale/wannier

max_forks_repo_name: VVitale/wannier


text:
export MarchenkoPastur

struct MarchenkoPastur <: ContinuousUnivariateDistribution
    gamma::Real
    MarchenkoPastur(gamma) = 0 < gamma <= 1 ? new(gamma) : error("Gamma must be in (0, 1]")
end

function minimum(d::MarchenkoPastur)
    (1 - sqrt(d.gamma))^2
end

function maximum(d::MarchenkoPastur)
    (1 + sqrt(d.gamma))^2
end

function pdf(d::MarchenkoPastur, x)
    lambdamin = minimum(d)
    lambdamax = maximum(d)

    if lambdamin < x < lambdamax
        return sqrt((lambdamax - x) * (x - lambdamin))/(d.gamma * x * 2 * pi)
    else
        return 0
    end
end

ext: jl


text:
# ---
# title: 1402. Reducing Dishes
# id: problem1402
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/reducing-dishes/description/>
# hidden: true
# ---
#
# A chef has collected data on the `satisfaction` level of his `n` dishes. Chef
# can cook any dish in 1 unit of time.
#
# _Like-time coefficient_  of a dish is defined as the time taken to cook that
# dish including previous dishes multiplied by its satisfaction level  i.e.
# `time[i]`*`satisfaction[i]`
#
# Return the maximum sum of  _Like-time coefficient_ that the chef can obtain
# after dishes preparation.
#
# Dishes can be prepared in **any** order and the chef can discard some dishes
# to get this maximum value.
#
#
#
# **Example 1:**
#
#
#
#     Input: satisfaction = [-1,-8,0,5,-9]
#     Output: 14
#     Explanation: After Removing the second and last dish, the maximum total _Like-time coefficient_ will be equal to (-1*1 + 0*2 + 5*3 = 14). Each dish is prepared in one unit of time.
#
# **Example 2:**
#
#
#
#     Input: satisfaction = [4,3,2]
#     Output: 20
#     Explanation: Dishes can be prepared in any order, (2*1 + 3*2 + 4*3 = 20)
#
#
# **Example 3:**
#
#
#
#     Input: satisfaction = [-1,-4,-5]
#     Output: 0
#     Explanation: People don't like the dishes. No dish is prepared.
#
#
# **Example 4:**
#
#
#
#     Input: satisfaction = [-2,5,-1,0,3,-3]
#     Output: 35
#
#
#
#
# **Constraints:**
#
#   * `n == satisfaction.length`
#   * `1 <= n <= 500`
#   * `-10^3 <= satisfaction[i] <= 10^3`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@test length(coalitions([1,2,3])) == 2^3
@test length(coalitions([1,2,3,4])) == 2^4
@test sum([w for (c,w) in coalitions([1,2,3])]) ≈ 1
@test sum([w for (c,w) in coalitions(collect(1:9))]) ≈ 1
@test sum([w for (c,w) in coalitions(collect(1:10))]) ≈ 1
it = ShapleyValues.statefuliterator(coalitions([1,2]))
@test next(it)[1] == Int64[]
@test next(it)[1] == [1,2]
@test next(it)[1] == [1]
@test next(it)[1] == [2]
@test done(it)
it = ShapleyValues.statefuliterator(coalitions([1,2,3]))
@test next(it)[1] == Int64[]
@test next(it)[1] == [1,2,3]
@test next(it)[1] == [1]
@test next(it)[1] == [2,3]
@test next(it)[1] == [2]
@test next(it)[1] == [1,3]
@test next(it)[1] == [3]
@test next(it)[1] == [1,2]
@test done(it)


text:
t#
# Benchmark example
#
using DirectConvolution
using DSP, BenchmarkTools,LinearAlgebra

# function bench_directconv(filter,signal)
#     wrapped_filter = LinearFilter(filter,0)
#     convolved = directConv(wrapped_filter,signal)

#     convolved
# end
function bench_directconv(filter,signal)
    convolved = similar(signal)
    n=length(convolved)

    directConv!(filter,0,-1,signal,convolved,1:n)

    convolved
end

function bench_dsp(filter,signal)
    convolved = conv(filter,signal)

    n=length(signal)
    convolved[1:n]
end
# function bench_check(;filter_length,signal_length)
#     r1 =
# end

function bench(;filter_length,signal_length)
    @assert filter_length>0 && signal_length>0

    filter = rand(filter_length)
    signal = rand(signal_length)

    t1 = @benchmark bench_directconv($filter,$signal)
    t2 = @benchmark bench_dsp($filter,$signal)

    t1_min=minimum(t1.times)
    t2_min=minimum(t2.times)

    println()
    println("filter = $filter_length, signal = $signal_length")
    println("DirectConvolution: $(t1_min)μs, $(t1.allocs) allocs")
    println("DSP              : $(t2_min)μs, $(t2.allocs) allocs")
    println("Ratio ",t2_min/t1_min)
end

t = bench(filter_length=5,signal_length=10)
t = bench(filter_length=5,signal_length=1000)
t = bench(filter_length=5,signal_length=10000)

t = bench(filter_length=15,signal_length=10)
t = bench(filter_length=15,signal_length=1000)
t = bench(filter_length=15,signal_length=10000)

t = bench(filter_length=50,signal_length=1000)
t = bench(filter_length=500,signal_length=1000)

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
# Finding greatest common divisor using Euclidean Algorithm
function gcdFunction(a, b)
    if a == 0
        return b
    end
    return gcdFunction(b % a, a)
end
function power(x, y, m)
    if y == 0
        return 1
    end
    p = power(x, y / 2, m) % m
    p = (p * p) % m

    if y % 2 == 0
        return p
    else
        return (x * p) % m
    end
end
function modInverse(a, m)
    gcd = gcdFunction(a, m)

    if gcd != 1
        return -1
    else
        return power(a, m - 2, m)
    end
end

a = 3
m = 11

ans = modInverse(a, m)

if ans == -1
    println("Inverse doesn't exist")
else
    println("Modular multiplicative inverse is $(ans-1)")
end

# Output

#Modular multiplicative inverse is 4

max_stars_repo_name: rajatenzyme/Coding-Journey-

max_issues_repo_name: rajatenzyme/Coding-Journey-

max_forks_repo_name: rajatenzyme/Coding-Journey-


text:
using SumOfSquares
using DynamicPolynomials

function gram_to_poly(n)
    model = MOIU.Model{Float64}()
    @polyvar x
    monos = monomials(x, 0:n)
    q, Q, con_Q = SumOfSquares.add_gram_matrix(model, MOI.PositiveSemidefiniteConeTriangle, monos, Float64)
    @time polynomial(q)
    return
end
gram_to_poly(1000)


text:
"""Module for marine sedimentation analysis."""
module Sediments

using ArchGDAL
using ArgParse
using Dates

include(".\\Utils\\Functions.jl")

export run_sediment

const agd = ArchGDAL

mutable struct Sediment <: Functions.AbstractAnalysisObject
 # Parameters
    dredged_mass::Float64                           # Sediment imput rate (kg/sec)
    time::Int64                                     # Final temporal step
    mean_depth::Float64                             # Depth (m)
    x_dispersion_coeff::Float64                     # Diffusion X component
    y_dispersion_coeff::Float64                     # Diffusion Y component
    x::Float64                                      # Target coordinate X
    y::Float64                                      # Target coordinate Y
    mean_flow_speed::Float64                        # Average current's speed
    direction::Float64                              # Direction fo the flow (°)
    mean_sedimentation_velocity::Float64            # Average sediment velocity
    time_intreval::Int64                            # Interval time of analysis, for integral discretization
    current_oscillatory_amplitude::Float64          # Amplitude of the oscillatory current
    tide::Int64                                     # Tidal cycle (h)
    ω::Float64
  # Computational results
    ew::Float64

    function Sediment(dredged_mass,time,mean_depth,x_dispersion_coeff,y_dispersion_coeff,x,y,mean_flow_speed,direction,mean_sedimentation_velocity,time_intreval,current_oscillatory_amplitude,tide)
        ω = (current_oscillatory_amplitude > 0) && (tide > 0) ? 2.0π / tide : 0.0
        return new(dredged_mass,time,mean_depth,x_dispersion_coeff,y_dispersion_coeff,x,y,mean_flow_speed,direction,mean_sedimentation_velocity,time_intreval,current_oscillatory_amplitude,tide,ω,0.0)
    end
end

function calc_e!( s::Sediment, i )
    s.ew = s.ω > 0 ? s.current_oscillatory_amplitude / ( s.ω * cos(deg2rad(s.ω)) - cos(deg2rad(s.ω * i * s.time_intreval)) ) : 0.0
    e1 = ℯ^(-(( s.x - s.mean_flow_speed * ( s.time - i * s.time_intreval) + s.ew ) / ( 4.0s.x_dispersion_coeff * (s.time - i * s.time_intreval) ) ))
    e2 = ℯ^(-( s.y^2 / ( 4.0s.y_dispersion_coeff * (s.time - i * s.time_intreval) ) ) - ( (s.mean_sedimentation_velocity * (s.time - i * s.time_intreval)) / s.mean_depth ) )
    return e1 * e2
end

function Functions.compute_concentration!( s::Sediment )
    if s.x > 0
        q = s.dredged_mass / ( 4.0π * s.mean_depth * √(s.x_dispersion_coeff * s.y_dispersion_coeff) )
        n = s.time ÷ s.time_intreval
        csum = 0.0
        @inbounds for i in 0:n-1
            csum += calc_e!(s, i) * ( 1 / ( s.time - ( i * s.time_intreval ) ) )
        end
        return q * csum * s.time_intreval
    end
    return 0.0
end

"""
    run_sediment(; dem_file::String, source_file::String, resolution::Float64, mean_flow_speed::Float64, mean_depth::Float64, x_dispersion_coeff::Float64,
                   y_dispersion_coeff::Float64, dredged_mass::Float64, tollerance::int64=2, flow_direction::Float64, mean_sedimentation_velocity::Float64,
                   time::Int64, time_intreval::Int64, current_oscillatory_amplitude::Float64=0.0, tide::Int64=0, output_path::String=".\\sediment_output_model.tiff" )

Create and save as `output_path` a raster containing the results of model of plumes of turbidity induced by dredging.

# Arguments
- `dem_file::String`: path to the raster of terrain.
- `source_file::String`: path to the shapefile containing the dredging source point.
- `resolution::Float64`: size of a cell in meters.
- `mean_flow_speed::Float64`: average current's speed.
- `mean_depth::Float64`: depth in meters.
- `x_dispersion_coeff::Float64`: coefficient of dispersion along the x axis.
- `y_dispersion_coeff::Float64,`: coefficient of dispersion along y axis.
- `contaminantCASNum::String`: CAS number identifier of a substance.
- `dredged_mass::Float64`: initial mass of the dredged substance.
- `tollerance::Int64=2`: value used to determine wether the concentration of pollutant in a cell is relevant.
    Specifically, a concentration value is considered relevant if its value is within "tollerance" orders of magnitute from the concentration on other cells.
- `flow_direction::Float64`: direction of the flow as an angle in degrees.
- `mean_sedimentation_velocity::Float64`: velocity of sedimentation.
- `time::Int64`: start time for the model.
- `time_intreval::Int64`: length of an epoch.
- `current_oscillatory_amplitude::Int64=0`: water oscillatory amplitude.
- `tide::Int64=0`: tidal cycle in hours.
- `output_path::String=".\\output_model_sediments.tiff"`: path of the resulting raster.
"""
function run_sediment(; dem_file::String, source_file::String, resolution::Float64, mean_flow_speed::Float64, mean_depth::Float64, x_dispersion_coeff::Float64,
                       y_dispersion_coeff::Float64, dredged_mass::Float64, tollerance::Int64=2, flow_direction::Int64, mean_sedimentation_velocity::Float64,
                       time::Int64, time_intreval::Int64, current_oscillatory_amplitude::Float64=0.0, tide::Int64=0, output_path::String=".\\sediment_output_model.tiff" )

 # messaggio+='ALGORITMO UTILIZZATO: Shao (Shao, Dongdong, et al. "Modeling dredging-induced turbidity plumes in the far field under oscillatory tidal currents." Journal of Waterway, Port, Coastal, and Ocean Engineering 143.3 (2016))\n\n'

    geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(geom) != 0
        throw(DomainError(source, "`source` must be a point."))
    end

    dem = agd.read(dem_file)
    refsys = agd.getproj(dem)
    geotransform = agd.getgeotransform(dem)

    if refsys != agd.toWKT(agd.getspatialref(geom))
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    x_source = agd.getx(geom, 0)
    y_source = agd.gety(geom, 0)
    r_source, c_source = Functions.toIndexes(dem, x_source, y_source)

    #   start_time = time.time()

    sediment = Sediment(dredged_mass, time, mean_depth, x_dispersion_coeff, y_dispersion_coeff, 0.0, 0.0, mean_flow_speed,
                        flow_direction, mean_sedimentation_velocity, time_intreval, current_oscillatory_amplitude, tide)
    points = Functions.expand(r_source, c_source, dredged_mass, tollerance, dem, sediment)

    minR = minimum( point -> point[1], points )
    minC = minimum( point -> point[2], points )
    maxR = maximum( point -> point[1], points )
    maxC = maximum( point -> point[2], points )

    geotransform[[1, 4]] = Functions.toCoords(dem, minR, minC)
    noData = Float32(agd.getnodatavalue(agd.getband(dem, 1)))
    data = fill(noData, maxR-minR+1, maxC-minC+1)
    for r in minR:maxR, c in minC:maxC
        match = findfirst( p -> p[1] == r && p[2] == c, points )
        if !isnothing(match)
            data[r-minR+1, c-minC+1] = points[match][3]
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, refsys, noData, output_path)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
function main()
    e = zeros(2, 500)
    e[2, 500] = 2
    e[1, 500] = 2.004

    for k in 498:-1:0
        e[2, k+1] = 1000 * ((499 - k) * e[2, k+2]/500 + 1)/(999 - k)
        e[1, k+1] = 1000 * ((499 - k) * e[1, k+2]/500 + e[2, k+1]/1000 + 1)/(999 - k)
    end

    println(round(e[1, 1], digits=8))
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
abstract type Adjoint_rep_Gaugefields_4D{NC,NumofBasis}  <: Adjoint_rep_Gaugefields{NC,4,NumofBasis}
end

include("./Adjoint_rep_gaugefields_4D_wing.jl")

function Base.size(U::Adjoint_rep_Gaugefields_4D{NC,NumofBasis}) where {NC,NumofBasis}
    return NumofBasis,NumofBasis,U.NX,U.NY,U.NZ,U.NT
end


text:
@testset "Utils" begin

    A = [missing 1.0; 2.0 missing]
    ClimateTools.replace_missing!(A)
    @test isnan(A[1,1])
    @test isnan(A[2,2])

    x = 20.0:0.5:150

    ω = ClimateTools.weight(x, frac=0.25, power=1.0)
    @test maximum(ω) == 1.0
    @test minimum(ω) == 0.0
    @test (ω[34]+ω[33])/2 == 0.5
    @test findfirst(ω .== 1.0) == 66 # idx where ω = 1.0 and cutoff is set to 0.25

    ω = ClimateTools.weight(x, frac=0.25, power=2.0)
    @test maximum(ω) == 1.0
    @test minimum(ω) == 0.0
    @test round((ω[34]+ω[33])/2, digits=2) == 0.25
    @test findfirst(ω .== 1.0) == 66 # idx where ω = 1.0 and cutoff is set to 0.25

end


text:
# Improved Newton's Method

# Safe version - disclosed version

#par1 - lower bound
#par2 - upper bound
function solve_Newton_safely(f,f1,par,par1,par2;nlim = 1e4,ϵ=1e-3)
    for i in 1:nlim
        par3 = par - f(par)/f1(par)
        if ϵ>abs(par3 -par)
            return(par3)
        elseif par3<=par1
            par3 = (par+par1)/2
        elseif par3>=par2
            par3 = (par+par2)/2
        end

        par=par3
    end
    println("Warning! unconverged after" * string(nlim) *"Iterations.")
    return(par3)
end

#where f1 is a constant (avoid finding gradient)
function solve_constGrad(f,f1,par;nlim = 1e4,ϵ=1e-3)
    for i in 1:nlim
        par3 = par - f(par)/f1
        if ϵ>abs(par3 -par)
            return(par3)
        end
        par=par3
    end
    println("Warning! unconverged after" * string(nlim) *"Iterations.")
    return(par3)
end

# by cotangent value
# where
# par1 and par2 are 2 starting points( to give 1st gradient value)
function solve_cotangent(f,par1,par2;nlim = 1e4,ϵ=1e-3)
    for i in 1:nlim
        par3 = par2 - f(par2)*(par2 - par1)/(f(par2) - f(par1))
        if ϵ>abs(par3 -par2)
            return(par3)
        end
        par1 = par2;par2 = par3
    end
    println("Warning! unconverged after" * string(nlim) *"Iterations.")
    return(par3)
end

using Plots

x = 0:.1:4
y=sin.(x)
plot(x,y)
plot!([0,4],[0,0])
solve_cotangent(sin,3,4)

max_stars_repo_path: Improved Newton Solver.jl

max_issues_repo_path: Improved Newton Solver.jl

max_forks_repo_path: Improved Newton Solver.jl


text:
@testset "Bayes" begin
    m=NaiveBayesClassifier([:spam, :ham])
    TextAnalysis.fit!(m, "this is ham", :ham);
    TextAnalysis.fit!(m, "this is spam", :spam);
    r = TextAnalysis.predict(m, "is this spam?")
    @test r[:spam] > r[:ham]
end


text:
"""
    module PhysicalProblem
"""
module PhysicalProblem

    using DocStringExtensions

    using StaticArrays

    import WaveProp.Geometry: ambient_dimension, geometric_dimension

    export
        AbstractPDE,
        Laplace,
        Helmholtz,
        Elastostatic,
        Maxwell,
        default_density_eltype,
        default_kernel_eltype

    include("pde.jl")

end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Flux, Flux.Data.MNIST, Statistics, BSON, Random, Plots; pyplot()
using Flux: onehotbatch, onecold, crossentropy
Random.seed!(0)

epochs = 30
eta = 5e-3
batchSize = 1000
trainRange, validateRange = 1:5000, 5001:10000

function minibatch(x, y, indexRange)
    xBatch = Array{Float32}(undef, size(x[1])..., 1, length(indexRange))
    for i in 1:length(indexRange)
        xBatch[:, :, :, i] = Float32.(x[indexRange[i]])
    end
    return (xBatch, onehotbatch(y[indexRange], 0:9))
end

trainLabels = MNIST.labels()[trainRange]
trainImgs = MNIST.images()[trainRange]
mbIdxs = Iterators.partition(1:length(trainImgs), batchSize)
trainSet = [minibatch(trainImgs, trainLabels, bi) for bi in mbIdxs]

validateLabels = MNIST.labels()[validateRange]
validateImgs = MNIST.images()[validateRange]
validateSet = minibatch(validateImgs, validateLabels, 1:length(validateImgs))

model1= Chain(flatten, Dense(784, 200,relu),Dense(200, 100,tanh),
				Dense(100, 10,sigmoid), softmax)

model2= Chain(Conv((5, 5), 1=>8, relu), MaxPool((2,2)),
                Conv((3, 3), 8=>16, relu), MaxPool((2,2)),
                flatten, Dense(400, 10), softmax)

opt1 = ADAM(eta); opt2 = ADAM(eta)
accuracyPaths = [[],[]]
accuracy(x, y, model) = mean(onecold(model(x)) .== onecold(y))
loss(x, y, model) = crossentropy(model(x), y)
cbF1() = push!(accuracyPaths[1],accuracy(validateSet..., model1))
cbF2() = push!(accuracyPaths[2],accuracy(validateSet..., model2))

model1(trainSet[1][1]); model2(trainSet[1][1])
for _ in 1:epochs
    Flux.train!((x,y)->loss(x,y,model1), params(model1), trainSet, opt1, cb=cbF1)
    Flux.train!((x,y)->loss(x,y,model2), params(model2), trainSet, opt2, cb=cbF2)
	print(".")
end

println("\nModel1 (Dense) accuracy = ", accuracy(validateSet..., model1))
println("Model2 (Convolutional) accuracy = ", accuracy(validateSet..., model2))
# cd(@__DIR__)
# BSON.@save "../data/mnistConv.bson" modelParams=cpu.(params(model2))
plot(accuracyPaths,label = ["Dense" "Convolutional"],
	ylim=(0.7,1.0), xlabel="Batch number", ylabel = "Validation Accuracy")

max_stars_repo_name: yoninazarathy/JuliaMLForGRDC

max_issues_repo_name: yoninazarathy/JuliaMLForGRDC

max_forks_repo_name: yoninazarathy/JuliaMLForGRDC


text:
abstract type F81 <: NASM end

const _πA(mod::F81) = mod.πA
const _πC(mod::F81) = mod.πC
const _πG(mod::F81) = mod.πG
const _πT(mod::F81) = mod.πT

function P(mod::F81, t::Array{Float64})
  return [P(mod, i) for i in t]
end


text:
###############################################################################
#################### Submatrix from the proportional terms ####################
###############################################################################
function (method::SemiDiscretization)(A::CyclicVector{<:AbstractMatrix}, rst::AbstractResult{d}) where d
   fill(SubMX(subMxRange(0, d), exp(A[1] * method.Δt)), rst.n_steps)
end
function (method::SemiDiscretization)(A::Vector{<:AbstractMatrix}, rst::AbstractResult{d}) where d
    SubMX.((subMxRange(0, d),), exp.(A .* method.Δt))
end

###############################################################################
###################### Submatrix from the delayed terms #######################
###############################################################################
#τ discretisation
function (method::SemiDiscretization)(τ::Real, rst::AbstractResult{d}) where d
    r = rOfDelay(τ, method)
    τerr = τ - r * method.Δt
    ranges = [subMxRange(r - k, d) for k in 0:methodorder(method)]
    return (τerr, ranges)
end
function (method::SemiDiscretization)(τs::Vector{<:Real}, rst::AbstractResult{d}) where d
    rs = rOfDelay.(τs, Ref(method))
    τerrs = τs .- rs .* method.Δt
    rangess = [[subMxRange(r - k, d) for k in 0:methodorder(method)] for r in rs]
    return (τerrs, rangess)
end

###############################################################################
# To distinguish between B as a constant matrix vs B as a function matrix

function (method::SemiDiscretization{<:NumericSD})(τ::Real, B::DelayMX{d,dT,<:AbstractMatrix}, rst::AbstractResult) where {d,dT}
    (τerr, ranges) = method(τ, rst)
    MXs = [[QuadGK.quadgk((t -> (exp(rst.A_avgs[i] * (method.Δt - t)) * lagr_el0(methodorder(method), k, τerr, method.Δt, t))), 0.0, method.Δt)[1] * B.MX for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    # SubMX.((ranges,), MXs)
    # TODO: this is for 1.0
    SubMX.(Ref(ranges), MXs)
end
function (method::SemiDiscretization{<:NumericSD})(τs::Vector{<:Real}, B::DelayMX{d,dT,<:AbstractMatrix}, rst::AbstractResult) where {d,dT}
    (τerrs, rangess) = method(τs, rst)
    MXs = [[QuadGK.quadgk((t -> (exp(rst.A_avgs[i] * (method.Δt - t)) * lagr_el0(methodorder(method), k, τerrs[i], method.Δt, t))), 0.0, method.Δt)[1] * B.MX for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    SubMX.(rangess, MXs)
end
function (method::SemiDiscretization{<:NumericSD})(τ::Real, B::DelayMX{d,dT,<:Function}, rst::AbstractResult) where {d,dT}
    (τerr, ranges) = method(τ, rst)
    MXs = [[QuadGK.quadgk((t -> (exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * B(t) * lagr_el0(methodorder(method), k, τerr, method.Δt, t-rst.ts[i]))), rst.ts[i], rst.ts[i + 1])[1] for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    # SubMX.((ranges,), MXs)
    # TODO: this is for 1.0
    SubMX.(Ref(ranges), MXs)
end
function (method::SemiDiscretization{<:NumericSD})(τs::Vector{<:Real}, B::DelayMX{d,dT,<:Function}, rst::AbstractResult) where {d,dT}
    (τerrs, rangess) = method(τs, rst)
    MXs = [[QuadGK.quadgk((t -> (exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * B(t) * lagr_el0(methodorder(method), k, τerrs[i], method.Δt, t-rst.ts[i]))), rst.ts[i], rst.ts[i + 1])[1] for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    SubMX.(rangess, MXs)
end

###############################################################################
# To distinguish between τ as constant vs τ as function
function (method::SemiDiscretization{<:NumericSD})(B::DelayMX{d,<:Real,mT}, rst::AbstractResult) where {d,mT}
    method(B.τ.τ, B, rst)
    # rest=method(B.τ.τ, B, rst)
    # println(typeof(rest))
    # rest
end

function (method::SemiDiscretization{<:NumericSD})(B::DelayMX{d,<:Function,mT}, rst::AbstractResult) where {d,mT}
    τis = [quadgk(B.τ, rst.ts[i], rst.ts[i + 1])[1] / method.Δt for i in 1:rst.n_steps]
    method(τis, B, rst)
    # rest=method(τis, B, rst)
    # println(typeof(rest))
    # rest
end

###############################################################################
###################### Submatrix from the additive terms ######################
###############################################################################
function (method::SemiDiscretization)(c::Additive{d,<:AbstractArray{<:Real}}, rst::AbstractResult{d}) where d
    Vs = SubV.([SVector{d}(vec(quadgk(t -> exp(rst.A_avgs[i] * (method.Δt - t)), zero(method.Δt), method.Δt)[1] * c.V)) for i in 1:rst.n_steps])
end
function (method::SemiDiscretization)(c::Additive{d,<:Function}, rst::AbstractResult{d}) where d
    Vs = SubV.([SVector{d}(vec(quadgk(t -> exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * c(t), rst.ts[i], rst.ts[i + 1])[1])) for i in 1:rst.n_steps])
end
function (method::SemiDiscretization)(c::Additive{d,<:Array{<:Nothing}}, rst::AbstractResult) where d
    Vs = fill(SubV(@SArray zeros(d)),length(rst.ts)-1)
end


text:
# full size example
# there are 6 cases
# please write the data output (e.g. to csv) by yourself
using JuMP
using Gurobi

include("solveDEA.jl")
include("variable.jl")

# CRS input-oriented
function case1(data)
    # make up the data set
    X = data[:,1:6] # 6
    Z = data[:,7:10] # 4
    Y = data[:,11:15] # 5
    scale = size(data)[1]
    xNo = size(X)[2]
    zNo = size(Z)[2]
    yNo = size(Y)[2]

    for k = 1:scale
        # crs = Model(solver = GurobiSolver(Presolve=0, OutputFlag=0))
        model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.

        @solveDEA_varibleInit(model)
        @solveDEA_varible(model, Lambda[1:25000] >= 0, "nc") # non-critical group 1
        @solveDEA_varible(model, Gamma[1:25000] >= 0, "nc") # non-critical group 2
        @solveDEA_varible(model, middle[1:4] >= 0, "c") # critial
        @solveDEA_varible(model, Theta >= 0, "c") # critial

        @objective(model, Min, Theta)

        @constraint(model, inputConA[i=1:xNo], sum{Lambda[r]*X[r,i], r = 1:scale} <= Theta*X[k,i])
        @constraint(model, outputConA[j=1:zNo], sum{Lambda[r]*Z[r,j], r = 1:scale} >= middle[j])
        @constraint(model, inputConB[i=1:zNo], sum{Gamma[r]*Z[r,i], r = 1:scale} <= middle[i])
        @constraint(model, outputConB[j=1:yNo], sum{Gamma[r]*Y[r,j], r = 1:scale} >= Y[k,j])

        #solve(model)
        solveDEA(model)

        #println("lambdas: $(getvalue(Lambda)))")

        # retrieve results
        println("Objective value: ", getobjectivevalue(model))

    end
end

println("begins")
data = readcsv("C:/Users/mb516/Desktop/solveDEA_local_0307/15-7-25k_01.csv")
case1(data)
# case2(data)
# case3(data)
# case4(data)
# case5(data)
# case6(data)

println("done!")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ADIOS2"
version = v"2.7.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ornladios/ADIOS2/archive/refs/tags/v2.7.1.tar.gz", "c8e237fd51f49d8a62a0660db12b72ea5067512aa7970f3fcf80b70e3f87ca3e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ADIOS2-2.7.1
# Don't define clock_gettime on macOS
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/clock_gettime.patch
# See <https://github.com/ornladios/ADIOS2/issues/2705>
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/gettid.patch
# PR <https://github.com/ornladios/ADIOS2/pull/2712>
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ndims.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/shlwapi.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/sockaddr_in.patch

mkdir build
cd build
archopts=
if [[ "$target" == *-apple-* ]]; then
    # MPICH's pkgconfig file "mpich.pc" lists these options:
    #     Libs:     -framework OpenCL -Wl,-flat_namespace -Wl,-commons,use_dylibs -L${libdir} -lmpi -lpmpi -lm    -lpthread
    #     Cflags:   -I${includedir}
    # cmake doesn't know how to handle the "-framework OpenCL" option
    # and wants to use "-framework" as a stand-alone option. This fails
    # gloriously, and cmake concludes that MPI is not available.
    archopts="-DMPI_C_ADDITIONAL_INCLUDE_DIRS='' -DMPI_C_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi' -DMPI_CXX_ADDITIONAL_INCLUDE_DIRS='' -DMPI_CXX_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi'"
elif [[ "$target" == x86_64-w64-mingw32 ]]; then
    # - The MSMPI Fortran bindings are missing a function; see
    #   <https://github.com/microsoft/Microsoft-MPI/issues/7>
    echo 'void __guard_check_icall_fptr(unsigned long ptr) {}' >cfg_stub.c
    gcc -c cfg_stub.c
    ar -crs libcfg_stub.a cfg_stub.o
    cp libcfg_stub.a $prefix/lib
    # - cmake's auto-detection for MPI doesn't work on Windows.
    # - The SST and Table ADIOS2 components don't build on Windows
    #   (reported in <https://github.com/ornladios/ADIOS2/issues/2705>)
    export FFLAGS="-I$prefix/src -I$prefix/include -fno-range-check"
    archopts="-DMPI_GUESS_LIBRARY_NAME=MSMPI -DMPI_C_LIBRARIES=msmpi64 -DMPI_CXX_LIBRARIES=msmpi64 -DMPI_Fortran_LIBRARIES='msmpifec64;msmpi64;cfg_stub' -DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
elif [[ "$target" == *-mingw* ]]; then
    archopts="-DMPI_GUESS_LIBRARY_NAME=MSMPI -DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
fi
# Fortran is not supported with Clang
cmake \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DADIOS2_USE_Blosc=ON \
    -DADIOS2_USE_Fortran=OFF \
    -DADIOS2_USE_MPI=ON \
    -DADIOS2_USE_PNG=ON \
    -DADIOS2_USE_ZFP=ON \
    -DADIOS2_USE_ZeroMQ=ON \
    -DADIOS2_BUILD_EXAMPLES=OFF \
    -DBUILD_TESTING=OFF \
    -DMPI_HOME=$prefix \
    ${archopts} \
    -DADIOS2_INSTALL_GENERATE_CONFIG=OFF \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    ..
cmake --build . --config Release --parallel $nproc
cmake --build . --config Release --parallel $nproc --target install
install_license ../Copyright.txt ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# 32-bit architectures are not supported; see
# <https://github.com/ornladios/ADIOS2/issues/2704>
platforms = filter(p -> nbits(p) ≠ 32, platforms)
platforms = expand_cxxstring_abis(platforms)
# Windows doesn't build with libcxx="cxx03"
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libadios2_c", :libadios2_c),
    LibraryProduct("libadios2_c_mpi", :libadios2_c_mpi),
    LibraryProduct("libadios2_core", :libadios2_core),
    LibraryProduct("libadios2_core_mpi", :libadios2_core_mpi),
    LibraryProduct("libadios2_cxx11", :libadios2_cxx11),
    LibraryProduct("libadios2_cxx11_mpi", :libadios2_cxx11_mpi),

    # Missing on Apple:
    # LibraryProduct("libadios2_taustubs", :libadios2_taustubs),

    # Missing on Windows:
    # LibraryProduct("libadios2_atl", :libadios2_atl),
    # LibraryProduct("libadios2_dill", :libadios2_dill),
    # LibraryProduct("libadios2_evpath", :libadios2_evpath),
    # LibraryProduct("libadios2_ffs", :libadios2_ffs),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Blosc_jll")),
    # We don't want to use Bzip2 because this would lock us into Julia ≥1.6
    # Dependency(PackageSpec(name="Bzip2_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="MPICH_jll")),
    Dependency(PackageSpec(name="MicrosoftMPI_jll")),
    Dependency(PackageSpec(name="ZeroMQ_jll")),
    Dependency(PackageSpec(name="libpng_jll")),
    Dependency(PackageSpec(name="zfp_jll")),
    # We cannot use HDF5 because we need a HDF5 configuration with MPI support
    # Dependency(PackageSpec(name="HDF5_jll")),
]

# Build the tarballs, and possibly a `build.jl` as well.
# GCC 4 is too old for Windows; it doesn't have <regex.h>
# GCC 5 is too old for FreeBSD; it doesn't have `std::to_string`
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")


text: println("Hello World!")


text:
let
    gw = LegacyGridWorld(sx=2, sy=2, rs=[GridWorldState(1,1)], rv=[10.0])

    pvec = fill(GridWorldAction(:left), 5)

    solver = VectorSolver(pvec)

    p = solve(solver, gw)

    io = IOBuffer()
    d = TextDisplay(io)
    display(d, p)
    @test String(take!(io)) == "VectorPolicy{GridWorldState,Symbol}:\n GridWorldState(1, 1, false) -> :left\n GridWorldState(2, 1, false) -> :left\n GridWorldState(1, 2, false) -> :left\n GridWorldState(2, 2, false) -> :left\n GridWorldState(0, 0, true) -> :left"

    for s1 in states(gw)
        @test action(p, s1) == GridWorldAction(:left)
    end

    p2 = VectorPolicy(gw, pvec)
    for s2 in states(gw)
        @test action(p2, s2) == GridWorldAction(:left)
    end

    p3 = ValuePolicy(gw)
    for s2 in states(gw)
        @inferred(action(p3, s2)) isa GridWorldAction
    end

    io = IOBuffer()
    d = TextDisplay(io)
    display(d, p3)
    @test String(take!(io)) == "ValuePolicy{LegacyGridWorld,Array{Float64,2},Symbol}:\n GridWorldState(1, 1, false) -> :up\n GridWorldState(2, 1, false) -> :up\n GridWorldState(1, 2, false) -> :up\n GridWorldState(2, 2, false) -> :up\n GridWorldState(0, 0, true) -> :up"
end


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

@doc raw"""PodDNSConfig defines the DNS parameters of a pod in addition to those generated from DNSPolicy.

    IoK8sApiCoreV1PodDNSConfig(;
        nameservers=nothing,
        options=nothing,
        searches=nothing,
    )

    - nameservers::Vector{String} : A list of DNS name server IP addresses. This will be appended to the base nameservers generated from DNSPolicy. Duplicated nameservers will be removed.
    - options::Vector{IoK8sApiCoreV1PodDNSConfigOption} : A list of DNS resolver options. This will be merged with the base options generated from DNSPolicy. Duplicated entries will be removed. Resolution options given in Options will override those that appear in the base DNSPolicy.
    - searches::Vector{String} : A list of DNS search domains for host-name lookup. This will be appended to the base search paths generated from DNSPolicy. Duplicated search paths will be removed.
"""
mutable struct IoK8sApiCoreV1PodDNSConfig <: SwaggerModel
    nameservers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nameservers
    options::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodDNSConfigOption} } # spec name: options
    searches::Any # spec type: Union{ Nothing, Vector{String} } # spec name: searches

    function IoK8sApiCoreV1PodDNSConfig(;nameservers=nothing, options=nothing, searches=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("nameservers"), nameservers)
        setfield!(o, Symbol("nameservers"), nameservers)
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("searches"), searches)
        setfield!(o, Symbol("searches"), searches)
        o
    end
end # type IoK8sApiCoreV1PodDNSConfig

const _property_map_IoK8sApiCoreV1PodDNSConfig = Dict{Symbol,Symbol}(Symbol("nameservers")=>Symbol("nameservers"), Symbol("options")=>Symbol("options"), Symbol("searches")=>Symbol("searches"))
const _property_types_IoK8sApiCoreV1PodDNSConfig = Dict{Symbol,String}(Symbol("nameservers")=>"Vector{String}", Symbol("options")=>"Vector{IoK8sApiCoreV1PodDNSConfigOption}", Symbol("searches")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodDNSConfig }) = collect(keys(_property_map_IoK8sApiCoreV1PodDNSConfig))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodDNSConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodDNSConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodDNSConfig }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodDNSConfig[property_name]

function check_required(o::IoK8sApiCoreV1PodDNSConfig)
    true
end

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


text:
using Unitful

"""
    shift_prefix(u, e)

Shift the SI prefix of the unit `u` by `e` orders of magnitude (which
needs to be a multiple of 3).
"""
function shift_prefix(u::Unitful.FreeUnits{N,D,A},e) where {N,D,A}
    rem(e,3) == 0 || throw(ArgumentError("Cannot shift prefix by $(e), only multiples of 3 allowed"))
    U = first(N)
    Ũ = typeof(U)(U.tens + e, U.power)
    Unitful.FreeUnits{(Ũ,N[2:end]...),D,A}()
end

"""
    SI_prefix_convert(q)

Shift the SI prefix of the quantity `q` such that no zeros appear
before the decimal point, and at most three figures do. Note that this
conversion is not necessarily numerically accurate, it is mostly
intended for display purposes.

# Examples

```jldoctest
julia> SolverTraces.SI_prefix_convert(1.0u"nm")
1.0 nm

julia> SolverTraces.SI_prefix_convert(0.10u"nm")
100.0 pm

julia> SolverTraces.SI_prefix_convert(100.0u"nm")
100.0 nm

julia> SolverTraces.SI_prefix_convert(999.0u"nm")
999.0 nm
```
"""
function SI_prefix_convert(q::Quantity)
    v = ustrip(q)
    v == 0.0 && return q
    e = 3*(abs(v) > 1 ? floor(Int, log10(abs(v))/3) : -ceil(Int,log10(1/abs(v))/3))
    (v/10.0^e)*shift_prefix(unit(q), e)
end


text:
export shortcut, adaptive_shortcut, adaptive_shortcut!
export smooth_solution!, time_discretize_solution!, time_space_solution!

### ADAPTIVE-SHORTCUT (Hsu 2000)

function shortcut(path::Path, CC::CollisionChecker)
    N = length(path)
    if N == 2
        return path
    end
    if is_free_motion(path[1], path[end], CC)
        return path[[1,end]]
    end
    mid = ceil(Int, N/2)
    return [shortcut(path[1:mid], CC)[1:end-1]; shortcut(path[mid:end], CC)]
end

function cut_corner(v1::State, v2::State, v3::State, CC::CollisionChecker)
    m1 = (v1 + v2)/2
    m2 = (v3 + v2)/2
    while !is_free_motion(m1, m2, CC)
        m1 = (m1 + v2)/2
        m2 = (m2 + v2)/2
    end
    return typeof(v1)[v1, m1, m2, v3]
end

function adaptive_shortcut(path::Path, CC::CollisionChecker, iterations::Int = 10)
    while (short_path = shortcut(path, CC)) != path
        path = short_path
    end
    for i in 1:iterations
        path = [path[1:1]; vcat([cut_corner(path[j-1:j+1]..., CC)[2:3] for j in 2:length(path)-1]...); path[end:end]]
        while (short_path = shortcut(path, CC)) != path
            path = short_path
        end
    end
    return path, cumsum([0; map(norm, diff(path))])
end

function adaptive_shortcut!(P::MPProblem, iterations::Int = 10)
    P.status == :solved || error("Cannot post-process unsolved problem! (adaptive-shortcut)")
    isa(P.SS.dist, Euclidean) || error("Adaptive-shortcut requires Euclidean SS")
    S = P.solution
    smoothed_path, smoothed_cumcost = adaptive_shortcut(P.V.V[S.metadata["path"]], P.CC, iterations)
    S.metadata["smoothed_path"] = smoothed_path
    S.metadata["smoothed_cumcost"] = smoothed_cumcost
    S.metadata["smoothed_cost"] = smoothed_cumcost[end]
    smoothed_cumcost[end]
end

### Smoothing

function smooth_solution!(P::MPProblem)
    P.status != :solved && error("Cannot post-process unsolved problem! (adaptive-shortcut)")
    isa(P.SS.dist, Euclidean) && return adaptive_shortcut!(P)
end

### Path discretization

function time_discretize_solution!(P::MPProblem, dt::Real, usesmoothed = true)
    if usesmoothed && haskey(P.solution.metadata, "smoothed_path")
        path = P.solution.metadata["smoothed_path"]
        costs = P.solution.metadata["smoothed_cumcost"]
    else
        path = P.V[P.solution.metadata["path"]]
        costs = P.solution.metadata["cumcost"]
    end
    controlseq = steering_control(P.SS, path...)
    P.solution.metadata["discretized_path"] = propagate(P.SS, path[1], controlseq, 0:dt:duration(controlseq))
end

function time_space_solution!(P::MPProblem, n::Integer, usesmoothed = true)
    if usesmoothed && haskey(P.solution.metadata, "smoothed_path")
        path = P.solution.metadata["smoothed_path"]
        costs = P.solution.metadata["smoothed_cumcost"]
    else
        path = P.V[P.solution.metadata["path"]]
        costs = P.solution.metadata["cumcost"]
    end
    controlseq = steering_control(P.SS, path...)
    P.solution.metadata["discretized_path"] = propagate(P.SS, path[1], controlseq, linspace(0., duration(controlseq), n))
end


text:
export ContinuousScalar
export ContinuousEuclid
export Circular, Circle

Base.convert(::Type{<:Tuple}, ::InstanceType{Manifolds.Euclidean{Tuple{N}, ℝ}} ) where N = tuple([:Euclid for i in 1:N]...)
Base.convert(::Type{<:Tuple}, ::InstanceType{Manifolds.Circle{ℝ}})  = error("#FIXME")#(:Circular,)
Base.convert(::Type{<:Tuple}, ::InstanceType{Manifolds.RealCircleGroup})  = (:Circular,)

## Euclid 1

"""
$(TYPEDEF)

Most basic continuous scalar variable in a `::DFG.AbstractDFG` object.

DevNotes
- TODO Consolidate with ContinuousEuclid{1}
"""
@defVariable ContinuousScalar TranslationGroup(1) [0.0;]

"""
    ContinuousEuclid{N}
Continuous Euclidean variable of dimension `N`.
"""
struct ContinuousEuclid{N} <: InferenceVariable end

ContinuousEuclid(x::Int) = ContinuousEuclid{x}()

# not sure if these overloads are necessary since DFG 775?
DFG.getManifold(::InstanceType{ContinuousEuclid{N}}) where N = TranslationGroup(N)
DFG.getDimension(val::InstanceType{ContinuousEuclid{N}}) where N = manifold_dimension(getManifold(val))

DFG.getPointType(::Type{ContinuousEuclid{N}}) where N = Vector{Float64}
DFG.getPointIdentity(M_::Type{ContinuousEuclid{N}}) where N = zeros(N) # identity_element(getManifold(M_), zeros(N))

Base.convert(::Type{<:ManifoldsBase.AbstractManifold}, ::InstanceType{ContinuousEuclid{N}}) where N = TranslationGroup(N)

## Circular

"""
$(TYPEDEF)

Circular is a `Manifolds.Circle{ℝ}` mechanization of one rotation, with `theta in [-pi,pi)`.
"""
@defVariable Circular RealCircleGroup() [0.0;]

#


text:
using PyCall
const PACKAGES = ["scipy"]

try
    pip = pyimport("pip")
catch
    get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
    download("https://bootstrap.pypa.io/get-pip.py", get_pip)
    run(`$(PyCall.python) $get_pip --user`)
end

run(`$(PyCall.python) -m pip install --user --upgrade pip setuptools`)
run(`$(PyCall.python) -m pip install --user $(PACKAGES)`)


text:
export lyapunovs, lyapunov
#######################################################################################
#                                      Discrete                                       #
#######################################################################################
"""
```julia
lyapunovs(ds::DynamicalSystem, N; kwargs...) -> [λ1, λ2, ..., λD]
```
Calculate the spectrum of lyapunov exponents [1] of `ds` by applying the
QR-decomposition method `N` times (see method "H2" of [2]).
Returns a vector with the *final*
values of the lyapunov exponents in descending order.
# Keyword Arguments:
* `Ttr` : Extra "transient" time to evolve the system before application of the
  algorithm. Should be `Int` for discrete systems. Defaults are
  system type dependent.
* `dt = 1.0` : (only for continuous) Time of individual evolutions
  between sucessive orthonormalization steps.
* `diff_eq_kwargs = Dict()` : (only for continuous)
  Keyword arguments passed into the solvers of the
  `DifferentialEquations` package (see `evolve` or `timeseries` for more info).

[1] : A. M. Lyapunov, *The General Problem of the Stability of Motion*,
Taylor & Francis (1992)

[2] : K. Geist *et al*, Progr. Theor. Phys. **83**, pp 875 (1990)
"""
function lyapunovs(ds::DiscreteDS, N::Real; Ttr::Real = 100)

  u = deepcopy(ds.state)
  D = length(u)
  eom = ds.eom
  jac = ds.jacob
  # Transient iterations
  for i in 1:Ttr
    u = eom(u)
  end

  # Initialization
  λ = zeros(eltype(u), D)
  Q = eye(eltype(u), D)
  K = copy(Q)
  # Main algorithm
  for i in 1:N
    u = eom(u)
    A_mul_B!(K, jac(u), Q)

    Q, R = qr(K)
    for i in 1:D
      λ[i] += log(abs(R[i,i]))
    end
  end
  λ./N
end# this works. add tests and good to go.

"""
```julia
lyapunov(ds::DynamicalSystem, Τ; kwargs...) -> λ
```
Calculate the maximum lyapunov exponent using a method due to Benettin [1], which simply
evolves two neighboring trajectories while constantly rescaling one of the two.
`T`  denotes the total time of evolution (should be `Int` for discrete systems).

# Keyword Arguments:
* `Ttr` : Extra "transient" time to evolve the system before application of the
  algorithm. Should be `Int` for discrete systems. Defaults are
  system type dependent.
* `d0 = 1e-9` : Initial & rescaling distance between two neighboring trajectories.
* `threshold = 10^3*d0` : Threshold to rescale the test trajectory.
* `diff_eq_kwargs = Dict()` : (only for continuous)
  Keyword arguments passed into the solvers of the
  `DifferentialEquations` package (see `evolve` or `timeseries` for more info).
* `dt = 0.1` : (only for continuous) Time of evolution between each check of
  distance e0xceeding the `threshold`.

*Warning*: Default values have been choosen to give accurate & fast results for
maximum lyapunov exponent expected between 0.1 to 1.0. Be sure to adjust
them properly for your system. Specifically for the continuous systems,
be sure that `exp(λ*dt) < threshold/d0`.

[1] : G. Benettin *et al.*, Phys. Rev. A **14**, pp 2338 (1976)
"""
function lyapunov(ds::DiscreteDS, N::Real = 100000; Ttr::Int = 100,
  d0=1e-9*one(eltype(ds.state)), threshold=10^4*d0)

  threshold <= d0 && throw(ArgumentError("Threshold must be bigger than d0!"))
  eom = ds.eom
  st1 = deepcopy(ds.state)

  # transient
  for i in 1:Ttr
    st1 = eom(st1)
  end

  st2 = st1 + d0
  dist = d0
  λ = zero(eltype(st1))
  i = 0
  while i < N
    #evolve until rescaling:
    while dist < threshold
      st1 = eom(st1)
      st2 = eom(st2)
      # if hasnan
      dist = norm(st1 - st2)
      i+=1
      i>=N && break
    end
    a = dist/d0
    λ += log(a)
    #rescale:
    st2 = st1 + (st2 - st1)/a
    dist = d0
  end
  λ /= i
end

function lyapunovs(ds::DiscreteDS1D, N::Real = 10000; Ttr = 100)

  eom = ds.eom
  der = ds.deriv
  x = deepcopy(ds.state)

  #transient
  for i in 1:Ttr
    x = eom(x)
  end

  λ = log(abs(der(x)))
  for i in 1:N
    x = eom(x)
    λ += log(abs(der(x)))
  end
  λ/N
end
lyapunov(ds::DiscreteDS1D, N::Int=10000; Ttr::Int = 100) = lyapunovs(ds, N, Ttr=Ttr)

#######################################################################################
#                                    Continuous                                       #
#######################################################################################
function lyapunov(ds::ContinuousDS, T::Real = 10000.0; Ttr = 10.0,
  d0=1e-9*one(eltype(ds.state)), threshold=10^4*d0, dt = 0.1,
  diff_eq_kwargs = Dict())

  const dict = Dict()

  threshold <= d0 && throw(ArgumentError("Threshold must be bigger than d0!"))

  # Transient
  if Ttr != 0
    evolve!(ds, Ttr; diff_eq_kwargs = diff_eq_kwargs)
  end
  # initialize:
  st1 = ds.state
  st2 = st1 + d0
  prob1 = ODEProblem(ds, dt)
  prob2 = ODEProblem(ds, dt)
  prob2.u0 = st2
  dist = d0
  λ = zero(eltype(st1))
  t = zero(T)
  i = 0
  # start evolution and rescaling:
  while t < T
    #evolve until rescaling:
    while dist < threshold
      # evolve
      if diff_eq_kwargs == dict
        st1 = evolve!(prob1)
        st2 = evolve!(prob2)
      else
        st1 = evolve!(prob1, diff_eq_kwargs)
        st2 = evolve!(prob2, diff_eq_kwargs)
      end
      dist = norm(st1 - st2)
      t += dt
      i += 1
      t ≥ T && break
    end
    # add computed scale to accumulator:
    a = dist/d0
    if a > 1e^4 && i <= 1
      warnstr = "Distance between test and original trajectory exceeded threshold "
      warnstr*= "after just 1 evolution step. "
      warnstr*= "Please decrease `dt`, increase `threshold` or decrease `d0`."
      warn(warnstr)
      errorstr = "Parameters choosen for `lyapunov` with "
      errorstr*= "`ContinuousDS` are not accurate."
      throw(ArgumentError(errorstr))
    end
    λ += log(a); i+= 1
    #rescale:
    st2 = st1 + (st2 - st1)/a
    prob2.u0 = st2
    dist = d0; i = 0
  end
  λ /= t
end

function lyapunovs(ds::ContinuousDS, N::Real;
  Ttr::Real = 1.0, diff_eq_kwargs::Dict = Dict(), dt::Real = 1.0)

  D = dimension(ds)
  T = eltype(ds.state)
  jac = ds.jacob
  const dict = Dict()
  # Transient
  if Ttr != 0
    evolve!(ds, Ttr; diff_eq_kwargs = diff_eq_kwargs)
  end
  # Initialization
  tbprob = tangentbundle_setup(ds, dt)
  λ = zeros(T, D)
  Q = eye(eltype(ds.state), D)

  # Main algorithm
  for i in 1:N
    tbprob.u0[:, 2:end] .= Q #update tangent dynamics state
    # evolve system and tangent dynamics:
    if diff_eq_kwargs == dict
      evolve!(tbprob)
    else
      evolve!(tbprob, diff_eq_kwargs)
    end

    # Perform QR:
    Q, R = qr(view(tbprob.u0, :, 2:D+1))
    for j in 1:D
      λ[j] += log(abs(R[j,j]))
    end
  end
  λ./(N*dt) #return spectrum
end


text:
function VectorialStorelessStatistic(arg0::jint, arg1::StorelessUnivariateStatistic)
    return VectorialStorelessStatistic((jint, StorelessUnivariateStatistic), arg0, arg1)
end

function clear(obj::VectorialStorelessStatistic)
    return jcall(obj, "clear", void, ())
end

function equals(obj::VectorialStorelessStatistic, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_dimension(obj::VectorialStorelessStatistic)
    return jcall(obj, "getDimension", jint, ())
end

function get_n(obj::VectorialStorelessStatistic)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::VectorialStorelessStatistic)
    return jcall(obj, "getResult", Vector{jdouble}, ())
end

function hash_code(obj::VectorialStorelessStatistic)
    return jcall(obj, "hashCode", jint, ())
end

function increment(obj::VectorialStorelessStatistic, arg0::Vector{jdouble})
    return jcall(obj, "increment", void, (Vector{jdouble},), arg0)
end


text:
# using DataFrames, CSV

# include("../src/LinearRegression.jl")

@testset "Cook's Distance" begin
        st_df = DataFrame(
                Y=[6.4, 7.4, 10.4, 15.1, 12.3 , 11.4],
                XA=[1.5, 6.5, 11.5, 19.9, 17.0, 15.5],
                XB=[1.8, 7.8, 11.8, 20.5, 17.3, 15.8],
                XC=[3., 13., 23., 39.8, 34., 31.],
                # values from SAS proc reg
                CooksD_base=[1.4068501943, 0.176809102, 0.0026655177, 1.0704009915, 0.0875726457, 0.1331183932],
                CooksD_multi=[1.7122291956, 18.983407026, 0.000118078, 0.8470797843, 0.0715921999, 0.1105843157],
                )

        t_lm_base = regress(@formula(Y ~ 1+ XA), st_df)
        results = predict_in_sample(t_lm_base, st_df, req_stats=["all"])
        @test isapprox(st_df.CooksD_base, results.cooksd)

        t_lm_multi = regress(@formula(Y ~ 1+ XA + XB), st_df)
        results = predict_in_sample(t_lm_multi, st_df, req_stats=["all"])
        @test isapprox(st_df.CooksD_multi, results.cooksd)
end


text:
# Geary's c test of Global Spatial Autocorrelation
struct GlobalGeary <: AbstractGlobalSpatialAutocorrelation
    n::Int
    C::Float64
    EC::Float64
    p::Float64
    Cperms::Vector{Float64}
    Cpermsmean::Float64
    Cpermsstd::Float64
    z::Float64
end

"""
    geary(x, W)
Compute the Geary's c test of global spatial autocorrelation.

# Optional Arguments
- `permutations=9999`: number of permutations for the randomization test.
- `rng=default_rng()`: random number generator for the randomization test.
"""
function geary(x::AbstractVector{T} where T<:Real, W::SpatialWeights; permutations::Int = 9999,
    rng::AbstractRNG = default_rng())::GlobalGeary

    n = length(x)
    z = x
    S0 = sum(sum.(W.weights::Vector{Vector{Float64}}))
    denominator = sum( (z .- mean(z)).^2 )

    # Auxiliar function to calculate Geary's C
    function geary_calc(z::AbstractVector, W::SpatialWeights)::Float64
        numerator = 0.0
        for i in 1:n
            numerator += sum( weights(W, i) .* (z[i] .- z[neighbors(W, i)]).^2 )
        end
        return (n - 1)/(2 * S0) * numerator / denominator
    end

    # Gearys's C
    C = geary_calc(z, W)
    EC = 1.0

    # Build permutations array
    Zi = Array{Float64}(undef, n, permutations)
    for i in 1:permutations
        Zi[:,i] = shuffle(rng, z)
    end

    # Calculate Geary's c for all the permutations
    Cperms = zeros(permutations)
    Threads.@threads for i in 1:permutations
        Cperms[i] = geary_calc(view(Zi, :, i), W)
    end

    larger = sum(Cperms .>= C)
    if (permutations - larger) < larger
        larger = permutations - larger
    end
    p = (larger + 1) / (permutations + 1)

    Cpermsstd = std(Cperms, corrected = false)
    Cpermsmean = mean(Cperms)
    z = (C - Cpermsmean)  / Cpermsstd

    return GlobalGeary(n, C, EC, p, Cperms, Cpermsmean, Cpermsstd, z)

end

score(x::GlobalGeary) = x.C;

scoreperms(x::GlobalGeary) = x.Cperms;

expected(x::GlobalGeary) = x.EC;

mean(x::GlobalGeary) = x.Cpermsmean;

std(x::GlobalGeary) = x.Cpermsstd;

zscore(x::GlobalGeary) = x.z;

pvalue(x::GlobalGeary) = x.p;

testname(x::GlobalGeary) = "Geary's c";


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle(zeros(1006), [fill(0.0001, 400); zeros(606)]),
    "info" => Dict(
        "name" => "fom",
        "dim" => 1006,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.TimeSteppers: time_step!
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using OceanTurbulenceParameterEstimation.Observations: FieldTimeSeriesCollector, IdentityNormalization, SyntheticObservations
using OceanTurbulenceParameterEstimation.InverseProblems: ConcatenatedOutputMap, transform_time_series

Ex, Ey = (2, 2)
sz = ColumnEnsembleSize(Nz = 8, ensemble = (Ex, Ey))
halo = ColumnEnsembleSize(Nz = sz.Nz)

grid = RectilinearGrid(size = sz, halo = halo, z = (-128, 0), topology = (Flat, Flat, Bounded))

closure = [ConvectiveAdjustmentVerticalDiffusivity() for i = 1:Ex, j = 1:Ey]

Qᵇ = [+1e-8 for i = 1:Ex, j = 1:Ey]
Qᵘ = [-1e-4 for i = 1:Ex, j = 1:Ey]
Qᵛ = [0.0 for i = 1:Ex, j = 1:Ey]

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

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

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

Δt = 10.0
stop_iteration = 3
times = [Δt * (n - 1) for n = 1:stop_iteration]

field_names = (:u,)

simulation = Simulation(model; Δt = Δt, stop_iteration = stop_iteration)

simulation_fields = fields(simulation.model)

fields_to_collect = NamedTuple(name => simulation_fields[name] for name in field_names)

time_series_collector = FieldTimeSeriesCollector(fields_to_collect, times)

function initialize!(simulation, time_series_collector)

    model_fields = fields(simulation.model)
    for name in keys(model_fields)
        model_field = model_fields[name]
        set!(model_field, 0.0)
    end

    set!(simulation.model, b = bᵢ)

    simulation.model.clock.time = 0.0
    simulation.model.clock.iteration = 0

    for time_series in time_series_collector.field_time_serieses
        time_series.data .= 0
    end

    simulation.callbacks[:data_collector] = Callback(time_series_collector, SpecifiedTimes(times...))

    simulation.stop_time = times[end]

    return nothing
end

function forward_map(simulation, time_series_collector)

    initialize!(simulation, time_series_collector)
    run!(simulation)

    map = ConcatenatedOutputMap()

    normalization = Dict(name => IdentityNormalization() for name in field_names)

    transposed_output = SyntheticObservations(time_series_collector.field_time_serieses,
        time_series_collector.grid,
        time_series_collector.times,
        nothing,
        nothing,
        normalization)

    output = transform_time_series(map, transposed_output)

    return output
end

x = forward_map(simulation, time_series_collector)[:, 1:1]

y = forward_map(simulation, time_series_collector)[:, 1:1]

@show x == y


text:
"""
Package for likelihood-free inference. Includes Bayesian methods `PMC`, `AdaptiveSMC`, `KernelABC` and `K2ABC`, and different point estimators.
"""
module LikelihoodfreeInference
using Distributions, Random, LinearAlgebra, StatsBase, Dates, PDMats, Roots, Requires, NLopt, Printf, DiffResults, ForwardDiff

export PMC, run!, EpsilonExponentialDecay, EpsilonLinearDecay, MultivariateUniform, AdaptiveSMC, MultiParticleNormal, TruncatedMultivariateNormal, defaultproposal, weights, particles, mean, PointEstimator, QDLoss, KernelLoss, K2ABC, Kernel, Gaussian, ModifiedGaussian, MedianHeuristic, ScottsHeuristic, Bandwidth, Smoothed, Scheduled, StandardMMD, KernelABC, KernelRecursiveABC, LogGaussianKernel, euclidean, energydistance, EnergyLoss, MMDLoss, KLLoss, kldistance, mmd

include("distributions.jl")
include("utils.jl")
include("kernels.jl")
include("smc.jl")
include("pointestimator.jl")
include("kernelabc.jl")
include("krabc.jl")
include("optimization.jl")

function __init__()
    include(joinpath(@__DIR__, "optional_optimization.jl"))
    @require StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" begin
        function StatsPlots.corrplot(x::AbstractABC;
                                     N = nothing, kwargs...)
            p = particles(x)
            w = weights(x)
            N = N === nothing ? min(1000, 10*length(p)) : N
            M = vcat([hcat(resample(p, w)...)' for _ in 1:div(N, length(p))]...)
            StatsPlots.corrplot(M; kwargs...)
        end
        function StatsPlots.histogram(x::AbstractABC;
                                      layout = nothing,
                                      kwargs...)
            d = length(particles(x)[1])
            n = ceil(Int, sqrt(d))
            layout = layout === nothing ? (n, ceil(Int,d/n)) : layout
            StatsPlots.histogram(hcat(particles(x)...)', weights = weights(x);
                                 layout = layout, legend = false, kwargs...)

        end
    end
end

end # module


text:
using ModelingToolkit, DiffEqBase, DiffEqJump, Test, LinearAlgebra
MT = ModelingToolkit

# basic MT SIR model with tweaks
@parameters β γ t
@variables S(t) I(t) R(t)
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I+t
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = VariableRateJump(rate₂,affect₂)
js      = JumpSystem([j₁,j₂], t, [S,I,R], [β,γ])
statetoid = Dict(MT.value(state) => i for (i,state) in enumerate(states(js)))
mtjump1 = MT.assemble_crj(js, j₁, statetoid)
mtjump2 = MT.assemble_vrj(js, j₂, statetoid)

# doc version
rate1(u,p,t) = (0.1/1000.0)*u[1]*u[2]
function affect1!(integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump1 = ConstantRateJump(rate1,affect1!)
rate2(u,p,t) = 0.01u[2]+t
function affect2!(integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = VariableRateJump(rate2,affect2!)

# test crjs
u = [100, 9, 5]
p = (0.1/1000,0.01)
tf = 1.0
mutable struct TestInt{U,V,T}
    u::U
    p::V
    t::T
end
mtintegrator = TestInt(u,p,tf)
integrator   = TestInt(u,p,tf)
@test abs(mtjump1.rate(u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
@test abs(mtjump2.rate(u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
mtjump1.affect!(mtintegrator)
jump1.affect!(integrator)
@test all(integrator.u .== mtintegrator.u)
mtintegrator.u .= u; integrator.u .= u
mtjump2.affect!(mtintegrator)
jump2.affect!(integrator)
@test all(integrator.u .== mtintegrator.u)

# test MT can make and solve a jump problem
rate₃   = γ*I
affect₃ = [I ~ I - 1, R ~ R + 1]
j₃ = ConstantRateJump(rate₃,affect₃)
js2 = JumpSystem([j₁,j₃], t, [S,I,R], [β,γ])
u₀ = [999,1,0]; p = (0.1/1000,0.01); tspan = (0.,250.)
u₀map = [S => 999, I => 1, R => 0]
parammap = [β => .1/1000, γ => .01]
dprob = DiscreteProblem(js2, u₀map, tspan, parammap)
jprob = JumpProblem(js2, dprob, Direct(), save_positions=(false,false))
Nsims = 30000
function getmean(jprob,Nsims)
  m = 0.0
  for i = 1:Nsims
    sol = solve(jprob, SSAStepper())
    m += sol[end,end]
  end
  m/Nsims
end
m = getmean(jprob,Nsims)

# test save_positions is working
jprob = JumpProblem(js2, dprob, Direct(), save_positions=(false,false))
sol = solve(jprob, SSAStepper(), saveat=1.0)
@test all((sol.t) .== collect(0.:tspan[2]))

#test the MT JumpProblem rates/affects are correct
rate2(u,p,t) = 0.01u[2]
jump2 = ConstantRateJump(rate2,affect2!)
mtjumps = jprob.discrete_jump_aggregation
@test abs(mtjumps.rates[1](u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
@test abs(mtjumps.rates[2](u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
mtjumps.affects![1](mtintegrator)
jump1.affect!(integrator)
@test all(integrator.u .== mtintegrator.u)
mtintegrator.u .= u; integrator.u .= u
mtjumps.affects![2](mtintegrator)
jump2.affect!(integrator)
@test all(integrator.u .== mtintegrator.u)

# direct vers
p = (0.1/1000,0.01)
prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
r1(u,p,t) = (0.1/1000.0)*u[1]*u[2]
function a1!(integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
j1 = ConstantRateJump(r1,a1!)
r2(u,p,t) = 0.01u[2]
function a2!(integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
j2 = ConstantRateJump(r2,a2!)
jset = JumpSet((),(j1,j2),nothing,nothing)
jprob = JumpProblem(prob,Direct(),jset, save_positions=(false,false))
m2 = getmean(jprob,Nsims)

# test JumpSystem solution agrees with direct version
@test abs(m-m2)/m < .01

# mass action jump tests for SIR model
maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
js3   = JumpSystem([maj1,maj2], t, [S,I,R], [β,γ])
statetoid = Dict(MT.value(state) => i for (i,state) in enumerate(states(js)))
ptoid     = Dict(MT.value(par) => i for (i,par) in enumerate(parameters(js)))
dprob = DiscreteProblem(js3, u₀map, tspan, parammap)
jprob = JumpProblem(js3, dprob, Direct())
m3 = getmean(jprob,Nsims)
@test abs(m-m3)/m < .01

# maj jump test with various dep graphs
js3b = JumpSystem([maj1,maj2], t, [S,I,R], [β,γ])
jprobb = JumpProblem(js3b, dprob, NRM())
m4 = getmean(jprobb,Nsims)
@test abs(m-m4)/m < .01
jprobc = JumpProblem(js3b, dprob, RSSA())
m4 = getmean(jprobc,Nsims)
@test abs(m-m4)/m < .01

# mass action jump tests for other reaction types (zero order, decay)
maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
maj2 = MassActionJump(γ, [S => 1], [S => -1])
js4   = JumpSystem([maj1,maj2], t, [S], [β,γ])
statetoid = Dict(MT.value(state) => i for (i,state) in enumerate(states(js)))
ptoid     = Dict(MT.value(par) => i for (i,par) in enumerate(parameters(js)))
dprob = DiscreteProblem(js4, [S => 999], (0,1000.), [β => 100.,γ => .01])
jprob = JumpProblem(js4, dprob, Direct())
m4 = getmean(jprob,Nsims)
@test abs(m4 - 2.0/.01)*.01/2.0 < .01

# test second order rx runs
maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
maj2 = MassActionJump(γ, [S => 2], [S => -1])
js4   = JumpSystem([maj1,maj2], t, [S], [β,γ])
statetoid = Dict(MT.value(state) => i for (i,state) in enumerate(states(js)))
ptoid     = Dict(MT.value(par) => i for (i,par) in enumerate(parameters(js)))
dprob = DiscreteProblem(js4, [S => 999], (0,1000.), [β => 100.,γ => .01])
jprob = JumpProblem(js4, dprob, Direct())
sol = solve(jprob, SSAStepper());

# issue #819
@testset "Combined system name collisions" begin
  sys1 = JumpSystem([maj1, maj2], t, [S], [β, γ], name = :sys1)
  sys2 = JumpSystem([maj1, maj2], t, [S], [β, γ], name = :sys1)
  @test_throws ArgumentError JumpSystem([sys1.γ ~ sys2.γ], t, [], [], systems = [sys1, sys2])
end


text:
#!/usr/bin/env julia

using HypertextLiteral
using NarrativeTest

subs = NarrativeTest.common_subs()

# Ignore the difference in the output of `print(Int)` between 32-bit and 64-bit platforms.
push!(subs, r"Int64" => s"Int(32|64)")

# Normalize printing of vector types.
if VERSION < v"1.6.0-DEV"
    Base.show_datatype(io::IO, x::Type{Vector{T}}) where {T} = print(io, "Vector{$T}")
end

# Normalize printing of type parameters.
if VERSION < v"1.6.0-DEV"
    function Base.show_datatype(io::IO, x::DataType)
        istuple = x.name === Tuple.name
        if (!isempty(x.parameters) || istuple) && x !== Tuple
            n = length(x.parameters)::Int
            if istuple && n > 3 && all(i -> (x.parameters[1] === i), x.parameters)
                print(io, "NTuple{", n, ", ", x.parameters[1], "}")
            else
                Base.show_type_name(io, x.name)
                print(io, '{')
                for (i, p) in enumerate(x.parameters)
                    show(io, p)
                    i < n && print(io, ", ")
                end
                print(io, '}')
            end
        else
            Base.show_type_name(io, x.name)
        end
    end
end

# Ignore line ending differences for Windows targets.
push!(subs, r"\r\n" => "\n")

# Set the width to 72 so that MD->PDF via pandoc fits the page.
ENV["COLUMNS"] = "72"

package_path(x) = relpath(joinpath(dirname(abspath(PROGRAM_FILE)), "..", x))
default = package_path.(["README.md", "docs/src"])
runtests(; default=default, subs=subs)


text:
using PortAudio, SampledSignals, DSP, Colors, Interpolations
const FFT_SCALE_DEFAULT = 1000.0
const fft_scale = Ref(FFT_SCALE_DEFAULT)
const fft_rescale_up_counter = Ref(0)
const fft_rescale_down_counter = Ref(0)
const FLOAT32_MAX = 3.4028235f38
const INT_32_MAX = 2^32-1

mutable struct AudioAnalysis
    spec_bufs::Array{AFArray{Complex{Float64}},1}
    audio_buffers::Array{AFArray{Float64},1}
    spec_buf_order::Array{Int,1}
    audio_buffer_order::Array{Int,1}
    delta_buffers::Array{AFArray{Float64},1}
    spec_obj::SpectrumBuf
    function AudioAnalysis(cpuAudio::SampleBuf, numBufs::Int)
        audio_buffers = [AFArray{Float64}(convert.(Float64, cpuAudio.data[:,1])) for i in 1:numBufs]
        spec_bufs = [fft(audio_buffers[1]) for j in 1:numBufs]
        audio_buf_order = [numBufs - (k - 1) for k in 1:numBufs]
        spec_buf_order = [numBufs - (l - 1) for l in 1:numBufs]
        delta_buffers = [abs(spec_bufs[1]) - abs(spec_bufs[end]) for m in 1:numBufs]
        spec_obj = SpectrumBuf(Array(spec_bufs[spec_buf_order[1]]), nframes(cpuAudio)/samplerate(cpuAudio))
        return new(spec_bufs, audio_buffers, spec_buf_order, audio_buf_order, delta_buffers, spec_obj)
    end
end

function rotate_array_clockwise!(arr::Array)
    if length(arr) > 1
        tmp1 = arr[1]
        tmp2 = arr[2]
        for i in eachindex(arr)
            if i != length(arr)
                tmp2 = arr[i+1]
                setindex!(arr, tmp1, i+1)
                tmp1 = tmp2
            else
                arr[1] = tmp1
            end
        end
    end
end

function push_audio!(ana::AudioAnalysis, cpuAudio::Array{Float32,1})
    ana.audio_buffers[ana.audio_buffer_order[end]] = AFArray(convert(Array{Float64}, cpuAudio))
    rotate_array_clockwise!(ana.audio_buffer_order)
end

function process_audio!(ana::AudioAnalysis, cpuAudio::SampledSignals.SampleBuf)
    push_audio!(ana, cpuAudio.data[:,1])

    # For the moment, the deltas will just share order with the spectrums, but that could be added in the future if needed.
    ana.spec_bufs[ana.spec_buf_order[end]] = fft(ana.audio_buffers[ana.audio_buffer_order[1]])
    ana.delta_buffers[ana.spec_buf_order[end]] = abs(ana.spec_bufs[ana.spec_buf_order[end]]) - abs(ana.spec_bufs[ana.spec_buf_order[1]])
    rotate_array_clockwise!(ana.spec_buf_order)

    ana.spec_obj.data .= ana.spec_bufs[ana.spec_buf_order[1]]
    ana.spec_obj.samplerate = nframes(cpuAudio)/samplerate(cpuAudio)

    handle_scaling(ana)
end

function binFFT(ana::AudioAnalysis, nbins::Int)
    rawspec = ana.spec_bufs[ana.spec_buf_order[1]]
    nbin = nbins + 1
    f(x) = 2.0^x
    dom = domain(ana.spec_obj)
    range = linspace(log(dom[2])/log(2), log(dom[end-1])/log(2), nbin)[1:nbins]
    bands = collect(abs.(f.(range))*length(dom)/dom[end])
    spec = abs(approx1(ana.spec_bufs[ana.spec_buf_order[1]], AFArray{Float64}(bands), AF_INTERP_CUBIC_SPLINE, 0.0f0))
    maxspec = maximum(spec)
    if length(spec) > nbins
        return Array(spec)[1:nbins], maxspec
    else
        return Array(spec), maxspec
    end
end

function binFFT(ana::AudioAnalysis, nbins::Int, cutoff::Int)
    rawspec = ana.spec_bufs[ana.spec_buf_order[1]][1:cutoff]
    nbin = nbins + 1
    f(x) = 2.0^x
    dom = domain(ana.spec_obj)
    range = linspace(log(dom[2])/log(2), log(dom[end-1])/log(2), nbin)[1:nbins]
    bands = collect(abs.(f.(range))*length(dom)/dom[end])
    spec = abs(approx1(ana.spec_bufs[ana.spec_buf_order[1]], AFArray{Float64}(bands), AF_INTERP_CUBIC_SPLINE, 0.0f0))
    maxspec = maximum(spec)
    if length(spec) > nbins
        return Array(spec)[1:nbins], maxspec
    else
        return Array(spec), maxspec
    end
end

@inline function hsl_to_rgb(color::ColorTypes.HSL)
    return hsl_to_rgb(color.h/360, color.s, color.l)
end

@inline function hsl_to_rgb(h::Real,s::Real,l::Real)::Array{UInt8}
    r = 0x00
    g = 0x00
    b = 0x00
    if s == 0
        r = round(UInt8, l*255)
        g = round(UInt8, l*255)
        b = round(UInt8, l*255)
    else
        q = l < 0.5 ? l * (1 + s) : l + s - l * s
        p = 2 * l - q
        r = hue2rgb(p, q, h + 1/3)*255
        g = hue2rgb(p, q, h)*255
        b = hue2rgb(p, q, h - 1/3)*255
        if isnan(r);r=0;end;
        if isnan(g);g=0;end;
        if isnan(b);b=0;end;
    end
    return round.(UInt8, [r, g, b])
end
function hue2rgb(p::Real, q::Real, t::Real)::Real
    if t < 0; t+=1; end;
    if t > 1; t-=1; end;
    if t < 1/6; return p + (q - p) * 6 * t; end;
    if t < 1/2; return q; end;
    if t < 2/3; return p + (q - p) * (2/3 - t) * 6; end;
    return p
end

function gfft(buf::SampledSignals.SampleBuf)
    return SpectrumBuf(Array(fft(ArrayFire.AFArray(buf.data))), nframes(buf)/samplerate(buf))
end

function gfft!(buf::SampledSignals.SampleBuf, gpu_buf1::AFArray, spec::SampledSignals.SpectrumBuf, gpu_buf2::AFArray)
    gpu_buf1 .= buf.data
    gpu_buf2 = fft(gpu_buf1)
    spec.data .= gpu_buf2
    spec.samplerate = nframes(buf)/samplerate(buf)
end

function handle_scaling(ana::AudioAnalysis)
    spec = abs(ana.spec_bufs[ana.spec_buf_order[1]])
    maxspec = maximum(spec)
    mean_delta = mean(abs(ana.delta_buffers[ana.spec_buf_order[1]]))
    spec_coef = (mean(spec)+maxspec+mean_delta)/3
    if maxspec > fft_scale[]*1.0
        spec.*=fft_scale[]/maxspec
        fft_scale[] = maxspec * 1.2
        fft_rescale_up_counter[] += 4
        if fft_rescale_up_counter[] > 35
            fft_scale[] *= 1.5
            fft_rescale_up_counter[] = 0
        end
    elseif spec_coef > fft_scale[]*0.50
        fft_rescale_up_counter[] += 2
        if fft_rescale_up_counter[] > 20
            fft_scale[] *= 1.25
        end
    elseif spec_coef > fft_scale[]*0.30
        fft_rescale_up_counter[] += 1
        if fft_rescale_up_counter[] > 15
            fft_scale[] *= 1.05
        end
    elseif spec_coef < fft_scale[] * 0.1 && maxspec > fft_scale[]/100
        fft_scale[] /= 1.05
        spec.*= !isnan(fft_scale[]/spec_coef * 0.25) ? fft_scale[]/spec_coef * 0.25 : fft_scale[]/0.1
        fft_rescale_down_counter[] += 4
        if fft_rescale_down_counter[] > 35
            fft_scale[] /= 1.5
            fft_rescale_down_counter[] = 0
        end
    elseif spec_coef < fft_scale[] * 0.20
        fft_rescale_down_counter[] += 1
        if fft_rescale_down_counter[] > 15
            fft_scale[] /= 1.05
        end
    end
    if maxspec < FFT_SCALE_DEFAULT/100 && fft_scale[] > FFT_SCALE_DEFAULT
        spec = zeros(AFArray{Float64}, size(spec))
        fft_scale[] = FFT_SCALE_DEFAULT
    end
    if maxspec > fft_scale[]
        spec = max(spec-fft_scale[], 0.0)+fft_scale[].*(spec.>fft_scale[])
    end
    return spec, maxspec
end

macro bounded(val, minVal, maxVal)
    return esc(:(min(max($val, $minVal), $maxVal)))
end