text:
using StanSample, MCMCChains

# Define the Stan language model

m8_2s = "
data{
    int N;
    vector[N] y;
}
parameters{
    real mu;
    real<lower=0> sigma;
}
model{
    y ~ normal( mu , sigma );
}
";

# Define the Stanmodel.

m_8_2s = SampleModel("m8.2s", m8_2s);

# Input data for cmdstan

m8_2_data = Dict("N" => 2, "y" => [-1, 1]);
m8_2_init = Dict("mu" => 0.0, "sigma" => 1.0);

# Sample using cmdstan

rc = stan_sample(m_8_2s; data=m8_2_data, init=m8_2_init);

# Describe the draws
if success(rc)
  chn = read_samples(m_8_2s; output_format=:mcmcchains)
  describe(chn)
end

ext: jl


text:
@doc """
	boundary_coef_mat(F::BoundaryFilter) -> Matrix

The boundary coefficients collected in a matrix where the `i`'th row
contains the coefficients of the `i`'th boundary scaling function.
"""->
function boundary_coef_mat(F::BoundaryFilter)
	vm = van_moment(F)
	coef_mat = zeros(Float64, vm, vm)

	for row in 1:vm
		coef = bfilter(F, row-1)
		@inbounds coef_mat[row,:] = sqrt2*coef[1:vm]
	end

	return coef_mat
end

@doc """
	DaubScaling(p::Integer, side::Char, R::Integer) -> x, Y

Compute the boundary scaling function at resolution `R` on `side` with `p` vanishing moments.
"""->
function DaubScaling(p::Integer, side::Char, R::Integer)
	B = bfilter(p, side)
	IF = ifilter(p, true)

	x = dyadic_rationals( support(B), R )
	Y = DaubScaling(B, IF, R)

	return x, Y'
end

#=
	DaubScaling(B::BoundaryFilter) -> Vector

Compute the boundary scaling function values at 0.
=#
function DaubScaling(B::BoundaryFilter)
	boundary_mat = boundary_coef_mat(B)
	E = eigval1(boundary_mat)
	# TODO: How to scale E?

	return E
end

#=
	DaubScaling(B, I) -> Matrix

Compute the boundary scaling function defined by boundary filter `B` and
interior filter `I` at the non-zero integers in their support.

The ouput is a matrix where the `k`'th row are the functions values of the `k-1` scaling function.
=#
function DaubScaling(B::BoundaryFilter, IF::InteriorFilter)
	internal = DaubScaling(IF)
	IS = support(IF)
	BS = support(B)

	# Function values at 0
	vm = van_moment(B)
	xvals = integers(B)
	Y = zeros(Float64, vm, length(xvals)+1)
	Y[:, x2index(0,BS)] = DaubScaling(B)

	xvals = integers(B)
	# The translations of the interior scaling function differ for the two sides
	interior_start = ( side(B) == 'L' ? vm : -vm-1 )
	interior_iter = ( side(B) == 'L' ? 1 : -1 )

	for x in xvals
		xindex = x2index(x, BS)
		doublex = 2*x
		doublex_index = x2index(doublex, BS)

		for k in 1:vm
			filterk = bfilter(B, k-1)

			# Boundary contribution
			if isinside(doublex, BS)
				for l in 1:vm
					Y[k, xindex] += sqrt2 * filterk[l] * Y[l, doublex_index]
				end
			end

			# Interior contribution
			fidx = vm
			interior_arg = doublex - interior_start
			flength = vm + 2*(k-1) #length(filterk)
			while fidx <= flength
				fidx += 1
				if isinside(interior_arg, IS)
					Y[k,xindex] += sqrt2 * filterk[fidx] * internal[ x2index(interior_arg,IS) ]
				end
				interior_arg -= interior_iter
			end
		end
	end

	return Y
end

@doc """
	DaubScaling(B, I, R) -> Matrix

Compute the boundary scaling function defined by boundary filter `B` and interior filter `I` at the dyadic rationals up to resolution `R` in their support.

The ouput is a matrix where the `k`'th row are the functions values of the `k-1` scaling function.
"""->
function DaubScaling(B::BoundaryFilter, IF::InteriorFilter, R::Int)
	R >= 0 || throw(DomainError())

	internal = DaubScaling(IF,R)
	Ny = length(internal)
	vm = van_moment(B)
	Y = zeros(Float64, vm, Ny)

	IS = support(IF)
	BS = support(B)

	# Base level
	cur_idx = dyadic_rationals(BS, R, 0)
	Y[:,cur_idx] = DaubScaling(B, IF)

	interior_start = ( side(B) == 'L' ? vm : -vm-1 )
	interior_iter = ( side(B) == 'L' ? 1 : -1 )

	# Recursion: Fill remaining levels
	xvals = dyadic_rationals(BS, R)
	for L in 1:R
		# Indices of x values on scale L
		cur_idx = dyadic_rationals(BS, R, L)

		for xindex in cur_idx
			doublex = 2*xvals[xindex]
			doublex_index = x2index(doublex, BS, R)

			for k in 1:vm
				filterk = bfilter(B, k-1)

				# Boundary contribution
				if isinside(doublex, BS)
					for l in 1:vm
						Y[k,xindex] += sqrt2 * filterk[l] * Y[l,doublex_index]
					end
				end

				# Interior contribution
				fidx = vm
				interior_arg = doublex - interior_start
				flength = vm + 2*(k-1) #length(filterk)
				while fidx <= flength
					fidx += 1
					if isinside(interior_arg, IS)
						Y[k,xindex] += sqrt2 * filterk[fidx] * internal[ x2index(interior_arg,IS,R) ]
					end
					interior_arg -= interior_iter
				end
			end
		end
	end

	# Hack until the right 0 values are computed
	# TODO: Fix this
	if side(B) == 'L'
		Y[:,1] = Y[:,2]
	else
		Y[:,end] = Y[:,end-1]
	end

	return Y
end


text:
# 1 NaN chunk, no NaNs at edges
x  = [1.0, 2.0, NaN, NaN, 2.3, 5.6]
fs = findall_nan_chunks(x)

@test length(fs) == 1
@test fs[1] == (3, 4)

# > 1 NaN chunk, no NaNs at edges

x  = [1.0, 2.0, NaN, NaN, 2.3, 5.6, NaN, NaN, NaN, 2.0]
fs = findall_nan_chunks(x)
@test length(findall_nan_chunks(x)) == 2
@test fs[1] == (3, 4)
@test fs[2] == (7, 9)

# NaN chunk at left boundary
x = [NaN, 1.0, 2.0, NaN, NaN, 2.3, 5.6, NaN, NaN, NaN, 2.0]
fs = findall_nan_chunks(x)
fs[1] == (1, 1)
fs[2] == (4, 5)
fs[3] == (8, 10)

# NaN chunk at right boundary
x = [2.3, 5.6, NaN, NaN, NaN]
fs = findall_nan_chunks(x)
fs[1] == (3, 5)

x  = [NaN, 1.0, 2.0, NaN, NaN, 2.3, NaN, 5.6, 8.7, NaN]
fc = findall_nan_chunks(x)
@test fc[1] == (1, 1)
@test fc[2] == (4, 5)
@test fc[3] == (7, 7)
@test fc[4] == (10, 10)
@test length(fc) == 4

# The number of NaNs in the ranges matches the total number of NaNs
x = rand(10000)
x[rand(1:1000, 1000)] .= NaN
tupdiff(t::Tuple{Int, Int}) = t[2] - t[1]
sum(tupdiff.(findall_nan_chunks(x)) .+ 1) == length(findall(isnan.(x)))

# In-place method matches regular method
v = zeros(Bool, length(x));
findall_nan_chunks!(v, x) == findall_nan_chunks(x)


text:
# v0.6.0

cp = deepcopy(obj)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Mongo, LibBSON
using Caesar

include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
# cloudGraph, addrdict = standardcloudgraphsetup()

## interactive operation
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

# Get the keys (and possibly the data)
# cursor = find(cloudGraph.mongo.cgBindataCollection, query())
# for o in cursor
#     # println( ", oid: ", o["_id"])
#     println("id: ", o["_id"]);
# end

# Get Roxana's encoded parameters
#REF: https://github.com/dehann/Caesar.jl/blob/master/examples/database/python/mongo_interaction.py
# mongoClient = cloudGraph.mongo.client
# mongoTestCollection = MongoCollection(mongoClient, "tester", "testing_collection")
# cursor = find(mongoTestCollection, query())
# for o in cursor
#     println("Key: ", o["key"], ", oid: ", o["_id"])  # Or simply,
#     #println(o)
# end

# example of string array
# myKeyToFind = BSONOID("58b0a34b5d76256bd2b8183c") # all NaNs
myKeyToFind = BSONOID("58b0a3525d76256bd2b81860") # some valid numbers

findsomthing = find(cloudGraph.mongo.cgBindataCollection, ("_id" => eq(myKeyToFind)))
myFavouriteKey = first( findsomthing );
mfkv = myFavouriteKey["val"];
# @show myFavouriteKey["val"]

ptrT = pointer(mfkv)
# Ptr{UInt8} @0x00000000080f9740
#
ptrTf = convert(Ptr{Float32}, ptrT)
# Ptr{Float32} @0x00000000080f9740

len = length(mfkv)
arr = Vector{Float32}(round(Int,len/4));

unsafe_copy!(pointer(arr),ptrTf,round(Int,len/4))

# fetch a binary png image
myKeyToFind = BSONOID("58af67255d7625647859fa71") # good example of opencv binary encoded png
findsomthing = find(cloudGraph.mongo.cgBindataCollection, ("_id" => eq(myKeyToFind)))
# @show myFavouriteKey["val"]
myFavouriteKey = first( findsomthing );
gf
file = open("test.png", "w")
write(file, myFavouriteKey["val"])
close(file)

run(`eog test.png`)

data = myFavouriteKey["val"]
using ImageMagick
img = ImageMagick.readblob(data);
# using Images
# using FileIO
# using Images, ImageView
# img = load("test.png")
# imshow(img)
# run(`eog test.png`)
# run(`rm test.png`)

img[1,1].r
imgA = zeros(480,640,3);
@time for i in 1:480, j in 1:640
  imgA[i,j,1] = img[i,j].r
  imgA[i,j,2] = img[i,j].g
  imgA[i,j,3] = img[i,j].b
end

# limitations in LibBSON.jl prevent Binary unpacking at this time, using PyCall to pymongo interim
# key = BSONOID("58c1a1bdb2ca1534d30b1c96")
# obj = find(cloudGraph.mongo.cgBindataCollection, ( "_id" => eq(key) ) )

# fetching data from Mongo using getimages.py
# IDs = getPoseExVertexNeoIDs(fg.cg.neo4j.connection, session=attrdict["session"], reqbackendset=false);
# mongk = getmongokeys(fg, :x1, IDs)
# mongo_key = bson.ObjectId(mongk["keyframe_rgb"])
# arr = gi.getbinarray(dbcoll, mongo_key)

myKeyToFind = BSONOID("58e302546689bc617c1994c1") # good example of opencv binary encoded png
findsomthing = find(cloudGraph.mongo.cgBindataCollection, ("_id" => eq(myKeyToFind)))
# @show myFavouriteKey["val"]
mkf = first( findsomthing );

typeof(mkf["val"])

buf = IOBuffer(mkf["val"])

