text:
################################
### crossvalidation routines ###
################################

# container object for XV results
immutable OMPCrossvalidationResults{T <: Float}
    mses :: Vector{T}
    path :: Vector{Int}
    b    :: Vector{T}
    bidx :: Vector{Int}
    k    :: Int
end

# subroutine to refit preditors after crossvalidation
function refit_omp{T <: Float}(
    x     :: DenseMatrix{T},
    y     :: DenseVector{T},
    k     :: Int;
    quiet :: Bool = true,
)
    # first use OMP to extract model
    betas = omp(x, y, k, quiet=quiet)

    # which components of β are nonzero?
    # cannot use binary indices here since we need to return Int indices
    bidx = find(betas[:,end])

    # allocate the submatrix of x corresponding to the inferred model
    x_inferred = x[:,bidx]

    # now estimate β with the ordinary least squares estimator β = inv(x'x)x'y
    # return it with the vector of MSEs
    xty = BLAS.gemv('T', x_inferred, sdata(y)) :: Vector{T}
    xtx = At_mul_B(x_inferred, x_inferred) :: Matrix{T}
    b   = zeros(T, length(bidx))
    try
        b = (xtx \ xty) :: Vector{T}
    catch e
        warn("caught error: ", e, "\nSetting returned values of b to -Inf")
        fill!(b, -Inf)
    end

    return b, bidx
end

"""
    one_fold(x, y, k, folds, fold) -> Vector

For a regularization path given by the integer `k`,
this function performs orthogonal matching pursuit on `x` and `y` and computes an out-of-sample error based on the indices given in `folds`.

Arguments:
- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `k` is an integer that specifies the maximum desired sparsity level to test.
- `folds` is an `Int` vector indicating which component of `y` goes to which fold, e.g. `folds = IHT.cv_get_folds(n,nfolds)`
- `fold` is the current fold to compute.

Optional Arguments:
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).

Output:
- `errors` is a vector of out-of-sample errors (MSEs) for the current fold.
"""
function one_fold{T <: Float}(
    x     :: DenseMatrix{T},
    y     :: DenseVector{T},
    k     :: Int,
    folds :: DenseVector{Int},
    fold  :: Int;
    quiet :: Bool  = true,
)

    # make vector of indices for folds
    test_idx = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx

    # allocate the arrays for the training set
    x_train = x[train_idx,:]
    y_train = y[train_idx]

    # compute the regularization path on the training set
    betas = omp(x_train, y_train, k, quiet=quiet)

    # compute the mean out-of-sample error for the TEST set
    #errors = vec(sumabs2(broadcast(-, y[test_idx], x[test_idx,:] * betas), 1)) ./ (2*test_size)
    xb = view(x, test_idx, :) * betas
    r  = broadcast(-, y[test_idx], xb)
    errors = vec(sumabs2(r, 1)) ./ (2*test_size)

    return errors :: Vector{T}
end

