text:
#Struct to hold model state
include("HGrid.jl")
include("CGrid.jl")
include("VGrid.jl")
include("UGrid.jl")
include("SigmaGrid.jl")

struct Fields{T <: Real, N <: Real}
    gh::HGrid{T,N}
    gu::UGrid{T,N}
    gv::VGrid{T,N}
    gc::CGrid{T,N}
    g3d::SigmaGrid{T,N}
    wu::UWavelets{T,N}
    wv::VWavelets{T,N}
end

"""
    setup_fields(grid, initial_conditions, solver_params, params, bed_array)

Acts as a constructor for the fields (no explicit constructor as fields `only ever called when setting up a model)
"""

function setup_fields(grid, initial_conditions, solver_params, params, bed_array)
    #Define masks for points on h-, u-, v- and c-grids that lie in model domain.
    h_mask = grid.h_mask
    u_mask = get_u_mask(h_mask)
    v_mask = get_v_mask(h_mask)
    c_mask = get_c_mask(h_mask)

    #Remove all points on u- and v-grids with homogenous Dirichlet conditions.
    u_mask[grid.u_iszero].=false
    v_mask[grid.v_iszero].=false

    #h-grid
    #gh=HGrid(grid, params) #uncomment if using the explicit constructor method
    h =  deepcopy(initial_conditions.initial_thickness)
    ηav = deepcopy(initial_conditions.initial_viscosity[:,:,1]) #set to the viscosity on the first level for now
    gh=HGrid(
    nxh=grid.nx,
    nyh=grid.ny,
    mask=h_mask,
    b = bed_array,
    h = h,
    ηav = ηav,
    )

    #u-grid
    gu=UGrid(
    nxu=grid.nx+1,
    nyu=grid.ny,
    dx=grid.dx,
    dy=grid.dy,
    mask=u_mask,
    levels=solver_params.levels
    )

    #v-grid
    gv=VGrid(
    nxv=grid.nx,
    nyv=grid.ny+1,
    dx=grid.dx,
    dy=grid.dy,
    mask=v_mask,
    levels=solver_params.levels
    )

    #c-grid
    gc=CGrid(
    nxc=grid.nx-1,
    nyc=grid.ny-1,
    mask=c_mask
    )

    #3D-grid
    η = deepcopy(initial_conditions.initial_viscosity)
    θ = deepcopy(initial_conditions.initial_temperature)
    Φ = deepcopy(initial_conditions.initial_damage)
    g3d=SigmaGrid(
    nxs=grid.nx,
    nys=grid.ny,
    nσs=grid.nσ,
    σ =grid.σ,
    η = η,
    θ = θ,
    Φ = Φ,
    glen_b = glen_b.(θ,Φ,params.glen_a_ref, params.glen_n, params.glen_a_activation_energy, params.glen_temperature_ref, params.gas_const),
    quadrature_weights = grid.quadrature_weights
    )

    #Wavelet-grid, u-component.
    wu=UWavelets(nxuw=grid.nx+1,nyuw=grid.ny,levels=solver_params.levels)

    #Wavelet-grid, v-component.
    wv=VWavelets(nxvw=grid.nx,nyvw=grid.ny+1,levels=solver_params.levels)
    return Fields(gh,gu,gv,gc,g3d,wu,wv)
end


text:
refVals = []
refPrecs = []