st = String(take!(buf))
bb = BSONObject(st)
typeof(bb)
bba = convert(Array, bb["pointcloud"][1])
convert(bb["pointcloud"]

# temporary helper function to read binary BSON data via pymongo
using PyCall
@pyimport pymongo

client = pymongo.MongoClient(addrdict["mongo addr"])
db = client[:CloudGraphs]
collection = "bindata"

IDs = getPoseExVertexNeoIDs(fg.cg.neo4j.connection, session=attrdict["session"], reqbackendset=false);
mongk = getmongokeys(fg, :x1, IDs)
arr = getbinarraymongo(db[collection], mongk["arrnameinneo4j"] )
# arr = getbinarraymongo(db[collection], "58c1a1bdb2ca1534d30b1c96" )

#


text:
#  A delay eigenvalue problem with one eigenvalue equal to one.
function dep1()
    A0=([1 2 3 ; 4 5 6; 1 -1 3]);
    A1=((-A0+[1 0 3;0 0 -1;0 0 10])*exp(1));
    Q=[1 0 3; 1 1 -4; 2 3 1];
    A0=Q\(A0*Q);
    A1=Q\(A1*Q);
    nep=DEP([A0,A1],[0,1])
    return nep;
end

# A symmetric delay eigenvalue problem with double eigenvalues
# Examle from H. Voss and M. M. Betcke, Restarting iterative projection methods for Hermitian nonlinear eigenvalue problems with minmax property, Numer. Math., 2017
function dep_symm_double(n::Integer=100)
    LL=-sparse(1:n,1:n,2*ones(n))+sparse(2:n,1:n-1,ones(n-1),n,n)+sparse(1:n-1,2:n,ones(n-1),n,n)

    x = range(0, stop = pi, length = n)
    h=x[2]-x[1];
    LL=LL/(h^2)
    LL=kron(LL,LL)

    b=broadcast((x,y)->-100*abs(sin(x+y)),x,transpose(x))
    a=broadcast((x,y)->8*sin(x)*sin(y),x,transpose(x))
    B=sparse(1:n^2,1:n^2,b[:])
    A=LL+sparse(1:n^2,1:n^2,a[:])

    nep=DEP([A,B],[0,2.0])
    return nep
end

# A delay eigenvalue problem with a double non-semisimple eigenvalue in λ=3πi
# Examle from E. Jarlebring, Convergence factors of Newton methods for nonlinear eigenvalue problems, LAA, 2012
function dep_double()
    denom = 8+5*pi;
    a1 = 2/5 *(65*pi + 32)/(denom);
    a2 = 9*pi^2*(13+5*pi)/(denom);
    a3 = 324/5 *pi^2*(5*pi+4)/(denom);
    b1 = (260*pi + 128 + 225*pi^2)/(10*denom);
    b2 = 45*pi^2/denom;
    b3 = 81*pi^2*(40*pi + 32 + 25*pi^2)/(10*denom);
    A0 = [ 0    1    0;  0    0    1;  -a3  -a2  -a1];
    A1 = [ 0    0    0;  0    0    0;  -b3  -b2  -b1];

    tau=1.0;
    nep=DEP([A0,A1],[0,tau])
    return nep
end

#  A quadratic problem with real eigenvalues
# Four smallest eigenvalues of the problem:
# 1.0e+03 * [ -2.051741417993845  -0.182101627437811  -0.039344930222838  -0.004039879577113]
function real_quadratic()
    A0 = [ 4     0     1     1;
           0     2     1     1;
           1     1     6    -2;
           1     1    -2     3];
    A1 = [ 167  -140    95  -131;
          -140   327    54    85;
            95    54   235   -81;
          -131    85   -81   181];
    A2 = [ 2     1    -1    -1;
           1     5    -3     2;
          -1    -3     3     0;
          -1     2     0     3];
    A=[A0,A1,A2]
    nep=PEP(A)
    return nep
end

# A quadratic delay eigenvalue problem in "The infinite Bi-Lanczos method for nonlinear eigenvalue problems",  Sarah W. Gaaf and Elias Jarlebring
function qdep0()
    qdepbase=joinpath(dirname(@__FILE__()),
                      "converted_misc" ,"qdep_infbilanczos_")
    A0=read_sparse_matrix(qdepbase * "A0.txt", Int64)
    A1=read_sparse_matrix(qdepbase * "A1.txt", Int64)
    tau::Float64 = 1
    quadfun = S -> S^2
    constfun = S -> one(S)
    expfun = S -> exp(-tau*S)

    AA = [-one(A0), A0, A1]
    fi = [quadfun, constfun, expfun]
    return SPMF_NEP(AA, fi)
end

# A quadratic delay eigenvalue problem in "A linear eigenvalue algorithm for the  nonlinear eigenvalue problem",  Elias Jarlebring, Wim Michiels, Karl Meerbergen \\
function qdep1()
    A0=[ 0.3000   -0.6000         0    0.4000
        -0.3000    0.4000   -0.8000    1.9000
         0.1000   -1.6000   -1.3000         0
        -1.4000   -0.9000    0.2000    0.9000];
    A1=[ 0.8000    0.2000   -1.3000   -0.3000
        -1.1000    0.9000    1.2000    0.5000
         0.5000    0.2000   -1.6000   -1.3000
         0.7000    0.4000   -0.4000         0];
    return SPMF_NEP([one(A0), A0, A1], [λ -> -λ^2, λ -> one(λ), λ -> exp(-λ)])
end

# A quadratic eigenvalue problem with chosen eigenvalues
function qep_fixed_eig(n::Integer=5, E::AbstractVecOrMat=NaN*ones(2*n))
    if any(isnan.(E))
        msws_rng = MSWS_RNG()
        E=vec(gen_rng_mat(msws_rng,2*n,1))
    end
    A1 = diagm(0 => E[1:n])
    A2 = diagm(0 => E[n+1:2*n])
    K = one(A1)
    nep=PEP([A1*A2,-A1-A2,K])
    return nep
end

# This problem stems from
# L. P. Shayer and S. A. Campbell.  Stability, bifurcation and multistability in a system of two coupled neurons with multiple time delays. SIAM J. Applied Mathematics , 61(2):673–700, 2000
# It is also a benchmark example in DDE-BIFTOOL
function neuron0()
    pars = [1/2; -1; 1; 2.34; 0.2; 0.2 ; 1.5] .+ 0im
    kappa = pars[1]
    beta = pars[2]
    A = [0 pars[3]; pars[4] 0]

    x = [0; 0]  # The zero (trivial) stationary solution
    # A non-trivial stationary solution
    #x=[3.201081590416643561697725111745656884148241428177442574927999582405266342752249e-01
    #   5.096324796647208606096018689631125587762848405395086474417800152349531876959548e-01]

    tauv = [0;0.2;0.2;1.5]
    A0 = -kappa * Matrix(1.0I, 2, 2)
    A1 = A[2,1] * [0 0; (1-tanh(x[2])^2) 0]
    A2 = A[1,2] * [0 (1-tanh(x[1])^2); 0 0]
    A3 = beta * diagm(0 => [(1-tanh(x[1])^2), (1-tanh(x[2])^2)])
    return DEP([A0, A1, A2, A3], tauv)
end

# DEP modelling a beam
function beam(n::Integer=100)
    h=1/n;
    ee = ones(n);
    A0 = spdiagm(-1 => ee[1:n-1], 0 => -2*ee, 1 => ee[1:n-1]);
    A0[end,end]=1/h;
    A0[end,end-1]=-1/h;
    A1=sparse([n],[n],[1.0]); # A1=en*en'
    tau=1.0;
    return DEP([A0,A1],[0,tau]);
end

function sine_nep()
    data_dir=joinpath(dirname(@__FILE__()),  "converted_sine")
    A0=read_sparse_matrix(joinpath(data_dir,"sine_A0.txt"), Int64);
    A1=read_sparse_matrix(joinpath(data_dir,"sine_A1.txt"), Int64);
    A2=read_sparse_matrix(joinpath(data_dir,"sine_A2.txt"), Int64);
    V=Matrix(read_sparse_matrix(joinpath(data_dir,"sine_V.txt"), Int64));
    Q=Matrix(read_sparse_matrix(joinpath(data_dir,"sine_Q.txt"), Int64));

    n=size(A0,1);
    Z=spzeros(n,n);
    pep=PEP([A0,A1,Z,Z,A2]);
    # Matrix  sine function. Note that the Term is rank two which is not exploited here
    sin_nep=SPMF_NEP([V*Q'], [S-> sin(S)]);
    nep=SPMFSumNEP(pep,sin_nep) # Note: nep has a low-rank term
    return nep;
end

# From the tutorial on NEP-PACK online user's manual (except we use LowRankNEP)
function schrodinger_movebc(n=1000,L0=1,L1=8,α=25*pi/2,V0=10.0)

    # Discreatization matrices
    xv=Vector(range(0,stop=L0,length=n))
    h=xv[2]-xv[1];
    n=size(xv,1);
    V=x->1+sin(α*x);
    Dn=spdiagm(-1 => [ones(n-2);0]/h^2, 0 => -2*ones(n-1)/h^2, 1 => ones(n-1)/h^2)
    Vn=spdiagm(0 => [V.(xv[1:end-1]);0]);
    In=spdiagm(0 => [ones(n-1);0])
    # Construct functions
    f1=S->one(S);
    f2=S->-S;
    hh=S-> sqrt(S+V0*one(S))
    g=S-> cosh((L1-L0)*hh(S))
    f=S-> inv(hh(S))*sinh((L1-L0)*hh(S))

    # First NEP
    nep1=SPMF_NEP([Dn-Vn,In],[f1,f2]);

    # Create the Low-Rank NEP
    Lv1=sparse([zeros(n-1,1);1.0]);
    Lv2=sparse([zeros(n-1,1);1.0]);
    Uv1=sparse([zeros(n-1,1);1.0]);
    Uv2=sparse([zeros(n-3,1); [1/(2*h), -2/h, 3/(2*h)]])
    nep2=LowRankFactorizedNEP([Lv1,Lv2],[Uv1,Uv2],[g,f]);

    # Create a sum of the two
    nep=SumNEP(nep1,nep2);
    return nep
end


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

# asinh, acosh, and atanh are heavily based on FDLIBM code:
# e_sinh.c, e_sinhf, e_cosh.c, e_coshf, s_tanh.c, s_tanhf.c, s_asinh.c,
# s_asinhf.c, e_acosh.c, e_coshf.c, e_atanh.c, and e_atanhf.c
# that are made available under the following licence:

# ====================================================
# Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
#
# Developed at SunSoft, a Sun Microsystems, Inc. business.
# Permission to use, copy, modify, and distribute this
# software is freely granted, provided that this notice
# is preserved.
# ====================================================

@inline function exthorner(x, p::Tuple)
	# polynomial evaluation using compensated summation.
	# much more accurate, especially when lo can be combined with other rounding errors
    hi, lo = p[end], zero(x)
    for i in length(p)-1:-1:1
        pi = p[i]
        prod = hi*x
        err = fma(hi, x, -prod)
        hi = pi+prod
        lo = fma(lo, x, prod - (hi - pi) + err)
    end
    return hi, lo
end

# Hyperbolic functions
# sinh methods
H_SMALL_X(::Type{Float64}) = 2.0^-28
H_MEDIUM_X(::Type{Float64}) = 22.0

H_SMALL_X(::Type{Float32}) = 2f-12
H_MEDIUM_X(::Type{Float32}) = 9f0

H_LARGE_X(::Type{Float64}) = 709.7822265633563 # nextfloat(709.7822265633562)

H_LARGE_X(::Type{Float32}) = 88.72283f0

SINH_SMALL_X(::Type{Float64}) = 2.1
SINH_SMALL_X(::Type{Float32}) = 3.0f0

# For Float64, use DoubleFloat scheme for extra accuracy
function sinh_kernel(x::Float64)
    x2 = x*x
    x2lo = fma(x,x,-x2)
    hi_order = evalpoly(x2, (8.333333333336817e-3, 1.9841269840165435e-4,
                             2.7557319381151335e-6, 2.5052096530035283e-8,
                             1.6059550718903307e-10, 7.634842144412119e-13,
                             2.9696954760355812e-15))
    hi,lo = exthorner(x2, (1.0, 0.16666666666666635, hi_order))
    return muladd(x, hi, muladd(x, lo, x*x2lo*0.16666666666666635))
end
# For Float32, using Float64 is simpler, faster, and doesn't require FMA
function sinh_kernel(x::Float32)
    x=Float64(x)
    res = evalpoly(x*x, (1.0, 0.1666666779967941, 0.008333336726447933,
                         0.00019841001151414065, 2.7555538207080807e-6,
                         2.5143389765825282e-8, 1.6260094552031644e-10))
    return Float32(res*x)
end

function sinh(x::T) where T<:Union{Float32,Float64}
    # Method
    # mathematically sinh(x) is defined to be (exp(x)-exp(-x))/2
    #    1. Sometimes replace x by |x| (sinh(-x) = -sinh(x)).
    #    2. Find the branch and the expression to calculate and return it
    #      a)   0 <= x < SINH_SMALL_X
    #               approximate sinh(x) with a  minimax polynomial
    #      b)   SINH_SMALL_X <= x < H_LARGE_X
    #               return sinh(x) = (exp(x) - exp(-x))/2
    #      d)   H_LARGE_X  <= x
    #               return sinh(x) = exp(x/2)/2 * exp(x/2)
    #               Note that this branch automatically deals with Infs and NaNs

    absx = abs(x)
    if absx <= SINH_SMALL_X(T)
        return sinh_kernel(x)
    elseif absx >= H_LARGE_X(T)
        E = exp(T(.5)*absx)
        return copysign(T(.5)*E*E, x)
    end
    E = exp(absx)
    return copysign(T(.5)*(E - 1/E),x)
end

COSH_SMALL_X(::Type{T}) where T= one(T)

function cosh_kernel(x2::Float32)
    return evalpoly(x2, (1.0f0, 0.49999997f0, 0.041666888f0, 0.0013882756f0, 2.549933f-5))
end

function cosh_kernel(x2::Float64)
    return evalpoly(x2, (1.0, 0.5000000000000002, 0.04166666666666269,
                         1.3888888889206764e-3, 2.4801587176784207e-5,
                         2.7557345825742837e-7, 2.0873617441235094e-9,
                         1.1663435515945578e-11))
end

function cosh(x::T) where T<:Union{Float32,Float64}
    # Method
    # mathematically cosh(x) is defined to be (exp(x)+exp(-x))/2
    #    1. Replace x by |x| (cosh(x) = cosh(-x)).
    #    2. Find the branch and the expression to calculate and return it
    #      a)   x <= COSH_SMALL_X
    #               approximate sinh(x) with a minimax polynomial
    #      b)   COSH_SMALL_X <= x < H_LARGE_X
    #               return cosh(x) = = (exp(x) + exp(-x))/2
    #      e)   H_LARGE_X  <= x
    #               return cosh(x) = exp(x/2)/2 * exp(x/2)
    #      			Note that this branch automatically deals with Infs and NaNs

    absx = abs(x)
    if absx <= COSH_SMALL_X(T)
        return cosh_kernel(x*x)
    elseif absx >= H_LARGE_X(T)
        E = exp(T(.5)*absx)
        return T(.5)*E*E
    end
    E = exp(absx)
    return T(.5)*(E + 1/E)
end

# tanh methods
TANH_LARGE_X(::Type{Float64}) = 44.0
TANH_LARGE_X(::Type{Float32}) = 18.0f0
TANH_SMALL_X(::Type{Float64}) = 1.0
TANH_SMALL_X(::Type{Float32}) = 1.3862944f0       #2*log(2)
@inline function tanh_kernel(x::Float64)
    return evalpoly(x, (1.0, -0.33333333333332904, 0.13333333333267555,
                        -0.05396825393066753, 0.02186948742242217,
                        -0.008863215974794633, 0.003591910693118715,
                        -0.0014542587440487815, 0.0005825521659411748,
                        -0.00021647574085351332, 5.5752458452673005e-5))
end
@inline function tanh_kernel(x::Float32)
    return evalpoly(x, (1.0f0, -0.3333312f0, 0.13328037f0,
                        -0.05350336f0, 0.019975215f0, -0.0050525228f0))
end
function tanh(x::T) where T<:Union{Float32, Float64}
    # Method
    # mathematically tanh(x) is defined to be (exp(x)-exp(-x))/(exp(x)+exp(-x))
    #    1. reduce x to non-negative by tanh(-x) = -tanh(x).
    #    2. Find the branch and the expression to calculate and return it
    #      a) 0 <= x < H_SMALL_X
    #             Use a minimax polynomial over the range
    #      b) H_SMALL_X <= x < TANH_LARGE_X
    #           1 - 2/(exp(2x) + 1)
    #      c) TANH_LARGE_X <= x
    #            return 1
    abs2x = abs(2x)
    abs2x >= TANH_LARGE_X(T) && return copysign(one(T), x)
    abs2x <= TANH_SMALL_X(T) && return x*tanh_kernel(x*x)
    k = exp(abs2x)
    return copysign(1 - 2/(k+1), x)
end

# Inverse hyperbolic functions
AH_LN2(::Type{Float64}) = 6.93147180559945286227e-01
AH_LN2(::Type{Float32}) = 6.9314718246f-01
# asinh methods
function asinh(x::T) where T <: Union{Float32, Float64}
    # Method
    # mathematically asinh(x) = sign(x)*log(|x| + sqrt(x*x + 1))
    # is the principle value of the inverse hyperbolic sine
    # 1. Find the branch and the expression to calculate and return it
    #    a) |x| < 2^-28
    #        return x
    #    b) |x| < 2
    #        return sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
    #    c) 2 <= |x| < 2^28
    #        return sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1)))
    #    d) |x| >= 2^28
    #        return sign(x)*(log(x)+ln2))
    if isnan(x) || isinf(x)
        return x
    end
    absx = abs(x)
    if absx < T(2)
        # in a)
        if absx < T(2)^-28
            return x
        end
        # in b)
        t = x*x
        w = log1p(absx + t/(T(1) + sqrt(T(1) + t)))
    elseif absx < T(2)^28
        # in c)
        t = absx
        w = log(T(2)*t + T(1)/(sqrt(x*x + T(1)) + t))
    else
        # in d)
        w = log(absx) + AH_LN2(T)
    end
    return copysign(w, x)
end

# acosh methods
@noinline acosh_domain_error(x) = throw(DomainError(x, "acosh(x) is only defined for x ≥ 1."))
function acosh(x::T) where T <: Union{Float32, Float64}
    # Method
    # mathematically acosh(x) if defined to be log(x + sqrt(x*x-1))
    # 1. Find the branch and the expression to calculate and return it
    #     a) x = 1
    #         return log1p(t+sqrt(2.0*t+t*t)) where t=x-1.
    #     b) 1 < x < 2
    #         return log1p(t+sqrt(2.0*t+t*t)) where t=x-1.
    #     c) 2 <= x <
    #         return log(2x-1/(sqrt(x*x-1)+x))
    #     d) x >= 2^28
    #         return log(x)+ln2
    # Special cases:
    #     if x < 1 throw DomainError

    isnan(x) && return x

    if x < T(1)
        return acosh_domain_error(x)
    elseif x == T(1)
        # in a)
        return T(0)
    elseif x < T(2)
        # in b)
        t = x - T(1)
        return log1p(t + sqrt(T(2)*t + t*t))
    elseif x < T(2)^28
        # in c)
        t = x*x
        return log(T(2)*x - T(1)/(x+sqrt(t - T(1))))
    else
        # in d)
        return log(x) + AH_LN2(T)
    end
end

# atanh methods
@noinline atanh_domain_error(x) = throw(DomainError(x, "atanh(x) is only defined for |x| ≤ 1."))
function atanh(x::T) where T <: Union{Float32, Float64}
    # Method
    # 1.Reduced x to positive by atanh(-x) = -atanh(x)
    # 2. Find the branch and the expression to calculate and return it
    #     a) 0 <= x < 0.5
    #         return 0.5*log1p(2x/(1-x))
    #     b) 0.5 <= x <= 1
    #         return 0.5*log((x+1)/(1-x))
    # Special cases:
    #    if |x| > 1 throw DomainError
    isnan(x) && return x

    absx = abs(x)

    if absx > 1
        atanh_domain_error(x)
    end
    if absx < T(0.5)
        # in a)
        t = log1p(T(2)*absx/(T(1)-absx))
    else
        # in b)
        t = log((T(1)+absx)/(T(1)-absx))
    end
    return T(0.5)*copysign(t, x)
end

max_stars_repo_name: ranocha/julia

max_issues_repo_name: ranocha/julia

max_forks_repo_name: ranocha/julia


text:
@testset "no linesearch" begin

function f_nolin!(fvec, x)
    fvec[1] = (x[1]+3)*(x[2]^3-7)+18
    fvec[2] = sin(x[2]*exp(x[1])-1)
end

r = nlsolve(f_nolin!, [ -0.5; 1.4], autodiff = :forward, method = :newton, linesearch = LineSearches.Static())
@test converged(r)

end


text:
using Gtk, Gtk.GLib

@testset "glib" begin

ccall((:gtk_init, Gtk.libgtk), Nothing,(Ptr{Nothing},Ptr{Nothing}),C_NULL,C_NULL)

hnd = ccall((:gtk_window_new, Gtk.libgtk),Ptr{GObject},(Cint,),0)

wrap = convert(GObject, hnd)
#detect type and create wrapper (if neccessary) at runtime
@test isa(wrap, GtkWindowLeaf)

wrap2 = convert(GObject, hnd)
@test wrap == wrap2

repr = Base.print_to_string(wrap) #should display properties
@test endswith(repr,')')
@test occursin("name=\"\"",repr)
@test occursin("visible=FALSE",repr)
@test occursin("title=NULL",repr)
@test occursin("type=GTK_WINDOW_TOPLEVEL",repr)

@test isa(convert(Gtk.GLib.GBoxedUnkown, Gtk.GLib.GBoxedUnkown(C_NULL)), Gtk.GLib.GBoxedUnkown)

x = Ref{Int}(1)

function g_timeout_add_cb()
    x[] = 2
    false
end

g_idle_add(g_timeout_add_cb)
sleep(0.5)
@test x[] == 2

x[] = 1 #reset
g_timeout_add(g_timeout_add_cb, 1)
sleep(0.5)
@test x[] == 2

# do syntax

x[] = 1 #reset
g_idle_add() do
  x[] = 2
  return false # only call once
end
sleep(0.5)
@test x[] == 2

x[] = 1 #reset
g_timeout_add(1) do
  x[] = 2
  return false # only call once
end
sleep(0.5)
@test x[] == 2

# macro syntax

x[] = 1 #reset
@idle_add begin
  x[] = 2
end
sleep(0.5)
@test x[] == 2

# deprecated

function g_timeout_add_cb(user_data)
    user_data[] = 2
    false
end

x[] = 1 #reset
g_idle_add(()->g_timeout_add_cb(x))
sleep(0.5)
@test x[] == 2

x[] = 1 #reset
g_timeout_add(()->g_timeout_add_cb(x), 1)
sleep(0.5)
@test x[] == 2

end

# TODO
# module Test2   ### fails inside @testset
#     import Gtk
#     using Gtk.GLib
#     const suffix = :Test2
#     @Gtype GtkWidget Gtk.libgtk gtk_widget
# end
#
# @test Test2.GtkWidgetTest2 != Gtk.GtkWidgetLeaf
# @test Test2.GtkWidget == Gtk.GtkWidget
# nothing

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using ARDESPOT
using Test

using POMDPs
using POMDPModels
using POMDPSimulators
using Random
using POMDPModelTools

include("memorizing_rng.jl")

pomdp = BabyPOMDP()

K = 10
rng = MersenneTwister(14)
rs = MemorizingSource(K, 50)
Random.seed!(rs, 10)
b_0 = initialstate_distribution(pomdp)
scenarios = [i=>rand(rng, b_0) for i in 1:K]
b = ScenarioBelief(scenarios, rs, 0, false)
pol = FeedWhenCrying()
r1 = ARDESPOT.branching_sim(pomdp, pol, b, 10)
r2 = ARDESPOT.branching_sim(pomdp, pol, b, 10)
@test r1 == r2

scenarios = [1=>rand(rng, b_0)]
b = ScenarioBelief(scenarios, rs, 0, false)
pol = FeedWhenCrying()
r1 = ARDESPOT.rollout(pomdp, pol, b, 10)
r2 = ARDESPOT.rollout(pomdp, pol, b, 10)
@test r1 == r2

# constant bounds
bounds = (reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTSolver(bounds=bounds)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=2)
@time hist = simulate(hr, pomdp, planner)

# policy lower bound
bounds = IndependentBounds(DefaultPolicyLB(FeedWhenCrying()), 0.0)
solver = DESPOTSolver(bounds=bounds)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=2)
@time hist = simulate(hr, pomdp, planner)

# Type stability
pomdp = BabyPOMDP()
bounds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTSolver(epsilon_0=0.1,
                      bounds=bounds,
                      rng=MersenneTwister(4)
                     )
p = solve(solver, pomdp)

b0 = initialstate_distribution(pomdp)
D = @inferred ARDESPOT.build_despot(p, b0)
@inferred ARDESPOT.explore!(D, 1, p)
@inferred ARDESPOT.expand!(D, length(D.children), p)
@inferred ARDESPOT.prune!(D, 1, p)
@inferred ARDESPOT.find_blocker(D, length(D.children), p)
@inferred ARDESPOT.make_default!(D, length(D.children))
@inferred ARDESPOT.backup!(D, 1, p)
@inferred ARDESPOT.next_best(D, 1, p)
@inferred ARDESPOT.excess_uncertainty(D, 1, p)
@inferred action(p, b0)

bounds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
rng = MersenneTwister(4)
solver = DESPOTSolver(epsilon_0=0.1,
                      bounds=bounds,
                      rng=rng,
                      random_source=MemorizingSource(500, 90, rng),
                      tree_in_info=true
                     )
p = solve(solver, pomdp)
a = action(p, initialstate_distribution(pomdp))

include("random_2.jl")

# visualization
show(stdout, MIME("text/plain"), D)
a, info = action_info(p, initialstate_distribution(pomdp))
show(stdout, MIME("text/plain"), info[:tree])

# from README:
using POMDPs, POMDPModels, POMDPSimulators, ARDESPOT

pomdp = TigerPOMDP()

solver = DESPOTSolver(bounds=(-20.0, 0.0))
planner = solve(solver, pomdp)

for (s, a, o) in stepthrough(pomdp, planner, "sao", max_steps=10)
    println("State was $s,")
    println("action $a was taken,")
    println("and observation $o was received.\n")
end


text:
export NewtonRaphson

abstract type AbstractNonLinearSolver end

"""
    NewtonRaphson <: AbstractNonLinearSolver

Newton-Raphson algorithm. Used to solve the non-linear equation
``g(x, u) = 0``, at a fixed control ``u``.

### Attributes
- `maxiter::Int` (default 20): maximum number of iterations
- `tol::Float64` (default `1e-8`): tolerance of the algorithm
- `verbose::Int` (default `NONE`): verbosity level

"""
struct NewtonRaphson <: AbstractNonLinearSolver
    maxiter::Int
    tol::Float64
    verbose::Int
end
NewtonRaphson(; maxiter=20, tol=1e-8, verbose=VERBOSE_LEVEL_NONE) = NewtonRaphson(maxiter, tol, verbose)

"""
    ConvergenceStatus

Convergence status returned by a non-linear algorithm.

### Attributes
- `has_converged::Bool`: states whether the algorithm has converged.
- `n_iterations::Int`: total number of iterations of the non-linear algorithm.
- `norm_residuals::Float64`: final residual.
- `n_linear_solves::Int`: number of linear systems ``Ax = b`` resolved during the run.

"""
struct ConvergenceStatus
    has_converged::Bool
    n_iterations::Int
    norm_residuals::Float64
    n_linear_solves::Int
end

# Sparse utilities
mutable struct Spmat{VTI<:AbstractVector, VTF<:AbstractVector}
    colptr::VTI
    rowval::VTI
    nzval::VTF

    # create 2 Spmats from complex matrix
    function Spmat{VTI, VTF}(mat::SparseMatrixCSC{Complex{Float64}, Int}) where {VTI, VTF}
        matreal = new(VTI(mat.colptr), VTI(mat.rowval), VTF(real.(mat.nzval)))
        matimag = new(VTI(mat.colptr), VTI(mat.rowval), VTF(imag.(mat.nzval)))
        return matreal, matimag
    end
    # copy constructor
    function Spmat{VTI, VTF}(mat) where {VTI, VTF}
        return new(VTI(mat.colptr), VTI(mat.rowval), VTF(mat.nzval))
    end
end

function _copy_csc!(J_dest, J_src, shift)
    @inbounds for i in 1:size(J_src, 2)
        for j in J_src.colptr[i]:J_src.colptr[i+1]-1
            row = J_src.rowval[j]
            @inbounds J_dest[row+shift, i] = J_src.nzval[j]
        end
    end
end

function _transfer_sparse!(J_dest::SparseMatrixCSC, J_src::SparseMatrixCSC, shift)
    _copy_csc!(J_dest, J_src, shift)
end

function _transfer_sparse!(J_dest::CUSPARSE.CuSparseMatrixCSR, J_src::CUSPARSE.CuSparseMatrixCSR, shift)
    nnz_start = J_dest.rowPtr[shift+1]
    nnz_ = nnz(J_src)
    Jnnz = @view J_dest.nzVal[nnz_start:nnz_start+nnz_-1]
    Jnnz .= J_src.nzVal
end

# projection operator
function project!(w::VT, u::VT, u♭::VT, u♯::VT) where VT<:AbstractArray
    w .= max.(min.(u, u♯), u♭)
end

# Utils function to solve transposed linear system  A' x = y
# Source code taken from:
# https://github.com/JuliaGPU/CUDA.jl/blob/master/lib/cusolver/wrappers.jl#L78L111
function csclsvqr!(A::CUSPARSE.CuSparseMatrixCSC{Float64},
                    b::CUDA.CuVector{Float64},
                    x::CUDA.CuVector{Float64},
                    tol::Float64,
                    reorder::Cint,
                    inda::Char)
    n = size(A,1)
    desca = CUSPARSE.CuMatrixDescriptor(
        CUSPARSE.CUSPARSE_MATRIX_TYPE_GENERAL,
        CUSPARSE.CUSPARSE_FILL_MODE_LOWER,
        CUSPARSE.CUSPARSE_DIAG_TYPE_NON_UNIT, inda)
    singularity = Ref{Cint}(1)
    CUSOLVER.cusolverSpDcsrlsvqr(CUSOLVER.sparse_handle(), n, A.nnz, desca, A.nzVal, A.colPtr, A.rowVal, b, tol, reorder, x, singularity)

    if singularity[] != -1
        throw(SingularException(singularity[]))
    end

    x
end


text:
# low-level interface, meaning works with vectors and matrices, not data frames and estimation objects

"Structure to store test results"
struct BootTestResult{T}
  stat::T; stattype::String
  p::T; padj::T
  reps::Int64; repsfeas::Int64
  nbootclust::Int64
  dof::Int64; dof_r::T
  plot::Union{Nothing, NamedTuple{(:X, :p), Tuple{Tuple{Vararg{Vector{T}, N} where N},Vector{T}}}}
  peak::Union{Nothing, NamedTuple{(:X, :p), Tuple{Vector{T}, T}}}
  ci::Union{Nothing, Matrix{T}}
  dist::Matrix{T}
  b::Vector{T}
  V::Matrix{T}
  auxweights::Union{Nothing,Matrix{T}}
  # M::StrBootTest
end

"Return test statistic"
teststat(o::BootTestResult) = o.stat

"Return numerator of test statistic"
statnumer(o::BootTestResult) = o.b

"Return denominator of test statistic"
statvar(o::BootTestResult) = o.V

"""Return type of test statistic: "t", "z", "F", or "χ²" """
stattype(o::BootTestResult) = o.stattype

"Return p value"
p(o::BootTestResult) = o.p

"Return p value after multiple-hypothesis adjustment, if any"
padj(o::BootTestResult) = o.padj

"Return requested number of replications"
reps(o::BootTestResult) = o.reps

"Return actual number of replications, subject to enumeration of Rademacher draws"
repsfeas(o::BootTestResult) = o.repsfeas

"Return number of bootstrapping clusters in test"
nbootclust(o::BootTestResult) = o.nbootclust

"Return degrees of freedom of test"
dof(o::BootTestResult) = o.dof

"Return residual degrees of freedom of test"
dof_r(o::BootTestResult) = o.dof_r

"""
Return data for confidence plot of test.
Return value is a 2-tuple with named entries `X` and `p` holding
the confidence sampling locations and p values respectively. `X` is in turn
a 1- or 2-tuple of vectors of sampling coordinates for each
dimension of the tested hypothesis.
"""
plotpoints(o::BootTestResult) = o.plot

"Return parameter value with peak p value in test"
peak(o::BootTestResult) = o.peak

"Return confidence interval matrix from test, one row per disjoint piece"
ci(o::BootTestResult) = o.ci

"Return bootstrap distribution of statistic or statistic numerator in bootstrap test"
dist(o::BootTestResult) = o.dist

"Return auxilliary weight matrix for wild bootstrap"
auxweights(o::BootTestResult) = o.auxweights

strint(x) = iszero(mod(x,1)) ? "$(Int64(x))" : "$x"
function Base.show(io::IO, o::BootTestResult{T}) where T
	s = stattype(o) * ( iszero(dof_r(o)) ? isone(dof(o)) ? " " : "(" * strint(dof(o)) * ")" :                                                     # z, χ²
	                                       isone(dof(o)) ? "(" * strint(dof_r(o)) * ")" : "(" * strint(dof(o)) * ", " * strint(dof_r(o)) * ")" )  # t, F
	Printf.@printf(io, "%s = %6.4f\n", s, teststat(o))
	Printf.@printf(io, "p%s = %6.4f\n", repeat(" ", length(s)-1), p(o))
	isdefined(o, :ci) && !isnothing(o.ci) && length(o.ci)>0 && print(io, "CI" * repeat(" ", length(s)-2) * " = $(round.(ci(o); sigdigits=4))\n")
end

# single entry point with arguments already converted to standardized types, to allow a smaller set of precompile() calls(?)
function __wildboottest(
	R::Matrix{T},
	r::Vector{T};
	resp::VecOrMat{T},
	predexog::Matrix{T},
	predendog::Matrix{T},
	inst::Matrix{T},
	R1::Matrix{T},
	r1::Vector{T},
	clustid::Matrix{Int64},
	nbootclustvar::Int64,
	nerrclustvar::Int64,
	issorted::Bool,
	hetrobust::Bool,
	nfe::Int64,
	feid::VecOrMat{Int64},
	fedfadj::Int64,
	obswt::VecOrMat{T},
	fweights::Bool,
	maxmatsize::Float16,
	ptype::Symbol,
	bootstrapc::Bool,
	liml::Bool,
	fuller::T,
	kappa::T,
	arubin::Bool,
	small::Bool,
	clusteradj::Bool,
	clustermin::Bool,
	scorebs::Bool,
	reps::Int64,
	imposenull::Bool,
	auxwttype::Symbol,
	rng::AbstractRNG,
	level::T,
	rtol::T,
	madjtype::Symbol,
	nH0::Int16,
	ml::Bool,
	scores::Matrix{T},
	beta::Vector{T},
	A::Symmetric{T,Matrix{T}},
	gridmin::VecOrMat{T},
	gridmax::VecOrMat{T},
	gridpoints::VecOrMat{T},
	diststat::Symbol,
	getci::Bool,
	getplot::Bool,
	getauxweights::Bool) where T

	M = StrBootTest{T}(R, r, R1, r1, resp, predexog, predendog, inst, obswt, fweights, liml, fuller, kappa, arubin,
	                   reps, auxwttype, rng, maxmatsize, ptype, imposenull, scorebs, !bootstrapc, clustid, nbootclustvar, nerrclustvar, issorted, hetrobust, small, clusteradj, clustermin,
	                   nfe, feid, fedfadj, level, rtol, madjtype, nH0, ml, beta, A, scores, getplot,
	                   gridmin, gridmax, gridpoints)

	if getplot || (level<1 && getci)
		plot!(M)
		plot = getplot & isdefined(M, :plotX) ? (X=Tuple(M.plotX), p=M.plotY) : nothing
		peak = M.peak
		ci = level<1 & getci ? M.ci : nothing
	else
		ci = plot = peak = nothing
	end

	padj = getp(M)  # trigger main (re)computation

	BootTestResult{T}(getstat(M),
	                  isone(nrows(R)) ? (small ? "t" : "z") : (small ? "F" : "χ²"),
	                  M.p, padj, M.B, M.BFeas, M.N✻, M.dof, M.dof_r, plot, peak, ci,
	                  getdist(M, diststat),
	                  getb(M), getV(M),
	                  getauxweights && reps>0 ? getv(M) : nothing #=, M=#)
end

vecconvert(T::DataType, X) = Vector(isa(X, AbstractArray) ? vec(    eltype(X)==T ? X : T.(X)) : X)
matconvert(T::DataType, X) = Matrix(isa(X, AbstractArray) ? reshape(eltype(X)==T ? X : T.(X), size(X,1), size(X,2)) : X)

function _wildboottest(T::DataType,
					  R::AbstractVecOrMat,
						r::AbstractVecOrMat;
					  resp::AbstractVecOrMat{<:Real},
					  predexog::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  predendog::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  inst::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  R1::AbstractVecOrMat=zeros(T,0,0),
						r1::AbstractVecOrMat=zeros(T,0),
					  hetrobust::Bool=true,
					  clustid::AbstractVecOrMat{<:Integer}=zeros(Int,0,0),  # bootstrap-only clust vars, then boot&err clust vars, then err-only clust vars
					  nbootclustvar::Integer=ncols(clustid),
					  nerrclustvar::Integer=nbootclustvar,
						issorted::Bool=false,
						nfe::Integer=0,
					  feid::AbstractVecOrMat{<:Integer}=Int8[],
					  fedfadj::Integer=length(feid)>0 ? -1 : 0,
					  obswt::AbstractVecOrMat{<:Real}=T[],
					  fweights::Bool=false,
					  maxmatsize::Number=0,
					  ptype::Symbol=:symmetric,
					  bootstrapc::Bool=false,
					  liml::Bool=false,
					  fuller::Number=0,
					  kappa::Number=NaN,
					  arubin::Bool=false,
					  small::Bool=true,
						clusteradj::Bool=small,
						clustermin::Bool=false,
					  scorebs::Bool=false,
					  reps::Integer=999,
					  imposenull::Bool=true,
					  auxwttype::Symbol=:rademacher,
					  rng::AbstractRNG=MersenneTwister(),
					  level::Number=.95,
					  rtol::Number=1e-6,
					  madjtype::Symbol=:none,
					  nH0::Integer=1,
					  ml::Bool=false,
					  scores::AbstractVecOrMat=Matrix{Float32}(undef,0,0),
					  beta::AbstractVecOrMat=T[],
					  A::AbstractMatrix=zeros(T,0,0),
					  gridmin::Union{VecOrMat{S},VecOrMat{Union{S,Missing}}} where S<:Number = T[],
					  gridmax::Union{VecOrMat{S},VecOrMat{Union{S,Missing}}} where S<:Number = T[],
					  gridpoints::Union{VecOrMat{S},VecOrMat{Union{S,Missing}}} where S<:Number = Int64[],
					  diststat::Symbol=:none,
					  getci::Bool=true,
					  getplot::Bool=getci,
					  getauxweights::Bool=false)

	nrows(R)>2 && (getplot = getci = false)

	@assert any(auxwttype .== (:rademacher, :mammen, :webb, :gamma, :normal)) "auxwttype shoud be :rademacher, :mammen, :webb, :gamma, or :normal"
	@assert any(ptype .==(:symmetric, :equaltail, :lower, :upper)) "ptype should be :symmetric, :equaltail, :lower, or :upper"
	@assert any(madjtype .== (:none, :bonferroni, :sidak)) "madjtype should be :none, :bonferroni, or :sidak"
	@assert any(diststat .== (:none, :t, :numer))
	@assert ml || ncols(resp)==1 "resp should have one column"
  @assert (length(predexog)==0 || nrows(predexog)==nrows(resp)) &&
	        (length(predendog)==0 || nrows(predendog)==nrows(resp)) &&
					(length(inst)==0 || nrows(inst)==nrows(resp)) "All data vectors/matrices must have same height"
	@assert ncols(inst) >= ncols(predendog) "Model has fewer instruments than instrumented variables"
  @assert length(feid)==0 || nrows(feid)==nrows(resp) "feid vector must have same height as data matrices"
	@assert ncols(feid)≤1 "feid should have one column"
  @assert length(clustid)==0 || nrows(clustid)==nrows(resp) "clustid must have same height as data matrices"
  @assert nrows(obswt)==0 || nrows(obswt)==nrows(resp) "obswt must have same height as data matrices"
	@assert ncols(obswt)≤1 "obswt must have one column"
  @assert nrows(R)==nrows(r) "R and r must have same height"
  @assert (ncols(R) == (ml ? nrows(beta) : ncols(predexog)+ncols(predendog)) && isone(ncols(r))) "Wrong number of columns in null specification"
  @assert nrows(R1)==nrows(r1) "R₁ and r₁ must have same height"
  @assert length(R1)==0 || ncols(R1)==ncols(predexog)+ncols(predendog) "Wrong number of columns in model constraint specification"
	@assert ncols(r)==1 "r should have one column"
	@assert length(R1)==0 || ncols(r1)==1 "r1 should have one column"
  @assert nbootclustvar ≤ ncols(clustid) "nbootclustvar > width of clustid"
  @assert nerrclustvar ≤ ncols(clustid) "nerrclustvar > width of clustid"
  @assert reps ≥ 0 "reps < 0"
  @assert level ≥ 0. && level≤1. "level must be in the range [0,1]"
  @assert rtol > 0. "rtol ≤ 0"
  @assert nH0 > 0 "nH0 ≤ 0"
	@assert !liml || (ncols(predendog)>0 && ncols(inst)>0) "For liml, non-empty predendog and inst arguments are needed"
	@assert fuller==0 || (ncols(predendog)>0 && ncols(inst)>0) "For Fuller liml, non-empty predendog and inst arguments are needed"
	@assert iszero(ncols(predendog)) || ncols(inst)>0 "predendog provided without inst"
	@assert !arubin || ncols(predendog)>0 "Anderson-Rubin test requested but predendog not provided"

	if getplot || getci
		@assert iszero(length(gridmin   )) || length(gridmin   )==nrows(R) "Length of gridmin doesn't match number of hypotheses being jointly tested"
		@assert iszero(length(gridmax   )) || length(gridmax   )==nrows(R) "Length of gridmax doesn't match number of hypotheses being jointly tested"
		@assert iszero(length(gridpoints)) || length(gridpoints)==nrows(R) "Length of gridpoints doesn't match number of hypotheses being jointly tested"
		@assert iszero(length(gridmin   )) || ncols(gridmin)   ==1 "gridmin should have one column"
		@assert iszero(length(gridmax   )) || ncols(gridmax)   ==1 "gridmax should have one column"
		@assert iszero(length(gridpoints)) || ncols(gridpoints)==1 "gridpoints should have one column"
	end

	_gridmin = Vector{T}(undef, length(gridmin))
	_gridmax = Vector{T}(undef, length(gridmax))
	_gridpoints = Vector{T}(undef, length(gridpoints))
	for i ∈ 1:length(gridmin)  # cumbersome loops because map() and list comprehensions mess up type inference(?!)
		_gridmin[i] = T(ismissing(gridmin[i]) ? NaN : gridmin[i])
	end
	for i ∈ 1:length(gridmax)
		_gridmax[i] = T(ismissing(gridmax[i]) ? NaN : gridmax[i])
	end
	for i ∈ 1:length(gridpoints)
		_gridpoints[i] = T(ismissing(gridpoints[i]) ? NaN : gridpoints[i])
	end

	__wildboottest(
		matconvert(T,R),
		vecconvert(T,r);
		resp=vecconvert(T,resp),
		predexog=matconvert(T,predexog),
		predendog=matconvert(T,predendog),
		inst=matconvert(T,inst),
		R1=matconvert(T,R1),
		r1=vecconvert(T,r1),
		clustid=matconvert(Int64,clustid),
		nbootclustvar=Int64(nbootclustvar),
		nerrclustvar=Int64(nerrclustvar),
		issorted,
		hetrobust,
		nfe,
		feid=vecconvert(Int64,feid),
		fedfadj,
		obswt=vecconvert(T,obswt),
		fweights,
		maxmatsize=Float16(maxmatsize),
		ptype,
		bootstrapc,
		liml,
		fuller=T(fuller),
		kappa=T(kappa),
		arubin,
		small,
		clusteradj,
		clustermin,
		scorebs,
		reps=Int64(reps),
		imposenull,
		auxwttype,
		rng,
		level=T(level),
		rtol=T(rtol),
		madjtype,
		nH0=Int16(nH0),
		ml,
		scores=matconvert(T,scores),
		beta=vecconvert(T,beta),
		A=Symmetric(matconvert(T,A)),
		gridmin=_gridmin,
		gridmax=_gridmax,
		gridpoints=_gridpoints,
		diststat,
		getci,
		getplot,
		getauxweights)
end

_wildboottest(T::DataType, R, r::Number; kwargs...) = _wildboottest(T, R, [r]; kwargs...)
_wildboottest(T::DataType, R::UniformScaling{Bool}, r; kwargs...) = _wildboottest(T, diagm(fill(T(R.λ),nrows(r))), r; kwargs...)
_wildboottest(R::Union{UniformScaling{Bool},AbstractVecOrMat}, r::Union{Number,AbstractVecOrMat}; kwargs...) = _wildboottest(Float64, R, r; kwargs...)

"""
wildboottest([T::DataType=Float64,] R::AbstractMatrix, r::AbstractVector;
             resp, <optional keyword arguments>) -> WildBootTests.BootTestResult

Function to perform wild-bootstrap-based hypothesis test

# Positional arguments
* `T::DataType`: data type for inputs, results, and computations: Float32 or Float64 (default)
* `R::AbstractMatrix` and `r::AbstractVector`: required matrix and vector expressing the null Rβ=r; see notes below

# Required keyword argument
* `resp::AbstractVector`: response/dependent variable (y or y₁ in Roodman et al. (2019))

# Optional keyword arguments
* `predexog::AbstractVecOrMat`: exogenous predictors, including constant term, if any (X/X₁)
* `predendog::AbstractVecOrMat`: endogenous predictors (Y₂)
* `inst::AbstractVecOrMat`: instruments (X₂)
* `R1::AbstractMatrix` and `r1::AbstractVector`: model constraints; same format as for `R` and `r`
* `clustid::AbstractVecOrMat{<:Integer}`: data vector/matrix of error and bootstrapping cluster identifiers; see Notes
* `nbootclustvar::Integer=1`: number of bootstrap-clustering variables
* `nerrclustvar::Integer=nbootclustvar`: number of error-clustering variables
* `issorted:Bool=false`: time-saving flag: data matrices are already sort by column types 2, then 3, then 1 (see notes)
* `hetrobust::Bool=true`: true unless errors are treated as iid
* `nfe::Integer=0`: number of fixed-effect groups; if 0 yet `feid` is provided, will be computed
* `feid::AbstractVector{<:Integer}`: data vector for one-way fixed effect group identifier
* `fedfadj::Integer=nfe`: degrees of freedom that fixed effects (if any) consume
* `obswt::AbstractVector=[]`: observation weight vector; default is equal weighting
* `fweights::Bool=false`: true for frequency weights
* `maxmatsize::Number`: maximum size of auxilliary weight matrix (v), in gigabytes
* `ptype::Symbol=:symmetric`: p value type (`:symmetric`, `:equaltail`, `:lower`, `:upper`)
* `bootstrapc::Bool=false`: true to request bootstrap-c instead of bootstrap-t
* `liml::Bool=false`: true for LIML or Fuller LIML
* `fuller::Number`: Fuller LIML factor
* `kappa::Number`: fixed κ for _k_-class estimation
* `arubin::Bool=false`: true for Anderson-Rubin test
* `small::Bool=true`: true to multiply test statistics by G/(G-1) × N/(N-k), where G, N, k are number of clusters, observations, and predictors
* `clusteradj::Bool=true`: false to drop G/(G-1) factor
* `clustermin::Bool=false``: for multiway clustering, true to base G/(G-1) factor for all clusterings ]on the smallest G across clusterings
* `scorebs::Bool=false`: true for score bootstrap instead of wild bootstrap
* `reps::Integer=999`: number of bootstrap replications; `reps` = 0 requests classical Rao (or Wald) test if `imposenull` = `true` (or `false`)
* `imposenull::Bool=true`: true to impose null
* `auxwttype::Symbol=:rademacher`: auxilliary weight type (`:rademacher`, `:mammen`, `:webb`, `:normal`, `:gamma`)
* `rng::AbstractRNG=MersenneTwister()`: randon number generator
* `level::Number=.95`: significance level (0-1)
* `rtol::Number=1e-6`: tolerance for ci bound determination
* `madjtype::Symbol=:none`: multiple hypothesis adjustment (`none`, `:bonferroni`, `:sidak`)
* `nH0::Integer=1`: number of hypotheses tested, including one being tested now
* `ml::Bool=false`: true for (nonlinear) ML estimation
* `scores::AbstractVecOrMat`: for ML, pre-computed scores
* `beta::AbstractVector`: for ML, parameter estimates
* `A::AbstractMatrix`: for ML, covariance estimates
* `gridmin`: vector of graph lower bounds; max length 2, `missing`/`NaN` entries ask wildboottest() to choose
* `gridmax`: vector of graph upper bounds; `missing`/`NaN` entries ask wildboottest() to choose
* `gridpoints`: vector of number of sampling points; `missing`/`NaN` entries ask wildboottest() to choose
* `diststat::Symbole=:none`: `:t` to save bootstrap distribution of t/z/F/χ² statistics; `:numer` to save numerators thereof
* `getci::Bool=true`: whether to return confidence interval
* `getplot::Bool=getci`: whether to generate plot data
* `getauxweights::Bool=false`: whether to save auxilliary weight matrix (v)

# Notes
`T`, `ptype`, `auxwttype`, `madjtype`, and `diststat` may also be strings. Examples: `"Float32"` and `"webb"`.

The columns of `R` in the statement of the null should correspond to those of the matrix [`predexog` `predendog`],
where `predendog` is non-empty only in regressions with instruments.

Order the columns of `clustid` this way:
1. Variables only used to define bootstrapping clusters, as in the subcluster bootstrap.
2. Variables used to define both bootstrapping and error clusters.
3. Variables only used to define error clusters.
`nbootclustvar` is then the number of columns of type 1 or 2; `nerrclustvar` is the number of columns of type 2 or 3. Typically `clustid` is a single column of type 2.

`wildboottest()` does not handle missing data values: all data and identifier matrices must
be restricted to the estimation sample.

"""
wildboottest(   R, r; kwargs...) = _wildboottest(                                                  R, r; Dict(a.first => isa(a.second, AbstractString) ? Meta.parse(a.second) : a.second for a ∈ kwargs)...)
wildboottest(T, R, r; kwargs...) = _wildboottest(isa(T, AbstractString) ? eval(Meta.parse(T)) : T, R, r; Dict(a.first => isa(a.second, AbstractString) ? Meta.parse(a.second) : a.second for a ∈ kwargs)...)


text:
export at, atend, atrow, setat!, fst, snd, third, last
export part, values, ckeys, cvalues, vec, rowpart, trimmedpart, take, takelast, takewhile
export drop, dropat, droplast, dropwhile, cut
export every
export partition, partsoflen
export getindex
export extract, extractvec, extractnested
export fieldvalues, dict
export isnil
export czip

#######################################
##  at, atend

@inline at{T,N}(a::NTuple{T,N},i) = a[i]
@inline at(a::AbstractArray, ind::Tuple) = a[ind...]
@inline at{T}(a::AbstractArray{T},i::AbstractArray) =
    len(i) == 1 ? (size(i,1) == 1 ? at(a, i[1]) : a[subtoind(i,a)]) : error("index has len>1")
@inline at{T}(a::AbstractArray{T,1},i::Number) = a[i]
#at{T,N}(a::AbstractArray{T,N},i) = slicedim(a,N,i)
@inline at{T}(a::AbstractArray{T,2},i::Number) = col(a[:,i])
@inline at{T}(a::AbstractArray{T,3},i::Number) = a[:,:,i]
@inline at{T}(a::AbstractArray{T,4},i::Number) = a[:,:,:,i]
@inline at{T}(a::AbstractArray{T,5},i::Number) = a[:,:,:,:,i]
@inline at{T}(a::AbstractArray{T,6},i::Number) = a[:,:,:,:,:,i]
@inline at{T}(a::AbstractArray{T,7},i::Number) = a[:,:,:,:,:,:,i]
@inline at{T}(a::AbstractArray{T,8},i::Number) = a[:,:,:,:,:,:,:,i]
@inline at{T,N}(a::AbstractArray{T,N},i::Number,args...) = at(at(a,i), args...)
@inline at(a::Dict,i) = a[i]
@inline at(a::Dict,ind...) = at(at(a, ind[1]), ind[2:end]...)
@inline at(a,i) = a[i]
atend(a, i) = at(a, len(a)-i+1)
atrow(a,i) = a[i,:]

#######################################
##  setat!

@inline setat!{T}(a::AbstractArray{T,1},i::Number,v) = (a[i] = v; a)
@inline setat!{T}(a::AbstractArray{T,2},i::Number,v) = (a[:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,3},i::Number,v) = (a[:,:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,4},i::Number,v) = (a[:,:,:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,5},i::Number,v) = (a[:,:,:,:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,6},i::Number,v) = (a[:,:,:,:,:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,7},i::Number,v) = (a[:,:,:,:,:,:,i] = v; a)
@inline setat!{T}(a::AbstractArray{T,8},i::Number,v) = (a[:,:,:,:,:,:,:,i] = v; a)
@inline setat!(a,i,v) = (a[i] = v; a)

@inline fst(a) = at(a,1)
@inline snd(a) = at(a,2)
@inline third(a) = at(a,3)

import Base.last
@inline last(a::Union{AbstractArray,AbstractString}) = at(a,len(a))
@inline last(a::Union{AbstractArray,AbstractString}, n) = trimmedpart(a,(-n+1:0)+len(a))

#######################################
##  part

part(a::AbstractArray, i::Real) = part(a,[i])
part{T}(a::Vector, i::AbstractArray{T,1}) = a[i]
part{T}(a::AbstractString, i::AbstractArray{T,1}) = string(a[i])
part{T}(a::NTuple{T},i::Int) = a[i]
part{T,T2,N}(a::AbstractArray{T2,N}, i::AbstractArray{T,1}) = slicedim(a,max(2,ndims(a)),i)
part{T1,T2}(a::AbstractArray{T1,1}, i::AbstractArray{T2,1}) = a[i]
dictpart(a, inds) = Dict(map(filter(collect(keys(a)), x->in(x,inds)), x->Pair(x,at(a,x))))
part(a::Dict, inds::AbstractVector) = dictpart(a, inds)
part(a::Dict, inds...) = dictpart(a, inds)
part{T<:Real}(a::AbstractArray,i::DenseArray{T,2}) = map(i, x->at(a,x))
part(a::AbstractArray,i::Base.ValueIterator) = part(a,typed(collect(i)))

import Base.values
values(a::Dict, ind, inds...) = values(a, [ind; inds...])
values(a::Dict, inds::AbstractArray) = mapvec(inds,x->at(a,x))
ckeys(a::Dict) = collect(keys(a))
cvalues(a::Dict) = collect(values(a))
values(a, inds...) = Any[getfield(a,x) for x in inds]

import Base.vec
vec(a::Dict) = [Pair(k,a[k]) for k in keys(a)]

rowpart(a::Matrix, i) = a[i, :]

trimmedpart(a, i::Int) = trimmedpart(a, [i])
trimmedpart(a, i::UnitRange) = part(a, max(1, minimum(i)):min(len(a),maximum(i)))
trimmedpart(a, i::AbstractArray) = part(a, i[(i .>= 1) & (i .<= len(a))])

import Base.take
take(a::Union{Array, UnitRange, AbstractString}, n::Int) = part(a, 1:min(n, len(a)))
takelast(a, n::Int = 1) = part(a, max(1,len(a)-n+1):len(a))
function takewhile(a, f)
    for i in 1:len(a)
        if !f(at(a,i))
            i == 1 && return []
            return take(a,i-1)
        end
    end
    a
end

import Base.drop
drop(a::AbstractString,i::Int) = part(a,i+1:len(a))
drop(a::AbstractArray,i::Int) = part(a,i+1:len(a))
dropat(a, ind) = part(a, setdiff(1:len(a), ind))

droplast(a) = isempty(a) ? a : part(a,1:(len(a)-1))
droplast(a,i) = isempty(a) ? a : part(a,1:(len(a)-i))
function dropwhile(a, f)
    for i in 1:len(a)
        if !f(at(a,i))
            return drop(a,i-1)
        end
    end
    []
end

cut(a, i) = (dropat(a,i),part(a,i))

every(a,n) = part(a,1:n:len(a))

function partition(a,n)
    n = min(n, len(a))
    ind = round(Int, linspace(1, len(a)+1, n+1))
    r = Array{Any}(n)
    for i = 1:n
        r[i] = part(a, ind[i]:ind[i+1]-1)
    end
    r
end

function partsoflen(a,n::Int)
    s = len(a)
    [part(a, i:floor(min(s,i+n-1))) for i in 1:n:s]
end

extract(a::Array, x::Any, default = nothing) = map(a, y->extract(y, x, default))
extract(a::Array, x::Symbol, default = nothing) = map(a, y->extract(y, x, default))
extractvec(a::Array, x::Any, default = nothing) = mapvec(a, y->extract(y, x, default))
extractvec(a::Array, x::Symbol, default = nothing) = mapvec(a, y->extract(y, x, default))
extract(a::Dict, x::Symbol, default = nothing) = get(a, x, default)
extract(a::Dict, x, default = nothing) = get(a, x, default)
extract(a, x::Symbol, default = nothing) = getfield(a,x)
extractnested(a::Array, args...) = map(a, x->at(x,args...))

fieldvalues(a) = [getfield(a,x) for x in sort(fieldnames(a))]
dict(a) = Dict([Pair(k,getfield(a,k)) for k in sort(fieldnames(a))])

isnil(a) = a == nothing || a == Void || (isa(a, Nullable) && isnull(a))

czip(a) = czip(a...)
function czip(a,b...)
    temp = Any[a, b...]
    mapvec(1:len(a), i->mapvec(temp, x->at(x,i)))
end


text:
@testset "schemas" begin

    using StatsModels: schema, apply_schema, FullRank

    f = @formula(y ~ 1 + a + b + c + b&c)
    df = (y = rand(9), a = 1:9, b = rand(9), c = repeat(["d","e","f"], 3))
    f = apply_schema(f, schema(f, df))
    @test f == apply_schema(f, schema(f, df))

end


text:
"""
ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля

РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке, следующий - за исключением двух последних клеток на Востоке, и т.д.
"""

function liners!(r::Robot)
    first = go_to_point!(r,West) #Робот двигается до западной границы
    second = go_to_point!(r,Sud) #Робот двигается южной границы

    putmarker!(r)
    side = Ost
    putmarkers!(r,side)
    side = inverse(side)

    back_by_markers!(r,side)
    back_to_point!(r)

    go_to_point!(r,Nord,second)
    go_to_point!(r,Ost,first) #Возвращение в исходную точку

end

function go_to_point!(r::Robot,side::HorizonSide) #Считается кол-во шагов
    steps=0
    while isborder(r,side) == false
        move!(r,side)
        steps=steps+1
    end
    return steps #Возвращает кол-во шагов
end

function go_to_point!(r::Robot,side::HorizonSide,steps::Int) #Возвращается кол-во шагов
    for a in (1:steps)
        move!(r,side)
    end
end

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

function back_to_point!(r::Robot) #Возвращение в левый нижний угол
    while isborder(r,West)==false
        move!(r,West)
    end
    while isborder(r,Sud)==false
        move!(r,Sud)
    end
end

function back_by_markers!(r::Robot,side::HorizonSide) #Робот возвращается по маркерам и маркирует следующую строку
    while isborder(r,Nord) == false
        move!(r,Nord)
        putmarkers!(r,side)
        side = inverse(side)
        while ismarker(r) == true
            move!(r,side)
        end
        side = inverse(side)
        if isborder(r,side) == false
          move!(r,side)
        end
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))  #Функция задаёт направление, противоположное данному

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
import Dates
println(Dates.now())

using Distributed
include("imports.jl")  # Load on main node.

# NOTE: Change this.
if length(ARGS) == 4
  SIMDIR = ARGS[1]
  RESULTS_DIR_PREFIX = ARGS[2]
  AWS_BUCKET_PREFIX = ARGS[3]
  ISTEST = parse(Int, ARGS[4]) > 0  # this argument should be 0 or 1 (for test)
else
  println("Usage: julia parsim.jl <simdir> <results_dir_prefix> <aws_bucket_prefix> <istest>")
  exit(1)
end

function setnumcores(n::Int)
  children_procs = filter(w -> w > 1, workers())
  rmprocs(children_procs)
  addprocs(n)
end

# NOTE: read configs
include("$(SIMDIR)/settings.jl")

for setting in settings
  setting[:results_dir] = "$(RESULTS_DIR_PREFIX)/$(setting[:outdir_suffix])"
  setting[:aws_bucket] = "$(AWS_BUCKET_PREFIX)/$(setting[:outdir_suffix])"
  if ISTEST
    setting[:nsamps] = 10
    setting[:nburn] = 10
  end
end

ncores = 24
setnumcores(ncores)
@everywhere include("imports.jl")  # Load on main node.

println("Sourcing files ..."); flush(stdout)
@everywhere SIMDIR = $SIMDIR
@everywhere include("$(SIMDIR)/simfn.jl")

@everywhere function sim(setting)
  results_dir = setting[:results_dir]
  aws_bucket = setting[:aws_bucket]
  mkpath(results_dir)

  println("Running $(results_dir)")
  flush(stdout)

  # FIXME: Can't do `asyncmap` and redirect outputs.
  # CytofRepFAM.Model.redirect_stdout_to_file("$(results_dir)/log.txt") do
  #   simfn(setting)
  # end

  simfn(setting)

  # Send to S3.
  CytofRepFAM.Model.s3sync(from=results_dir, to=aws_bucket,
                           tags=`--exclude '*.nfs'`)

  # Remove results to save space.
  # rm(results_dir, recursive=true)
  return
end

# NOTE:
# - f::Function: A function which takes one argument of type Dict{Any}
# - settings::Vector{Dict}): A vector of settings
println("Starting jobs ..."); flush(stdout)
ntasks = 6  # 6 runs at a time, because each run does 4-core PT
@time _ = asyncmap(sim, settings, ntasks=ntasks);
println("DONE with all runs!")
println(Dates.now())


text:
"""Initializes a SCM calibration process."""

using ArgParse
using CalibrateEDMF
using CalibrateEDMF.Pipeline

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "Inverse problem config file"
    arg_type = String
    "--job_id"
    help = "Job identifier"
    arg_type = String
    default = "12345"
end
parsed_args = parse_args(ARGS, s)
include(parsed_args["config"])
config = get_config()

init_calibration(config; job_id = parsed_args["job_id"], config_path = parsed_args["config"])


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

JLLWrappers.@generate_main_file_header("nauty")
JLLWrappers.@generate_main_file("nauty", UUID("55c6dc9b-343a-50ca-8ff2-b71adb3733d5"))
end  # module nauty_jll


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 43179a7a-0dc3-4987-914b-5f4bbdb0c672
begin
	using ProgressMeter
	using PyCall
	np = pyimport("numpy")
	cel = pyimport("celluloid")
	using PyPlot
end

# ╔═╡ ea404548-ac63-49ba-bb15-e1b4b564c153
md"
# Random Walk "

# ╔═╡ 8902895a-f524-4152-93c6-1c175aeefbd0
md"## Unbiaised random walk"

# ╔═╡ 4df7f63d-8f57-4469-bc04-5b79b3c47ae6
function unbiaisedrw(A,B,vit,h,l)
	v = vit
    r = rand(Float64)*360

	A = A + v*cosd(r)
	B = B + v*sind(r)

	#boundaries conditions
	if A < 0
		A = 0
	end
	if A > l - abs(v*cosd(r))
		A = l
	end
	if B < 0
		B = 0
	end
	if B > h - abs(v*sind(r))
		B = h
	end
    return [A,B]
end;

# ╔═╡ 2c6e5885-6ad6-412b-8dcb-d98314109d5a
md" ## Biaised random walk"

# ╔═╡ a29b8a42-fe93-11eb-0544-fdfb80f03003
function biaisedrw(A,B,vit,h,l)
	v = vit
    r = rand(Float64)*360
    b = rand(Float64)

	A = A + v*cosd(r)
	if b > 0.98
		B = B + v*abs(sind(r))
	else
		B = B + v*sind(r)
	end

	# boundaries conditions
	if A < 0
		A = 0
	end
	if A > l - abs(v*cosd(r))
		A = l
	end
	if B < 0
		B = 0
	end
	if B > h - abs(v*sind(r))
		B = h
	end
    return [A,B]
end;

# ╔═╡ e5f498de-785e-45c6-a3d0-10f9d9369b5a
begin
	nbrsim = 8000;
	nbrcell = 20;
	haut = 150.0;
	long = 150.0;
	v = 1;

	list1 = [rand(0:0.001:long/3,nbrcell),rand(0:0.001:haut/3,nbrcell)]
	list2 = [rand(0:0.001:long/3,nbrcell),rand(haut/3:0.001:2*haut/3,nbrcell)]
	list3 = [rand(0:0.001:long/3,nbrcell),rand(2*haut/3:0.001:haut,nbrcell)]
	list4 = [rand(long/3:0.001:2*long/3,nbrcell),rand(0:0.001:haut/3,nbrcell)]
	list5 = [rand(long/3:0.001:2*long/3,nbrcell),rand(haut/3:0.001:2*haut/3,nbrcell)]
	list6 = [rand(long/3:0.001:2*long/3,nbrcell),rand(2*haut/3:0.001:haut,nbrcell)]
	list7 = [rand(2*long/3:0.001:long,nbrcell),rand(0:0.001:haut/3,nbrcell)]
	list8 = [rand(2*long/3:0.001:long,nbrcell),rand(haut/3:0.001:2*haut/3,nbrcell)]
	list9 = [rand(2*long/3:0.001:long,nbrcell),rand(2*haut/3:0.001:haut,nbrcell)]

	list0 = [[vcat(list1[1],list2[1],list3[1],list4[1],list5[1],list6[1],list7[1],list8[1],list9[1]),vcat(list1[2],list2[2],list3[2],list4[2],list5[2],list6[2],list7[2],list8[2],list9[2])]]

	list = list0
	for i = 1:nbrsim
		newlist = [[],[]]
		for i = 1:9*nbrcell
			X = list[end][1][i]
			Y = list[end][2][i]

			metro = biaisedrw(X,Y,v,haut,long) # replace by unbiaisedrw for unbiaised 												 random walk
			append!(newlist[1],[metro[1]])
			append!(newlist[2],[metro[2]])
		end
		append!(list,[newlist])
	end
end

# ╔═╡ dcedc848-5219-4614-9781-378defe5496a
md" ## Plot"

# ╔═╡ 13337591-67b3-4cc9-a086-ddcc6be1ed4a
begin
	a = 0.7
	f = 8000
	s = 100
	plt.figure(figsize=(6,6))
	plt.scatter(list[f][1][1:nbrcell],list[f][2][1:nbrcell],s=s,c="#4477AA",alpha = 0.5)
	plt.scatter(list[f][1][nbrcell+1:2*nbrcell],
		list[f][2][nbrcell+1:2*nbrcell],s=s,c="#DEA73A",alpha = a)
	plt.scatter(list[f][1][2*nbrcell+1:3*nbrcell],
		list[f][2][2*nbrcell+1:3*nbrcell],s=s,c="#44AA77",alpha = a)
	plt.scatter(list[f][1][3*nbrcell+1:4*nbrcell],
		list[f][2][3*nbrcell+1:4*nbrcell],s=s,c="#AAAA44",alpha = a)
	plt.scatter(list[f][1][4*nbrcell+1:5*nbrcell],
		list[f][2][4*nbrcell+1:5*nbrcell],s=s,c="gray",alpha = a)
	plt.scatter(list[f][1][5*nbrcell+1:6*nbrcell],
		list[f][2][5*nbrcell+1:6*nbrcell],s=s,c="#DD77AA",alpha = a)
	plt.scatter(list[f][1][6*nbrcell+1:7*nbrcell],
		list[f][2][6*nbrcell+1:7*nbrcell],s=s,c="#781C81",alpha = a)
	plt.scatter(list[f][1][7*nbrcell+1:8*nbrcell],
		list[f][2][7*nbrcell+1:8*nbrcell],s=s,c="#404096",alpha = a)
	plt.scatter(list[f][1][8*nbrcell+1:9*nbrcell],
		list[f][2][8*nbrcell+1:9*nbrcell],s=s,c="#D92120",alpha = a)

	plt.tick_params(left=false,
                bottom=false,
                labelleft=false,
                labelbottom=false)

	plt.gca().set_aspect("equal")
	plt.draw()
	plt.show()
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"

[compat]
ProgressMeter = "~1.7.1"
PyCall = "~1.92.3"
PyPlot = "~2.9.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[PyPlot]]
deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
git-tree-sha1 = "67dde2482fe1a72ef62ed93f8c239f947638e5a2"
uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
version = "2.9.0"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
"""

# ╔═╡ Cell order:
# ╟─ea404548-ac63-49ba-bb15-e1b4b564c153
# ╠═43179a7a-0dc3-4987-914b-5f4bbdb0c672
# ╟─8902895a-f524-4152-93c6-1c175aeefbd0
# ╠═4df7f63d-8f57-4469-bc04-5b79b3c47ae6
# ╟─2c6e5885-6ad6-412b-8dcb-d98314109d5a
# ╠═a29b8a42-fe93-11eb-0544-fdfb80f03003
# ╠═e5f498de-785e-45c6-a3d0-10f9d9369b5a
# ╟─dcedc848-5219-4614-9781-378defe5496a
# ╠═13337591-67b3-4cc9-a086-ddcc6be1ed4a
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: Flamanjaune/PHD-Thesis

max_issues_repo_name: Flamanjaune/PHD-Thesis

max_forks_repo_name: Flamanjaune/PHD-Thesis


text:
@marginalrule typeof(+)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    return (in1 = prod(ProdAnalytical(), NormalMeanVariance(mean(m_out) - mean(m_in2), var(m_out)), m_in1), in2 = m_in2)
end

@marginalrule typeof(+)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_in1, W_in1 = weightedmean_precision(m_in1)
    xi_in2, W_in2 = weightedmean_precision(m_in2)

    return MvNormalWeightedMeanPrecision([ xi_in1 + xi_out; xi_in2 + xi_out ], [ W_in1 + W_out W_out; W_out W_in2 + W_out])
end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rasterpath, zipurl, zipname

@testset "SRTM" begin

    zip_url = URI(scheme = "https", host = "srtm.csi.cgiar.org", path = "/wp-content/uploads/files/srtm_5x5/TIFF/srtm_01_02.zip")
    tile_index1 = CartesianIndex(2, 1)   # [y, x] order
    tile_index2 = CartesianIndex(2, 2)   # [y, x] order
    @test zipurl(SRTM; tile_index = tile_index1) == zip_url
    @test zipurl(SRTM; tile_index = tile_index1) == zip_url
    @test zipname(SRTM; tile_index = tile_index1) == "srtm_01_02.zip"

    raster_path1 = joinpath(ENV["RASTERDATASOURCES_PATH"], "SRTM", "srtm_01_02.tif")
    raster_path2 = joinpath(ENV["RASTERDATASOURCES_PATH"], "SRTM", "srtm_02_02.tif")
    @test rasterpath(SRTM; tile_index = tile_index1) == raster_path1
    @test getraster(SRTM; tile_index = tile_index1) == raster_path1
    @test isfile(raster_path1)

    lon1, lat1 = -180, 55       # Coordinates of [0, 0] pixel of tile x=1, y=2
    @test RasterDataSources._wgs84_to_tile_index(lon1, lat1) == tile_index1
    @test getraster(SRTM; bounds = (lon1, lat1, lon1, lat1)) == reshape([raster_path1], 1, 1)
    lon2, lat2 = -172.5, 52.5   # Coordinates of [3000, 3000] pixel of tile x=2, y=2
    @test RasterDataSources._wgs84_to_tile_index(lon2, lat2) == tile_index2
    @test getraster(SRTM; bounds = (lon1, lat1, lon2, lat2)) == reshape([raster_path1, raster_path2], 1, 2)
end


text:
using Rita
using Test
using QuickCheck
using LinearAlgebra: norm

import QuickCheck.generator
generator(::Type{Rita.Dir}, size) = Rita.Dir(rand(0:2*pi), rand(-pi:pi))

# Some simple 2d snellius in order to test some simple cases
function snellius_2d(in_angle::Float64, n1::Float64, n2::Float64)::Float64
    asin(sin(in_angle)*n1/n2)
end

@test Rita.refract_dir(Rita.Dir([0,-1, 0]), Rita.Dir([0, 0, -1]), 1., 1.) ≈ Rita.Dir([0,-1, 0])
@test Rita.refract_dir(Rita.Dir([0, 0, 1]), Rita.Dir([0, 0, -1]), 1., 1.) ≈ Rita.Dir([0, 0, 1])
# no total reflection
condproperty((dir::Rita.Dir, n::Rita.Dir, n1::Float64, n2::Float64) -> norm(Rita.refract_dir(dir, n, n1, n2)) ≈ 1.,
             100, 1000, (dir::Rita.Dir, n::Rita.Dir, n1::Float64, n2::Float64) -> 1.0 <= n1 <= n2 <= 5.0)


text:
### ℝ^{D-1} to D-dim simplex.

# Takes x ∈ (0,1)^D to somewhere within the interior of the D-dim simplex,
#   then we map it to the boundary of the (D+1)-dim simplex.
function unitcubetosimplexinterior(x::Vector{T})::Vector{T} where T
    D = length(x)

    s = Vector{T}(undef,D)
    for i = 1:D-1
        s[i] = x[i]*prod( one(T)-x[j] for j = i+1:D)
    end
    s[D] = x[D]

    return s
end

function unitcubetosimplex(x::Vector{T})::Vector{T} where T
    D = length(x)

    s = unitcubetosimplexinterior(x)
    push!(s, one(T) - sum(s[i] for i = 1:D))

    return s
end

function Euclideantounitcube(r::Vector{T})::Vector{T} where T
    D = length(r)

    x = Vector{T}(undef,D)
    for i = 1:D
        x[i] = realtounitinterval(r[i])
    end

    return x
end

function realtounitintervalV(x::T)::T where T
    return (one(T) + x/sqrt(one(T)+x^2))/2-0.5
end

function realtounitinterval2(x::T, a::T)::T where T
    return (one(T) + a*x/sqrt(one(T)+(a*x)^2))/2
end

function realtounitinterval(x::T)::T where T
    return (one(T) + x/sqrt(one(T)+x^2))/2
end
# function realtounitinterval(x::T)::T where T
#     return logicsticfunc(x)
# end

function derivativeofrealtounitinterval(x::T)::T where T
    return one(T)/(2*sqrt(one(T)+x^2)^3)
end

# package up.
function Euclideantosimplex(r::Vector{T})::Vector{T} where T
    x = Euclideantounitcube(r) # to the unit cube.
    s = unitcubetosimplex(x)

    return s
end

function Euclideantosimplexinterior(r::Vector{T})::Vector{T} where T
    x = Euclideantounitcube(r) # to the unit cube.
    s = unitcubetosimplexinterior(x)

    return s
end

### other direction: D-dim simplex to ℝ^{D-1}.
function simplextounitcube(s::Vector{T})::Vector{T} where T
    D = length(s) -1

    x = Vector{T}(undef,D)
    for i = 1:D-1
        x[i] = s[i]/(one(T) - sum( s[j] for j = i+1:D))
    end
    x[D] = s[D]

    return x
end

function simplexinteriortounitcube(s::Vector{T})::Vector{T} where T
    D = length(s)

    x = Vector{T}(undef,D)
    for i = 1:D-1
        x[i] = s[i]/(one(T) - sum( s[j] for j = i+1:D))
    end
    x[D] = s[D]

    return x
end

function unitcubetoEuclidean(x::Vector{T})::Vector{T} where T
    D = length(x)

    r = Vector{T}(undef,D)
    for i = 1:D
        r[i] = unitintervaltoreal(x[i])
    end

    return r
end

# inverse of realtounitinterval()
function unitintervaltoreal(y::T)::T where T
    u = 2*y-1
    return u/sqrt(one(T)-u^2)
end
# # test code.
# x = randn()
# h = realtounitinterval(x)
# x_rec = unitintervaltoreal(h)
# println("x = ", x)
# println("x_rec = ", x_rec)
# println("discrepancy = ", abs(x-x_rec))
# function unitintervaltoreal(y::T)::T where T
#     return invlogicsticfunc(y)
# end

# package up.
function simplextoEuclidean(s::Vector{T})::Vector{T} where T
    x = simplextounitcube(s) # to the unit cube.
    r = unitcubetoEuclidean(x)

    return r
end

function simplexinteriortoEuclidean(s::Vector{T})::Vector{T} where T
    x = simplexinteriortounitcube(s) # to the unit cube.
    r = unitcubetoEuclidean(x)

    return r
end

###### Jacobian of Euclideantosimplexinterior()

# derivative of realtounitinterval().
function ∂realtounitinterval(x::T)::T where T <: Real
    tmp = 2*sqrt(1+x^2)*(1+x^2)
    return one(T)/tmp
end

function ln∂realtounitinterval(x::T)::T where T <: Real
    tmp = log(2) + 3/2*log((1+x^2))
    return -tmp
end

# derivative of simplexinteriortounitcube().
function ∂simplexinteriortounitcube(u::Vector{T}, a::Int, b::Int) where T <: Real
    D = length(u)
    if a != b
        return -u[a]*prod( one(T)-u[j] for j = a+1:D if j != b)
    end

    if a == D
        return one(T)
    end

    return prod( one(T)-u[j] for j = a+1:D)
end

# ℝ^D to boundary of D-simplex.
function computeJacobiandiagforW( y::Vector{T}) where T <: Real
    D = length(y)

    # pre-compute.
    one_minus_u = collect( one(T) - realtounitinterval(y[i]) for i = 1:D )

    # compute Jacobian.
    J_diag = Vector{T}(undef,D)
    for i = 1:D-1

        factor1 = ∂realtounitinterval(y[i])
        factor2 = prod( one_minus_u[j] for j = i+1:D)

        J_diag[i] = factor1*factor2
    end
    J_diag[D] = ∂realtounitinterval(y[D])

    return J_diag
end

function computelogdetJacobianforW( y::Vector{T})::T where T <: Real
    D = length(y)

    # pre-compute.
    ln_one_minus_u = collect( log(one(T) - realtounitinterval(y[i])) for i = 1:D )

    # compute Jacobian.
    ln_det_J = zero(T)
    for i = 1:D-1

        ln_factor1 = ln∂realtounitinterval(y[i])
        ln_factor2 = sum( ln_one_minus_u[j] for j = i+1:D)

        update = ln_factor1 + ln_factor2
        ln_det_J += update
    end
    ln_det_J += ln∂realtounitinterval(y[D])

    return ln_det_J
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
using JuMP

import GLPK
lp_solver = with_optimizer(GLPK.Optimizer)


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 = "libfyaml"
version = v"0.7.12"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/pantoniou/libfyaml/releases/download/v$(version)/libfyaml-$(version).tar.gz",
                  "485342c6920e9fdc2addfe75e5c3e0381793f18b339ab7393c1b6edf78bf8ca8"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libfyaml-*

# Replace `alloca.h` with `stdlib.h`
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/alloca.patch

./bootstrap.sh
./configure \
    --build=${MACHTYPE} \
    --host=${target} \
    --prefix=${prefix}
make -j${nproc}
make -j${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = filter(!Sys.iswindows, platforms)

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

    ExecutableProduct("fy-compose", :fy_compose),
    ExecutableProduct("fy-dump", :fy_dump),
    ExecutableProduct("fy-filter", :fy_filter),
    ExecutableProduct("fy-join", :fy_join),
    ExecutableProduct("fy-testsuite", :fy_testsuite),
    ExecutableProduct("fy-tool", :fy_tool),
    ExecutableProduct("fy-ypath", :fy_ypath),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
abstract type AbstractGeodesicPoint{T} end

@with_kw struct GeodesicPoint{T} <: AbstractGeodesicPoint{T}
    retcode::Symbol
    t::T
    u::AbstractVector{T}
    v::AbstractVector{T}
    # we don't store the problem parameters
    # and can create a specialistion for the carter method
    # then provide dispatched accessor methods
    # p::P
end

function geodesic_point_type(m::AbstractMetricParams{T}) where {T}
    GeodesicPoint{T}
end

# TODO: GeodesicPath structure for the full geodesic path
# do we want to support this?

function unpack_solution(sol::SciMLBase.AbstractODESolution{T,N,S}) where {T,N,S}
    u = sol.u
    p = sol.prob.p
    t = sol.t
    (u, t, p)
end

function unpack_solution(simsol::SciMLBase.AbstractEnsembleSolution{T,N,V}) where {T,N,V}
    map(unpack_solution, simsol)
end

function get_endpoint(
    m::AbstractMetricParams{T},
    sol::SciMLBase.AbstractODESolution{T,N,S},
) where {T,N,S}
    us, ts, _ = unpack_solution(sol)
    u = us[end].x[2]
    v = us[end].x[1]
    t = ts[end]
    GeodesicPoint(sol.retcode, t, u, v)
end

function get_endpoint(
    m::AbstractMetricParams{T},
    simsol::SciMLBase.AbstractEnsembleSolution{T,N,S},
) where {T,N,S}
    map(sol -> get_endpoint(m, sol), simsol)
end

export get_endpoint, GeodesicPoint


text:
function bubbleSort!(a)
    n = length(a)
    for i in 1:n-1
        for j in 1:n-i
            if a[j] > a[j+1]
                a[j], a[j+1] = a[j+1], a[j]
            end
        end
    end
    return a
end

data = [65, 51, 32, 12, 23, 84, 68, 1]
bubbleSort!(data)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using Statistics
using LinearAlgebra
using StochasticDiffEq

# Estimate the parameters of a Lorenz SDE

u0 = [1.0,1.0]
tspan = (0.0,10.0)
function multiplicative_noise!(du,u,p,t)
  x,y = u
  du[1] = p[5]*x
  du[2] = p[6]*y
end
p = ptrue = [1.5,1.0,3.0,1.0,0.1,0.1]

function lotka_volterra!(du,u,p,t)
  x,y = u
  α,β,γ,δ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = δ*x*y - γ*y
end

# define log.(u) via Ito

ξ0 = log.(u0)
function g!(du,u,p,t)
  x,y = u
  du[1] = p[5]
  du[2] = p[6]
  nothing
end

function f!(du,u,p,t)
  x,y = u
  α,β,γ,δ = p
  du[1] = dx = α - β*exp(y)
  du[2] = dy = δ*exp(x) - γ
  nothing
end
function f(u, p, t)
  du = zeros(2)
  f!(du, u, p, t)
  du
end

# f is linear in the parameters with Jacobian
function f_jac(J,u,p,t)
  x,y = u
  J .= false
  J[1,1] = 1.0
  J[1,2] = -exp(y)
  J[2,3] = -1.0
  J[2,4] = exp(x)
  nothing
end
# and intercept
function ϕ0(du,u,p,t)
  du .= false
  nothing
end

using Plots

prob_sde = SDEProblem(lotka_volterra!,multiplicative_noise!,u0,tspan,p)
ensembleprob = EnsembleProblem(prob_sde)
@time data = solve(ensembleprob,SOSRI(),saveat=0.1,trajectories=1000)
plot(EnsembleSummary(data))

prob_sde2 = SDEProblem(f!,g!,ξ0,tspan,p)
ensembleprob2 = EnsembleProblem(prob_sde2)
data2 = solve(ensembleprob2,SOSRI(),saveat=0.1,trajectories=1000)

@time data2exp = solve(ensembleprob2,SOSRI(),saveat=0.1,trajectories=1000)
for u in data2exp.u
  u.u .= [exp.(x) for x in u.u]
end

plot(EnsembleSummary(data2exp))

quvar(u) = sum((du - f(u, p, t)*dt).^2 for (t,u, dt, du) in zip(u.t, u.u, diff(u.t), diff(u.u)))

p = (0.5 .+ 1.5rand(length(ptrue))).*ptrue
ps = [copy(p)]
for k in 1:5
  global p
  # Solve the regression problem

  sdekernel = MSDE.SDEKernel(f!,g!,tspan,p)

  ϕprototype = zeros((2,4)) # prototypes for vectors
  yprototype = zeros((2,))
  R = MSDE.Regression!(sdekernel,yprototype,paramjac_prototype=ϕprototype,paramjac=f_jac,intercept=ϕ0)

  G = MSDE.conjugate(R, data2, 0.1I(4))

  # Estimate
  p̂ = mean(G)
  se = sqrt.(diag(cov(G)))
  display(map((p̂, se, p) -> "$(round(p̂, digits=2)) ± $(round(se, digits=2)) (true: $p)", p̂, se, ptrue))

  p[1:4] = rand(convert(Gaussian{(:μ, :Σ)}, G))

  # Estimate covariance
  T = sum(u.t[end] - u.t[1] for u in data2)

  Q = sum(quvar(u) for u in data2)

  p[end-1:end] = sqrt.(Q/T)
  push!(ps, copy(p))
end
ps
its = eachindex(ps)
cols = [k for j in its, k in eachindex(ptrue)]
pl = scatter(reduce(hcat, ps)', c=cols, labels=["α" "β" "γ" "δ" "s1" "s2"])
plot!(reduce(hcat, fill(ptrue, length(ps)))', c=cols, labels=["αtrue" "βtrue" "γtrue" "δtrue" "s1true" "s2true"])


text:
"""
pickett(args...)

pickett is plot recipe which plots the Pickett Plot for Water Saturation

Picket Plot is a graphical representation of the Archie Equation for estimating Sw

by rearranging the Archie equation

``{ Sw }^{ n }=\\frac { a\\quad { Rw } }{ { \\Phi  }^{ m }\\quad Rt } \\\\ \\\\ nLog\\left( Sw \\right) =Log\\left( a\\quad Rw \\right) -mLog\\left( \\phi  \\right) -Log(Rt)\\\\ \\\\ When\\quad Sw=1\\\\ \\\\ Log\\left( Rt \\right) =-mLog\\left( \\phi  \\right) +Log\\left( a\\quad Rw \\right)  ``

Where a plot of scale Log-Log produces a straight line representing Sw=1 with a slope of ``-m`` and intercept at Phi=100% of ``a*Rw``

Then, the desired Sw lines are plotted together with Rt and Phie from Logs.

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Rt|Mandatory|--|Array{Number,1}|Deep resistivity from Logs|
|Phie|Mandatory|--|Array{Number,1}|Effective Porosity from Logs|
|a|Optional|a=1|a=Number|a constant from archie Equation|
|m|Optional|m=2|a=Number|m cementation constant from archie Equation|
|n|Optional|n=2|a=Number|n exponent from archie Equation|
|Sw|Optional|Sw=range(0.2,stop=1,length=5)|Sw=range(Initial Sw,stop=Final Sw,length=Number)|Range of Sw lines|
|WellName|Optional|WellName=false|WellName="string"|Plot name|
"""
@userplot pickett

@recipe function f(h::pickett;
                    a=1,M=2,n=2,      #Default values for Archie equation
                    Sw=range(0.2,stop=1,length=5),  #Default Sw range to plot
                    WellName=false)

Rt, Phie, Rw=h.args

PhieRange=range(0.01,stop=1,length=5)        #Porosity range to map Sw lines
PhieRT1=(1/Rw).^(1/-M)                       #Porosiy at Rt=1 to estimate others Sw lines
rts=map(x->(a.*Rw)/(PhieRT1 .^M .* x.^n),Sw) # Rt at Sw Range with porosity PhieRT1
ResPhie1=map(x->(x.*PhieRT1.^M)./a,rts)        # Intercept of every SW to plot each lines at Phie=1
lines=zeros(size(PhieRange,1),size(Sw,1))   #Array of zeros of lines Array
Txpos=zeros(size(Sw,1))
Typos=0.9-0.1*size(Sw,1):0.1:0.9            #Array of Saturation Text positions

j=0
for i in ResPhie1                              #Calculate Sw Lines
    j=j+1
    lines[:,j]=map(x->a.*i.*x.^-M,PhieRange)
    Txpos[j]=a.*i.*Typos[j].^-M
end

tck=map(string,round.(Sw,digits=2))

## Picket Plot

xaxis := :log
yaxis := :log
xlim :=[0.1,20000]
ylim := [0.01,1]
xminorgrid := true
yminorgrid := true
title_location := :left
xlabel := "Rt[Ohm m]"
ylabel := "Phie[]"
xformatter --> :plain
yformatter --> :plain
yticks --> [0.01,0.1,1]
legend --> :bottomleft

c=range(RGB(0.0,1.0,0.0),stop=RGB(0.0,0.0,1.0),length=length(Sw))

cg=ColorGradient(c)
legendtitle := "Water Saturations[]"
xmirror := true

@series begin
    seriestype := :path
    line_z := (1:length(Sw))'
    linecolor := cg
    linewidth := 3
    lab := map(string,round.(Sw,digits=2))
    lines, PhieRange
end

@series begin
    seriestype := :scatter
    lab := "LogData"
    Rt, Phie
end

end #End Function


text:
using Documenter, LsqFit

makedocs(
    format = :html,
    sitename = "LsqFit.jl",
    doctest = false,
    strict = false,
    pages = Any[
            "Home" => "index.md",
            "Getting Started" => "getting_started.md",
            "Tutorial" => "tutorial.md",
            "API References" => "api.md",
            ],
    # Use clean URLs, unless built as a "local" build
    html_prettyurls = !("local" in ARGS),
    html_canonical = "https://julianlsolvers.github.io/LineSearches.jl/stable/"
    )

deploydocs(
    repo = "github.com/JuliaNLSolvers/LsqFit.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
export combineLatestUpdates

import Base: show

"""
    combineLatestUpdates(sources...; strategy = PushEach())
    combineLatestUpdates(sources::S, strategy::G = PushEach()) where { S <: Tuple, U }

`combineLatestUpdates` is a more effiecient version of `combineLatest(sources) + map_to(sources)` operators chain.

# Arguments
- `sources`: input sources
- `strategy`: optional update strategy for batching new values together

See also: [`Subscribable`](@ref), [`subscribe!`](@ref), [`PushEach`](@ref), [`PushEachBut`](@ref), [`PushNew`](@ref), [`PushNewBut`](@ref), [`PushStrategy`](@ref)
"""
function combineLatestUpdates end

combineLatestUpdates(; strategy = PushEach())                                       = error("combineLatestUpdates operator expects at least one inner observable on input")
combineLatestUpdates(args...; strategy = PushEach())                                = combineLatestUpdates(args, strategy)
combineLatestUpdates(sources::S, strategy::G = PushEach()) where { S <: Tuple, G }  = CombineLatestUpdatesObservable{S, G}(sources, strategy)

##

struct CombineLatestUpdatesInnerActor{L, W} <: Actor{L}
    index   :: Int
    wrapper :: W
end

Base.show(io::IO, ::CombineLatestUpdatesInnerActor{L, W}) where { L, W } = print(io, "CombineLatestUpdatesInnerActor($L, $I)")

on_next!(actor::CombineLatestUpdatesInnerActor{L, W}, data::L) where { L, W } = next_received!(actor.wrapper, data, actor.index)
on_error!(actor::CombineLatestUpdatesInnerActor{L, W}, err)    where { L, W } = error_received!(actor.wrapper, err, actor.index)
on_complete!(actor::CombineLatestUpdatesInnerActor{L, W})      where { L, W } = complete_received!(actor.wrapper, actor.index)

##

struct CombineLatestUpdatesActorWrapper{S, A, G, U}
    sources       :: S
    actor         :: A
    nsize         :: Int
    strategy      :: G # Push update strategy
    updates       :: U # Updates
    subscriptions :: Vector{Teardown}
end

function CombineLatestUpdatesActorWrapper(sources::S, actor::A, strategy::G) where { S, A, G }
    updates       = getustorage(S)
    nsize         = length(sources)
    subscriptions = fill!(Vector{Teardown}(undef, nsize), voidTeardown)
    return CombineLatestUpdatesActorWrapper(sources, actor, nsize, strategy, updates, subscriptions)
end

push_update!(wrapper::CombineLatestUpdatesActorWrapper) = push_update!(wrapper.nsize, wrapper.updates, wrapper.strategy)

dispose(wrapper::CombineLatestUpdatesActorWrapper) = begin fill_cstatus!(wrapper.updates, true); foreach(s -> unsubscribe!(s), wrapper.subscriptions) end

function next_received!(wrapper::CombineLatestUpdatesActorWrapper, data, index::Int)
    vstatus!(wrapper.updates, index, true)
    ustatus!(wrapper.updates, index, true)
    if all_vstatus(wrapper.updates) && !all_cstatus(wrapper.updates)
        push_update!(wrapper)
        next!(wrapper.actor, wrapper.sources)
    end
end

function error_received!(wrapper::CombineLatestUpdatesActorWrapper, err, index)
    if !(cstatus(wrapper.updates, index))
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::CombineLatestUpdatesActorWrapper, index::Int)
    if !all_cstatus(wrapper.updates)
        cstatus!(wrapper.updates, index, true)
        if ustatus(wrapper.updates, index)
            vstatus!(wrapper.updates, index, true)
        end
        if all_cstatus(wrapper.updates) || (!(vstatus(wrapper.updates, index)))
            dispose(wrapper)
            complete!(wrapper.actor)
        end
    end
end

##

@subscribable struct CombineLatestUpdatesObservable{S, G} <: Subscribable{S}
    sources  :: S
    strategy :: G
end

function on_subscribe!(observable::CombineLatestUpdatesObservable{S, G}, actor::A) where { S, G, A }
    wrapper = CombineLatestUpdatesActorWrapper(observable.sources, actor, observable.strategy)
    W       = typeof(wrapper)

    for (index, source) in enumerate(observable.sources)
        @inbounds wrapper.subscriptions[index] = subscribe!(source, CombineLatestUpdatesInnerActor{eltype(source), W}(index, wrapper))
        if cstatus(wrapper.updates, index) === true && vstatus(wrapper.updates, index) === false
            dispose(wrapper)
            break
        end
    end

    if all_cstatus(wrapper.updates)
        dispose(wrapper)
    end

    return CombineLatestUpdatesSubscription(wrapper)
end

##

struct CombineLatestUpdatesSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: CombineLatestUpdatesSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::CombineLatestUpdatesSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::CombineLatestUpdatesObservable{D}) where D  = print(io, "CombineLatestUpdatesObservable($D)")
Base.show(io::IO, ::CombineLatestUpdatesSubscription)           = print(io, "CombineLatestUpdatesSubscription()")

##


text:
using Pope.NoiseAnalysis
using ARMA
using HDF5
using Random
import ReferenceMicrocalFiles
using Pope
using Test

# generate some fake data to make a basis from
x=1:400
npulses = 3000
data = zeros(Float64,length(x),npulses)
# components
a=x[:]
b=20*(sinpi.(x/10) .+ 1)
c=10*(cospi.(x/10) .+ 1)
Random.seed!(0)
for i in 1:size(data,2)
    data[:,i] .+= a .* rand(1:10)
    data[:,i] .+= b .* rand(1:10)
    data[:,i] .+= c .* rand(1:10)
end
# data=round.(data)

model = ARMA.ARMAModel([2.41644, -4.25252, 1.85175], [1, -1.90467, .905899])
covar = ARMA.model_covariance(model, length(x))

noise_result = NoiseResult(ARMA.model_covariance(model, length(x)),
    ARMA.model_psd(model, length(x)),0,0,"",model)
frac_keep = 0.8
n_loop = 5
n_pulses_for_train=npulses
n_basis = 3
n_presamples = 0
tsvd_method_string="TSVD"
pulse_file = "dummy filename"
tsvd_basis, tsvd_basisinfo = Pope.create_basis_one_channel(data,noise_result,
    frac_keep, n_loop,
    n_pulses_for_train, n_basis,tsvd_method_string, n_presamples,
    pulse_file,-1)

h5open("artifacts/dummy_model.h5","w") do h5 Pope.hdf5save(h5,tsvd_basisinfo) end
loaded_basisinfo = h5open("artifacts/dummy_model.h5","r") do h5 Pope.hdf5load(Pope.SVDBasisWithCreationInfo,h5) end
loaded_basis = h5open("artifacts/dummy_model.h5","r") do h5 Pope.analyzer_from_preknowledge(h5) end

@testset "TSVD basis" begin
    @test tsvd_basisinfo.pulse_file == pulse_file
    @test loaded_basisinfo.svdbasis.basis==loaded_basis.basis==tsvd_basisinfo.svdbasis.basis

    # Test that inverted pulses don't get zeroed out
    invpulse = 200 .- [0, 0, 0, 10, 100, 80, 65, 55, 50, 46, 43, 41, 40.]
    noiseautocorr = 0.0*invpulse
    noiseautocorr[1] = 1.0
    # Make sure the following doen't error or return NaN singular values (it used to):
    a, _ = Pope.TSVD_tsvd_mass3(hcat(invpulse, 2*invpulse, 3*invpulse), 3, 3, noiseautocorr)
    @test !any(isnan.(a))
end

rmfdir = splitdir(pathof(ReferenceMicrocalFiles))[1]
noisepath = normpath("$(rmfdir)/../ljh/20150707_C_chan13.noi")
ljhpath = normpath("$(rmfdir)/../ljh/20150707_D_chan13.ljh")
noise_result_path = "artifacts/noise_result.h5"
model_path = "artifacts/model.h5"

# take a constant, arrival time, and average pulse from data
# train on residuals

tsvd_method_string="TSVDmass3"
pulse_file = "dummy filename"
n_basis=3 # test the case where the only 3 elements return are the mass3 (dc value, average pulse, pulse derivative)
mass3_basis, mass3_basisinfo = Pope.create_basis_one_channel(data,noise_result,
    frac_keep, n_loop,
    n_pulses_for_train, n_basis,tsvd_method_string, n_presamples,
    pulse_file,-1)

n_basis=5
mass3_basis, mass3_basisinfo = Pope.create_basis_one_channel(data,noise_result,
    frac_keep, n_loop,
    n_pulses_for_train, n_basis,tsvd_method_string, n_presamples,
    pulse_file,-1)

@testset "scripts with SVDBasis" begin
    # invoke scripts through bash to use of env variables
    @test success(run(`../scripts/noise_analysis.jl $noisepath -o $noise_result_path`))
    @test success(run(`../scripts/basis_create.jl $ljhpath $noise_result_path -o $model_path`))
    if !haskey(ENV,"POPE_NOMATPLOTLIB")
        @test success(run(`../scripts/basis_plots.jl $model_path`))
        @test success(run(`../scripts/noise_plots.jl $noise_result_path`))
    end
end