"""
    pfold(x, y, k, folds, q [, pids=procs()]) -> Vector

This function is the parallel execution kernel in `cv_omp`. It is not meant to be called outside of `cv_omp`.
For floating point data `x` and `y` and an integer `k`, `pfold` will distribute `q` crossvalidation folds across the processes supplied by the optional argument `pids`.
It calls `one_fold` for each fold, then collects the vectors of MSEs each process, applys a reduction, and finally returns the average MSE vector across all folds.
"""
function pfold{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    k        :: Int,
    folds    :: DenseVector{Int},
    q        :: Int;
    pids     :: Vector{Int} = procs(),
    quiet    :: Bool  = true,
)
    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray(T, (k,q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                                one_fold(x, y, k, folds, current_fold, quiet=quiet)
                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

"""
    cv_omp(x,y) -> Vector

This function will perform `q`-fold cross validation for the ideal model size in OMP least squares regression using the `n` x `p` design matrix `x` and the response vector `y`.
Each path is asynchronously spawned using any available processor.
For each path, one fold is held out of the analysis for testing, while the rest of the data are used for training.
The function to compute each path, `one_fold()`, will return a vector of out-of-sample errors (MSEs).
Arguments:
- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
Optional Arguments:
- `q` is the number of folds to compute. Defaults to `max(3, min(CPU_CORES, 5))`, where `CPU_CORES`is the Julia variable to query the number of available CPU cores.
- `path` is an `Int` vector that specifies which model sizes to include in the path. Defaults to `path = collect(1:min(p,20))`.
- `folds` is the partition of the data. Defaults to `IHT.cv_get_folds(n,q)`.
- `pids`, a vector of process IDs. Defaults to `procs()`, which recruits all available processes.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).
   *NOTA BENE*: each processor outputs feed to the console without regard to the others,
   so setting `quiet = false` can yield very messy output!
Output:
An `OMPCrossvalidationResults` object with the following fields:
- `mses` is the averaged MSE over all folds.
- `k` is the best crossvalidated model size.
- `path` is the regularization path used in the crossvalidation.
- `b`, a vector of `k` floats
- `bidx`, a vector of `k` indices indicating the support of the best model.
"""
function cv_omp{T <: Float}(
    x     :: DenseMatrix{T},
    y     :: DenseVector{T};
    q     :: Int  = cv_get_num_folds(3,5),
    k     :: Int  = min(size(x,2), 20),
    folds :: DenseVector{Int} = cv_get_folds(sdata(y),q),
    pids  :: Vector{Int} = procs(),
    quiet :: Bool = true,
)
    # do not allow crossvalidation with fewer than 3 folds
    q > 2 || throw(ArgumentError("Number of folds q = $q must be at least 3."))

    # problem dimensions?
    n,p = size(x)

    # check for conformable arrays
    n == length(y) || throw(DimensionMismatch("Row dimension of x ($n) must match number of rows in y ($(length(y)))"))

    # want to compute a path for each fold
    # the folds are computed asynchronously over processes enumerated by pids
    # master process then reduces errors across folds and returns MSEs
    mses = pfold(x, y, k, folds, q, pids=pids, quiet=quiet)

    # what is the best model size?
    path = collect(1:k)
    kbest = path[indmin(mses)] :: Int

    # print results
    !quiet && print_cv_results(mses, path, kbest)

    # refit best model
    #b, bidx = refit_omp(x, y, kbest, quiet=quiet)
    betas = omp(x, y, kbest, quiet=quiet)
    bidx  = find(betas[:,end])
    b     = full(betas[bidx,end])

    return OMPCrossvalidationResults{T}(mses, path, b, bidx, kbest)
end


text:
function genData(m::SearchDiscovery,nCons::Int,nProdTot::Int,nProdCons::Int;
					fname::String = "" , weights::AbstractWeights = fweights(fill(1,nProdTot-1)))

	Random.seed!(m.seed)
	# Initial checks
	@argcheck round((nProdCons-1-m.nA0)/m.nd)==(nProdCons-1-m.nA0)/m.nd  "Number of products (per consumer) does not match m.nA0 and m.nd combo."

	# Search cost specification  (relevant for Weitzman)
	cfun	= 	getCfun(m.cfun)
	zdfun   =   getZdfun(m.zdfun)

	# Consumer index into full data
	seli = [ 1+ii*nProdCons:(ii+1)*nProdCons for ii = 0:nCons-1]

	# Products' positions (same for all consumers)
	pos = fill([zeros(Int64,1+m.nA0) ; repeat(collect(Int64,1:(nProdCons-1-m.nA0)/m.nd),inner=m.nd)],nCons)

	# Generate products

	pid,chars = genProducts(m::Model,seli,nProdTot,nProdCons,weights=weights)

	# Gather data in NamedTuple for convenience
	d = (seli = seli,pid = pid,chars = chars,pos = pos)

	# Discovery values (same for all)
	# NOTE: Distribution does not include outside option dummy
	nChar = length(m.dChars)
	μxb = mean(m.dChars)'*m.β[1:nChar] .+ mean(m.dV)
	σxb = m.β[1:nChar]'*cov(m.dChars)*m.β[1:nChar] + var(m.dV)
	Gb = Normal(0,sqrt(σxb))
	Ξ = calcZd(Gb,cfun(m.cs),m.cd,m.nd) # calculates based on mean zero (see paper for details)
	zd = μxb .+ zdfun.(Ξ,m.dze,unique(pos[1]))

	# Utility
	e   = rand(m.dE,size(d.chars,1))
	v   = rand(m.dV,size(d.chars,1))
	xβ  = d.chars * m.β
	u   = xβ .+ v .+ e

	# Search values
	ξ = calcZs(0.0,m.dE,cfun(m.cs))[1]
	zs = xβ .+ v .+ ξ

	# Generate search paths
	path,purch,stop,searched = genSearchPaths(m,d,nProdCons,zd,zs,u)

	# Gather again data
	path = [path[path[:,i].>0,i] for i =1:nCons] # reduces path to array of arrays
	searched = [searched[:,i] for i = 1:nCons]
	d = (seli = seli,pid = pid,chars = chars,pos = pos,path = path, purch = purch, stop = stop, searched = searched, u=u , zs = zs, zd = zd)

	if !isempty(fname)
		@info "Saved data in bson format to $fname."
		@save fname d
	end

	return d
end

###################################################
function genSearchPaths(m::SearchDiscovery,d,maxProd,zd,zs,u)
	n = length(d.seli)

	# Same for all consumers
	path = zeros(Int64,maxProd,n)
	purch = fill(0,n)
	stop  = fill(0,n)
	searched = fill(false,maxProd,n)

  @inbounds Threads.@threads for k = 1:Threads.nthreads()
		threadRange = getrange(n)
		locPathi 	= zeros(Int64,maxProd)
		locSearchi  = fill(false,maxProd)
		locO 		= fill(false,maxProd,2)

		for ii in threadRange # Loop over consumers allocated to specific thread
			# Select consumer
			r = d.seli[ii]
			nProd = length(d.seli[ii])
			# Reset locals
			locPathi .= 0
			locSearchi .= false
			locO .= false
			locO[1,2] = true
			if m.nA0 > 0
				locO[2:m.nA0+1,1] .= true # Initial awareness set
			end
			# generate search paths (updates locPathi)
			purchi,stopi = genSearchPathi!(m,locPathi,locSearchi,zs[r],u[r],locO,zd,d.pos[ii])
			# Fill in values
			path[:,ii] = locPathi
			purch[ii] = purchi
			stop[ii] = stopi
			searched[:,ii] = locSearchi
		end

	end

	return path,purch,stop,searched
end

###################################################
function genSearchPathi!(m::SearchDiscovery,path,searched,zs,u,O,zd,pos)

	# Initial status
	cc_pos = 0  # position counter
	ind_pos = 1+m.nA0+1 # position index
	cc_s = 1 # counter of how many searches
	ind_p = 1 # index of maximum utility (initially outside option)
	u_max = u[ind_p] # max utility
	maxPos = maximum(pos)
	# Generate search paths
	while true
		# Discover more (unless reached end)
		val_s,ind_s = findmax_subsel(zs,O[:,1])
		if cc_pos < maxPos && zd[cc_pos+1] >  max(val_s,u_max)
			cc_pos += 1
			O[ind_pos:ind_pos+m.nd-1,1] .= true # put in awareness set
			ind_pos += m.nd
		# Purchase and end search
		elseif u_max >= val_s
			# Note that only reaches this condition if zd < y
			# -Inf added since S otherwise might be empty
			purch = ind_p # d[i].pid[ind_p]
			# path = path[1:cc_s-1]
			stop  = cc_pos
			return purch,stop
			break # end with purchase
		# Search best from awareness set
		else
			# Update awareness and consideration set
			O[ind_s,1] = false
			O[ind_s,2] = true
			path[cc_s] = ind_s
			searched[ind_s] = true
			cc_s += 1
			if u[ind_s] >= u_max
				ind_p = ind_s
				u_max = u[ind_p]
			end
		end
	end
end

######################################################
function calcWelfareDemand(m::SearchDiscovery,d; e = [] , v = [] )

	# Draw random numbers if not provided
	if isempty(e)
		Random.seed!(m.seed+4123)
		e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	end
	if isempty(v)
		Random.seed!(m.seed+2333)
		v = rand(m.dV,size(d.chars,1),m.options.nDraws)
	end

	# Derived values
	nCons = length(d.seli)
	maxProd = maximum(vcat([length(d.seli[i]) for i = 1:length(d.seli)]))

	# Form utility
	xβv = d.chars * m.β .+ v
	u   = xβv .+ e

	# Search cost specification
	cfun	= 	getCfun(m.cfun)
	zdfun   =   getZdfun(m.zdfun)
	cs = cfun(m.cs) # Note, need this here as cfun(m.cs) does not work in threaded

	# Search values (same for all consumers)
	ξ   = calcZs(0.0,m.dE,cs)[1]
	zs  = xβv .+ ξ

	# Discovery values
	μx = mean(m.dChars)' * m.β[1:end-m.outDummy] .+ mean(m.dV)
	σx = sum(m.β[1:end-m.outDummy].^2 .* var(m.dChars)) .+ var(m.dV)
	G   = Normal(0, sqrt(σx))
	zd0 = calcZd(G,cs,m.cd,m.nd)
	zd  = μx .+ zdfun.(zd0,m.dze,0:length(d.seli[1]))

	# Pre allocate
	sumU = zeros(Threads.nthreads())
	sumC = zeros(Threads.nthreads())
	demand = zeros(Int64,maxProd,Threads.nthreads())
	@inbounds Threads.@threads for k = 1:Threads.nthreads()
		threadRange = getrange(length(d.seli))

		locPathi 	= zeros(Int64,maxProd)
		locSearchi  = fill(false,maxProd)
		locO 		= fill(false,maxProd,2)

		locU = 0.0
		locC = 0.0
		locD = zeros(Int64,maxProd)
		for ii in threadRange
			@views posi = d.pos[ii]
			for dd in 1:size(e,2)
				@views zsd = zs[d.seli[ii],dd]
				@views ud = u[d.seli[ii],dd]
				locPathi .= 0
				locO .= false
				locO[1,2] = true
				if m.nA0 > 0
					locO[2:m.nA0+1,1] .= true # Initial awareness set
				end
				purchi,stopi = genSearchPathi!(m,locPathi,locSearchi,
								zsd,ud,locO,zd,posi)
				locU += ud[purchi]
				locC += cs * (findfirst(isequal.(locPathi,0)) - 1) + m.cd * stopi
				locD[purchi] += 1
			end
		end
		sumU[k] = locU / size(e,2)
		sumC[k] = locC / size(e,2)
		demand[:,k] .= locD
	end
	return sum(sumU) - sum(sumC), sum(sumU), sum(sumC), sum(demand,dims=2) ./ size(e,2) ./ nCons
end

ext: jl


text:
using StatsFuns
using ForwardDiff: Dual

function check_rmath(fname, statsfun, params, aname, a, isprob, rtol)
    v = statsfun(params..., a)
    rv = statsfun(params..., Dual(a)).value
    if isprob
        rd = abs(v / rv - 1.0)
        if rd > rtol
            error("$fname deviates too much from Rmath at " *
                "params = $params, $aname = $a:\n" *
                "  v = $v (rv = $rv)\n  |v/rv - 1| = $rd > $rtol.")
        end
    else
        τ = (1.0 + abs(rv)) * rtol
        ad = abs(v - rv)
        if ad > τ
            error("$fname deviates too much from Rmath at " *
                "params = $params, $aname = $a:\n" *
                "  v = $v (rv = $rv)\n  |v - rv| = $ad > $τ.")
        end
    end
end

function genericcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
  pdf = string(basename, "pdf")
  logpdf = string(basename, "logpdf")
  stats_pdf = eval(Symbol(pdf))
  stats_logpdf = eval(Symbol(logpdf))
  for i = 1:length(X)
    x = X[i]
    check_rmath(pdf, stats_pdf, params, "x", x, true, rtol)
    check_rmath(logpdf, stats_logpdf, params, "x", x, false, rtol)
  end
end

function genericcomp_tests(basename, configs)
  println("\ttesting $basename ...")
  for (params, data) in configs
      genericcomp(basename, params, data)
  end
end

### Test cases

@testset "Generic" begin
    genericcomp_tests("beta", [
        ((1.0, 1.0), 0.01:0.01:0.99),
        ((2.0, 3.0), 0.01:0.01:0.99),
        ((10.0, 2.0), 0.01:0.01:0.99),
    ])

    genericcomp_tests("binom", [
        ((1, 0.5), 0.0:1.0),
        ((1, 0.7), 0.0:1.0),
        ((8, 0.6), 0.0:8.0),
        ((20, 0.1), 0.0:20.0),
        ((20, 0.9), 0.0:20.0),
        ((20, 0.9), 0:20),
    ])

    genericcomp_tests("chisq", [
        ((1,), 0.0:0.1:8.0),
        ((4,), 0.0:0.1:8.0),
        ((9,), 0.0:0.1:8.0),
    ])

    genericcomp_tests("fdist", [
        ((1, 1), (0.0:0.1:5.0)),
        ((2, 1), (0.0:0.1:5.0)),
        ((5, 2), (0.0:0.1:5.0)),
        ((10, 1), (0.0:0.1:5.0)),
        ((10, 3), (0.0:0.1:5.0)),
    ])

    genericcomp_tests("gamma", [
        ((1.0, 1.0), (0.05:0.05:12.0)),
        ((0.5, 1.0), (0.05:0.05:12.0)),
        ((3.0, 1.0), (0.05:0.05:12.0)),
        ((9.0, 1.0), (0.05:0.05:12.0)),
        ((2.0, 3.0), (0.05:0.05:12.0)),
    ])

    # genericcomp_tests("nbeta", [
    #     ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
    #     ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
    #     ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
    #     ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
    # ])

    # genericcomp_tests("nchisq", [
    #     ((2, 1), 0.0:0.2:8.0),
    #     ((2, 3), 0.0:0.2:8.0),
    #     ((4, 1), 0.0:0.2:8.0),
    #     ((4, 3), 0.0:0.2:8.0),
    # ])

    # genericcomp_tests("nfdist", [
    #     ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
    #     ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
    #     ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
    # ])

    genericcomp_tests("norm", [
        ((0.0, 1.0), -6.0:0.01:6.0),
        ((2.0, 1.0), -3.0:0.01:7.0),
        ((0.0, 0.5), -3.0:0.01:3.0),
        ((0, 1), -3.0:3.0),
        ((0, 1), -3.0:0.01:3.0)
    ])

    # genericcomp_tests("ntdist", [
    #     ((0, 1), -4.0:0.1:10.0),
    #     ((0, 4), -4.0:0.1:10.0),
    #     ((2, 1), -4.0:0.1:10.0),
    #     ((2, 4), -4.0:0.1:10.0),
    # ])

    genericcomp_tests("pois", [
        ((1.0,), 0:30),
        ((10.0,), 0:42)
    ])

    genericcomp_tests("tdist", [
        ((1,), -5.0:0.1:5.0),
        ((2,), -5.0:0.1:5.0),
        ((5,), -5.0:0.1:5.0),
    ])

    #genericcomp_tests("srdist", [
    #    ((1,2), (0.0:0.1:5.0)),
    #    ((2,2), (0.0:0.1:5.0)),
    #    ((5,3), (0.0:0.1:5.0)),
    #    ((10,2), (0.0:0.1:5.0)),
    #    ((10,5), (0.0:0.1:5.0))
    #])
end


text:
#=
    types that mirror the layout used in Arb
        for the mag, arf, and arb structs
    see Cstructs.txt for the C library versions
=#

type MagFloat <: AbstractFloat
    radius_exponentOf2::Int
    radius_significand::UInt   ## radius is unsigned (nonnegative), by definition

    function MagFloat()
         z = new(zero(Int), zero(UInt64))
         ccall(@libarb(mag_init), Void, (Ptr{MagFloat}, ), &z)
         finalizer(z, c_release_mag)
         return z
    end
end

function c_release_mag(x::MagFloat)
  ccall(@libarb(mag_clear), Void, (Ptr{MagFloat}, ), &x)
end

    #       P is the precision in bits as a parameter
    #
type ArfFloat{P} <: AbstractFloat
    exponentOf2 ::Int
    nwords_sign::Int           ## Int, as words is an offset; lsb holds sign of significand
    significand1::UInt         ## UInt, as each significand word is a subspan of significand
    significand2::UInt         ##   the significand words are unsigned (sign is in nwords_sign)

    function ArfFloat()
         z = new{P}(0,0%UInt,0%UInt,0%UInt)
         ccall(@libarb(arf_init), Void, (Ptr{ArfFloat{P}}, ), &z)
         finalizer(z, c_release_arf)
         return z
    end
end

function c_release_arf{P}(x::ArfFloat{P})
  ccall(@libarb(arf_clear), Void, (Ptr{ArfFloat{P}}, ), &x)
end

    #       P is the precision in bits as a parameter
    #
type ArbFloat{P} <: AbstractFloat
                               ##     ArfFloat{P}
    exponentOf2 ::Int          ##        exponentOf2
    nwords_sign::Int           ##        nwords_sign
    significand1::UInt         ##        significand1
    significand2::UInt         ##        significand2
                               ###    ArbMag{P}
    radius_exponentOf2::Int    ####      radius_exponentOf2
    radius_significand::UInt   ####      radius_significand

    function ArbFloat()
         z = new{P}(0,0%UInt,0%UInt,0%UInt,0,0%UInt)
         ccall(@libarb(arb_init), Void, (Ptr{ArbFloat{P}}, ), &z)
         finalizer(z, c_release_arb)
         return z
    end
end

function c_release_arb{P}(x::ArbFloat{P})
  ccall(@libarb(arb_clear), Void, (Ptr{ArbFloat{P}}, ), &x)
end


text:
import NaiveNASflux
using Flux

@testset "Utils" begin

    struct BogusLayer <: NaiveNASflux.FluxLayer end

    @testset "Sizes" begin

        @test nin(Dense(3,4)) == 3
        @test nout(Dense(3,4)) == 4

        @test nin(Conv((2,), 3=>4)) == 3
        @test nout(Conv((2,), 3=>4)) == 4
        @test nin(Conv((1,2), 3=>6)) == 3
        @test nout(Conv((1,2), 3=>6)) == 6
        @test nin(Conv((1,2,3), 4=>5)) == 4
        @test nout(Conv((1,2,3), 4=>5)) == 5

        @test nin(ConvTranspose((2,), 3=>4)) == 3
        @test nout(ConvTranspose((2,), 3=>4)) == 4
        @test nin(ConvTranspose((1,2), 3=>6)) == 3
        @test nout(ConvTranspose((1,2), 3=>6)) == 6
        @test nin(ConvTranspose((1,2,3), 4=>5)) == 4
        @test nout(ConvTranspose((1,2,3), 4=>5)) == 5

        @test nin(DepthwiseConv((2,), 3=>4*3)) == 3
        @test nout(DepthwiseConv((2,), 3=>4*3)) == 12
        @test nin(DepthwiseConv((1,2), 3=>6*3)) == 3
        @test nout(DepthwiseConv((1,2), 3=>6*3)) == 18
        @test nin(DepthwiseConv((1,2,3), 4=>5*4)) == 4
        @test nout(DepthwiseConv((1,2,3), 4=>5*4)) == 20

        @test nin(CrossCor((2,), 3=>4)) == 3
        @test nout(CrossCor((2,), 3=>4)) == 4
        @test nin(CrossCor((1,2), 3=>6)) == 3
        @test nout(CrossCor((1,2), 3=>6)) == 6
        @test nin(CrossCor((1,2,3), 4=>5)) == 4
        @test nout(CrossCor((1,2,3), 4=>5)) == 5

        @test nin(Flux.Diagonal(3)) == nout(Flux.Diagonal(3)) == 3

        @test nin(LayerNorm(3)) == nout(LayerNorm(3)) == 3
        @test nin(BatchNorm(3)) == nout(BatchNorm(3)) == 3
        @test nin(InstanceNorm(3)) == nout(InstanceNorm(3)) == 3
        @test nin(GroupNorm(3,1)) == nout(GroupNorm(3,1)) == 3

        @test nin(RNN(3,4)) == 3
        @test nout(RNN(3,4)) == 4
        @test nin(LSTM(3,4)) == 3
        @test nout(LSTM(3,4)) == 4
        @test nin(GRU(3,4)) == 3
        @test nout(GRU(3,4)) == 4

        @test_throws ArgumentError nin(BogusLayer(), "dummy")
        @test_throws ArgumentError nout(BogusLayer(), "dummy")
    end

    @testset "Dims" begin
        @test actdim(Dense(3,4)) == 1

        @test actdim(Conv((1,2), 3=>6)) == 3

        @test actdim(ConvTranspose((1,2), 3=>6)) == 3

        @test actdim(DepthwiseConv((1,2), 3=>6)) == 3

        @test actdim(CrossCor((1,2), 3=>6)) == 3

        @test actdim(Flux.Diagonal(1)) == indim(Flux.Diagonal(2)) == outdim(Flux.Diagonal(3)) == 1

        @test actdim(RNN(3,4)) == 1
        @test actdim(LSTM(3,4)) == 1
        @test actdim(GRU(3,4)) == 1

        @test_throws ArgumentError actdim(BogusLayer())
        @test_throws ArgumentError indim(BogusLayer())
        @test_throws ArgumentError outdim(BogusLayer())
    end
end


text:
include("path.jl")
include("model.jl")
using TrajectoryOptimization
using TrajOptPlots
import TrajectoryOptimization: set_state!
using Plots

function gen_car_prob(U0, path, N; s0=0.0, L=10.0, δ0=0.0, integration=RK3)

	# Discretization
	sf = L + s0
	ds = (sf - s0)/(N-1)

	s_traj = range(s0,sf,step=ds)
	@assert N == length(s_traj)
	# k_traj = k_itp.(s_traj)
	# d_path = Dict(s_traj.=>k_traj)

	# Model
	car = BicycleCar(path=path)
	n,m = size(car)
	p_c = 0.1 # probability of contingency plan

	x,u = zeros(car)
	car.μ = 0.9
	fiala(a) = fiala_tire_model(car.μ, car.Cαf, a, 0.0, 10000)
	plot(fiala.(range(-0.2,0.2, length=100)))

	# Objective
	Ux0 = 15.
	Ux_des = 15.
	x0 = @SVector [δ0,     # δ
	               0,      # fx
	               0,      # r
	               0,      # Uy
	               Ux0, # Ux
	               0,      # dpsi
	               0,      # e
	               0]      # t

	xd = @SVector [0,      # δ
	               0,      # fx
	               0,      # r
	               0,      # Uy
	               Ux_des, # Ux
	               0,      # dpsi
	               0,      # e
	               0]      # t

	Qd = @SVector [0.000,      # δ
	               0.000,      # fx
	               0.00,      # r
	               0.00,      # Uy
	               0.1,      # Ux
	               10.0,      # dpsi
	               10.0,      # e
	               0.0]      # t

	Qf = @SVector [0.000,      # δ
	               0.000,      # fx
	               0.000,      # r
	               0.000,      # Uy
	               100.0,      # Ux
	               100.0,      # dpsi
	               100.0,      # e
	               0.00]      # t

	Rd   = @SVector [0.01,  # δ_dot
	                 1.0]    # fx_dot
	Rd_c = @SVector [0.3282,  # δ_dot
	                 4e-8]    # fx_dot

	# Q = Diagonal([(1-p_c)*Qd; p_c*Qd])
	# R = Diagonal([Rd; Rd_c])

	# Qd = @SVector fill(0.1,n)
	# Rd = @SVector fill(0.01,m)
	Q = Diagonal(Qd)
	R = Diagonal(Rd)

	obj = LQRObjective(Q,R,Q,xd,N)

	# Bound Constraints
	δ_dot_bound = deg2rad(90)  # deg/s
	δ_bound = deg2rad(27)  # deg
	Fx_max = car.μ*car.mass*car.g / 1000  # friction cone?
	Ux_min = 1.0  # m/s
	e_bnd = 1.0 # m

	x_min = fill(-Inf,n)
	x_min[1] = -δ_bound
	x_min[5] = Ux_min
	x_min[7] = -e_bnd

	x_max = fill(Inf,n)
	x_max[1] = δ_bound
	x_max[2] = Fx_max
	x_max[7] = e_bnd

	u_min = [-δ_dot_bound; -Inf]
	u_max = [ δ_dot_bound;  Inf]

	bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min, u_min=u_min, u_max=u_max)

	# Brake constraint
	brake = BrakeForceConstraint(car)

	# Constraint Set
	conSet = ConstraintSet(n,m,N)
	add_constraint!(conSet, bnd, 1:N-1)
	add_constraint!(conSet, brake, 1:N)

	# Problem
	prob = Problem(car, obj, xd, sf, constraints=conSet, x0=x0, t0=s0,
		integration=integration)
	initial_controls!(prob, U0)

	return prob
end

# # Straight with circle
# t = range(0,20,length=901)
# radius = 20.
# X_curve = cos.(pi .- pi*t/30)*radius
# Y_curve = sin.(pi*t/30)*radius
# X_straight = fill(-radius, 100)
# Y_straight = collect(range(-5,0,length=100))
# X = [X_straight; X_curve]
# Y = [Y_straight; Y_curve]
#
# # Spiral
# t = range(0,20,length=901)
# radius = range(100,10,length=901)
# X = cos.(pi*t / 10) .* radius
# Y = sin.(pi*t / 10) .* radius
#
# s,ϕ,k = pathToLocal(X,Y)
# path = CarPath(s,ϕ,k,X,Y)
# plot(path.X,path.Y, aspect_ratio=:equal)
# s[end]

# Circle Path
line = StraightPath(30., pi/2)
arc = ArcPath(line, 20, 3pi/4)
path = DubinsPath([line, arc])
# path = line
plot(path, aspect_ratio=:equal)

N = 101
U0 = [@SVector [0,0] for k = 1:N-1]
prob = gen_car_prob(U0, path, N, s0=0.0, δ0=deg2rad(0), L=30.0, integration=RK3)

solver = iLQRSolver(prob)
solver = AugmentedLagrangianSolver(prob)
initial_controls!(solver, U_sol)
solver.opts.verbose = true
solver.opts.opts_uncon.verbose = true
solve!(solver)
plot(solver)
states(solver)[end][5]

cost(solver)
iterations(solver)
plot(controls(solver))
plot(states(solver), 1:2)
evaluate!(prob.constraints, get_trajectory(solver))
b = [v[1] for v in prob.constraints[2].vals]
plot(b)

plot(controls(prob))
plot(states(prob), 1:2)
plot(states(prob), 4:5)
plot(states(prob), 6:6)
plot(states(prob), 7:7)
plot(states(prob), 8:8)
states(prob)[end][5]
cost(prob)
TO.stage_cost(prob.obj.cost[10], states(prob)[10], controls(prob)[10])

Fx = states(prob.Z, 2)
Ux = states(prob.Z, 5)
dpsi = states(prob.Z, 6)
Fxf = [FWD_force_model(prob.model, fx*1e3)[1] for fx in Fx]
Fxr = [FWD_force_model(prob.model, fx*1e3)[2] for fx in Fx]
plot(Ux)
plot(Fxf)
plot(Fxr)

Uxdot = [dynamics(prob.model, z)[5] for z in prob.Z]
plot(Uxdot)

# Solve the initial problem
function initial_solve()
	N = 51
	U0 = [@SVector zeros(2) for k = 1:N-1]
	prob = gen_car_prob(U0, path, N, s0=0.0, δ0=deg2rad(0), L=30.0, integration=RK2)
	solver = AugmentedLagrangianSolver(prob)
	solver.opts.penalty_initial = 0.1
	solver.opts.opts_uncon.verbose = false
	solver.opts.verbose = true
	initial_controls!(solver, U0)
	rollout!(solver)
	solve!(solver)
end

function advance_problem!(solver)
	Z = get_trajectory(solver)

	# Advance the time and controls by one time step
	shift_fill!(Z)

	# Advance the duals
	solver.opts.reset_duals = false
	shift_fill!(get_constraints(solver))
	solver.opts.penalty_initial = 10.0
	solver.opts.penalty_scaling = 10.0
	solver.opts.cost_tolerance_intermediate = 1e-3
	solver.opts.constraint_tolerance = 1e-3

	# Set the new initial condition
	set_initial_state!(solver, state(Z[2]))

	nothing
end
solver = initial_solve()
TO.get_initial_state(solver)
states(solver)[end][5]
plot(solver)
plot(controls(solver))
findmax_violation(solver)

Nruns = 150
times = zeros(Nruns)
solver = initial_solve()
solver.opts.verbose = false
for i = 1:Nruns
	t = @elapsed begin
		advance_problem!(solver)
		solve!(solver)
	end
	s = get_trajectory(solver)[1].t
	println("run $i: iters=$(iterations(solver)), s=$s")
	times[i] = t
end
times
@show median(times)*1000
@show maximum(times)*1000
plot(solver)
plot(states(solver),7:7)

# advance_problem!(solver)
# U_sol = deepcopy(controls(solver))
#
# solver.opts.verbose = true
# solver.opts.reset_duals = true
# solver.opts.reset_penalties
# initial_controls!(solver, U_sol)
# solve!(solver)
# max_violation(solver)
# findmax_violation(solver)
# plot(solver)
# get_times(solver)[1]
# x0_sol = copy(TO.get_initial_state(solver))
#
# prob = gen_car_prob(U_sol, path, N, s0=get_times(solver)[1], δ0=deg2rad(0), L=30.0, integration=RK3)
# solver0 = AugmentedLagrangianSolver(prob)
# set_initial_state!(solver0, x0_sol)
# initial_controls!(solver0, U_sol)
# solver0.opts.verbose = true
# solve!(solver0)
#
# controls(solver) ≈ controls(solver0)
# TO.get_initial_state(solver) ≈ TO.get_initial_state(solver0)

#
# plot(path, aspect_ratio=:equal)
# e = [z.z[7] for z in Z]
# s = [z.t for z in Z]
# x,y = localToGlobal(path, e, s)

# # Simulate the next time step
# U_guess = U_sol[1:end]
# initial_controls!(prob, U_sol)
# solver2 = AugmentedLagrangianSolver(prob)
# rollout!(solver2)
# plot(solver2)
# solve!(solver)
# iterations(solver)
# plot(solver)
#
Z = get_trajectory(solver)
plot(states(Z),1:2)  # steering, accel
plot(states(Z),3:3)  # yaw rate
plot(states(Z),4:4)  # slip velocity
plot(states(Z),5:5)  # velocity
plot(states(Z),6:6)  # heading error
plot(states(Z),7:7)  # lateral error
plot(states(Z),8:8)  # time
plot(controls(Z))
conSet = get_constraints(solver)
b = [v[1] for v in conSet.constraints[2].vals]
plot(b)
# #
# # plot(path.X,path.Y, aspect_ratio=:equal)
# # plot!(x,y)
plot(solver)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -9.75097970893739 <= q <= 8.436779502208893)

                     add_NL_constraint(m, :(log(1 + exp(0.5456914178120993 + -0.994137029231358*log(1 + exp(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3]))) + 0.38451384631825425*log(1 + exp(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3]))) + 0.9673084611668079*log(1 + exp(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))))) + log(1 + exp(0.6143810270772536 + 0.32603207961941205*log(1 + exp(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3]))) + 0.7660966505536861*log(1 + exp(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3]))) + 0.30918160231413605*log(1 + exp(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))))) + log(1 + exp(-0.4079478457750749 + 0.06443443960646311*log(1 + exp(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3]))) + 0.09351478982242645*log(1 + exp(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3]))) + -0.7258007123919756*log(1 + exp(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@json struct PullRequestBranch
    default_merge_strategy::String
    merge_strategies::Vector{String}
    name::String
end
@json_struct PullRequestBranch

@json struct PullRequestCommit
    hash::String
end
@json_struct PullRequestCommit

@json struct PullRequestEndpoint
    branch::Branch
    commit::PullRequestCommit
    repository::Repo
end
@json_struct PullRequestEndpoint

"""
    RenderedPullRequestMarkup

title, description, and reason are Dicts like:
Dict("raw"=> "<string>", "markup"=> "<string>", "html"=> "<string>")
"""
@json struct RenderedPullRequestMarkup
    description::Dict{String, Any}
    reason::Dict{String, Any}
    title::Dict{String, Any}
end
@json_struct RenderedPullRequestMarkup

@json struct PullRequest
    author::User
    close_source_branch::Bool
    closed_by::Union{User, Nothing}
    comment_count::Int
    created_on::DateTime
    description::String
    destination::PullRequestEndpoint
    id::Int
    links::NamedTuple
    merge_commit::PullRequestCommit
    participants::Vector{User}
    reason::String
    rendered::RenderedPullRequestMarkup
    reviewers::Vector{User}
    source::PullRequestEndpoint
    state::String
    summary::Dict{String, Any}
    task_count::Int
    title::String
    updated_on::DateTime
end
@json_struct PullRequest

@not_implemented(::BitbucketAPI, ::typeof(get_pull_request), ::String, ::String, ::UUID)
@not_implemented(::BitbucketAPI, ::typeof(get_pull_request), ::UUID, ::UUID)

endpoint(::BitbucketAPI, ::typeof(get_pull_requests), workspace::String, repo::String) =
    Endpoint(:GET, "/repositories/$workspace/$repo/pullrequests")
@not_implemented(::BitbucketAPI, ::typeof(get_pull_requests), ::UUID)
into(::BitbucketAPI, ::typeof(get_pull_requests)) = PullRequest

endpoint(::BitbucketAPI, ::typeof(create_pull_request), workspace::String, repo::String) =
    Endpoint(:POST, "/repositories/$workspace/$repo/pullrequests")
@not_implemented(::BitbucketAPI, ::typeof(create_pull_request), ::UUID)
into(::BitbucketAPI, ::typeof(create_pull_request)) = PullRequest

endpoint(::BitbucketAPI, ::typeof(update_pull_request), workspace::String, repo::String, id::Int64) =
    Endpoint(:PUT, "/repositories/$workspace/$repo/pullrequests/$id")
@not_implemented(::BitbucketAPI, ::typeof(update_pull_request), ::UUID, ::Int64)
into(::BitbucketAPI, ::typeof(update_pull_request)) = PullRequest

@not_implemented(::BitbucketAPI, ::typeof(subscribe_to_pull_request), ::UUID, ::UUID)

@not_implemented(::BitbucketAPI, ::typeof(unsubscribe_from_pull_request), ::UUID, ::UUID)


text:
module RollingTimeWindows

import Dates
export RollingTimeWindow

struct RollingTimeWindow
    timestamps::AbstractArray{T} where {T<:Dates.AbstractDateTime}
    period::Dates.Period
end

struct RollingTimeWindowState
    index::Int
    start_timestamp::Union{Nothing, T} where {T<:Dates.AbstractDateTime}
end
RollingTimeWindowState() = RollingTimeWindowState(1, nothing)

function Base.iterate(
    rolling::RollingTimeWindow,
    state::RollingTimeWindowState=RollingTimeWindowState()
)
    max_index = length(rolling.timestamps)
    if state.index > max_index
        return nothing
    end
    index = state.index
    start_index = index
    start_timestamp = state.start_timestamp === nothing ? floor(rolling.timestamps[start_index], rolling.period) : state.start_timestamp
    end_timestamp = start_timestamp + rolling.period
    while index <= max_index && rolling.timestamps[index] < end_timestamp
        index += 1
    end
    end_index = index - 1
    next_state = RollingTimeWindowState(index, end_timestamp)
    return ((start_index:end_index, start_timestamp, end_timestamp), next_state)
end

Base.IteratorSize(IterType::RollingTimeWindow) = Base.SizeUnknown()
Base.IteratorEltype(IterType::RollingTimeWindow) = Base.HasEltype()
Base.eltype(IterType::RollingTimeWindow) = Tuple{UnitRange{Int}, T, T} where {T <: Dates.AbstractDateTime}

end


text:
function bits_for_digits(digits)
    Int64(ceil(digits * log2(10)))
end

function mp(A::UpperTriangular, digits)
    UpperTriangular(mp(A.data, digits))
end

function with_digits(fun, digits)
    MPFR_TLS.setprecision(bits_for_digits(digits)) do
        return fun();
    end
end

function mp(A::Array{T, N}, digits::Int64) where {T <: Complex, N}
    MPFR_TLS.setprecision(bits_for_digits(digits)) do
        return convert(Array{Complex{MPFR_TLS.BigFloat}, N}, A)
    end
end

function mp(x::T, digits::Int64) where T <: Complex
    MPFR_TLS.setprecision(bits_for_digits(digits)) do
        return convert(Complex{MPFR_TLS.BigFloat}, x)
    end
end

function mp(A::Array{T, N}, digits::Int64) where {T <: AbstractFloat, N}
    MPFR_TLS.setprecision(bits_for_digits(digits)) do
        return convert(Array{MPFR_TLS.BigFloat, N}, A)
    end
end

function mp(x::T, digits::Int64) where T <: AbstractFloat
    MPFR_TLS.setprecision(bits_for_digits(digits)) do
        return convert(MPFR_TLS.BigFloat, x)
    end
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
using ExtremeLearning
using Documenter

DocMeta.setdocmeta!(ExtremeLearning, :DocTestSetup, :(using ExtremeLearning); recursive=true)

makedocs(;
    modules=[ExtremeLearning],
    authors="Gabriel Baraldi",
    repo="https://github.com/gbaraldi/ExtremeLearning.jl/blob/{commit}{path}#{line}",
    sitename="ExtremeLearning.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gbaraldi.github.io/ExtremeLearning.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Functions" => "functions.md"
    ],
)

deploydocs(;
    repo="github.com/gbaraldi/ExtremeLearning.jl",
)


text:
using Plots, SNN; plotly(); reload("SNN")

S = SNN.Rate(;N = 200)
SS = SNN.FLSynapse(S, S; σ = 1.5, p = 1.0)
P, C = [S], [SS]

SNN.monitor(SS, [:f, :z])

A = 1.3 / 1.5; fr = 1 / 60ms
f(t) = (A /1.0) * sin(1π * fr * t) +  (A / 2.0) * sin(2π * fr * t) +
(A / 6.0) * sin(3π * fr * t) +   (A / 3.0) * sin(4π * fr * t)

for t in 0:0.1ms:1440ms
    SS.f = f(t)
    SNN.train!(P, C, 0.1ms)
end

for t in 1440ms:0.1ms:1700ms
    SS.f = f(t)
    SNN.sim!(P, C, 0.1ms)
end

plot([SNN.getrecord(SS, :f) SNN.getrecord(SS, :z)], label = ["f" "z"]);
vline!([1440ms/0.1ms], color = :cyan, label = "") |> gui


text:
export @require

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

@init begin
  push!(Base.package_callbacks, loadmod)
end

loaded(mod::Symbol) = getthing(Main, mod) != nothing

const modlisteners = Dict{Symbol, Vector{Function}}()

listenmod(f, mod::Symbol) =
  loaded(mod) ? f() :
    modlisteners[mod] = push!(get(modlisteners, mod, Function[]), f)

function loadmod(mod)
  fs = get(modlisteners, mod, Function[])
  delete!(modlisteners, mod)
  map(f->Base.invokelatest(f), fs)
end

importexpr(mod::Symbol) = Expr(:import, mod)
importexpr(mod::Expr) = Expr(:import, map(Symbol, split(string(mod), "."))...)

function withpath(f, path)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end

function err(f, listener, mod)
  try
    f()
  catch e
    warn("Error requiring $mod from $listener:")
    showerror(STDERR, e, catch_backtrace())
    println(STDERR)
  end
end

macro require(mod, expr)
  ex = quote
    listenmod($(QuoteNode(mod))) do
      withpath(@__FILE__) do
        err($(current_module()), $(string(mod))) do
          $(esc(:(eval($(Expr(:quote, Expr(:block,
                                           importexpr(mod),
                                           expr)))))))
        end
      end
    end
  end
  quote
    if isprecompiling()
      @init @guard $(ex)
    else
      $(ex)
    end
  end
end


text:
# ColorScore.m is run in ColorGraph.m Calculates the score for a graph and
# passes it back as GraphScore
#tic()
m=0
aScore=0
bScore=0
cScore=0
while m<Nodes
    #Count
    m=m+1
    #Determine originating Node
    origin=rainbow[m,1]
    if origin==0
        # Then act as if it is connect to itself
        # determine terminal node color
        termColor=rainbow[m,2]
        # Currently just summing properties, but could replace with more()
        # complicated operation
        aScore=aScore+A[7,termColor]
        bScore=bScore+B[7,termColor]
        cScore=cScore+C[7,termColor]
    else()
        orgColor=rainbow[rainbow[m,1],2]
        # determine terminal node color
        termColor=rainbow[m,2]
        # Currently just summing properties, but could replace with more()
        # complicated operation
        aScore=aScore+A[orgColor,termColor]
        bScore=bScore+B[orgColor,termColor]
        cScore=cScore+C[orgColor,termColor]
    end
end
COUNTRUNS=COUNTRUNS+1
# Record Graph's Score
GraphScore=[aScore,bScore,cScore]
clear m origin orgColor termColor aScore bScore cScore
#toc()

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
# Property (*) means that the step size α*d is small -> β = 0
struct CGRestart{T} end
struct FletcherReeves{T} end
# using an "intial vectors" function we can initialize s if necessary or nothing if not to save on vectors
function optimize(obj::OnceDiffed, x0, cg::FletcherReeves, linesearch = Backtracking(); maxiter=100)
    T=eltype(x0)
    fx, ∇fx = obj(x0, copy(x0))
    Tf = typeof(fx)
    fz = fx
    d = -copy(∇fx)
    ∇fz = copy(∇fx)
    k = 0
    x = copy(x0)
    while norm(∇fx, Inf) ≥ 1e-8 && k ≤ maxiter
        α, f_α, ls_success = NLSolvers.find_steplength(linesearch, obj, d, x, Tf(1.0), fx, ∇fx)
        z=x+α*d
        fz, ∇fz = obj(z, ∇fz)

        βkp1 = calc_β(cg, ∇fz, ∇fx)
        fx = fz
        x .= z
        ∇fx .= ∇fz

        d .= -∇fx + βkp1*d
        k = k+1
    end
    x, fx
end
cg_restart(::FletcherReeves{<:CGRestart{true}}) = true

function calc_β(cg, ∇fz, ∇fx)
    T = eltype(∇fz)
    βkp1 = dot(∇fz, ∇fz)/dot(∇fx, ∇fx)
    if cg_restart(cg)
        βkp1 = max(βkp1, T(0))
    end
    βkp1
end
# function calc_β(method, ∇fx, ∇fz)
#
# if method == :fr
# return dot(∇fz, ∇fz)/dot(∇fx, ∇fx)
# elseif method == :pr
# β = dot(∇fz, s)/dot(∇fx, ∇fx) # dot(∇fx, ∇fx) can be saved between iteratoins only s is needed
# return    max(β, 0) # We need this according to Powell 1986 Gilbert Powell 1980 show it is globally  convergenct for exact and inexact
# # pr+ max(beta, 0)
# elseif method ==:hs
# dot(∇fz, s)/dot(s,p)
# end


text:
"""
Base struct for all DiffEqEnvironments.jl IC samplers.
"""
struct ICSampler
    sample
end

(ics::ICSampler)() = ics.sample()

"""
Samples initial condition from Uniform distribution bounded by upper and lower bounds
"""
function UniformSampler(s0_lb::S, s0_ub::S; rng=Random.GLOBAL_RNG) where {S <:Union{Real,Vector{<:Real}}}
    length(s0_lb) == length(s0_ub) || throw(ArgumentError("Lengths $(length(s0_lb)) and $(length(s0_ub)) of bounds don't match."))
    return ICSampler(() -> rand(rng, length(s0_lb)) .* (s0_ub .- s0_lb) .+ s0_lb)
end


text:
module SeisJOLI

    using JOLI

    # convert shot record to common-midpoint offset
    include("operators/joSRtoCMO.jl")

    # units-dependent Radon transform
    include("operators/joRadon.jl")

end


text:
module AccuracyAtTopPrimal

using LinearAlgebra, Statistics, Flux, Roots, Parameters

using Zygote: @adjoint, @nograd

export PatMat, PatMatNP, TopPush, TopPushK, τFPL, TopMean, Grill, GrillNP
export fnr, fpr, hinge, quadratic, threshold

abstract type AbstractThreshold end

include("thresholds.jl")
include("utilities.jl")

end # module


text:
# module NeighbourOperators

# using Random: randperm
# using Distributions: NegativeBinomial
# using StatsBase: sample, ProbabilityWeights
# using Solutions: Solution, swap!, permu!

# export small_shift!,
#        swap_close_planes!,
#        swap_close_costly_planes!,
#        permutation_operator!

function small_shift!(sol::Solution; p=0.5)
    n = sol.inst.nb_planes
    shifted_plane = sample(1:n)
    new_position = find_shift_destination(shifted_plane, p, n)
    shift!(sol, shifted_plane, new_position)
end

function shift_costly_plane!(sol::Solution; p=0.5)
    shifted_plane = cost_weighted_sample(sol)
    new_position = find_shift_destination(shifted_plane, p,
                                          sol.inst.nb_planes)
    shift!(sol, shifted_plane, new_position)
end

"""Given a position, returns a random but valid destination.

The size of the shift is a random variable 1 + X,
where X follows a Negative Binomial of parameters (1, p)
"""
function find_shift_destination(current_position::Int, p::Float64, n_planes::Int)
    negative_binomial = NegativeBinomial(1, p)
    direction = sample([-1, 1])
    norm = 1 + rand(negative_binomial)
    new_position = current_position + direction * norm
    new_position = min(max(new_position, 1), n_planes)
    return new_position
end

"""Swaps a random plane and one of its neighbour."""
function swap_close_planes!(sol::Solution)
    n = sol.inst.nb_planes
    swapped_plane_1_id = sample(1:n)
    if swapped_plane_1_id == n
        swapped_plane_2_id = n-1
    elseif swapped_plane_1_id == 1
        swapped_plane_2_id = 2
    else
        swapped_plane_2_id = swapped_plane_1_id + sample([-1, 1])
    end
    swap!(sol, swapped_plane_1_id, swapped_plane_2_id)
end

"""Swaps a random plane according to the cost distribution and one of its neighbour."""
function swap_close_costly_planes!(sol::Solution)
    n = sol.inst.nb_planes
    swapped_plane_1_id = cost_weighted_sample(sol)
    if swapped_plane_1_id == n
        swapped_plane_2_id = n-1
    elseif swapped_plane_1_id == 1
        swapped_plane_2_id = 2
    else
        swapped_plane_2_id = swapped_plane_1_id + sample([-1, 1])
    end
    swap!(sol, swapped_plane_1_id, swapped_plane_2_id)
end

function permutation_operator!(sol::Solution; permuted_ratio=0.1)
    n_permuted_planes = Int(round(permuted_ratio * sol.inst.nb_planes))
    println(n_permuted_planes)
    if n_permuted_planes <= 2
        swap!(sol)
    else
        permuted_planes = sample(1:sol.inst.nb_planes, n_permuted_planes,
                                replace=false)
        permu!(sol, permuted_planes,
               permuted_planes[randperm(length(permuted_planes))])
    end
end

function cost_weighted_sample(sol::Solution)
    sample(1:sol.inst.nb_planes, ProbabilityWeights(sol.costs))
end

# end


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic

"""
    ShowDemoWindowPopups()
The properties of popups windows are:
- they block normal mouse hovering detection outside them. (*)
- unless modal, they can be closed by clicking anywhere outside them, or by pressing ESCAPE.
- their visibility state (~bool) is held internally by imgui instead of being held by the
    programmer as we are used to with regular Begin() calls.
User can manipulate the visibility state by calling OpenPopup().

(*) One can use IsItemHovered(ImGuiHoveredFlags_AllowWhenBlockedByPopup) to bypass it and detect hovering even when normally blocked by a popup.
Those three properties are connected. The library needs to hold their visibility state because it can close popups at any time.

Typical use for regular windows:
  bool my_tool_is_active = false; if (CImGui.Button("Open")) my_tool_is_active = true; [...] if (my_tool_is_active) Begin("My Tool", &my_tool_is_active) { [...] } End();
Typical use for popups:
  if (CImGui.Button("Open")) CImGui.OpenPopup("MyPopup"); if (CImGui.BeginPopup("MyPopup") { [...] EndPopup(); }

With popups we have to go through a library call (here OpenPopup) to manipulate the visibility state.
This may be a bit confusing at first but it should quickly make sense. Follow on the examples below.
"""
function ShowDemoWindowPopups()
    CImGui.CollapsingHeader("Popups & Modal windows") || return

    if CImGui.TreeNode("Popups")
        CImGui.TextWrapped("When a popup is active, it inhibits interacting with windows that are behind the popup. Clicking outside the popup closes it.")

        @cstatic selected_fish=Cint(-1) names=["Bream", "Haddock", "Mackerel", "Pollock", "Tilefish"] toggles=[true, false, false, false, false] begin
            # simple selection popup
            # (If you want to show the current selection inside the Button itself, you may want to build a string using the "###" operator to preserve a constant ID with a variable label)
            CImGui.Button("Select..") && CImGui.OpenPopup("my_select_popup")
            CImGui.SameLine()
            CImGui.TextUnformatted(selected_fish == -1 ? "<None>" : names[selected_fish])
            if CImGui.BeginPopup("my_select_popup")
                CImGui.Text("Aquarium")
                CImGui.Separator()
                for i = 1:length(names)
                    CImGui.Selectable(names[i]) && (selected_fish = i;)
                end
                CImGui.EndPopup()
            end

            # showing a menu with toggles
            CImGui.Button("Toggle..") && CImGui.OpenPopup("my_toggle_popup")
            if CImGui.BeginPopup("my_toggle_popup")
                for i = 1:length(names)
                    CImGui.MenuItem(names[i], "", pointer(toggles)+(i-1)*sizeof(Bool))
                end
                if CImGui.BeginMenu("Sub-menu")
                    CImGui.MenuItem("Click me") && @info "Trigger `Click me` | find me here: $(@__FILE__) at line $(@__LINE__)"
                    CImGui.EndMenu()
                end

                CImGui.Separator()
                CImGui.Text("Tooltip here")
                CImGui.IsItemHovered() && CImGui.SetTooltip("I am a tooltip over a popup")

                CImGui.Button("Stacked Popup") && CImGui.OpenPopup("another popup")
                if CImGui.BeginPopup("another popup")
                    for i = 1:length(names)
                        CImGui.MenuItem(names[i], "", pointer(toggles)+(i-1)*sizeof(Bool))
                    end
                    if CImGui.BeginMenu("Sub-menu")
                        CImGui.MenuItem("Click me") && @info "Trigger `Click me` | find me here: $(@__FILE__) at line $(@__LINE__)"
                        CImGui.EndMenu()
                    end
                    CImGui.EndPopup()
                end
                CImGui.EndPopup()
            end
        end # @cstatic
        # call the more complete ShowExampleMenuFile which we use in various places of this demo
        CImGui.Button("File Menu..") && CImGui.OpenPopup("my_file_popup")
        if CImGui.BeginPopup("my_file_popup")
            ShowExampleMenuFile()
            CImGui.EndPopup()
        end

        CImGui.TreePop()
    end

    if CImGui.TreeNode("Context menus")
        # BeginPopupContextItem() is a helper to provide common/simple popup behavior of essentially doing:
        #    if (IsItemHovered() && IsMouseReleased(0))
        #       OpenPopup(id);
        #    return BeginPopup(id);
        # For more advanced uses you may want to replicate and cuztomize this code. This the comments inside BeginPopupContextItem() implementation.
        @cstatic value = Cfloat(0.5) begin
            CImGui.Text(@sprintf("Value = %.3f (<-- right-click here)", value))
            if CImGui.BeginPopupContextItem("item context menu")
                CImGui.Selectable("Set to zero") && (value = 0.0;)
                CImGui.Selectable("Set to PI") && (value = 3.1415;)
                CImGui.PushItemWidth(-1)
                @c CImGui.DragFloat("##Value", &value, 0.1, 0.0, 0.0)
                CImGui.PopItemWidth()
                CImGui.EndPopup()
            end
        end

        # We can also use OpenPopupOnItemClick() which is the same as BeginPopupContextItem() but without the Begin call.
        # So here we will make it that clicking on the text field with the right mouse button (1) will toggle the visibility of the popup above.
        CImGui.Text("(You can also right-click me to the same popup as above.)")
        CImGui.OpenPopupOnItemClick("item context menu", 1)

        # When used after an item that has an ID (here the Button), we can skip providing an ID to BeginPopupContextItem().
        # BeginPopupContextItem() will use the last item ID as the popup ID.
        # In addition here, we want to include your editable label inside the button label. We use the ### operator to override the ID (read FAQ about ID for details)
        @cstatic name="Label1"*"\0"^26 begin
            buf = @sprintf("Button: %s###Button", strip(name, '\0')) # `###` operator override ID ignoring the preceding label
            CImGui.Button(buf)
            if CImGui.BeginPopupContextItem()
                CImGui.Text("Edit name:")
                CImGui.InputText("##edit", name, length(name))
                CImGui.Button("Close") && CImGui.CloseCurrentPopup()
                CImGui.EndPopup()
            end
        end
        CImGui.SameLine()
        CImGui.Text("(<-- right-click here)")

        CImGui.TreePop()
    end

    if CImGui.TreeNode("Modals")
        CImGui.TextWrapped("Modal windows are like popups but the user cannot close them by clicking outside the window.")

        CImGui.Button("Delete..") && CImGui.OpenPopup("Delete?")

        if CImGui.BeginPopupModal("Delete?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
            CImGui.Text("All those beautiful files will be deleted.\nThis operation cannot be undone!\n\n")
            CImGui.Separator()

            # @cstatic dummy_i=Cint(0) @c CImGui.Combo("Combo", &dummy_i, "Delete\0Delete harder\0")

            @cstatic dont_ask_me_next_time=false begin
                CImGui.PushStyleVar(CImGui.ImGuiStyleVar_FramePadding, (0, 0))
                @c CImGui.Checkbox("Don't ask me next time", &dont_ask_me_next_time)
                CImGui.PopStyleVar()
            end

            CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
            CImGui.SetItemDefaultFocus()
            CImGui.SameLine()
            CImGui.Button("Cancel",(120, 0)) && CImGui.CloseCurrentPopup()
            CImGui.EndPopup()
        end

        CImGui.Button("Stacked modals..") && CImGui.OpenPopup("Stacked 1")
        if CImGui.BeginPopupModal("Stacked 1", C_NULL, CImGui.ImGuiWindowFlags_MenuBar)
            if CImGui.BeginMenuBar()
                if CImGui.BeginMenu("File")
                    CImGui.MenuItem("Dummy menu item") && @info "Trigger `Dummy menu item` | find me here: $(@__FILE__) at line $(@__LINE__)"
                    CImGui.EndMenu()
                end
                CImGui.EndMenuBar()
            end
            CImGui.Text("Hello from Stacked The First\nUsing style.Colors[ImGuiCol_ModalWindowDimBg] behind it.");

            # testing behavior of widgets stacking their own regular popups over the modal.
            @cstatic item=Cint(1) color=Cfloat[0.4,0.7,0.0,0.5] begin
                @c CImGui.Combo("Combo", &item, "aaaa\0bbbb\0cccc\0dddd\0eeee\0\0")
                CImGui.ColorEdit4("color", color)
            end
            CImGui.Button("Add another modal..") && CImGui.OpenPopup("Stacked 2")

            # Also demonstrate passing a bool* to BeginPopupModal(), this will create a regular close button which will close the popup.
            # Note that the visibility state of popups is owned by imgui, so the input value of the bool actually doesn't matter here.
            dummy_open = true
            if @c CImGui.BeginPopupModal("Stacked 2", &dummy_open)
                CImGui.Text("Hello from Stacked The Second!")
                CImGui.Button("Close") && CImGui.CloseCurrentPopup()
                CImGui.EndPopup()
            end

            CImGui.Button("Close") && CImGui.CloseCurrentPopup()
            CImGui.EndPopup()
        end

        CImGui.TreePop()
    end

    if CImGui.TreeNode("Menus inside a regular window")
        CImGui.TextWrapped("Below we are testing adding menu items to a regular window. It's rather unusual but should work!")
        CImGui.Separator()
        # NB: As a quirk in this very specific example, we want to differentiate the parent of this menu from the parent of the various popup menus above.
        # To do so we are encloding the items in a PushID()/PopID() block to make them two different menusets. If we don't, opening any popup above and hovering our menu here
        # would open it. This is because once a menu is active, we allow to switch to a sibling menu by just hovering on it, which is the desired behavior for regular menus.
        CImGui.PushID("foo")
        CImGui.MenuItem("Menu item", "CTRL+M")
        if CImGui.BeginMenu("Menu inside a regular window")
            ShowExampleMenuFile()
            CImGui.EndMenu()
        end
        CImGui.PopID()
        CImGui.Separator()
        CImGui.TreePop()
    end
end


text:
function setup_feasibility()
    options = Dict(:verbose => false,
                   :eps_abs => 1e-06,
                   :eps_rel => 1e-06,
                   :max_iter => 5000)
    return options
end

tol = 1e-3

@testset "feasibility" begin

    @testset "feasibility_problem" begin
        n = 30
        m = 30
        P = spzeros(n, n)
        q = zeros(n)
        A = sprandn(m, n, 0.8)
        u = randn(m)
        l = u
        options = setup_feasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test isapprox(norm(A * results.x - u), 0., atol=tol)

    end

end


text:
import JLD2, FileIO

# With the following standard FileIO.save, FileIO.load, JLD2.@save, JLD2.@load should work
struct JLD2KnetArray{T,N}; array::Array{T,N}; end
JLD2.writeas(::Type{KnetArray{T,N}}) where {T,N} = JLD2KnetArray{T,N}
JLD2.wconvert(::Type{JLD2KnetArray{T,N}}, x::KnetArray{T,N}) where {T,N} = JLD2KnetArray(Array(x))
JLD2.rconvert(::Type{KnetArray{T,N}}, x::JLD2KnetArray{T,N}) where {T,N} = KnetArray(x.array)

# For people who use Knet.load, Knet.@save etc.
save(x...; o...) = FileIO.save(x...; o...)
load(x...; o...) = FileIO.load(x...; o...)
macro save(x...); a=:(JLD2.@save); append!(a.args,x); esc(a); end
macro load(x...); a=:(JLD2.@load); append!(a.args,x); esc(a); end

# To load a file saved with Knet-1.4.3 or earlier:
load143(x...; o...)=jld2serialize(FileIO.load(x...; o...))


text:
module LinkedLists

import Base: isempty, empty!, length, first, last, iterate, isdone
import Base: eltype, pushfirst!, popfirst!, deleteat!
import Base: show, println, indexin, lastindex, push!, pop!, insert!, splice!, keys
import Base: in , append!, prepend!
import Base: getindex, setindex!

abstract type AbstractList{T} end
abstract type AbstractNode{T} end

include("list.jl")

export LinkedList, ListNode

export AbstractList, AbstractNode

export SLinkedList, SListNode, indextoposition, positiontoindex

end # module


text:
using RayTracer
using Test

@testset "creating_a_canvas" begin
    c = RayTracer.canvas(10,20)
    @test c.width == 10
    @test c.height == 20
    for x=1:c.width, y=1:c.height
        @test c.grid[x,y] ≈ RayTracer.color(0,0,0)
    end
end

@testset "writing_pixels_to_canvas" begin
    c = RayTracer.canvas(10,20)
    red = RayTracer.color(1,0,0)
    RayTracer.write_pixel!(c,2,3,red)
    @test RayTracer.pixel_at(c,2,3) ≈ red
end

@testset "constructing_ppm_header" begin
    c = RayTracer.canvas(5,3)
    ppm = RayTracer.canvas_to_ppm(c)
    f = open(ppm, "r")
    header = readlines(f)[1:3]
    close(f)
    @test header == ["P3","5 3","255"]
end

@testset "constructing_ppm_pixel_data" begin
    c = RayTracer.canvas(5,3)
    c1 = RayTracer.color(1.5,0,0)
    c2 = RayTracer.color(0,0.5,0)
    c3 = RayTracer.color(-0.5,0,1)
    RayTracer.write_pixel!(c,1,1,c1)
    RayTracer.write_pixel!(c,3,2,c2)
    RayTracer.write_pixel!(c,5,3,c3)
    ppm = RayTracer.canvas_to_ppm(c)
    f = open(ppm, "r")
    body = readlines(f)[4:6]
    @test body == [
        "255 0 0 0 0 0 0 0 0 0 0 0 0 0 0",
        "0 0 0 0 0 0 0 128 0 0 0 0 0 0 0",
        "0 0 0 0 0 0 0 0 0 0 0 0 0 0 255"
    ]
end

@testset "constructing_ppm_pixel_data_large" begin
    c = RayTracer.canvas(10,2)
    c1 = RayTracer.color(1,0.8,0.6)
    for height in 1:c.height
        for width in 1:c.width
            RayTracer.write_pixel!(c,width,height,c1)
        end
    end
    ppm = RayTracer.canvas_to_ppm(c)
    f = open(ppm, "r")
    body = readlines(f)[4:7]
    @test body == [
        "255 204 153 255 204 153 255 204 153 255 204 153 255 204 153 255 204",
        "153 255 204 153 255 204 153 255 204 153 255 204 153",
        "255 204 153 255 204 153 255 204 153 255 204 153 255 204 153 255 204",
        "153 255 204 153 255 204 153 255 204 153 255 204 153"
    ]
end

@testset "ppm_files_ending_newline" begin
    c = RayTracer.canvas(5,3)
    ppm = RayTracer.canvas_to_ppm(c)
    f = open(ppm, "r")
    endchar = read(f, String)[end]
    @test endchar == '\n'
end


text:
using BinDeps
ver = "1.0.0"
@BinDeps.setup
libyepppdep = library_dependency("libyeppp", aliases=["yeppp"])
if Sys.ARCH != :x86_64
    warn("Your package is not built successfully. You can still use this package if downding Yeppp! and adding the binary library file corresponding to your platform to the search path.")
else

    provides(Binaries, URI("http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-$ver.tar.bz2"), libyepppdep, unpacked_dir = "yeppp-$ver/binaries/windows/amd64/", os = :Windows)
    provides(Binaries, URI("http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-$ver.tar.bz2"), libyepppdep, unpacked_dir = "yeppp-$ver/binaries/macosx/x86_64/", os = :Darwin)
    @static Sys.islinux() ? push!(BinDeps.defaults, Binaries) : nothing
    provides(Binaries, URI("http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-$ver.tar.bz2"), libyepppdep, unpacked_dir = "yeppp-$ver/binaries/linux/x86_64/", os = :Linux)
end
@BinDeps.install Dict(:libyeppp => :libyeppp)
@static Sys.islinux() ? pop!(BinDeps.defaults) : nothing


text:
# https://adventofcode.com/2020/day/3

# You start on the open square (.) in the top-left corner and need to reach the
# bottom (below the bottom-most row on your map).

# The toboggan can only follow a few specific slopes (you opted for a cheaper model
# that prefers rational numbers); start by counting all the trees you would encounter
# for the slope right 3, down 1:

# From your starting position at the top-left, check the position that is right 3
# and down 1. Then, check the position that is right 3 and down 1 from there, and
# so on until you go past the bottom of the map.

function part1(patch, right, down=1)
    trees_count = 0
    position = 1
    time_to_move = 1

    for row in split(patch, '\n', keepempty=false)
        if time_to_move == 1
            if row[position] == '#'
                trees_count += 1
            end
            position = mod1(position + right, length(row))
            time_to_move = down
        else
            time_to_move -= 1
        end
    end
    return trees_count
end

example = "
..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#
"

@assert part1(example, 3) == 7

# Determine the number of trees you would encounter if, for each of the following
# slopes, you start at the top-left corner and traverse the map all the way
# to the bottom:
#
#  Right 1, down 1.
#  Right 3, down 1. (This is the slope you already checked.)
#  Right 5, down 1.
#  Right 7, down 1.
#  Right 1, down 2.
#
# In the above example, these slopes would find 2, 7, 3, 4, and 2 tree(s) respectively;
# multiplied together, these produce the answer 336.

function part2(patch, right, down)
    result = 1
    for (r, d) in zip(right, down)
        result *= part1(patch, r, d)
    end
    return result
end

@assert part2(example, [1 3 5 7 1], [1 1 1 1 2]) == 336

test = read("data/day-03.txt", String)
println("Part 1: $(result1 = part1(test, 3))")
println("Part 2: $(result2 = part2(test, [1 3 5 7 1], [1 1 1 1 2]))")

@assert result1 == 203
@assert result2 == 3316272960

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
function widen_invokes!(b::CodeInfoTools.Builder)
    counter = 0
    for (v, st) in b
        st isa Expr || continue
        st.head == :call || continue
        st.args[1] == invoke || continue
        insert!(b, v + counter, Expr(:tuple, st.args[2], st.args[4]...))
        insert!(b, v + 1 + counter, Expr(:call, typeof, Core.SSAValue(v)))
        replace!(b, v + 2 + counter, Expr(:call, invoke,
                                          GlobalRef(Mixtape, :call),
                                          Core.SSAValue(v + 1),
                                          st.args[2], st.args[4]...))
        counter += 2
    end
    return b
end


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

#const DEFAULT_DATASET_PATH = "../../../datasets"

println("--> Dataset 1")
#config = "test_dataset1_config.txt"
config = ARGS[1]
generations = parse(Int, ARGS[2])
st1 = run_NEAT_epochs(config, generations, false)
print_best_solution_perfomance(st1)
print("* Calculating information plane points")
calculate_information_plane_points(st1)

#=
# st = get_simulation_state("checkpoint-dt1-limited-10000g.jld2")
# st = get_simulation_state("Resultados_23:07:2021/not-limited/checkpoint-dt1-not-limited-10000g.jld2")
# st = get_simulation_state("Resultados_23:07:2021/limited/checkpoint-dt1-limited-10000g.jld2")
plot_neat_fitness(st1, imagename="fitness-dt1.png")
plot_neat_IP(st1, imagename="IP-dt1.png", nbins=10)
plot_neat_species(st1, imagename="species-dt1.png")
plot_output_distribution(st1, imagename="distribution-dt1.png")
plot_complexity(st1, imagename="solution-complexity-dt1.png")
plot_accuracy(st1, imagename="accuracy-dt1.png")
plot_total_accuracy(st1, imagename="total-accuracy-dt1.png")
plot_output_distribution_test(st1, imagename="output-test-distribution-best-NEAT.png")
plot_confusion_matrix(st1; imagename="confusion-matrix-testset-dt1.png", set=:test)
=#
#pop = getPopulation("checkpoint-dt1.jld2")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
module DatasetFunctions
using Statistics
using Dates
using DataFrames
import ..MasslistFunctions

		export nansum, nanmean, nanmax, nanmin, masknans, matlab2datetime, matlab2unixtime, calculateMeanTraces

	function nansum(x::AbstractArray)
		if length(filter(!isnan, x)) > 0
			return Statistics.sum(filter(!isnan, x))
		elseif length(filter(!isnan, x)) == 0
			return NaN
		end
	end

	function nansum(x::AbstractArray,d)
		return mapslices(nansum,x;dims = d)
	end

	function nanmean(x::AbstractArray)
		if length(filter(!isnan, x)) > 0
			return Statistics.mean(filter(!isnan,x))
		elseif length(filter(!isnan, x)) == 0
			return NaN
		end
	end

	function nanmean(x::AbstractArray,d)
		return mapslices(nanmean,x;dims = d)
	end

	function nanmax(x::AbstractArray)
		if length(filter(!isnan, x)) > 0
			return maximum(filter(!isnan, x))
		elseif length(filter(!isnan, x)) == 0
			return NaN
		end
	end

	function nanmin(x::AbstractArray)
		if length(filter(!isnan, x)) > 0
			return minimum(filter(!isnan, x))
		elseif length(filter(!isnan, x)) == 0
			return NaN
		end
	end

	function nanmax(x::AbstractArray, d)
		return mapslices(nanmax,x; dims = d)
	end

	function nanmin(x::AbstractArray, d)
		return mapslices(nanmin,x; dims = d)
	end

	function matlab2datetime(X :: Array{Float64})
		dtarray = Dates.Day.(Int64.(floor.(X.-719529)))+Dates.DateTime(1970,1,1)+Dates.Second.(Int64.(round.((X.-floor.(X)).*(24*60*60))))
		return dtarray
	end

	function matlab2unixtime(X :: Array{Float64})
		utarray = (X.-719529).*(24*60*60)
		return utarray
	end

	function masknans(X :: AbstractArray)
		return findall(x -> isnan(x) != true, X)
	end

	function timeRange(st, et, X :: Array{DateTime})
		return findall(st .< X .< et)
	end

	function calculateMeanTraces(stagestimes::Array{DateTime,1}, measResult::Main.ResultFileFunctions.MeasurementResult)
		meanTraces = zeros(1, length(measResult.MasslistMasses))
		for i in range(1,stop = (length(stagestimes)-1))
			if (stagestimes[i+1]-Minute(10)) > stagestimes[i]
				stage = stagestimes[i+1]-Minute(10) .< measResult.Times .< stagestimes[i+1]
			else
				stage = stagestimes[i] .< measResult.Times .< stagestimes[i+1]
			end
			a = measResult.Traces[stage,:]
			if (size(a)[1] == 0)
				meanTraces = cat(dims=1, meanTraces, zeros(1,length(measResult.MasslistMasses)))
			else
				meanTraces = cat(dims=1, meanTraces, mean(a; dims=1))
			end
		end
		laststage = stagestimes[end] .< measResult.Times
		meanTraces = cat(dims=1, meanTraces, mean(measResult.Traces[laststage,:]; dims=1))
		meanTraces = meanTraces[2:end, :]
		legStrings = []
		for i = 1:length(measResult.MasslistMasses)
		  push!(legStrings,"$(MasslistFunctions.sumFormulaStringFromCompositionArray(measResult.MasslistCompositions[:,i]))")
		end
		meanTraces = DataFrame(meanTraces)
		rename!(meanTraces, [f => t for (f,t) = zip(names(meanTraces), legStrings)])
		return meanTraces
	end

end

#=

static void invertBooleanArray(boolean[][] arr) {
           for (int i = 0; i < arr.length; i++)
               for (int j = 0; j < arr[0].length; j++)
                   arr[i][j] = !arr[i][j];
       }

static void invertBooleanArray(boolean[] arr) {
           for (int i = 0; i < arr.length; i++)
                   arr[i] = !arr[i];
       }

static void invertBooleanArray(boolean[][][] arr) {
           for (int i = 0; i < arr.length; i++)
               for (int j = 0; j < arr[0].length; j++)
               	for (int k = 0; k < arr[0].length; k++)
                   arr[i][j][k] = !arr[i][j][k];
       }
=#


text:
@testset "Event Difference" begin
    testEvent = 3
    testEventTimes = [1, 2, 3]
    @testset "Individual" begin
        inferTests && @inferred HawkesProcesses.event_diff(testEvent, testEventTimes)
        @test HawkesProcesses.event_diff(3, testEventTimes) == [2, 1]
        @test HawkesProcesses.event_diff(2, testEventTimes) == [1]
        @test HawkesProcesses.event_diff(1, testEventTimes) == []
    end
    @testset "Full Array" begin
        testEventTimes = [1,2,3]
        @test HawkesProcesses.event_difference_list(testEventTimes) == [[], [1], [2, 1]]
    end
end


text:
module unit_cube_tet_examples
using FinEtools
using FinEtools.MeshExportModule
using LinearAlgebra
using Arpack

function unit_cube_esnice()
    println("""
    Vibration modes of unit cube  of almost incompressible material.

    Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    tetrahedral. International Journal for Numerical Methods in
    Engineering 67: 841-867.
    """)
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 16;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;

    MR = DeforModelRed3D
    fens,fes  = T4block(a,b,h, na,nb,nh)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    mode = 17
    scattersysvec!(u, v[:,mode])
    File =  "unit_cube_esnice.vtk"
    vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
    @async run(`"paraview.exe" $File`)
    true

end # unit_cube_esnice

function allrun()
    println("#####################################################")
    println("# unit_cube_esnice ")
    unit_cube_esnice()
    return true
end # function allrun

end # module unit_cube_tet_examples


text:
char(x) = convert(Char, x)
char(x::FloatingPoint) = char(iround(x))

integer(x::Char) = int(x)

## char promotions ##

promote_rule(::Type{Char}, ::Type{Int8})    = Int32
promote_rule(::Type{Char}, ::Type{Uint8})   = Int32
promote_rule(::Type{Char}, ::Type{Int16})   = Int32
promote_rule(::Type{Char}, ::Type{Uint16})  = Int32
promote_rule(::Type{Char}, ::Type{Int32})   = Int32
promote_rule(::Type{Char}, ::Type{Uint32})  = Uint32
promote_rule(::Type{Char}, ::Type{Int64})   = Int64
promote_rule(::Type{Char}, ::Type{Uint64})  = Uint64
promote_rule(::Type{Char}, ::Type{Int128})  = Int128
promote_rule(::Type{Char}, ::Type{Uint128}) = Uint128

## character operations & comparisons ##

# numeric operations
# TODO: this should be removed, but needs to be here as long as Char <: Integer
+(x::Char   , y::Char   ) = int(x)+int(y)

# ordinal operations
+(x::Char   , y::Integer) = char(int(x)+int(y))
+(x::Integer, y::Char   ) = y+x
-(x::Char   , y::Char   ) = int(x)-int(y)
-(x::Char   , y::Integer) = char(int(x)-int(y))

# bitwise operations
(~)(x::Char) = char(~uint32(x))
(&)(x::Char, y::Char) = char(uint32(x) & uint32(y))
(|)(x::Char, y::Char) = char(uint32(x) | uint32(y))
($)(x::Char, y::Char) = char(uint32(x) $ uint32(y))

bswap(x::Char) = char(bswap(uint32(x)))

<<(x::Char, y::Int32)  = uint32(x) << y
>>(x::Char, y::Int32)  = uint32(x) >>> y
>>>(x::Char, y::Int32) = uint32(x) >>> y

< (x::Char, y::Char) = uint32(x) <  uint32(y)
<=(x::Char, y::Char) = uint32(x) <= uint32(y)

## traits ##

sizeof(::Type{Char}) = 4

## printing & showing characters ##

print(io::IO, c::Char) = (write(io,c); nothing)
show(io::IO, c::Char) = (print(io,'\''); print_escaped(io,UTF32String(c),"'"); print(io,'\''))

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
typealias SubPaving{N,T} Vector{IntervalBox{N,T}}

type Paving{N,T}
    separator::Separator
    inner::SubPaving{N,T}
    boundary::SubPaving{N,T}
    ϵ::Float64
end

function setdiff{N,T}(x::IntervalBox{N,T}, subpaving::SubPaving{N,T})
    working = [x]
    new_working = IntervalBox{N,T}[]

    local have_split

    for y in subpaving
        for x in working
            have_split = false

            diff = setdiff(x, y)

            if diff != [x]
                have_split = true
                append!(new_working, diff)
            end
        end

        !have_split && push!(new_working, x)

        working = new_working
        new_working = IntervalBox{N,T}[]

    end

    return working

end

setdiff{N,T}(X::SubPaving{N,T}, Y::SubPaving{N,T}) = vcat([setdiff(x, Y) for x in X]...)

function setdiff{N,T}(x::IntervalBox{N,T}, paving::Paving{N,T})
    Y = setdiff(x, paving.inner)
    Z = setdiff(Y, paving.boundary)
    return Z
end

function show{N,T}(io::IO, p::Paving{N,T})
    print(io, """Paving:
                 - tolerance ϵ = $(p.ϵ)
                 - inner approx. of length $(length(p.inner))
                 - boundary approx. of length $(length(p.boundary))"""
              )
end


text:
"""
Calculate ethane critical properties
T --- temperature in Kelvin
p --- pressure in bar
ro --- density in mol/l
"""
function LJ2Critical()
	sub=newCritical()
	sub.T=1.32
	sub.p=0.13006	#MPa to bar
	sub.ro=0.31	#mol/l #322 kg/m^3
	return sub
end

"""
Calculate ethane properties at triple point
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function LJ2Triple()
	sub=newTriple()
	sub.T=90.368
	sub.p=1.14*10^-5	#Pa to bar
#	sub.lRo=999.793/18.015268	#kg/m^3 to mol/l
#	sub.vRo=0.00485458/18.015268	#kg/m^3s to mol/l
	return sub
end

"""
calculate water properties at saturation line
Equation of State for the Lennard-Jones Fluid
Article in Journal of Physical and Chemical Reference Data · June 2016
DOI: 10.1063/1.4945000
T --- temperature
p --- pressure
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function ethaneSaturate(T)
a=[-6.48647577, 1.47010078, -1.66261122, 3.57898378, -4.79105705]
b=[1.56138026, -0.381552776, 0.0785372040, 0.0370315089]
c=[-1.89879145, -3.65459262, 0.850562745, 0.363965487, -1.50005943, -2.26690389]
sub=newSaturate()
crit=ethaneCritical()
if (T<crit.T)
	v=1.0-T/crit.T
	sub.p=crit.p*exp(crit.T/T*(a[1]*v+a[2]*v^1.5+a[3]*v^2.5+a[4]*v^3.5+a[5]*v^4.0))
	sub.lRo=crit.ro*exp(b[1]*v^0.329+b[2]*v^(4.0/6.0)+b[3]*v^(8.0/6.0)+b[4]*v^(19.0/6.0))
	sub.vRo=crit.ro*exp(crit.T/T*(c[1]*v^(0.346)+c[2]*v^(5.0/6.0)+c[3]*v+c[4]*v^2.0+c[5]*v^3.0+c[6]*v^5))
	sub.T=T
else
	sub.p=-1
	sub.lRo=-1
	sub.vRo=-1
	sub.T=-1
end
return sub
end

"""

T --- temperature, Kelvin
p --- pressure, bar
ro --- density, mol/l
z ---
s --- entropy,
s0 --- ideal part of entropy,
ds --- nonideal part derived k_B, dimensionless
g --- Gibbs energy
g0 --- ideal part of Gibbs energy
dg --- nonideal part of Gibbs energy
f ---  Helmholtz free energy
f0 --- ideal part
df --- nonideal part, dimensionless
u --- internal energy
du --- nonideal part
h --- entalpy
dh --- nonideal part
w --- speed of sound,
cp ---
cv ---
"""
function LJ2Singlephase(T,ro)

sub=newSingle()
crit=LJ2Critical()
delta=ro/crit.ro
tau=crit.T/T
R=1 #8.3144598
mm=0.03006904

nOne=6
nTwo=6
nThree=11

c=[0.00, 0.00, 0.00, 0.00, 0.00,
   0.00, 1.00, 2.00, 1.00, 2.00,	#10
   2.00, 1.00, 0.00, 0.00, 0.00,
   0.00, 0.00, 0.00, 0.00, 0.00,	#20
   0.00, 0.00, 0.00]

d=[4.00, 1.00, 1.00, 2.00, 2.00,
   3.00, 5.00, 2.00, 2.00, 3.00,	#10
   1.00, 1.00, 1.00, 1.00, 2.00,
   3.00, 3.00, 2.00, 1.00, 2.00,	#20
   3.00, 1.00, 1.00]

t=[ 1.000, 0.320, 0.505, 0.672, 0.843,
	0.898, 1.294, 2.590, 1.786, 2.770,	#10
	1.786, 1.205, 2.830, 2.548, 4.650,
	1.385, 1.460, 1.351, 0.660, 1.496,	#20
	1.830, 1.616, 4.970]

n=[0.52080730e-2, 0.21862520e1,
    -0.21610160e1, 0.14527000e1,
    -0.20417920e1, 0.18695286e0,
    -0.90988445e-1, -0.49745610e0,
    0.10901431e0, -0.80055922e0,
    -0.56883900e0, -0.62086250e0,
    -0.14667177e1, 0.18914690e1,
    -0.13837010e0, -0.38696450e0,
    0.12657020e0, 0.60578100e0,
    0.11791890e1, -0.47732679e0,
    -0.99218575e1, -0.57479320e0,
    0.37729230e-2]

al=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 2.067, 1.522, 8.820,
	1.722, 0.679, 1.883, 3.925, 2.461,	#20
	28.20, 0.753, 0.820]

bt=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.625, 0.638, 3.910,
	0.156, 0.157, 0.153, 1.160, 1.730,	#20
	383.0, 0.112, 0.119]