#! format: off
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_update/test/Ocean/SplitExplicit/simple_box_2dt.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.65970172440775332046e-01,  1.61603740424529351838e-01,  4.57235954923597420763e-03,  2.21382850057524928344e-02 ],
 [  "oce Q_3D",   "u[2]", -2.16238066786232502325e-01,  2.44156974607258908661e-01, -2.57322116506966246802e-03,  2.49380331699317753236e-02 ],
 [  "oce Q_3D",      "η", -7.98711302045854054654e-01,  3.47810107484795683064e-01, -2.92595480887335743312e-04,  2.91360229800833869795e-01 ],
 [  "oce Q_3D",      "θ",  3.55385987992295690474e-03,  9.91576780951165659417e+00,  2.49883145933055006438e+00,  2.17846939512263038097e+00 ],
 [   "oce aux",      "w", -1.66960462888799841472e-04,  1.65914321067276468273e-04,  5.22671200400887584530e-07,  1.42481192436579197488e-05 ],
 [   "oce aux",   "pkin", -8.84676906676938301644e+00,  0.00000000000000000000e+00, -3.26734804711110893294e+00,  2.49685702267286213640e+00 ],
 [   "oce aux",    "wz0", -2.12828008736777837579e-05,  3.08227212254399767725e-05, -1.55080158053877423592e-10,  7.82426303256959697840e-06 ],
 [   "oce aux", "u_d[1]", -1.56747674354107358052e-01,  1.23337664405837849069e-01, -7.85428994855595028748e-05,  1.15205292893568473495e-02 ],
 [   "oce aux", "u_d[2]", -2.12928714468491042666e-01,  2.30595749843429731474e-01,  2.75559843466661209680e-05,  1.75135782183724539318e-02 ],
 [   "oce aux", "ΔGu[1]", -2.33711927426000212944e-06,  3.31423592619396043514e-06, -4.04240250953691765080e-08,  3.37695011150041160627e-07 ],
 [   "oce aux", "ΔGu[2]", -2.38751859804088027546e-06,  2.23171788862484417825e-06,  1.26953194368922374152e-06,  7.64696321260357341382e-07 ],
 [   "oce aux",      "y",  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573163901915703900e+06 ],
 [ "baro Q_2D",   "U[1]", -1.77749182159153633620e+01,  6.89466489062268976795e+01,  4.64951906527867731000e+00,  1.81233032109943827948e+01 ],
 [ "baro Q_2D",   "U[2]", -3.48869693974309527107e+01,  8.56944551257798678989e+01, -2.60002834745782607229e+00,  1.80715674809898771969e+01 ],
 [ "baro Q_2D",      "η", -7.98718706462308802863e-01,  3.47811046235604659493e-01, -2.92587313874264492736e-04,  2.91375274868600042666e-01 ],
 [  "baro aux",  "Gᵁ[1]", -3.31423592619396051306e-03,  2.33711927426000221075e-03,  4.04240250953691941898e-05,  3.37711728311069743838e-04 ],
 [  "baro aux",  "Gᵁ[2]", -2.23171788862484409693e-03,  2.38751859804088020431e-03, -1.26953194368922394480e-03,  7.64734176576895561574e-04 ],
 [  "baro aux", "U_c[1]", -1.77747948208394639380e+01,  6.89456621278419703458e+01,  4.64949958798823903550e+00,  1.81231282979409442646e+01 ],
 [  "baro aux", "U_c[2]", -3.48859937354993405734e+01,  8.56944690523459655651e+01, -2.59996531383233087098e+00,  1.80715328159406638520e+01 ],
 [  "baro aux",    "η_c", -7.98711302045854054654e-01,  3.47810107484795683064e-01, -2.92595480887348482688e-04,  2.91374653217579937525e-01 ],
 [  "baro aux", "U_s[1]", -1.77749182159153633620e+01,  6.89466489062268976795e+01,  4.64951906527867731000e+00,  1.81233032109943827948e+01 ],
 [  "baro aux", "U_s[2]", -3.48869693974309527107e+01,  8.56944551257798678989e+01, -2.60002834745782607229e+00,  1.80715674809898771969e+01 ],
 [  "baro aux",    "η_s", -7.98718706462308802863e-01,  3.47811046235604659493e-01, -2.92587313874264492736e-04,  2.91375274868600042666e-01 ],
 [  "baro aux",  "Δu[1]", -4.52063405565574774267e-04,  3.99219316002806577735e-04, -1.09219330063404869891e-05,  1.02341200736279461128e-04 ],
 [  "baro aux",  "Δu[2]", -2.42431689759263984492e-04,  3.83040345117153786629e-04,  5.18877254212526286387e-05,  7.47027376508314802477e-05 ],
 [  "baro aux", "η_diag", -7.98659935505348972384e-01,  3.47918275378131136577e-01, -2.91845106891554467807e-04,  2.91369719089321743688e-01 ],
 [  "baro aux",     "Δη", -5.74411984127720653959e-04,  7.42725686140199847785e-04, -7.50373995767764799686e-07,  1.62142594044870064609e-04 ],
 [  "baro aux",      "y",  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15578885204060329124e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",      "η",    12,    12,     8,    12 ],
 [  "oce Q_3D",      "θ",    12,    12,    12,    12 ],
 [   "oce aux",      "w",    12,    12,     8,    12 ],
 [   "oce aux",   "pkin",    12,    12,    12,    12 ],
 [   "oce aux",    "wz0",    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",      "y",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",      "η",    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",    "η_c",    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",    "η_s",    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux", "η_diag",    12,    12,     8,    12 ],
 [  "baro aux",     "Δη",     9,     9,     6,    10 ],
 [  "baro aux",      "y",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

ext: jl


text:
using Documenter, ThreadedIterables;

makedocs(
	format = Documenter.HTML(),
	sitename = "ThreadedIterables.jl",
	modules = [ThreadedIterables],
	authors = "Marek Dědič",
	pages = Any[
		"Home" => "index.md",
		"Getting started" => "getting_started.md",
		"Reference" => Any[
			"reference/foreach.md",
			"reference/map.md",
			"reference/reduce.md",
			"reference/mapreduce.md",
			"reference/mapfoldl.md",
			"reference/mapfoldr.md"
		]
	]
);


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function ticpos(deglen::T, pixlen::T, ticsize::T) where {T<:AbstractFloat}
    minpix = deglen*60/pixlen
    incr = minpix*ticsize

    if incr >= 30
        units = "Degrees"
    elseif incr > 0.5
        units = "Arc Minutes"
    else
        units = "Arc Seconds"
    end

    if incr >= 120
        incr = T(4)
    elseif incr >= 60
        incr = T(2)
    elseif incr >= 30
        incr = T(1)
    elseif incr > 15
        incr = T(30)
    elseif incr >= 10
        incr = T(15)
    elseif incr >= 5
        incr = T(10)
    elseif incr >= 2
        incr = T(5)
    elseif incr >= 1
        incr = T(2)
    elseif incr > 0.5
        incr = T(1)
    elseif incr >= 0.25
        incr = T(30)
    elseif incr >= 0.16
        incr = T(15)
    elseif incr >= 0.08
        incr = T(10)
    elseif incr >= 0.04
        incr = T(5)
    elseif incr >= 0.02
        incr = T(2)
    else
        incr = T(1)
    end

    if units == "Degrees"
        minpix = minpix/60
    elseif units == "Arc Seconds"
        minpix = minpix*60
    end

    ticsize = incr/minpix
    return ticsize, incr, units
end

"""
    ticpos(deglen, pixlen, ticsize) -> ticsize, incr, units

### Purpose ###

Specify distance between tic marks for astronomical coordinate overlays.

### Explanation ###

User inputs number an approximate distance between tic marks,
and the axis length in degrees. `ticpos` will return a distance
between tic marks such that the separation is a round multiple
in arc seconds, arc minutes, or degrees.

### Arguments ###

* `deglen`: length of axis in degrees, positive scalar
* `pixlen`: length of axis in plotting units (pixels), postive scalar
* `ticsize`: distance between tic marks (pixels).  This value will be
   adjusted by `ticpos` such that the distance corresponds to a round
   multiple in the astronomical coordinate.

### Output ###

The 3-tuple `(ticsize, incr, units)`:

* `ticsize`: distance between tic marks (pixels), positive scalar
* `incr`: incremental value for tic marks in round units given
   by the `units` parameter
* `units`: string giving units of ticsize, either 'Arc Seconds',
  'Arc Minutes', or 'Degrees'

### Example ###

Suppose a 512 x 512 image array corresponds to 0.2 x 0.2 degrees on the sky.
A tic mark is desired in round angular units, approximately every 75 pixels.
Then

```jldoctest
julia> using AstroLib

julia> ticpos(0.2, 512, 75)
(85.33333333333333, 2.0, "Arc Minutes")
```

i.e. a good tic mark spacing is every 2 arc minutes, corresponding
to 85.333 pixels.

### Notes ###

All the arguments taken as input are assumed to be positive in nature.

Code of this function is based on IDL Astronomy User's Library.
"""
ticpos(deglen::Real, pixlen::Real, ticsize::Real) =
    ticpos(promote(float(deglen), float(pixlen), float(ticsize))...)


text:
module GenericRefFEsTests
##

using Gridap
using Gridap.Helpers

# using Gridap.RefFEs.RaviartThomasRefFEs

using LinearAlgebra

using Test

using Gridap.GenericRefFEs

#####
# Create a general Ref FE constructor, which generalizes the current implementations
# for nodal (Lagrangian) and non-nodal (RaviartThomas) RefFEs.
# * A method that generates the geomap from ref n-face to n-face in polytope
# * A functor to be provided by the user (optionally) for every n-face dimension
#####
# F = Gridap.RefFEs.GenericRefFEs
#
# p = Polytope(1,1,1)
#
# order = 4
#
# ref1 = RTRefFE(p,Float64,3)
# ref2 = F.MyNedelecRefFE(p,3)
# ref1.shfbasis.changeofbasis == ref2.shfbasis.changeofbasis

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

function test_reffe(reffe,order)

  dofsb = dofbasis(reffe)

  p = polytope(reffe)

  shb = shfbasis(reffe)

  nfdofs = nfacedofs(reffe)

  b = dofsb
  shfs = reffe.shfbasis
  kk = evaluate(b,shfs)
  @test kk ≈ Matrix(1.0I, size(kk))

  fun(x) = VectorValue(x[1],2*x[1])
  D = 2
  T = VectorValue{2,Float64}
  field = AnalyticalField(fun,D)

  isa(dofsb,DOFBasis{2,T})
  isa(field,Field{2,T})

  evaluate(b,field)
end
##
p = Polytope(1,1)

for order in 1:4
  reffe = NedelecRefFE(p,Float64,order)
  test_reffe(reffe,order)
end

for order in 1:4
  reffe = RTRefFE(p,Float64,order)
  test_reffe(reffe,order)
end

##
end # module


text:
include("../utils.jl")

mutable struct Heap{A, H}
    array::A
    size::H
end

function Heap(A)
    return Heap(A, length(A))
end

size(heap) = heap.size

parent(i) = floor(Int, i/2)
left(i) = 2i
right(i) = 2i + 1

function max_heapify!(heap, i)
    l = left(i)
    r = right(i)

    largest = i
    A = heap.array

    if l ≤ size(heap) && A[l] > A[largest]
        largest = l
    end

    if r ≤ size(heap) && A[r] > A[largest]
        largest = r
    end

    if largest != i
        exchange!(A, i, largest)
        max_heapify!(heap, largest)
    end
end

function build_max_heap!(A)
    heap = Heap(A)

    start = floor(Int, size(heap) / 2)
    for i in range(start, 1, step = -1)
        max_heapify!(heap, i)
    end

    return heap
end

function heap_sort!(A)
    heap = build_max_heap!(A)
    for i in range(length(A), 2, step = -1)
        exchange!(heap.array, 1, i)
        heap.size -= 1
        max_heapify!(heap, 1)
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    @test my_f(2,1)==5
end


text:
@testset "Particle" begin

    a = Acoustic(1.0, 1.0, 2)
    a2 = Acoustic(1.0, 2.0, 2)

    @testset "Types" begin

        circle1 = Sphere((0.0,0.0),1.0)
        circle2 = Circle((0.0,5.0),2.0)
        a = Acoustic(1.0,1.0,2)
        homog_particles = [Particle(a,circle1), Particle(a,circle2)]

        # Check types comparisons work as user would expect
        @test typeof(homog_particles) <: AbstractParticles
        @test typeof(homog_particles) <: AbstractParticles{Float64}
        @test typeof(homog_particles) <: AbstractParticles{Float64,2}

        circle = Sphere((0.0,0.0),1.0)
        rect = Box((2.0,2.0),(3.0,2.0))
        diff_shape_particles = [Particle(a,circle), Particle(a,rect)]

        # Check types comparisons work as user would expect
        @test typeof(diff_shape_particles) <: AbstractParticles
        @test typeof(diff_shape_particles) <: AbstractParticles{Float64}
        @test typeof(diff_shape_particles) <: AbstractParticles{Float64,2}

        a2 = Acoustic(1.0,1.0,2)
        a3 = Acoustic(1.0,1.0,3)
        sphere = Sphere((0.0,0.0,0.0),1.0)

        circular_particle = Particle(a2,circle)
        @test typeof(circular_particle) <: Particle{Float64,2}

        spherical_particle = Particle(a3,sphere)
        @test typeof(spherical_particle) <: AbstractParticle{Float64,3}

        # Dimension mismatch throws error
        @test_throws MethodError Particle(a3,circle)
        @test_throws MethodError Particle(a2,sphere)

        # This is a valid vector of valid particles, but not of type Particles
        # because the dimensions don't match
        invalid_particles = [circular_particle, spherical_particle]
        @test_throws TypeError invalid_particles::(Vector{Pt} where {Dim, Pt <: AbstractParticle{Float64,Dim}})

        # does not throw an error
        diff_shape_particles::(Vector{Pt} where {Dim, Pt <: AbstractParticle{Float64,Dim}})
        @test true

    end

    @testset "Comparisons" begin

        circle = Sphere((1.0,3.0),2.0)
        circle_identical = Sphere((1.0,3.0),2.0)
        circle_congruent = Sphere((4.0,7.0),2.0)
        rect = Box((1.0,2.0),(2.0,4.0))

        # Construct three particles, with two the same
        p = Particle(a2,circle)
        p_reference = p
        p_identical = Particle(a2,circle_identical)
        p_different = Particle(a2,rect)
        p_congruent = Particle(a2,circle_congruent)

        # Test comparison operators
        @test p == p_identical
        @test p != p_different
        @test !(p == p_different)
        @test iscongruent(p, p_congruent)
        @test !iscongruent(p, p_different)

        # Check that Julia fallback === works
        @test p === p_reference
        # Particles are immutable, so are compared at the bit levelS
        @test p === p_identical

        # Construct two almost identical particles
        p1 = Particle(a, Sphere((0.0,0.0),1.0))
        p1b = Particle(a, Sphere((-0.0,0.0),1.0))

        # isequal is strict about +/-0.0 but IEEE compliant == is not strict
        @test p1 == p1b
        @test !isequal(p1,p1b)

    end

    @testset "Plot" begin
        # Just run it to see if we have any errors (yes thats a very low bar)

        # Single particle
        particle = Particle(a, Sphere((1.0,0.0), 2.0))
        plot(particle)

        # Vector of particles of same shape
        particles = [
            Particle(a, Sphere((1.0,0.0), 2.0)),
            Particle(a, Sphere((3.0,-1.0), 1.0))
        ]
        plot(particles)

        # Vector of particles of different shape
        particles = [
            CapsuleParticle(
                Particle(Acoustic(2.0,2.0,2),Sphere((5.0,3.0),1.0)),
                Particle(Acoustic(1.0,1.,2),Sphere((5.0,3.0),2.0))
            ),
            Particle(a, Box((1.0,0.0),(2.0, 3.0))),
            Particle(a, Sphere((3.0,-1.0), 1.0))
        ]
        plot(particles)

        @test true
    end

end

@testset "Capsule Particle" begin

    circle_in = Sphere((0.0,0.0),1.0)
    circle_out = Sphere((0.0,0.0),2.0)
    a_out = Acoustic(3.0,3.0,2)
    a = Acoustic(1.0,1.0,2)
    concen_particles = [ Particle(a_out,circle_out),Particle(a,circle_in)]
    @test typeof(CapsuleParticle(concen_particles...)) <: AbstractParticle{Float64,2}

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    RSOCtoSOCBridge{T,F,G} <: Bridges.Constraint.AbstractBridge

`RSOCtoSOCBridge` implements the following reformulation:

  * ``||x||_2^2 \\le 2tu`` into ``||v|| \\le w ``, with the substitution rules
    ``t = \\frac{w}{\\sqrt 2} + \\frac{v_1}{\\sqrt 2}``,
    ``u = \\frac{w}{\\sqrt 2} - \\frac{v_1}{\\sqrt 2}``, and
    ``x = (v_2,\\ldots,v_N)``.

## Source node

`RSOCtoSOCBridge` supports:

  * `G` in [`MOI.RotatedSecondOrderCone`](@ref)

## Target node

`RSOCtoSOCBridge` creates:

  * `F` in [`MOI.SecondOrderCone`](@ref)
"""
struct RSOCtoSOCBridge{T,F,G} <:
       SetMapBridge{T,MOI.SecondOrderCone,MOI.RotatedSecondOrderCone,F,G}
    constraint::MOI.ConstraintIndex{F,MOI.SecondOrderCone}
end

const RSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{RSOCtoSOCBridge{T},OT}

function concrete_bridge_type(
    ::Type{<:RSOCtoSOCBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RotatedSecondOrderCone},
) where {T}
    S = MOI.Utilities.promote_operation(/, T, MOI.Utilities.scalar_type(G), T)
    Y = MOI.Utilities.promote_operation(-, T, S, S)
    Z = MOI.Utilities.promote_operation(+, T, S, S)
    F = MOI.Utilities.promote_operation(vcat, T, Z, Y, G)
    return RSOCtoSOCBridge{T,F,G}
end

function MOI.Bridges.map_set(
    ::Type{<:RSOCtoSOCBridge},
    set::MOI.RotatedSecondOrderCone,
)
    return MOI.SecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:RSOCtoSOCBridge},
    set::MOI.SecondOrderCone,
)
    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.map_function(::Type{<:RSOCtoSOCBridge{T}}, func) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    t, u, x = scalars[1], scalars[2], scalars[3:end]
    ts = MOI.Utilities.operate!(/, T, t, sqrt(T(2)))
    us = MOI.Utilities.operate!(/, T, u, sqrt(T(2)))
    return MOI.Utilities.operate(vcat, T, ts + us, ts - us, x)
end

# The map is an involution
function MOI.Bridges.inverse_map_function(BT::Type{<:RSOCtoSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is symmetric
function MOI.Bridges.adjoint_map_function(BT::Type{<:RSOCtoSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is a symmetric involution
function MOI.Bridges.inverse_adjoint_map_function(
    BT::Type{<:RSOCtoSOCBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end


text:
using Test
using Sounds

@testset "Sound" begin
    sound = Sound(joinpath("..", "examples", "piano.wav"))
    @show sound

    @test duration(sound)   == 32e6
    @test nchannels(sound)  == 2
    @test samplerate(sound) == 44100
    @test nsamples(sound)   == 2822400
    @test timepos(sound)    == 0
    @test getstatus(sound)  == :stopped

    @test islooping(sound) == false
    loop!(sound, true)
    @test islooping(sound) == true

    @test volume(sound) == 100
    volume!(sound, 50)
    @test volume(sound) == 50

    @test pitch(sound) == 1.0
    pitch!(sound, 2.0)
    @test pitch(sound) == 2
end

@testset "Music" begin
    music = Music(joinpath("..", "examples", "piano.wav"))
    @show music

    @test duration(music)   == 32e6
    @test nchannels(music)  == 2
    @test samplerate(music) == 44100
    @test timepos(music)    == 0
    @test getstatus(music)  == :stopped

    @test islooping(music) == false
    loop!(music, true)
    @test islooping(music) == true

    @test volume(music) == 100
    volume!(music, 50)
    @test volume(music) == 50

    @test pitch(music) == 1.0
    pitch!(music, 2.0)
    @test pitch(music) == 2
end


text:
using DEBStandardModel
using Test

@testset "AmPtests" begin
    include("scaledDEBtest.jl") end

max_stars_repo_name: HenriLaurie/DEBStandardModel

max_issues_repo_name: HenriLaurie/DEBStandardModel

max_forks_repo_name: HenriLaurie/DEBStandardModel


text:
module ECDSA

using BitConverter
using Secp256k1: Point, KeyPair, Signature, N, G
export KeyPair

KeyPair{:ECDSA}(𝑑::Integer) = 𝑑 ∉ 1:N-1 ? throw(NotInField()) : KeyPair{:ECDSA}(𝑑, 𝑑 * G)

"""
    ECDSA.sign(kp::KeyPair{:ECDSA}, 𝑧::Integer) -> Signature{:ECDSA}

Returns a Signature{:ECDSA} for a given `KeyPair{:ECDSA}` and data `𝑧` and in
which 𝑠 = (𝑧 + 𝑟𝑑) / 𝑘, 𝑘 being a random integer.
"""
function sign(kp::KeyPair{:ECDSA}, 𝑧::Integer)
    𝑘 = rand(big.(0:N))
    𝑟 = (𝑘 * G).𝑥.𝑛
    𝑘⁻¹ = powermod(𝑘, N - 2, N)
    𝑠 = mod((𝑧 + 𝑟 * kp.𝑑) * 𝑘⁻¹, N)
    if 𝑠 > N / 2
        𝑠 = N - 𝑠
    end
    return Signature{:ECDSA}(𝑟, 𝑠)
end

"""
    verify(𝑄::Point, 𝑧::Integer, sig::Signature{:ECDSA}) -> Bool

Returns true if Signature{:ECDSA} is valid for 𝑧 given 𝑄, false if not
"""
function verify(𝑄::Point, 𝑧::Integer, sig::Signature{:ECDSA})
    𝑠⁻¹ = powermod(sig.𝑠, N - 2, N)
    𝑢 = mod(𝑧 * 𝑠⁻¹, N)
    𝑣 = mod(sig.𝑟 * 𝑠⁻¹, N)
    𝑅 = 𝑢 * G + 𝑣 * 𝑄
    return 𝑅.𝑥.𝑛 == sig.𝑟
end

end  # module ECDSA


text:
struct BPFState{T} <: AbstractFilterState
    ensemble::WeightedParticleEnsemble{T}
    function BPFState(ens::WeightedParticleEnsemble{T}) where T
        new{T}(ens)
    end
end

struct BPF{TS<:HiddenStateModel, TO<:ObservationModel, T<:Real} <: AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime}
    state_model::TS
    obs_model::TO
    N::Int
    alpha::T
end

#####################
###### METHODS ######
#####################

initial_condition(filter::BPF) = filter.state_model.init
no_of_particles(filter::BPF) = filter.N
no_of_particles(st::BPFState) = no_of_particles(st.ensemble)
state_model(filter::BPF) = filter.state_model
obs_model(filter::BPF) = filter.obs_model

mean(st::BPFState) = mean(st.ensemble)
cov(st::BPFState)  = cov(st.ensemble)
var(st::BPFState)  = var(st.ensemble)

function Base.show(io::IO, ::MIME"text/plain", filter::BPF)
    print(io, "Bootstrap particle filter algorithm
    hidden state model:                              ", state_model(filter),"
    observation model:                               ", obs_model(filter),"
    ensemble size:                                   ", no_of_particles(filter),"
    threshold for effective number of particles (%): ", 100*filter.alpha)
end

function Base.show(io::IO, filter::BPF)
    print(io, "BPF with ", no_of_particles(filter)," particles ")
end

######################
### MAIN ALGORITHM ###
######################

function initialize(filter::BPF)
    return BPFState(filter)
end

function propagate!(filter_state::BPFState, filt_algo::BPF, dt)
    propagate!(filter_state.ensemble, state_model(filt_algo), dt)
end

function assimilate!(observation, filter_state::BPFState, filt_algo::BPF, dt)
    ensemble  = filter_state.ensemble

    update_weights!(ensemble, obs_model(filt_algo), observation, dt)

    if eff_no_of_particles(ensemble) < filt_algo.alpha * no_of_particles(ensemble)
        resample!(ensemble)
    end

    return ensemble.positions
end

function update!(filter_state::BPFState, filt_algo::BPF, dY, dt)
    propagate!(filter_state, filt_algo, dt)
    assimilate!(dY, filter_state, filt_algo, dt)
end

########################
### HELPER FUNCTIONS ###
########################

function update_weights!(ensemble::WeightedParticleEnsemble, model::DiffusionObservationModel, dY, dt)
    H         = mapslices(observation_function(model), ensemble.positions, dims=1)
    mean_H    = StatsBase.mean(H, ensemble.weights, dims=2)

    @inbounds for i in 1:no_of_particles(ensemble), j in 1:obs_dim(model)
        ensemble.weights[i] += ensemble.weights[i] * ( H[j,i] - mean_H[j,1] ) * (dY[j] - mean_H[j,1] * dt)
    end
end

function update_weights!(ensemble::WeightedParticleEnsemble, model::CountingObservationModel, dN, dt)
    H         = mapslices(observation_function(model), ensemble.positions, dims=1)
    mean_H    = StatsBase.mean(H, ensemble.weights, dims=2)

    @inbounds for i in 1:no_of_particles(ensemble), j in 1:obs_dim(model)
        ensemble.weights[i] += ensemble.weights[i] * ( H[j,i] - mean_H[j,1] ) * (dN[j] - mean_H[j,1] * dt) / mean_H[j,1]
    end
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

function BPFState(filter::BPF)
    N   = no_of_particles(filter)
    ens = WeightedParticleEnsemble(hcat([initialize(filter.state_model) for i in 1:N]...), StatsBase.ProbabilityWeights(fill(1/N, N)))
    return BPFState(ens)
end

function BPFState(filt_prob::AbstractFilteringProblem, N)
    return BPFState(state_model(filt_prob), obs_model(filt_prob), N)
end

function BPFState(state_model::HiddenStateModel, obs_model::ObservationModel, N)
    ens = WeightedParticleEnsemble(hcat([initialize(state_model) for i in 1:N]...), StatsBase.ProbabilityWeights(fill(1/N, N)))
    return BPFState(ens)
end

function BPF(filt_prob::AbstractFilteringProblem, N, alpha)
    st_mod = state_model(filt_prob)
    ob_mod = obs_model(filt_prob)

    return BPF(st_mod, ob_mod, N, alpha)
end


text:
include("../btg.jl")
include("btgBayesopt.jl")
#include("../computation/finitedifference.jl")
using Plots
using Random
Random.seed!(8);

#####
##### Define optimization problem, initialize GP with burn-in points, initialize btg parameters
#####
Himmelblau(x) = (x[1]^2 + x[2] -11)^2 + (x[1]+x[2]^2-7)^2 #function to optimize
y, x = sample_points(Himmelblau, [-5, -5], [5, 5]; num = 10) #burn-in points which are used as training data in GP
x = vcat(x, [5 5])
y = vcat(y, 890)
x = vcat(x, [-4.999 4.999])
y = vcat(y, 530)
x = vcat(x, [4.999 -4.999])
y = vcat(y, 610)
x = vcat(x, [-4.999 -4.999])
y = vcat(y, 249.6042059640019)
x = vcat(x, [5 2.55573])
y = vcat(y, 294.57103540062496)
x = vcat(x, [5.0 2.3475172404355718])
y = vcat(y, 279.5672977301108)

cur_min = minimum(mapslices(Himmelblau, x, dims = 2))

@info "y-vals", y[1:5]
### Set BTG parameters and get function handles for pdf, cdf, dpdf, etc.
Fx = linear_polynomial_basis(x)
train = extensible_trainingData(x, Fx, y) #training data triple: location, covariates, labels
#rangeθ = reshape(select_single_theta_range(x), 1, 2) #rangetheta is auto-selected
#rangeθ = [0.0001 1000]
#rangeθ = [10.0 300]
rangeθ  = [50.0 1000.0]
@info "rangeθ", rangeθ
rangeλ = [-1.0 1]
btg1 = btg(train, rangeθ, rangeλ);
 #get function handles for pdf, cdf, derivtives
lx = [1, -5, -5]; ux = [3000.0, 5, 5]
#(_, _, plt1, pol) = checkDerivative(cdf_fixed, cdf_gradient_fixed, vstars[1],  nothing, 4, 8)

"""
Takes n BO steps using a pre-initialized btg1 object, which already contains some points in its kernel system
Returns log of points chosen from optimizing UCB optimization problem

INPUTS: f, expensive function we are optimizing with BO
 - min: minimum value of f on points already inside BO kernel system (we assume it's initialized with burn-in points)

NOTE: assumes linear polynomial basis, e.g. covariates for point [x, y] is [1, x, y]

OUTPUTS:
    -  x_hist: history of places what the GP chooses to sample (via optimizing the acquisition function)
    -  f_hist: history of actual function evaluations at x_hist
    -  gp_hist: history of values the gp think is at 25th quantile
    -  min_hist: keeps track of minimum over all previous f_evals, including burn-in
    -  init_hist: history of initializations for acqusition function optimization routines
    -  quant_hist: history of cdf_evals at chosen sample points (should be 0.25)
"""
function BO(btg1::btg, f, n, lx, ux, min) #BO loop
    x_hist = [] #keep track of selected points
    f_hist = []
    gp_hist = [] #compare what the gp thinks the objective looks like to what it actually is (f_hist)
    min_hist = [min]
    init_hist = []
    quant_hist = []
    #debug histories
    training_size_hist = []
    argmin_hist = [x[argmin(getLabel(btg1.trainingData)):argmin(getLabel(btg1.trainingData)), :]]

    for i = 1:n #take 10 BO steps
        println("========================ITERATION: ", i)
        (_, cur_cdf, _, cdf_gradient) = solve(btg1; derivatives = true);
        #(vstar, initval, quant_eval) = optimize_acqusition(cdf, cdf_gradient, lx, ux; maxiter = 300, initial = argmin_hist[end]);
        (vstar, initval, quant_eval) = optimize_acqusition(cur_cdf, cdf_gradient, lx, ux; maxiter = 300, initial = [], quant = 0.5);
        s_star = reshape(vstar[2:end], 1, length(vstar[2:end]))
        push!(quant_hist, quant_eval)
        push!(init_hist, initval[2:end])
        push!(gp_hist, vstar[1])
        push!(x_hist, s_star) #save
        push!(training_size_hist, btg1.trainingData.n)
        f_star = f(s_star)
        push!(f_hist, f_star)
        if f_star < min_hist[end]
            push!(min_hist, f_star)
            push!(argmin_hist, s_star)
        else
            push!(min_hist, min_hist[end])
            push!(argmin_hist, argmin_hist[end])
        end
        #@info "f_star", f_star
        updateBTG!(btg1, s_star, linear_polynomial_basis(s_star), [f_star]) #update with chosen point
    end
    return x_hist, f_hist, min_hist, gp_hist, init_hist, quant_hist, training_size_hist
end
(x_hist, f_hist, min_hist, gp_hist, init_hist, quant_hist, training_size_hist) = BO(btg1, Himmelblau, 5, lx, ux, cur_min);

include("analyze.jl")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
if scene_mode == "data"
    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = CrowdNavControlParameter(model_dir, env_config, policy_config,
                                         policy_name, ego_pos_goal_vec,
                                         tcalc, dtr, target_speed);
else
    @error @error "scene_mode: $(scene_mode) is not supported!"
end


text:
function triangle(n)
    n < 0 && throw(DomainError("Invalid number"))
    n == 0 && return []
    rows = [[1]]
    for i in 2:n
        previous_row = [rows[i-1]..., 0]
        push!(rows, previous_row + reverse(previous_row))
    end
    rows
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using BinDeps

@BinDeps.setup

libfdt = library_dependency("libfdt")

builddir = joinpath(BinDeps.depsdir(libfdt), "build")
src = joinpath(BinDeps.depsdir(libfdt), "src")
prefix = joinpath(BinDeps.depsdir(libfdt), "usr")
bindir = joinpath(prefix, "lib")
binpath = joinpath(bindir, "libfdt."*BinDeps.shlib_ext)

if isdir(builddir)
    rm(builddir, recursive = true)
end

if isdir(prefix)
    rm(prefix, recursive = true)
end

for path in [prefix, builddir, src, bindir]
    !isdir(path) && mkdir(path)
end

if !isfile(binpath)
    @unix_only begin
        cd(src)
        for mkcmd in (:gnumake, :gmake, :make)
            try
                if success(`$mkcmd`)
                    cp(joinpath(src, "libfdt.so"), binpath)
                    break
                end
            catch
                continue
            end
        end
    end

    @windows_only begin
        Base.warn("No build available on Windows yet.")
    end
end

provides(Binaries, bindir, libfdt)
@BinDeps.install [ :libfdt => :libfdt]

max_stars_repo_name: tanmaykm/Chrono

max_issues_repo_name: tanmaykm/Chrono

max_forks_repo_name: tanmaykm/Chrono


text:
## Aircraft analysis case
using AeroMDAO

## Wing
wing = Wing(foils     = Foil.(fill(naca4((0,0,1,2)), 3)),
            chords    = [1.0, 0.6, 0.2],
            twists    = [2.0, 0.0, -2.0],
            spans     = [4.0, 0.2],
            dihedrals = [5., 30.],
            sweep_LEs = [5., 30.]);
print_info(wing, "Wing")

# Horizontal tail
htail = Wing(foils     = Foil.(fill(naca4((0,0,1,2)), 2)),
             chords    = [0.7, 0.42],
             twists    = [0.0, 0.0],
             spans     = [1.25],
             dihedrals = [0.],
             sweep_LEs = [6.39])
print_info(htail, "Horizontal Tail")

# Vertical tail
vtail = HalfWing(foils     = Foil.(fill(naca4((0,0,0,9)), 2)),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.0],
                 dihedrals = [0.],
                 sweep_LEs = [7.97])
print_info(vtail, "Vertical Tail")

## Assembly
wing_panels  = panel_wing(wing,                 # Wing or HalfWing type
                          [20, 3],              # Number of spanwise panels for half-wing
                          10;                   # Chordwise panels
                        #   spacing = "cosine"  # Spacing distribution: Default works well for symmetric
                         )

htail_panels = panel_wing(htail, [6], 6;
                          position	= [4., 0, 0],
                          angle 	= deg2rad(-2.),
                          axis 	  	= [0., 1., 0.],
                          spacing   = "uniform"
                         )
vtail_panels = panel_wing(vtail, [6], 5;
                          position 	= [4., 0, 0],
                          angle 	= π/2,
                          axis 	 	= [1., 0., 0.],
                          spacing   = "uniform"
                         )

aircraft = Dict("Wing" 			  	=> wing_panels,
                "Horizontal Tail" 	=> htail_panels,
                "Vertical Tail"   	=> vtail_panels)

S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing);

## Case
ac_name = "My Aircraft"
ρ       = 1.225
ref     = [0.25c, 0., 0.]
V, α, β = 1.0, 1.0, 0.0
Ω       = [0.0, 0.0, 0.0]
fs      = Freestream(V, α, β, Ω)

@time data =
    solve_case(aircraft, fs;
               rho_ref     = ρ, 		# Reference density
               r_ref       = ref, 		# Reference point for moments
               area_ref    = S, 		# Reference area
               span_ref    = b, 		# Reference span
               chord_ref   = c, 		# Reference chord
               name        = ac_name,	# Aircraft name
               print       = true,		# Prints the results for the entire aircraft
            #    print_components = true,	# Prints the results for each component
              );

## Data collection
comp_names = (collect ∘ keys)(data) # Gets aircraft component names from analysis
comp  = comp_names[1]			   # Pick your component
nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes, Γs = data[comp]; #  Get the nearfield, farfield, force and moment coefficients, and other data for post-processing
print_coefficients(nf_coeffs, ff_coeffs, comp)

## Stability case
@time dv_data =
    solve_stability_case(aircraft, fs;
                         rho_ref     = ρ,
                         r_ref       = ref,
                         area_ref    = S,
                         span_ref    = b,
                         chord_ref   = c,
                         name        = ac_name,
                         print       = true,
                         print_components = true,
                        );

## Data collection
names = (collect ∘ keys)(dv_data)
comp  = names[1]
nf, ff, dvs = dv_data[comp];
print_coefficients(nf, ff, comp)
print_derivatives(dvs, comp)

## Streamlines

# Chordwise distribution
# using StaticArrays

# num_points = 100
# max_z = 2
# y = span(wing) / 10
# seed = SVector.(fill(-0.1, num_points), fill(y, num_points), range(-max_z, stop = max_z, length = num_points))

# Spanwise distribution
span_points = 10
init        = chop_leading_edge(wing, span_points)
dx, dy, dz  = 0, 0, 1e-3
seed        = [ init .+ Ref([dx, dy, dz]) ;
                init .+ Ref([dx, dy,-dz]) ];

distance = 8
num_stream_points = 200
streams = plot_streams(fs, seed, horseshoes, Γs, distance, num_stream_points);

##
using Plots
gr(size = (200, 100), dpi = 300)

##
horseshoe_coords = plot_panels(horseshoe_panels[:]);

z_limit = b
plot(xaxis = "x", yaxis = "y", zaxis = "z",
     camera = (30, 60),
    #  xlim = (-z_limit/2, z_limit/2),
     aspect_ratio = 1,
     zlim = (-z_limit/2, z_limit/2),
     size = (1280, 720)
    )
plot!.(horseshoe_coords, color = :gray, label = :none)
plot!.(streams, color = :green, label = :none)
plot!()

## Exaggerated CF distribution for plot, only works with GR and not Plotly

# Forces
wind_CFs = body_to_wind_axes.(CFs, fs.alpha, fs.beta)
CDis     = @. getindex(wind_CFs, 1)
CYs	     = @. getindex(wind_CFs, 2)
CLs      = @. getindex(wind_CFs, 3)

hs_pts = bound_leg_center.(horseshoes)
hs_xs  = getindex.(hs_pts, 1)
hs_ys  = getindex.(hs_pts, 2)
hs_zs  = getindex.(hs_pts, 3)

quiver!(hs_xs[:], hs_ys[:], hs_zs[:], quiver=(CDis[:], CYs[:], CLs[:]) .* 500)
plot!()

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# ---
# title: 71. Simplify Path
# id: problem71
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/simplify-path/description/>
# hidden: true
# ---
#
# Given an **absolute path** for a file (Unix-style), simplify it. Or in other
# words, convert it to the **canonical path**.
#
# In a UNIX-style file system, a period `'.'` refers to the current directory.
# Furthermore, a double period `'..'` moves the directory up a level.
#
# Note that the returned canonical path must always begin with a slash `'/'`,
# and there must be only a single slash `'/'` between two directory names. The
# last directory name (if it exists) **must not**  end with a trailing `'/'`.
# Also, the canonical path must be the **shortest** string representing the
# absolute path.
#
#
#
# **Example 1:**
#
#
#
#     Input: path = "/home/"
#     Output: "/home"
#     Explanation: Note that there is no trailing slash after the last directory name.
#
#
# **Example 2:**
#
#
#
#     Input: path = "/../"
#     Output: "/"
#     Explanation: Going one level up from the root directory is a no-op, as the root level is the highest level you can go.
#
#
# **Example 3:**
#
#
#
#     Input: path = "/home//foo/"
#     Output: "/home/foo"
#     Explanation: In the canonical path, multiple consecutive slashes are replaced by a single one.
#
#
# **Example 4:**
#
#
#
#     Input: path = "/a/./b/../../c/"
#     Output: "/c"
#
#
#
#
# **Constraints:**
#
#   * `1 <= path.length <= 3000`
#   * `path` consists of English letters, digits, period `'.'`, slash `'/'` or `'_'`.
#   * `path` is a valid Unix path.
#
#
## @lc code=start
using LeetCode

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


text:
"""
    smc!(model::SMCModel, smcio::SMCIO)
Run the SMC algorithm for the given model and input/output arguments.

If ```smcio.nthreads == 1``` the algorithm will run in serial.
"""
function smc!(model::SMCModel, smcio::SMCIO)
  if smcio.nthreads > 1
    _smcParallel!(model, smcio)
  else
    _smcSerial!(model, smcio)
  end
end

"""
    csmc!(model::SMCModel, smcio::SMCIO, ref::Vector{Particle}, refout::Vector{Particle})

Run the conditional SMC algorithm for the given model, input/output arguments,
reference path and output path.

It is permitted for ref and refout to be the same. If ```smcio.nthreads == 1```
the algorithm will run in serial.
"""
function csmc!(model::SMCModel, smcio::SMCIO{Particle}, ref::Vector{Particle},
  refout::Vector{Particle}) where Particle
  if smcio.nthreads > 1
    _smcParallel!(model, smcio, ref, refout)
  else
    _smcSerial!(model, smcio, ref, refout)
  end
end

@inline function _getZetas(smcio::SMCIO{Particle}, hat::Bool, p::Int64) where
  Particle
  @assert 1 <= p <= smcio.n
  p == smcio.n && return smcio.zetas
  @assert smcio.fullOutput "full SMC output required for this operation"
  return smcio.allZetas[p]
end

@inline function _getWs(smcio::SMCIO, hat::Bool, p::Int64)
  @assert 1 <= p <= smcio.n
  if hat
    p == smcio.n && return smcio.ws
    @assert smcio.fullOutput "full SMC output required for this operation"
    return smcio.allWs[p]
  else
    if p == 1 || smcio.resample[p-1] == true
      return smcio.internal.vecOnes
    else
      @assert smcio.fullOutput "full SMC output required for this operation"
      return smcio.allWs[p-1]
    end
  end
end

"""
    eta(smcio::SMCIO{Particle}, f::F, hat::Bool, p::Int64) where {Particle, F<:Function}

Compute:
  - ```!hat```: ``\\eta^N_p(f)``
  - ```hat```:  ``\\hat{\\eta}_p^N(f)``
"""
function eta(smcio::SMCIO{Particle}, f::F, hat::Bool, p::Int64) where
  {Particle, F<:Function}
  zetas::Vector{Particle} = _getZetas(smcio, hat, p)
  ws::Vector{Float64} = _getWs(smcio, hat, p)
  v = f(zetas[1]) * ws[1]
  for i = 2:smcio.N
    @inbounds v += f(zetas[i]) * ws[i]
  end
  return v / smcio.N
end

"""
    allEtas(smcio::SMCIO, f::F, hat::Bool) where F<:Function

Compute ```eta(smcio::SMCIO, f::F, hat::Bool, p)``` for p in {1, …, smcio.n}
"""
function allEtas(smcio::SMCIO, f::F, hat::Bool) where F<:Function
  T = typeof(f(smcio.zetas[1]) / smcio.N)
  result::Vector{T} = Vector{T}(undef, smcio.n)
  for p = 1:smcio.n
    @inbounds result[p] = eta(smcio, f, hat, p)
  end
  return result
end

"""
    slgamma(smcio::SMCIO, f::F, hat::Bool, p::Int64) where {Particle, F<:Function}

Compute:
- ```!hat```: ``(\\eta^N_p(f) \\geq 0, \\log |\\gamma^N_p(f)|)``
- ```hat```:  ``(\\hat{\\eta}^N_p(f) \\geq 0, \\log |\\hat{\\gamma}_p^N(f)|)``
The result is returned as a ```Tuple{Bool, Float64}```: the first component
represents whether the returned value is non-negative, the second is the
logarithm of the absolute value of the approximation.
"""
function slgamma(smcio::SMCIO, f::F, hat::Bool, p::Int64) where
  F<:Function
  @assert 1 <= p <= smcio.n "p must be between 1 and smcio.n"
  idx::Int64 = p - 1 + hat
  logval = idx == 0 ? 0.0 : smcio.logZhats[idx]
  v::Float64 = eta(smcio, f, hat, p)
  logval += log(abs(v))
  return (v >= 0, logval)
end

"""
    allGammas(smcio::SMCIO, f::F, hat::Bool) where F<:Function

Compute ```slgamma(smcio::SMCIO, f::F, hat::Bool, p)``` for p in {1, …, smcio.n}
"""
function allGammas(smcio::SMCIO, f::F, hat::Bool) where F<:Function
  result::Vector{Tuple{Bool, Float64}} =
    Vector{Tuple{Bool, Float64}}(undef, smcio.n)
  for p = 1:smcio.n
    @inbounds result[p] = slgamma(smcio, f, hat, p)
  end
  return result
end

function _centreFunction(smcio::SMCIO, f::F, hat::Bool,
  p::Int64) where F<:Function
  mu = eta(smcio, f, hat, p)
  function cf(x)
    return f(x) - mu
  end
  return cf
end

"""
    V(smcio::SMCIO{Particle}, f::F, hat::Bool, centred::Bool, p::Int64) where
      {Particle, F<:Function}
Compute:
- ```!hat & !centred```: ``V^N_p(f)``
- ```!hat & centred```: ``V^N_p(f-\\eta_p^N(f))``
- ```hat & !centred```: ``\\hat{V}_p^N(f)``
- ```hat & centred```:  ``\\hat{V}_p^N(f-\\hat{\\eta}_p^N(f))``
"""
function V(smcio::SMCIO{Particle}, f::F, hat::Bool, centred::Bool,
  p::Int64) where {Particle, F<:Function}

  centred && return V(smcio, _centreFunction(smcio, f, hat, p), hat, false, p)

  zetas::Vector{Particle} = _getZetas(smcio, hat, p)
  ws::Vector{Float64} = _getWs(smcio, hat, p)
  if p == smcio.n
    eves = smcio.eves
  else
    eves = smcio.allEves[p]
  end

  etafSq = eta(smcio, f, hat, p)^2
  N::Int64 = smcio.N
  i::Int64 = 1
  v::Float64 = 0.0
  while i <= N
    @inbounds e = eves[i]
    tmp::Float64 = 0.0
    @inbounds while i <= N && eves[i] == e
      @inbounds tmp += ws[i] * f(zetas[i])
      i += 1
    end
    tmp /= N
    v += tmp*tmp
  end
  q::Int64 = 1
  for i = 1:p-1
    q += smcio.resample[i]
  end
  return etafSq - (etafSq-v)*(N/(N-1))^q
end

## currently single-threaded only
"""
    vpns(smcio::SMCIO, f::F, hat::Bool, centred::Bool, n::Int64) where F<:Function

Compute a vector of the values of, for p = 1, …, n,
- ```!hat & !centred```: ``v^N_{p,n}(f)``
- ```!hat & centred```:  ``v^N_{p,n}(f-\\eta_n^N(f))``
- ```hat & !centred```:  ``\\hat{v}_{p,n}^N(f)``
- ```hat & centred```:   ``\\hat{v}_{p,n}^N(f-\\hat{\\eta}_n^N(f))``

Note: if ```essThreshold <= 1.0```, and resampling did not occur at every time,
the length of the output will be less than n.
"""
function vpns(smcio::SMCIO, f::F, hat::Bool, centred::Bool,
  n::Int64) where F<:Function

  centred && return vpns(smcio, _centreFunction(smcio, f, hat, n), hat, false, n)

  N = smcio.N
  as = smcio.allAs
  eves = smcio.allEves
  ws = smcio.allWs
  zetas = smcio.allZetas[n]

  etafSq = eta(smcio, f, hat, n)^2

  Slocal = smcio.internal.scratch1
  if hat
    for i = 1:N
      @inbounds Slocal[i] = f(zetas[i]) * ws[n][i] / N
    end
  else
    if n > 1 && !smcio.resample[n-1]
      for i = 1:N
        @inbounds Slocal[i] = f(zetas[i]) * ws[n-1][i] / N
      end
    else
      for i = 1:N
        @inbounds Slocal[i] = f(zetas[i]) / N
      end
    end
  end
  Stmp = smcio.internal.scratch2
  Glocal = smcio.internal.scratch2
  fill!(Glocal, 0.0)
  if n > 1
    for i = 1:N
      @inbounds Glocal[eves[n-1][i]] += ws[n-1][i] / N
    end
  end

  result = Vector{Float64}(undef, n)
  factor0::Float64 = (N/(N-1))^n
  factor1::Float64 = (N-1)*factor0
  vp::Float64 = 0.0
  for i = 1:N
    @inbounds vp += Slocal[i]*Slocal[i] * (1.0 - Glocal[eves[n][i]])
  end
  result[n] = factor1 * vp

  for p = n-1:-1:1
    # Slocal = S[p+1,:]
    # Glocal = G[p-1,:]
    fill!(Glocal, 0.0)
    if p != 1
      for i = 1:N
        @inbounds Glocal[eves[p-1][i]] += ws[p-1][i] / N
      end
    end
    i::Int64 = 1
    vp = 0.0
    while i <= N
      @inbounds a::Int64 = as[p][i]
      R1::Float64 = 0.0
      R2::Float64 = 0.0
      @inbounds while i <= N && as[p][i] == a
        @inbounds R1 += Slocal[i]
        @inbounds R2 += Slocal[i]*Slocal[i]
        i += 1
      end
      @inbounds vp += (R1*R1 - R2) * (1.0 - Glocal[eves[p][a]])
    end
    @inbounds result[p] = factor1 * vp
    copyto!(Stmp, Slocal)
    fill!(Slocal, 0.0)
    for i = 1:N
      @inbounds Slocal[as[p][i]] += Stmp[i]
    end
  end

  vp = 0.0
  for i = 1:N
    @inbounds vp += Slocal[i]*Slocal[i]
  end
  mstar::Float64 = factor0 * (etafSq - vp)
  result .-= mstar

  actualResult = Vector{Float64}(undef, sum(smcio.resample[1:n-1]) + 1)
  j::Int64 = 0
  for i = 1:length(actualResult)-1
    j += 1
    while !smcio.resample[j]
      j += 1
    end
    actualResult[i] = result[j]
  end
  actualResult[length(actualResult)] = result[n]
  return actualResult
end

"""
    v(smcio::SMCIO, f::F, hat::Bool, centred::Bool, n::Int64) where F<:Function

Compute:
- ```!hat & !centred```: ``v^N_n(f)``
- ```!hat & centred```: ``v^N_n(f-\\eta_n^N(f))``
- ```hat & !centred```: ``\\hat{v}_n^N(f)``
- ```hat & centred```:  ``\\hat{v}_n^N(f-\\hat{\\eta}_n^N(f))``
"""
function v(smcio::SMCIO, f::F, hat::Bool, centred::Bool,
  n::Int64) where F<:Function
  return sum(vpns(smcio, f, hat, centred, n))
end


text:
using ParserCombinator.Parsers.GML
using ZipFile

println("10k-49963")
s = open(readall, "gml/10k-49963.gml")
parse_raw(s; debug=false)
@time parse_raw(s; debug=false)

#println("10k-49963 raw")
#io = open("gml/10k-49963.gml")
#parse_dict(io; debug=false)
#io = open("gml/10k-49963.gml")
#@time parse_dict(io; debug=false)

# initial (without :!)
#  223.622 seconds      (284 M allocations: 66954 MB, 18.56% gc time)
# with :!
#  235.137 seconds      (274 M allocations: 66469 MB, 6.03% gc time)
# with push!
#   87.233 seconds      (273 M allocations: 11556 MB, 5.63% gc time)
#   99.867 seconds      (273 M allocations: 11556 MB, 5.25% gc time)
# with :! on space
#   83.584 seconds      (235 M allocations: 9603 MB, 5.39% gc time)
#   85.376 seconds      (235 M allocations: 9603 MB, 5.44% gc time)
# slurping spaces with [, ]
#   96.082 seconds      (274 M allocations: 11600 MB, 5.46% gc time)
# string input
#   74.171 seconds      (235 M allocations: 9617 MB, 5.67% gc time)
# parse_one instead of parse_try (with string)
#   OOM
# with fixed string
#   74.171 seconds      (238 M allocations: 9788 MB, 5.71% gc time)


text:
using Documenter, ADSeismic

makedocs(sitename="ADSeismic.jl",
pages = Any["index.md",
        "ad.md",
        "forward_simulation.md",
        "backward_inversion.md",
        "NNFWI.md",
        "contents.md",
        "api.md"])

deploydocs(
    repo = "github.com/kailaix/ADSeismic.jl.git",
)


text:
module Sobol

using StaticArrays
using Sobol: SobolSeq, next!
using LinearAlgebra: det, norm
using NBodyIPs: lengths_and_angles

"""
* N : body-order
* D : dictionary
* r0, r1 : upper and lower bound on the bond-lengths
"""
function filtered_sobol(x0, x1, filter; npoints=nothing,
                                        nfailed=nothing)
   @assert npoints isa Integer && nfailed isa Integer
   x0, x1 = min.(x0, x1), max.(x0, x1)
   return _sobol_inner(SVector(x0...), SVector(x1...), npoints, nfailed, filter)
end

function _sobol_inner(x0::SVector{DIM, T}, x1::SVector{DIM, T},
                      npoints, nfailed, filter ) where {DIM, T}
   # Sobol sequence in the [x0, x1]^dim hypercube
   s = SobolSeq(DIM, x0, x1)
   # temporary storage
   t = zero(MVector{DIM, T})
   # output storage
   X = SVector{DIM, T}[]
   # generate points
   npts = 0
   nfl = 0
   while (npts < npoints) && (nfl < nfailed)
      next!(s, t)
      if filter(t)
         push!(X, SVector(t))
         npts += 1
      else
         nfl += 1
      end
   end
   @show npts, nfl
   return X
end

# function _cartesian_seq(npoint, N, D, r0, r1; inv_t=inv_t)
#    dim = (N*(N-1))÷2
#    ndim = ceil(Int, npoint^(1/dim))
#    x0 = D.transform(r0)
#    x1 = D.transform(r1)
#    xx = linspace(x0, x1, ndim) |> collect
#    oo = ones(ndim)
#    if N == 2
#       return xx
#    end
#    if N == 3
#       Xmat = [kron(xx, oo, oo) kron(oo, xx, oo) kron(oo, oo, xx)]'
#       X = vecs(Xmat)
#       X1 = X[ [ is_simplex(inv_t.(x)) for x in X ] ]
#       return X1
#    end
#    error("`_cartesian_seq` is only implemented for N = 2, 3")
# end

function _z2cart!(z, R)
   R[1] = SVector(0.0, 0.0, 0.0)
   R[2] = SVector(z[1], 0.0, 0.0)
   R[3] = SVector(z[2], z[3], 0.0)
   for n = 4:length(R)
      # n = 4 :        4, 5, 6
      # n = 5 :        7, 8, 9
      R[n] = SVector(z[3*(n-3)+1], z[3*(n-3)+2], z[3*(n-3)+3])
   end
   return R
end

Cart2BA(N) = let J = SVector(collect(1:N-1)...)
      return (R -> lengths_and_angles(R[2:end], J))
end

# function _z2cart_filter!(r0, r1, z, R)
#    R = _z2cart!(z, R)
#    r = norm.(R)
#    return (r0 <= minimum(r)) && (maximum(r) <= r1)
# end

"""
Create a Sobol grid in cartesian space which is then converted into a
?Sobol?-like grid in bond-length or bond-angle space.
"""
function filtered_cart_sobol(r1, N, converter, filter; npoints=nothing,
                                                           nfailed=nothing)
   R = [ (@SVector zeros(3)) for n = 1:N ]
   DescType = typeof(vcat(converter(R)...))
   X = Vector{DescType}()
   dim = 3*N-6
   z1 = r1 * SVector(ones(dim)...)
   X = _filtered_cart_sobol_inner(z1, npoints, nfailed, X,
                              z -> converter(_z2cart!(z, R)), filter)
   return X
end

function _filtered_cart_sobol_inner(z1::SVector{DIM, T},
                                    npoints, nfailed, Xout,
                                    transform, filter) where {DIM, T}
   # Sobol sequence in the [x0, x1]^dim hypercube
   s = SobolSeq(DIM, -z1, z1)
   # temporary storage
   t = zero(MVector{DIM, T})
   # generate points
   npts = 0
   nfl = 0
   while (npts < npoints) && (nfl < nfailed)
      next!(s, t)
      desc = transform(t)
      if filter(desc)
         push!(Xout, vcat(desc...))
         npts += 1
      else
         nfl += 1
      end
   end
   @show npts, nfl
   return Xout
end

# ------------------------------------------------------------------
#               FILTERS
# ------------------------------------------------------------------

function cayley_menger(r::SVector{6, T}) where {T}
   A = @SMatrix( T[
      #       r12     r13      r14
      0.0     r[1]^2  r[2]^2   r[3]^2   1.0;
      # r12           r23      r24
      r[1]^2  0.0     r[4]^2   r[5]^2   1.0;
      # r13    r23             r34
      r[2]^2  r[4]^2  0.0      r[6]^2   1.0;
      # r14    r24    r34
      r[3]^2  r[5]^2  r[6]^2   0.0      1.0;
      1.0     1.0     1.0      1.0      0.0 ])
   return det(A)
end

function cayley_menger(r::SVector{10, T}) where {T}
   A = @SMatrix( T[
      0.0     1.0     1.0     1.0      1.0      1.0;
      #                r12     r13      r14     r15
      1.0     0.0     r[1]^2  r[2]^2   r[3]^2   r[4]^2;
      #       r12             r23      r24      r25
      1.0     r[1]^2  0.0     r[5]^2   r[6]^2   r[7]^2;
      #       r13     r23              r34      r35
      1.0     r[2]^2  r[5]^2  0.0      r[8]^2   r[9]^2;
      #       r14     r24     r34               r45
      1.0     r[3]^2  r[6]^2  r[6]^2   0.0      r[10]^2;
      #       r15     r25     r35      r45
      1.0     r[4]^2  r[7]^2  r[6]^2   r[10]^2  1.0 ])
   return det(A)
end

# bond-length filters
bl_is_simplex(r::MVector) = bl_is_simplex(SVector(r))

bl_is_simplex(r::SVector{1}) = true

bl_is_simplex(r::SVector{3}) = (abs(r[2]-r[1]) <= r[3] <= r[2]+r[1])

bl_is_simplex(r::SVector{6}) =
   bl_is_simplex(SVector(r[1], r[2], r[4])) &&   # r12, r13, r23
   bl_is_simplex(SVector(r[1], r[3], r[5])) &&   # r12, r14, r2
   bl_is_simplex(SVector(r[2], r[3], r[6])) &&   # r13, r14, r34
   bl_is_simplex(SVector(r[4], r[5], r[6])) &&    # r23, r24, r34
   (cayley_menger(r) >= 0)

# r3^2 = r1^2 + r2^2 - 2 ψ r1 r2
# cosangle3(r1, r2, r3) = (r1^2 + r2^2 - r3^2) / (2 * r1 * r2)

ba_is_simplex(r::SVector{1}) = true

ba_is_simplex(rθ::SVector{3}) = true
ba_is_simplex(r::SVector{2}, θ::SVector{1}) = true

ba_is_simplex(rψ::SVector{6}) = ba_is_simplex(SVector(rψ[1], rψ[2], rψ[3]),
                                              SVector(rψ[4], rψ[5], rψ[6]) )

# sqrt(2 - 2*w[1]) = √(2 - 2 * R̂₁ ⋅ R̂₂) = √|R̂₁ - R̂₂|^2 = |R̂₁ - R̂₂|
# but we want |R₁-R₂|^2 = r₁^2 - 2 R₁⋅R₂ + r₂^2 + r₂^2
#                       = r₁^2 + r₁^2 - 2 r₁ r₂ w
_w2r(r1, r2, w) = sqrt(r1^2 + r2^2 - 2*r1*r2*w)

ba_is_simplex(r::SVector{3}, w::SVector{3}) =
   bl_is_simplex(SVector(r[1], r[2], r[3],
                         _w2r(r[1], r[2], w[1]),
                         _w2r(r[1], r[3], w[2]),
                         _w2r(r[2], r[3], w[3])))

"""
This 5B ba_is_simplex only implements a necessary condition for a 5B cluster
to be a simplex. I have no idea how close this is to being sharp. Basically,
this just checks that all 5 4B sub-clusters are simplices.
"""
function ba_is_simplex(r::SVector{4}, w::SVector{6})
   # r = r1, r2, r3, r4
   # w = w12, w13, w14, w23, w24, w34
   r1, w1 = SVector(r[1], r[2], r[3]), SVector(w[1], w[2], w[4])
   r2, w2 = SVector(r[1], r[2], r[4]), SVector(w[1], w[3], w[5])
   r3, w3 = SVector(r[1], r[3], r[4]), SVector(w[2], w[3], w[6])
   r4, w4 = SVector(r[2], r[3], r[4]), SVector(w[4], w[5], w[6])
   # r12, r13, r14, r23, r24, r34 => these make one more tetrahedron
   r5 = SVector( _w2r(r[1], r[2], w[1]),
                 _w2r(r[1], r[3], w[2]),
                 _w2r(r[1], r[4], w[3]),
                 _w2r(r[2], r[3], w[4]),
                 _w2r(r[2], r[4], w[5]),
                 _w2r(r[3], r[4], w[6]) )
   return ba_is_simplex(r1, w1) &&
          ba_is_simplex(r2, w2) &&
          ba_is_simplex(r3, w3) &&
          ba_is_simplex(r4, w4) &&
          bl_is_simplex(r5)
end

# WHAT IS THIS SUPPOSED TO BE?!?!
# ba_is_simplex(r::SVector{3}, ψ::SVector{3}) =
#    bl_is_simplex(vcat(SVector(1.0,1.0,1.0), sqrt.(2 - 2*ψ)))

"""
an auxiliary function to help transform between real and transformed
   coordinates where the Sobol sequence is constructed ...
"""
function inv_transform(x::T, r0::T, r1::T, D)::T where {T}
   TOL = 1e-6
   # Secant Bisection Method (transform is monotone)
   # Roots.jl is too slow (why?!?)
   t0 = transform(D, r0) - x
   t1 = transform(D, r1) - x
   r = 0.5*(r0+r1)
   while abs(r1 - r0) > TOL
      r = (r1 * t0 - r0 * t1) / (t0 - t1)
      t = transform(D, r) - x
      if abs(t) < 1e-7
         return r
      end
      if t*t0 > 0
         r0 = r
         t0 = t
      else
         r1 = r
         t1 = t
      end
   end
   return r
end

# using BenchmarkTools, StaticArrays
# r = @SVector rand(6)
# cayley_menger(r)
# @btime cayley_menger($r)
# @btime bl_is_simplex($r)
# @btime bl_is_simplex_2($r)
#
# # cnt = 0
# # for n = 1:200
# #    r = @SVector rand(6)
# #    if (bl_is_simplex(r) && cayley_menger(r) < -0.001) ||
# #             (!bl_is_simplex(r) && cayley_menger(r) > 0.001)
# #       r = round.(r, 4)
# #       cnt += 1
# #       @show r
# #       @show cayley_menger(r), bl_is_simplex(r)
# #    end
# # end
#
# ##
# tic()
# for n = 1:100_000
#    # r = vcat((@SVector rand(3)), 2 * ((@SVector rand(3)) - 0.5))
#    r = 1.0 + (@SVector rand(6))
#    if bl_is_simplex(r) != bl_is_simplex_2(r)
#       r = round.(r, 4)
#       @show r
#       @show cayley_menger(r)
#       @show bl_is_simplex(r), bl_is_simplex_2(r)
#    end
# end
# toc()

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
"""
`Blocking`

Module that contains functions performing the Flyvbjerg-Petersen
(J. Chem. Phys. 91, 461 (1989)) blocking analysis for evaluating
the standard error on a correlated data set. A "M-test" is also
implemented based on Jonsson (Phys. Rev. E 98, 043304, (2018)).

Note: This module is deprecated and will be removed soon. Use [`StatsTools`](@ref) instead.
"""
module Blocking

using DataFrames, Statistics

export autocovariance, covariance
export blocker, blocking, blockingErrorEstimation, mtest
export autoblock, blockAndMTest
export growthWitness, gW, smoothen

# """
# Calculate the variance of the dataset v
# """
# function variance(v::Vector)
#     n = 0::Int
#     sum = 0.0::Float64
#     sumsq = 0.0::Float64
#     for x in v
#         n += 1
#         sum += x
#         sumsq += x^2
#     end
#     return (sumsq-sum^2/n)/(n-1)
# end
#
# """
# Calculate the standard deviation of the dataset v
# """
# sd(v::Vector) = sqrt(variance(v))
#
"""
    se(v::Vector;corrected::Bool=true)
Calculate the standard error of the dataset `v`. If `corrected` is `true`
(the default) then the sum in `std` is scaled with `n-1`, whereas the sum
is scaled with `n` if corrected is `false` where `n = length(v)`.
"""
se(v::Vector;corrected::Bool=true) = std(v;corrected=corrected)/sqrt(length(v))

"""
    blocker(v::Vector) -> new_v::Vector
Reblock the data by successively taking the mean of two adjacent data points to
form a new vector with a half of the `length(v)`. The last data point will be
discarded if `length(v)` is odd.
"""
function blocker(v::Vector{T}) where T
    P = promote_type(T, Float64)
    new_v = Array{P}(undef,(length(v)÷2))
    for i  in 1:length(v)÷2
        new_v[i] = 0.5*(v[2i-1]+v[2i])
    end
    return new_v
end

"""
    blocking(v::Vector; typos = nothing) -> df
Perform a blocking analysis according to Flyvberg and Peterson
[JCP (1989)](http://aip.scitation.org/doi/10.1063/1.457480)
for single data set and return a `DataFrame` with
statistical data for each blocking step. M-test data according to Jonsson
[PRE (2018)](https://link.aps.org/doi/10.1103/PhysRevE.98.043304) is also
provided.

Keyword argument `typos`
* `typos = nothing` - correct all presumed typos.
* `typos = :FP` - use Flyvberg and Peterson (correct) standard error and Jonsson formul for M.
* `typos = :Jonsson` - calculate `M` and standard error as written in Jonsson.
"""
function blocking_old(v::Vector; typos = nothing)
    df = DataFrame(blocks = Int[], mean = Float64[], stdev = Float64[],
                    std_err = Float64[], std_err_err = Float64[], gamma = Float64[], M = Float64[])
    while length(v) >= 2
        n = length(v)
        mean = sum(v)/n
        var = v .- mean
        sigmasq = sum(var.^2)/n # uncorrected sample variance
        if typos === nothing
            gamma = sum(var[1:n-1].*var[2:n])/(n-1)
        else # typos ∈ {:FP, :Jonsson}
            gamma = sum(var[1:n-1].*var[2:n])/n
            # sample covariance ŷ(1) Eq. (6) [Jonsson]
            # but why is the denominator n and not (n-1)????
        end
        mj = n*((n-1)*sigmasq/(n^2)+gamma)^2/(sigmasq^2)
        stddev = sqrt(sigmasq)
        if typos === nothing || typos == :FP
            stderr = stddev/sqrt(n-1) # [F&P] Eq. (28)
        else
            stderr = stddev/sqrt(n) # [Jonsson] Fig. 2
        end
        stderrerr = stderr*1/sqrt(2*(n-1)) # [F&P] Eq. (28)
        v = blocker(v)
        #println(n, mean, stddev, stderr)
        push!(df,(n, mean, stddev, stderr, stderrerr, gamma, mj))
    end
    return df
end

"""
    autocovariance(v::Vector,h::Int; corrected::Bool=true)
``\\hat{\\gamma}(h) =\\frac{1}{n}\\sum_{t=1}^{n-h}(v_{t+h}-\\bar{v})(v_t-\\bar{v})^*``
Calculate the autocovariance of dataset `v` with a delay `h`. If `corrected` is `true`
(the default) then the sum is scaled with `n-h`, whereas the sum
is scaled with `n` if corrected is `false` where `n = length(v)`.
"""
function autocovariance(v::Vector,h::Int; corrected::Bool=true)
    n = length(v)
    mean_v = mean(v)
    covsum = zero(mean_v)
    for i in 1:n-h
        covsum += (v[i]-mean_v)*conj(v[i+h]-mean_v)
    end
    gamma = covsum/ifelse(corrected, (n-h), n)
    return gamma
end

"""
    blocking(v::Vector; corrected::Bool=true) -> df
Perform a blocking analysis according to Flyvberg and Peterson
[JCP (1989)](http://aip.scitation.org/doi/10.1063/1.457480)
for single data set and return a `DataFrame` with
statistical data for each blocking step. M-test data according to Jonsson
[PRE (2018)](https://link.aps.org/doi/10.1103/PhysRevE.98.043304) is also
provided.
If `corrected` is `true` (the default) then the sum in `var` is scaled
with `n-1` and in `autocovariance` is scaled with `n-h`, whereas the sum
is scaled with `n` for both if corrected is `false` where `n = length(v)`.
"""
function blocking(v::Vector{T}; corrected::Bool=true) where T
    P = promote_type(T, Float64)
    df = DataFrame(blocks = Int[], mean = P[], stdev = Float64[],
                    std_err = Float64[], std_err_err = Float64[], gamma = P[], M = P[])
    while length(v) >= 2
        n = length(v) # size of current dataset
        mean_v = mean(v)
        variance = var(v; corrected=corrected) # variance
        gamma = autocovariance(v,1; corrected=corrected) # sample covariance ŷ(1) Eq. (6) [Jonsson]
        mj = n*((n-1)*variance/(n^2)+gamma)^2/(variance^2) # the M value Eq. (12) [Jonsson]
        stddev = sqrt(variance) # standard deviation
        stderr = stddev/sqrt(n) # standard error
        stderrerr = stderr/sqrt(2*(n-1)) # error on standard error Eq. (28) [F&P]
        v = blocker(v) # re-blocking the dataset
        push!(df,(n, mean_v, stddev, stderr, stderrerr, gamma, mj))
    end
    return df
end

"""
    covariance(x::Vector,y::Vector; corrected::Bool=true)
Calculate the covariance between the two data sets `x` and `y` with equal length.
If `corrected` is `true` (the default) then the sum is scaled with
 `n-1`, whereas the sum is scaled with `n` if corrected is `false`
 where `n = length(x) = length(y)`.
"""
function covariance(vi::Vector,vj::Vector; corrected::Bool=true)
    # if length(vi) != length(vj)
    #     @warn "Two data sets with non-equal length! Truncating the longer one."
    #     if length(vi) > length(vj)
    #         vi = vi[1:length(vj)]
    #     else
    #         vj = vj[1:length(vi)]
    #     end
    # end
    n = length(vi)
    meani = mean(vi)
    meanj = mean(vj)
    covsum = zero(meani)
    for i in 1:n
        covsum += (vi[i]-meani)*conj(vj[i]-meanj)
    end
    cov = covsum/ifelse(corrected, (n-1), n)
    return cov
end

"""
    combination_division(x::Vector,y::Vector; corrected::Bool=true)
Find the standard error on the quotient of means `x̄/ȳ` from two data sets,
note that the standard errors are different on ``(x̄/ȳ) \\neq \\bar{(\\frac{x}{y})}``.
If `corrected` is `true` (the default) then the sums in both variance and covariance
are scaled with `n-1`, whereas the sums are scaled with `n` if corrected is `false`
 where `n = length(x) = length(y)`.
"""
function combination_division(vi::Vector,vj::Vector; corrected::Bool=true)
    # if length(vi) != length(vj)
    #     @warn "Two data sets with non-equal length! Truncating the longer one."
    #     if length(vi) > length(vj)
    #         vi = vi[1:length(vj)]
    #     else
    #         vj = vj[1:length(vi)]
    #     end
    # end
    n = length(vi)
    meani = mean(vi)
    meanj = mean(vj)
    meanf = meani/meanj
    sei = se(vi;corrected=corrected)
    sej = se(vj;corrected=corrected)
    cov = covariance(vi,vj;corrected=corrected)
    sef = abs(meanf*sqrt((sei/meani)^2 + (sej/meanj)^2 - 2.0*cov/(n*meani*meanj)))
    return sef
end

"""
    blocking(x::Vector,y::Vector) -> df::DataFrame
Perform a blocking analysis for the quotient of means `x̄/ȳ` from two data sets.
If `corrected` is `true` (the default) then the sums in both variance and covariance
are scaled with `n-1`, whereas the sums are scaled with `n` if corrected is `false`
 where `n = length(x) = length(y)`.
Entries in returned dataframe:
* `blocks` = number of blocks in current blocking step;
* `mean_x`, `SD_x`, `SE_x`, `SE_SE_x` = the mean, standard deviation, standard error and error on standard error estimated for dataset `x`;
* `mean_y`, `SD_y`, `SE_y`, `SE_SE_y` = ditto. for dataset `y`;
* `Covariance` = the covariance between data in `x` and `y`;
* `mean_f` = `x̄/ȳ`;
* `SE_f` = standard error estimated for `x̄/ȳ`.
"""
function blocking(vi::Vector{T1},vj::Vector{T2}; corrected::Bool=true) where {T1, T2}
    P = promote_type(T1, T2, Float64)
    df = DataFrame(blocks=Int[], mean_x=P[], SD_x=Float64[], SE_x=Float64[], SE_SE_x=Float64[],
            mean_y=P[], SD_y=Float64[], SE_y=Float64[], SE_SE_y=Float64[], Covariance=P[],
            mean_f=P[], SE_f=Float64[])
    # if length(vi) != length(vj)
    #     @warn "Two data sets with non-equal length! Truncating the longer one."
    #     if length(vi) > length(vj)
    #         vi = vi[1:length(vj)]
    #     else
    #         vj = vj[1:length(vi)]
    #     end
    # end
    while length(vi) >= 2
        n = length(vi)
        meani = mean(vi)
        meanj = mean(vj)
        meanf = meani/meanj
        sdi = std(vi;corrected=corrected)
        sdj = std(vj;corrected=corrected)
        sei = sdi/sqrt(n)
        sej = sdj/sqrt(n)
        sesei = sei*1/sqrt(2*(n-1))
        sesej = sej*1/sqrt(2*(n-1))
        cov = covariance(vi,vj;corrected=corrected)
        #sef = sei/sej
        sef = combination_division(vi,vj;corrected=corrected)
        vi = blocker(vi)
        vj = blocker(vj)
        #println(n, mean, stddev, stderr)
        push!(df,(n, meani, sdi, sei, sesei, meanj, sdj, sej, sesej, cov, meanf, sef))
    end
    return df
end

# no longer needed, using the M test now
# """
# estimating stnadard error from blocking analysis based on the overlapping of
# error bars, if all the error bars (or more than 3 on a roll) behind current
# one are overlapping with it, return the current standard error with error bar.
# """
# function blockingErrorEstimation(df::DataFrame)
#     e = df.std_err[1:end-1] # ignoring the last data point
#     ee = df.std_err_err[1:end-1] # ignoring the last data point
#     n = length(e)
#     ind = collect(1:length(e))
#     e_upper = map(x->e[x]+ee[x],ind) # upper bounds
#     e_lower = map(x->e[x]-ee[x],ind) # lower bounds
#     i = 1 # start from the first data point
#     plateau = false
#     while i < n
#         count = 0 # set up a counter for checking overlapped error bars
#         for j in (i+1):n # j : all data points after i
#             if e_lower[i] >= e_lower[j] && e_upper[i] <= e_upper[j]
#                 count += 1
#                 #println("i: ",i," j: ",j," c: ",count)
#                 # some tolerance, say if there are 3 overlaps on a roll could be a plateau
#                 if count > 3 && (i + count) == j
#                     plateau = true
#                     println("\x1b[32mplateau detected\x1b[0m")
#                     return e[i], ee[i], plateau
#                 end
#             end
#         end # for
#         if count == (n-i)
#             println("\x1b[32mNO plateau is detected, take the best estimation\x1b[0m")
#             return e[i], ee[i], plateau
#         else
#             i += 1 # move on to next point
#         end
#     end # while
#     println("\x1b[32mNO plateau, NO error bar overlap, take the second last point\x1b[0m")
#     return e[i], ee[i], plateau # return the last ponit
# end

"""
    mtest(df::DataFrame; warn = true) -> k
The "M test" based on Jonsson, M. Physical Review E, 98(4), 043304, (2018).
Expects `df` to be output of a blocking analysis with column `df.M` containing
relevant M_j values, which are compared to a χ^2 distribution.
Returns the row number `k` where the M-test is passed.
If the M-test has failed `mtest()` returns the value `-1` and optionally prints
a warning message.
"""
function mtest(df::DataFrame; warn = true)
    # the χ^2 99 percentiles
    q = [6.634897,  9.210340,  11.344867, 13.276704, 15.086272,
        16.811894, 18.475307, 20.090235, 21.665994, 23.209251,
        24.724970, 26.216967, 27.688250, 29.141238, 30.577914,
        31.999927, 33.408664, 34.805306, 36.190869, 37.566235,
        38.932173, 40.289360, 41.638398, 42.979820, 44.314105,
        45.641683, 46.962942, 48.278236, 49.587884, 50.892181]
    Mj = df.M
    M = reverse(cumsum(reverse(Mj)))
    #println(M)
    k = 1
    while k <= length(M)-1
       if M[k] < q[k]
           # if info
           #     stder = round(df.std_err[k],digits=3)
           #     stderer = round(df.std_err_err[k],digits=3)
           #     println("\x1b[32mM test passed, the smallest k is $k\x1b[0m")
           #     println("\x1b[32mStandard error estimation: $stder ± $stderer\x1b[0m")
           # end
           return k
       else
           k += 1
       end
    end
    if warn
        @warn "M test failed, more data needed"
    end
    return -1 # indicating the the M-test has failed
end

# using Statistics, StatsBase
#
# function mm(x::Vector, x̄)
#     n = length(x)
#     σ2 = varm(x,x̄, corrected = false)
#     γ = cov(x[1:n-1],x[2:n], corrected = false)
#     m = n * ((n - 1) * σ2 + γ)^2 / σ2
#     return m
# end
#
# """
#     autoblock(x,y)
# Perform automated blocking analysis for `x̄/ȳ`.
# """
# function block2(x::Vector,y::Vector)
#     n = length(x)
#     @assert length(y)==n "Vectors do not have the same length"
#     ms = Vector{Float64}(undef,trunc(Int,log2(n)))
#     blocking_step = 1
#     while n ≥ 2
#         x̄ = mean(x)
#         ȳ = mean(y)
#         f̄ = x̄/ȳ
#         σ2x = varm(x,x̄, corrected = false)
#         γx = cov(x[1:n-1],x[2:n], corrected = false)
#         mx = mm(x, x̄)
#         my = mm(y, ȳ)
#         sef = combination_division(x,y)
#         x = blocker(x)
#         y = blocker(y)
#         n = length(x)
#         blocking_step += 1
#     end
# end

"""
    v̄, σ, σσ, k, df = blockAndMTest(v::Vector)
Perform a blocking analysis and M-test on `v` returning the mean `v̄`,
standard error `σ`, its error `σσ`, the number of blocking steps `k`, and
the `DataFrame` `df` with blocking data.
"""
function blockAndMTest(v::Vector; corrected::Bool=true)
    df = blocking(v;corrected=corrected)
    k = mtest(df, warn=false)
    v̄ = df.mean[1]
    if k>0
        σ = df.std_err[k]
        σσ = df.std_err_err[k]
    else
        @warn "M test failed, more data needed"
        σ = maximum(df.std_err)
        σσ = maximum(df.std_err_err)
    end
    return v̄, σ, σσ, k, df
end

function blockTestShiftAndProjected(df::DataFrame; start = 1, stop = size(df)[1], corrected::Bool=true)
    s̄, σs, σσs, ks, dfs = blockAndMTest(df.shift[start:stop];corrected=corrected)
    v̄, σv, σσv, kv, dfv = blockAndMTest(df.vproj[start:stop];corrected=corrected)
    h̄, σh, σσh, kh, dfh = blockAndMTest(df.hproj[start:stop];corrected=corrected)
    dfp = blocking(df.hproj[start:stop], df.vproj[start:stop];corrected=corrected)
    k = max(ks, kv, kh)
    @show ks, kv, kh
    ks==kv==kh || @warn "k values are not the same."
    return s̄, σs, dfp.mean_f[1], dfp.SE_f[k], ks, kv, kh
end

"""
    autoblock(df::DataFrame; start = 1, stop = size(df)[1])
    -> s̄, σs, ē, σe, k
Determine mean shift `s̄` and projected energy `ē` with respective standard
errors `σs` and `σe` by blocking analsis from the `DataFrame` `df` returned
from `fciqmc!()`. The number `k` of blocking
steps and decorrelation time `2^k` are obtained from the M-test for the
shift and also applied to the projected energy, assuming that the projected
quantities decorrelate on the same time scale. Only the real part of the shift
is considered. Returns a named tuple.
"""
function autoblock(df::DataFrame; start = 1, stop = size(df)[1], corrected::Bool=true)
    s̄, σs, σσs, ks, dfs = blockAndMTest(real.(df.shift[start:stop]);corrected=corrected) # shift
    if hasproperty(df, :hproj)
        dfp = blocking(df.hproj[start:stop], df.vproj[start:stop];corrected=corrected) # projected
        return (s̄ = s̄, σs = σs, ē = dfp.mean_f[1], σe = dfp.SE_f[ks], k = ks)
    else
        return (s̄ = s̄, σs = σs, ē = missing, σe = missing, k = ks)
    end
end

# call signature for chaining
autoblock(nt::NamedTuple{(:df,:eqsteps)}) = autoblock(nt.df; start = nt.eqsteps)

# version for replica run
"""
    autoblock(dftup::Tuple; start = 1, stop = size(dftup[1])[1])
    -> s̄1, σs1, s̄2, σs2, ē1, σe1, ē2, σe2, ēH, σeH, k
Replica version. `dftup` is the tuple of `DataFrame`s returned from replica
`fciqmc!()`. Returns a named tuple with shifts and three variational energy
estimators and respective errors obtained from blocking analysis. The larger
of the `k` values from M-tests on the two shift time series is used.
"""
function autoblock(dftup::Tuple; start = 1, stop = size(dftup[1])[1], corrected::Bool=true)
    (df_mix, (df_1, df_2)) = dftup # unpack the three DataFrames
    s̄1, σs1, σσs1, ks1, dfs1 = blockAndMTest(real.(df1.shift[start:stop]);corrected=corrected) # shift 1
    s̄2, σs2, σσs2, ks2, dfs2 = blockAndMTest(real.(df2.shift[start:stop]);corrected=corrected) # shift 2
    xdy = df_mix.xdy[start:stop]
    s1_xdy = real.(dfs1.shift[start:stop]).*xdy
    s2_xdy = real.(dfs2.shift[start:stop]).*xdy
    xHy = df_mix.xHy[start:stop]
    df_var_1 = blocking(s1_xdy, xdy;corrected=corrected)
    df_var_2 = blocking(s2_xdy, xdy;corrected=corrected)
    df_var_H = blocking(xHy, xdy;corrected=corrected)
    dfp = blocking(df.hproj[start:stop], df.vproj[start:stop];corrected=corrected)
    ks = max(ks1, ks2)
    return (s̄1=s̄1, σs1=σs1, s̄2=s̄2, σs2=σs2,
        ē1 = df_var_1.mean_f[1], σe1 = df_var_1.SE_f[ks],
        ē2 = df_var_2.mean_f[1], σe2 = df_var_2.SE_f[ks],
        ēH = df_var_H.mean_f[1], σeH = df_var_H.SE_f[ks], k = ks)

end

#G_b^{(n)} = \\bar{S}^{(n)} - \\frac{\\log\\vert\\mathbf{c}^{(n+b)}\\vert - \\log\\vert\\mathbf{c}^{(n)}\\vert}{b d\\tau},
"""
    growthWitness(norm::AbstractArray, shift::AbstractArray, dt; b = 30, pad = :true) -> g
    growthWitness(df::DataFrame; b = 30, pad = :true) -> g
Compute the growth witness
```math
G_b^{(n)} = S̄^{(n)} - \\frac{\\log\\vert\\mathbf{c}^{(n+b)}\\vert - \\log\\vert\\mathbf{c}^{(n)}\\vert}{b d\\tau},
```
where `S̄` is an average of the `shift` over `b` time steps and \$\\vert\\mathbf{c}^{(n)}\\vert ==\$ `norm[n]`.
The parameter `b ≥ 1` averages the derivative quantity over `b` time steps and helps suppress noise.

If `pad` is set to `:false` then the returned array `g` has the length `length(norm) - b`.
If set to `:true` then `g` will be padded up to the same length as `norm` and `shift`.
"""
function growthWitness(norm::AbstractArray, shift::AbstractArray, dt; b = 30, pad = :true)
    l = length(norm)
    @assert length(shift) == l "`norm` and `shift` arrays need to have the same length."
    n = l - b
    g = Vector{eltype(shift)}(undef, pad ? l : n)
    offset = pad ? b÷2 : 0 # use offset only if pad == :true
    for i in 1:n
        g[i + offset] = -(1/(b*dt) * log(norm[i+b]/norm[i]) - 1/(b+1) * sum(shift[i:i+b]))
    end
    if pad # pad the vector g at both ends
        g[1:offset] .= g[offset]
        g[offset+n+1 : end] .= g[offset+n]
    end
    return g
end
growthWitness(df::DataFrame; b = 30, pad = :true) = growthWitness(df.norm, df.shift, df.dτ[1]; b=b, pad=pad)

"""
    smoothen(noisy::AbstractVector, b; pad = :true)
Smoothen the array `noisy` by averaging over a sliding window of length `b`.
Pad to `length(noisy)` if `pad == true`. Otherwise, the returned array will have
the length `length(noisy) - b`.
"""
function smoothen(noisy::AbstractVector, b; pad = :true)
    l = length(noisy)
    n = l - b
    smooth = Vector{promote_type(eltype(noisy),Float64)}(undef, pad ? l : n)
    offset = pad ? b÷2 : 0 # use offset only if pad == :true
    for i in 1:n
        smooth[i + offset] = 1/b * sum(noisy[i:i+b-1])
    end
    if pad # pad the vector g at both ends
        smooth[1:offset] .= smooth[offset+1]
        smooth[offset+n+1 : end] .= smooth[offset+n]
    end
    return smooth
end

"""
    gW(norm::AbstractArray, shift::AbstractArray, dt, [b]; pad = :true) -> g
    gW(df::DataFrame, [b]; pad = :true) -> g
Compute the growth witness
```math
G^{(n)} = S^{(n)} - \\frac{\\vert\\mathbf{c}^{(n+1)}\\vert - \\vert\\mathbf{c}^{(n)}\\vert}{\\vert\\mathbf{c}^{(n)}\\vert d\\tau},
```
where `S` is the `shift` and \$\\vert\\mathbf{c}^{(n)}\\vert ==\$ `norm[n, 1]`.
Setting `b ≥ 1` a sliding average over `b` time steps is computed.

If `pad` is set to `:false` then the returned array `g` has the length `length(norm) - b`.
If set to `:true` then `g` will be padded up to the same length as `norm` and `shift`.
"""
function gW(norm::AbstractArray, shift::AbstractArray, dt)
    l = length(norm)
    @assert length(shift) == l "`norm` and `shift` arrays need to have the same length."
    n = l - 1
    g = Vector{eltype(shift)}(undef, l)
    for i in 1:n
        g[i] = shift[i] - (norm[i+1] - norm[i])/(dt*norm[i])
    end
    # pad the vector g at the end
    g[n+1] = g[n]
    return g
end
function gW(norm::AbstractArray, shift::AbstractArray, dt, b; pad = :true)
    g_raw = gW(norm, shift, dt)
    return smoothen(g_raw, b; pad)
end
gW(df::DataFrame) = gW(df.norm, df.shift, df.dτ[1])
gW(df::DataFrame, b; pad = :true) = gW(df.norm, df.shift, df.dτ[1], b; pad=pad)

end # module Blocking


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 = "KaHyPar"
version = v"1.2.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/kahypar/kahypar.git", "971c38a4c7b7ed1ecf001d0263c1099cbc72886e"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd kahypar/
git submodule update --init --recursive --depth=1
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make install.library
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "musl"),
    Platform("x86_64", "macos"; ),
    Platform("x86_64", "freebsd"; )
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libkahypar", :libkahypar)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Boost breaks ABI in every single version because they embed the full version number in
    # the SONAME, so we're compatible with one and only one version at a time.
    Dependency(PackageSpec(name="boost_jll", uuid="28df3c45-c428-5900-9ff8-a3135698ca75"); compat="=1.71.0")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    encode(s::AbstractString)

Encode string to base64.
"""
function encode(s::AbstractString)
    io = IOBuffer()
    io64_encode = Base64EncodePipe(io)
    for char in s
        write(io64_encode, char)
    end
    close(io64_encode)
    String(take!(io))
end

"""
    get(api::AbstractString)

Send GET request to api.
"""
function get(api::AbstractString)
    get(api, default())
end

"""
    get(api::AbstractString, a::Authorized)

Send GET request to api.
"""
function get(api::AbstractString, a::Authorized)
    resp = HTTP.request("GET", OATH_URL*api,
        ["Authorization" => "bearer "*a.token,
        "User-Agent" => a.useragent])
    sleep(1) # wait for 1 second to ensure no more than 60 requests per minute
    String(resp.body)
end

"""
    post(api::AbstractString,
         data::AbstractString)

Send POST request to API, using default credentials.
"""
function post(api::AbstractString, body::AbstractString)
    post(api, body, default())
end

"""
    post(api::AbstractString,
         data::AbstractString,
         a::AbstractCredentials)

Send POST request to API.
"""
function post(api::AbstractString, body::AbstractString, a::Authorized)
    resp = HTTP.request("POST", OATH_URL*api,
                        ["Authorization" => "bearer "*a.token,
                         "User-Agent" => a.useragent],
                         body)
    sleep(1) # wait for 1 second to ensure no more than 60 requests per minute
    String(resp.body)
end


text:
# I think it makes more sense to let the array into the function.
# That way the coefficient conversions happen once.
#=
function cauchy{S,T}(f::Fun{S,T},z::Array,s...)
    ret=Array(Complex{Float64},size(z)...)
    for k=1:size(z,1),j=1:size(z,2)
        @inbounds ret[k,j]=cauchy(f,z[k,j],s...)
    end
    ret
end
=#

for op in (:(stieltjes),:(cauchy),:(logkernel),:(stieltjesintegral),:(cauchyintegral))
    @eval begin
        $op{F<:Fun}(v::Vector{F},z)=mapreduce(f->$op(f,z),+,v)
        $op{F<:Fun}(v::Vector{F})=map($op,v)
        $op(v::Vector{Any},z)=mapreduce(f->$op(f,z),+,v)
        $op(S::PiecewiseSpace,v,z)=$op(pieces(Fun(v,S)),z)
        $op{S<:PiecewiseSpace,T}(f::Fun{S,T}) = (v = $op(pieces(f)); Fun(vec(coefficientmatrix(v).'),ApproxFun.SumSpace(map(space,v))))
        $op(S::PiecewiseSpace,v)=depiece($op(pieces(Fun(v,S))))

        $op{F<:Fun}(v::Vector{F},z,s::Bool)=mapreduce(f->(z in domain(f))?$op(f,z,s):$op(f,z),+,v)
        $op(v::Vector{Any},z,s::Bool)=mapreduce(f->(z in domain(f))?$op(f,z,s):$op(f,z),+,v)
        $op(S::PiecewiseSpace,v,z,s::Bool)=$op(pieces(Fun(v,S)),z,s)
    end
end

hilbert{F<:Union{Fun,Any}}(v::Vector{F},x)=mapreduce(f->(x in domain(f))?hilbert(f,x):-stieltjes(f,x)/π,+,v)
hilbert(S::PiecewiseSpace,v,z)=hilbert(pieces(Fun(v,S)),z)

function cauchy{S<:ArraySpace,T}(v::Fun{S,T},z::Number,s...)
    m=mat(v)
    Complex{Float64}[cauchy(m[k,j],z,s...) for k=1:size(m,1),j=1:size(m,2)]
end


text:
module AlgebraicContracts

# Submodules
include("StaticContracts.jl")
include("ContractedMachine.jl")

end # module


text:
x = 10
y = 20
sum(abs, Int64[x,y])

max_stars_repo_path: My First Julia Scribt.jl

max_stars_repo_name: Ashmaawyy/My-First-Julia-Project

max_issues_repo_path: My First Julia Scribt.jl

max_issues_repo_name: Ashmaawyy/My-First-Julia-Project

max_forks_repo_path: My First Julia Scribt.jl

max_forks_repo_name: Ashmaawyy/My-First-Julia-Project


text:
using Pages,Blink

using Compat; import Compat.readstring

Pages.start()

Endpoint("/examples/pages") do request::Request
    readstring(joinpath(dirname(@__FILE__),"PagesJL.html"))
end

function example_pages()
    nwin = 3
    for i = 1:nwin
        w = Window(Dict(:title => "Pages: Window #$i", :url => "http://localhost:8000/examples/pages"))
        tools(w)
        wait(Pages.conditions["connected"])
    end

    Pages.broadcast("say","Hello everyone!")
    for (sid,s) in Pages.sessions
        c = s.client
        Pages.message(c,"say","You are connection #$(c.id).")
    end
end
example_pages()


text:
using ManifoldLearning
import StatsBase

import NearestNeighbors
struct KDTree <: ManifoldLearning.AbstractNearestNeighbors
    k::Integer
    fitted::NearestNeighbors.KDTree
end
Base.show(io::IO, NN::KDTree) = print(io, "KDTree(k=$(NN.k))")
StatsBase.fit(::Type{KDTree}, X::AbstractMatrix{T}, k::Integer) where {T<:Real} = KDTree(k, NearestNeighbors.KDTree(X))
function ManifoldLearning.knn(NN::KDTree, X::AbstractVecOrMat{T}; self=false) where {T<:Real}
    m, n = size(X)
    k = NN.k
    @assert n > k "Number of observations must be more then $(k)"

    idxs, dist = NearestNeighbors.knn(NN.fitted, X, k+1, true)
    D = Array{T}(undef, k, n)
    E = Array{Int32}(undef, k, n)
    for i in eachindex(idxs)
        E[:, i] = idxs[i][2:end]
        D[:, i] = dist[i][2:end]
    end
    return D, E
end

import FLANN
struct FLANNTree{T <: Real} <: ManifoldLearning.AbstractNearestNeighbors
    k::Integer
    index::FLANN.FLANNIndex{T}
end
Base.show(io::IO, NN::FLANNTree) = print(io, "FLANNTree(k=$(NN.k))")
function StatsBase.fit(::Type{FLANNTree}, X::AbstractMatrix{T}, k::Integer) where {T<:Real}
    params = FLANNParameters()
    idx = FLANN.flann(X, params)
    FLANNTree(k, idx)
end
function ManifoldLearning.knn(NN::FLANNTree, X::AbstractVecOrMat{T}; self=false) where {T<:Real}
    E, D = FLANN.knn(NN.index, X, NN.k+1)
    sqrt.(@view D[2:end, :]), @view E[2:end, :]
end