gm=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.710, 0.860, 1.940,
	1.480, 1.490, 1.945, 3.020, 1.110,	#20
	1.170, 1.330, 0.240]

ep=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.2053, 0.4090, 0.6000,
	1.2030, 1.8290, 1.3970, 1.3900, 0.5390,	#20
	0.9340, 2.3690, 2.4300]

phidelta=0.0
phir=0.0
phitau=0.0
phidelta2=0.0
phitau2=0.0
dphideltadtau=0.0

i=0
for ni=1:nOne
	i=i+1
	phidelta=phidelta+n[i]*d[i]*delta^(d[i]-1.0)*tau^t[i]
	phir=phir+n[i]*delta^d[i]*tau^t[i]
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1)
	phidelta2=phidelta2+n[i]*d[i]*(d[i]-1.0)*delta^(d[i]-2)*tau^t[i]
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2.0)
	dphideltadtau=dphideltadtau+n[i]*d[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)
end
for ni=1:nTwo
	i=i+1
	phidelta=phidelta+n[i]*exp(-delta^c[i])*(delta^(d[i]-1.0)*tau^t[i]*(d[i]-c[i]*delta^c[i]))
	phir=phir+n[i]*delta^(d[i])*tau^t[i]*exp(-delta^c[i])
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1.0)*exp(-delta^c[i]);
	phidelta2=phidelta2+n[i]*exp(-delta^c[i])*(delta^(d[i]-2)*tau^(t[i])*((d[i]-c[i]*delta^c[i])*(d[i]-1.0-c[i]*delta^c[i])-c[i]^2*delta^c[i]))
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2)*exp(-delta^c[i])
	dphideltadtau=dphideltadtau+n[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)*(d[i]-c[i]*delta^c[i])*exp(-delta^c[i])
end
for pn=1:nThree
	i=i+1
	phidelta=phidelta+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta-2.0*al[i]*(delta-ep[i]))
	phir=phir+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)
	phitau=phitau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(t[i]/tau-2*bt[i]*(tau-gm[i]))
	phidelta2=phidelta2+n[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(-2*al[i]*delta^d[i]+4*al[i]^2*delta^d[i]*(delta-ep[i])^2-4*d[i]*al[i]*delta^(d[i]-1)*(delta-ep[i])+d[i]*(d[i]-1)*delta^(d[i]-2))
	phitau2=phitau2+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*((t[i]/tau-2*bt[i]*(tau-gm[i]))^2-t[i]/tau^2-2*bt[i])
	dphideltadtau=dphideltadtau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta -2*al[i]*(delta-ep[i]))*(t[i]/tau-2*bt[i]*(tau-gm[i]))
end

sub.T=T
sub.ro=ro
sub.p=(1.0+phidelta*delta)*ro*T	#bar
sub.z=1.0+phidelta*delta	#
#ideal part
n0=[6.262265814,
-1.515151515,
1.5,
0.0,
0.0,
0.0,
0.0,
0.0]
gamma0=[0.0, 0.0, 0.0, 1.4091052332, 4.0099170712, 6.5967098342, 13.9798102659,1.0]

#entropy
phi0=log(delta)+n0[1]+n0[2]*tau+n0[3]*log(tau)

phi0tau=n0[2]+n0[3]/tau
phi0delta=1/delta
phi0delta2=-1/delta^2
phi0tau2=-n0[3]/tau^2
phi0taudelta=0

#Entropy
sub.s=(tau*(phi0tau+phitau)-phi0-phir)*R
sub.s0=(tau*(phi0tau)-phi0)*R
sub.ds=(tau*(phitau)-phir)	#devidet to R

#out.entrfrac=(tau*(phi0tau+phitau)-phi0-phir)/(tau*(phi0tau)-phi0)

dtaudtau=-tau/T
dphidt=dtaudtau*dphideltadtau

#sub.dpdt=(delta*cur_temp*dphidt+delta*phidelta+1)*R*T/100;
#phi0
#phir
#phidelta
#phidelta*delta
#(1+phi0+phir+phidelta*delta)
sub.g=(1+phi0+phir+phidelta*delta)*T	#R*cur_temp*
sub.g0=(1+phi0)*T
sub.dg=(phir+phidelta*delta)	#deriver RT
#
sub.f=(phi0+phir)*T
sub.f0=phi0*T
sub.df=phir

#Energy
sub.u=tau*(phi0tau+phitau)*T
sub.du=tau*(phitau)
#
sub.h=(1+tau*(phi0tau+phitau)+delta*phidelta)*R*T
sub.dh=1+tau*(phitau)+delta*phidelta;
#speed of sound
sub.w=sqrt(R/mm*T*(1+2*delta*phidelta+delta^2*phidelta2-(1+delta*phidelta-delta*tau*dphideltadtau)^2/(tau^2*(phi0tau2+phitau2))))	#0.018 перевод в килограмы
#
sub.cv=-tau^2*(phi0tau2+phitau2)*R
sub.cp=(-tau^2*(phi0tau2+phitau2)+(1+delta*phidelta-delta*tau*dphideltadtau)^2/(1+2*delta*phidelta+delta^2*phidelta2))*R
sub.dpdt=R/100.0*T*(1.0+2*delta*phidelta+delta^2*phidelta2)
#debug
# println("phi0 $(phi0)")
# println("phi0delta $(phi0delta)")
# println("phi0delta2 $(phi0delta2)")
# println("phi0tau $(phi0tau)")
# println("phi0tau2 $(phi0tau2)")
# println("phi0taudelta $(phi0taudelta)")
#
# println("phir $(phir)")
# println("phidelta $(phidelta)")
# println("phidelta2 $(phidelta2)")
# println("phitau $(phitau)")
# println("phitau2 $(phitau2)")
# println("dphideltadtau $(dphideltadtau)")
#enddebug
return sub
end


text:
using Tensors, BenchmarkTools

abstract type AbstractMaterialStateOld end

mutable struct VariableStateOld <: AbstractMaterialStateOld
     stress::Vector{Float64}
    dstress::Vector{Float64}
    strain::Vector{Float64}
    dstrain::Vector{Float64}
    backstress1::Vector{Float64}
    dbackstress1::Vector{Float64}
    backstress2::Vector{Float64}
    dbackstress2::Vector{Float64}
    plastic_strain::Vector{Float64}
    dplastic_strain::Vector{Float64}
    cumeq::Float64
    dcumeq::Float64
    R::Float64
    dR::Float64
end

function update!(state::VariableStateOld)
    state.stress[:] .+= state.dstress[:]
    state.strain[:] .+= state.dstrain[:]
    state.backstress1[:] .+= state.dbackstress1[:]
    state.backstress2[:] .+= state.dbackstress2[:]
    state.plastic_strain[:] .+= state.dplastic_strain[:]
    state.cumeq += state.dcumeq
    state.R += state.dR
    fill!(state.dstress, 0.0)
    fill!(state.dstrain, 0.0)
    fill!(state.dbackstress1, 0.0)
    fill!(state.dbackstress2, 0.0)
    fill!(state.dplastic_strain, 0.0)
    state.dcumeq = 0.0
    state.dR = 0.0
end

function bench_chaboche_style_oldstate(N)
    stress = zeros(6)
    dstress = zeros(6)
    strain = zeros(6)
    dstrain = zeros(6)
    backstress1 = zeros(6)
    dbackstress1 = zeros(6)
    backstress2 = zeros(6)
    dbackstress2 = zeros(6)
    plastic_strain = zeros(6)
    dplastic_strain = zeros(6)
    cumeq = 0.0
    dcumeq = 0.0
    R = 0.0
    dR = 0.0
    state = VariableStateOld(stress, dstress, strain, dstrain,
        backstress1, dbackstress1, backstress2, dbackstress2,
        plastic_strain, dplastic_strain, cumeq, dcumeq, R, dR)

    for i in 1:N
        state.dstress[:] .= randn(6)
        state.dstrain[:] .= randn(6)
        state.dbackstress1[:] .= randn(6)
        state.dbackstress2[:] .= randn(6)
        state.dplastic_strain[:] .= randn(6)
        state.dcumeq = norm(dplastic_strain)
        state.dR = randn()
        update!(state) # Update paremeters
    end
end
N = 1000
println("Benchmark Chaboche Old State")
@btime bench_chaboche_style_oldstate(N)


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

JLLWrappers.@generate_main_file_header("bliss")
JLLWrappers.@generate_main_file("bliss", UUID("508c9074-7a14-5c94-9582-3d4bc1871065"))
end  # module bliss_jll


text:
module Project2

using Base.Filesystem
using Flux
using Flux: @epochs, params, onecold, onehotbatch, MLUtils.DataLoader, logitcrossentropy, train!, glorot_uniform
using Random
using StableRNGs
using WAV

include("load.jl")
include("perceptron.jl")
include("utils.jl")

end # module

max_stars_repo_name: kpmrozowski/wut-mini-deep-learning

max_issues_repo_name: kpmrozowski/wut-mini-deep-learning

max_forks_repo_name: kpmrozowski/wut-mini-deep-learning


text:
using LayoutPointers, ArrayInterface, Aqua, Test

@testset "LayoutPointers.jl" begin
  Aqua.test_all(LayoutPointers)

  println("Grouped Strided Pointers")
  @time @testset "Grouped Strided Pointers" begin
    M, K, N = 4, 5, 6
    A = Matrix{Float64}(undef, M, K);
    B = Matrix{Float64}(undef, K, N);
    C = Matrix{Float64}(undef, M, N);
    struct SizedWrapper{M,N,T,AT<:AbstractMatrix{T}} <: AbstractMatrix{T} ; A::AT; end
    SizedWrapper{M,N}(A::AT) where {M,N,T,AT<:AbstractMatrix{T}} = SizedWrapper{M,N,T,AT}(A)
    Base.size(::SizedWrapper{M,N}) where {M,N} = (M,N);
    Base.getindex(A::SizedWrapper, i...) = getindex(parent(A), i...)
    Base.parent(dw::SizedWrapper) = dw.A
    ArrayInterface.parent_type(::Type{SizedWrapper{M,N,T,AT}}) where {M,N,T,AT} = AT
    LayoutPointers.memory_reference(dw::SizedWrapper) = LayoutPointers.memory_reference(parent(dw))
    ArrayInterface.contiguous_axis(dw::SizedWrapper) = LayoutPointers.contiguous_axis(parent(dw))
    ArrayInterface.contiguous_batch_size(dw::SizedWrapper) = LayoutPointers.contiguous_batch_size(parent(dw))
    LayoutPointers.val_stride_rank(dw::SizedWrapper) = LayoutPointers.val_stride_rank(parent(dw))
    function ArrayInterface.strides(dw::SizedWrapper{M,N,T}) where {M,N,T}
      if LayoutPointers.val_stride_rank(dw) === Val((1,2))
        return LayoutPointers.One(), LayoutPointers.StaticInt{M}()
      else#if LayoutPointers.val_stride_rank(dw) === Val((2,1))
        return LayoutPointers.StaticInt{N}(), LayoutPointers.One()
      end
    end
    ArrayInterface.offsets(dw::SizedWrapper) = LayoutPointers.offsets(parent(dw))
    LayoutPointers.val_dense_dims(dw::SizedWrapper{T,N}) where {T,N} = LayoutPointers.val_dense_dims(parent(dw))

    GC.@preserve A B C begin
      fs = (false,true)#[identity, adjoint]
      for ai ∈ fs, bi ∈ fs, ci ∈ fs
        At = ai ? A : (similar(A')');
        Bt = bi ? B : (similar(B')');
        Ct = ci ? C : (similar(C')');
        spdw = LayoutPointers.DensePointerWrapper{(true,true)}(LayoutPointers.stridedpointer(At))
        gsp, pres = @inferred(LayoutPointers.grouped_strided_pointer((spdw,Bt,Ct), Val{(((1,1),(3,1)),((1,2),(2,1)),((2,2),(3,2)))}()))
        if ai === ci
          @test sizeof(gsp.strides) == 2sizeof(Int)
        end
        # Test to confirm that redundant strides are not stored in the grouped strided pointer
        @test sizeof(gsp) == sizeof(Int) * (6 - (ai & ci) - ((!ai) & bi) - ((!bi) & (!ci)))
        @test sizeof(gsp.offsets) == 0
        pA, pB, pC = @inferred(LayoutPointers.stridedpointers(gsp))
        @test pA === stridedpointer(At)
        @test pB === stridedpointer(Bt)
        @test pC === stridedpointer(Ct)
        Btsw = SizedWrapper{K,N}(Bt)
        gsp2, pres2 = @inferred(LayoutPointers.grouped_strided_pointer((At,Btsw,Ct), Val{(((1,1),(3,1)),((1,2),(2,1)),((2,2),(3,2)))}()));
        @test sizeof(gsp2) == sizeof(Int) * (5 - (ai & ci) - ((!ai) & bi) - ((!bi) & (!ci)))

        pA2, pB2, pC2 = @inferred(LayoutPointers.stridedpointers(gsp2))
        @test pointer(pA2) == pointer(At)
        @test pointer(pB2) == pointer(Bt)
        @test pointer(pC2) == pointer(Ct)
        @test strides(pA2) == strides(pA)
        @test strides(pB2) == strides(pB)
        @test strides(pC2) == strides(pC)
      end
    end

    data_in_large = Array{Float64}(undef, 4, 4, 4, 4, 1);
    data_in = view(data_in_large, :, 1, :, :, 1);
    tmp1= Array{Float64}(undef, 4, 4, 4);
    sp_data_in, sp_tmp1 = LayoutPointers.stridedpointers(LayoutPointers.grouped_strided_pointer((data_in,tmp1), Val((((1,1),(2,1)),)))[1])
    @test sp_data_in === stridedpointer(data_in)
    @test sp_tmp1 === stridedpointer(tmp1)
  end

end


text:
## Reactive signal update performance: o, h, l, c vectors from 2D array
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")
ohlc_valarr = ohlc_BA.values

s_ohlc = Reactive.Input(ohlc_valarr[1,:])
s_close = lift(s -> s[4], s_ohlc, typ=Float64)

vsig = [s_close.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in 2:nmax
        push!(s_ohlc, ohlc_valarr[i,:])
        push!(vsig, s_close.value)
    end
end

println("OHLC updates o, h, l, c vector from 2D array rows:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: 0.135502571 seconds (17479632 bytes allocated, 28.15% gc time)
### no advantage of pre-separating values-array from TimeArray (good news)

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()