text:
# Activation function layers: have single input, same size output, can overwrite.

abstract Actf <: Op

ninputs(::Actf)=1
canoverwrite(::Actf)=true
back_reads_x(::Actf)=false
back_reads_y(::Actf)=true

for (ltype,lforw,lback) in
    ((:Sigm, :sigmforw, :sigmback),
     (:Tanh, :tanhforw, :tanhback),
     (:Relu, :reluforw, :reluback),
     (:Soft, :softforw, :softback),
     (:Soft73, :soft73forw, :soft73back),
     (:Copy, :copyforw, :copyback))
    @eval begin
        type $ltype <: Actf; $ltype(;o...)=new(); end
        forw(l::$ltype, x, y; o...)=($lforw(x,y;o...);gpusync();y)
        back(l::$ltype, dy, dx; y=nothing, o...)=(dx==nothing&&return;$lback(y,dy,dx;o...);gpusync();dx) # need o... for mask etc.
    end
end

function infersize(::Actf,xdims,ydims)
    if xdims==ydims==nothing
        nothing
    elseif xdims==nothing
        (ydims, ydims)
    elseif ydims==nothing
        (xdims, xdims)
    else
        @assert length(xdims) == length(ydims)
        dims = map(xdims, ydims) do x,y
            x == y ? x :
            x == 0 ? y :
            y == 0 ? x :
            throw(DimensionMismatch())
        end
        (dims, dims)
    end
end

### Implementations

@doc "@knet function copy(x) copies its input to its output" :copy
back_reads_y(::Copy)=false
copyforw(x,y;o...)=(x===y ? y : copysync!(y,x))
copyback(y,dy,dx;o...)=(dx===dy ? dx : copysync!(dx,dy))

@doc "@knet function sigm(x) computes the sigmoid activation function: 1/(1+exp(-x))" :sigm
sigmforw(x::Array,y::Array;o...)=(@inbounds for i=1:length(y); y[i]=(1/(1+exp(-x[i]))); end)
sigmback(y::Array,dy::Array,dx::Array;o...)=(@inbounds for i=1:length(dx); dx[i]=dy[i]*y[i]*(1-y[i]); end)
@gpu sigmforw(x::CudaArray,y::CudaArray;o...)=cudnnActivationForward(x,y; mode=CUDNN_ACTIVATION_SIGMOID)
@gpu sigmback(y::CudaArray,dy::CudaArray,dx::CudaArray;o...)=cudnnActivationBackward(y, dy, y, dx; mode=CUDNN_ACTIVATION_SIGMOID)

@doc "@knet function tanh(x) computes the hyperbolic tangent activation function." :tanh
tanhforw(x::Array,y::Array;o...)=(@inbounds for i=1:length(y); y[i]=tanh(x[i]); end)
tanhback(y::Array,dy::Array,dx::Array;o...)=(@inbounds for i=1:length(dx); dx[i]=dy[i]*(1+y[i])*(1-y[i]); end)
@gpu tanhforw(x::CudaArray,y::CudaArray;o...)=cudnnActivationForward(x,y; mode=CUDNN_ACTIVATION_TANH)
@gpu tanhback(y::CudaArray,dy::CudaArray,dx::CudaArray;o...)=cudnnActivationBackward(y, dy, y, dx; mode=CUDNN_ACTIVATION_TANH)

@doc "@knet function relu(x) computes the rectified linear activation function: (x<0 ? 0 : x)" :relu
# This is too slow, use broadcast:  TODO: find out why. same is not true for tanh.
# reluforw(x::Array,y::Array;o...)=(@inbounds for i=1:length(y); y[i]=(x[i]<0 ? 0 : x[i]) end)
# reluforw(x::Array,y::Array;o...)=broadcast!(max,y,x,0)
reluforw(x::Array,y::Array;o...)=(@inbounds for i=1:length(y); if x[i]>0; y[i]=x[i]; else; y[i]=0; end; end)
# reluback(y::Array,dy::Array,dx::Array;o...)=(@inbounds for i=1:length(dx); dx[i]=(y[i]==0 ? 0 : dy[i]) end)
# reluback(y::Array,dy::Array,dx::Array;o...)=(copy!(dx,dy); dx[y.==0]=0)
reluback(y::Array,dy::Array,dx::Array;o...)=(@inbounds for i=1:length(dx); if y[i]==0; dx[i]=0; else; dx[i]=dy[i]; end; end)
@gpu reluforw(x::CudaArray,y::CudaArray;o...)=cudnnActivationForward(x,y; mode=CUDNN_ACTIVATION_RELU)
@gpu reluback(y::CudaArray,dy::CudaArray,dx::CudaArray;o...)=cudnnActivationBackward(y, dy, y, dx; mode=CUDNN_ACTIVATION_RELU)

@doc "@knet function soft(x) computes the softmax activation function: exp(x[i,j])/sum(exp(x[:,j]))" :soft
# z = wx			;; z is the input to the soft layer
# qi = (exp zi) / (Σ exp zj)	;; q is the output of the soft layer
# ∂qi/∂zk = [(i=k)(exp zi)(Σ exp zj) - (exp zi)(exp zk)] / (Σ exp zj)^2
#         = (i=k) qi - qi qk
# ∂J/∂zk = Σ (∂J/∂qi)(∂qi/∂zk)	;; derivative wrt the input z
#        = Σ (1-pi/qi)((i=k) qi - qi qk)
#        = Σ ((i=k) qi - qi qk - (i=k) pi + pi qk)
#        = qk - pk - qk Σ (qi - pi)
#        = qk - pk
#
# Note that softback expects xgrad from softloss, not ygrad.
# See the softloss doc for an explanation.
back_reads_y(::Soft)=false
softback(y,dy,dx;o...)=(dx===dy ? dx : copysync!(dx,dy))
@gpu softforw(x::CudaArray,y::CudaArray;o...)=cudnnSoftmaxForward(x,y)
function softforw(x::Array,y::Array;o...)
    (st,nx) = size2(x)
    @inbounds for j=1:nx
        i1=(j-1)*st+1
        i2=j*st
        xmax = typemin(eltype(x))
        ysum = zero(Float64)
        for i=i1:i2; x[i] > xmax && (xmax = x[i]); end
        for i=i1:i2; ysum += (y[i]=exp(x[i] - xmax)); end
        for i=i1:i2; y[i] /= ysum; end
    end
    return y
end

# The exception of expecting xgrad instead of ygrad in the soft layer
# is ok if it is the last layer and is used with softloss, which is
# aware of this exception.  However, if soft is to be used internally
# in a model this does not work, the layer will receive ygrad.  For
# now we will define a soft73 layer to prevent code break and
# eventually retire the soft layer in favor of models that do not have
# a final actf and that use what we used to call xentloss.

# For the math let us use x for the input to softmax, y for output, J
# for loss.  The softforw function computes:
# ;; yi = exp(xi)/Σj exp(xj)
# Softback should compute:
# ;; ∂yi/∂xk = [(i=k)(exp xi)(Σ exp xj) - (exp xi)(exp xk)] / (Σ exp xj)^2
# ;; = (i=k) yi - yi yk
# ;; ∂J/∂xk = Σi (∂J/∂yi)(∂yi/∂xk)
# ;; = Σi (∂J/∂yi)((i=k) yk - yi yk)
# ;; = yk ((∂J/∂yk) - Σi yi (∂J/∂yi))

soft73forw(x,y;o...)=softforw(x,y;o...) # forward does not change

@gpu soft73back(y::CudaArray,dy::CudaArray,dx::CudaArray;o...)=(cudnnSoftmaxBackward(y, dy, dx); gpusync(); dx)

function soft73back(y::Array,dy::Array,dx::Array;o...)
    (st,nx) = size2(dy)
    for j=1:nx
        i1=(j-1)*st+1
        i2=j*st
        sumydy = zero(Float64)
        for i=i1:i2; sumydy += y[i] * dy[i]; end
        for i=i1:i2; dx[i] = y[i] * (dy[i] - sumydy); end
    end
    return dx
end

@doc "@knet function axpb(x;a=1,p=1,b=0) computes y=ax^p+b elementwise." :axpb
type Axpb <: Actf; a; p; b; Axpb(;a=1,p=1,b=0,o...)=new(a,p,b); end
back_reads_x(f::Axpb)=(f.p!=1)
back_reads_y(f::Axpb)=false

function forw(f::Axpb, x, y; o...)
    length(x)==length(y) || throw(DimensionMismatch())
    axpbforw(f.a,x,f.p,f.b,y)
    return y
end

function axpbforw{T}(a::Number,x::Array{T},p::Number,b::Number,y::Array{T})
    a = convert(T,a); p = convert(T, p); b = convert(T, b)
    @inbounds for i=1:length(y)
        yi = x[i]
        p!=1 && (yi^=p)
        a!=1 && (yi*=a)
        b!=0 && (yi+=b)
        y[i] = yi
    end
end

@gpu function axpbforw{T}(a::Number,x::CudaArray{T},p::Number,b::Number,y::CudaArray{T})
    T <: Float32 ? ccall((:axpbforw32,libknet),Void,(Cint,Cfloat, Ptr{Cfloat}, Cfloat, Cfloat, Ptr{Cfloat}), length(x), convert(Cfloat,a), x, convert(Cfloat,p), convert(Cfloat,b), y) :
    T <: Float64 ? ccall((:axpbforw64,libknet),Void,(Cint,Cdouble,Ptr{Cdouble},Cdouble,Cdouble,Ptr{Cdouble}),length(x), convert(Cdouble,a),x, convert(Cdouble,p),convert(Cdouble,b),y) :
    error("$T not supported")
    gpusync()
end

function back(f::Axpb, dy, dx; x=nothing, o...)
    dx==nothing && return
    issimilar(dx,dy) || throw(DimensionMismatch())
    if f.p != 1
        issimilar(x,dy) || throw(DimensionMismatch())
        axpbback(f.a,x,f.p,dy,dx)
    elseif f.a != 1             # dJ/dx=dJ/dy*dy/dx=dJ/dy*a
        dx===dy || copysync!(dx,dy)
        scale!(f.a,dx)
    else
        dx===dy || copysync!(dx,dy)
    end
end

function axpbback{T}(a::Number, x::Array{T}, p::Number, dy::Array{T}, dx::Array{T})
    ap = convert(T,a*p)
    p1 = convert(T,p-1)
    @inbounds for i=1:length(dx)
        dx[i] = dy[i] * ap * x[i]^p1
    end
end

@gpu function axpbback{T}(a::Number, x::CudaArray{T}, p::Number, dy::CudaArray{T}, dx::CudaArray{T})
    T <: Float32 ? ccall((:axpbback32,libknet),Void,(Cint,Cfloat, Ptr{Cfloat}, Cfloat, Ptr{Cfloat}, Ptr{Cfloat}),  length(x), convert(Cfloat,a),  x, convert(Cfloat,p),  dy, dx) :
    T <: Float64 ? ccall((:axpbback64,libknet),Void,(Cint,Cdouble,Ptr{Cdouble},Cdouble,Ptr{Cdouble},Ptr{Cdouble}), length(x), convert(Cdouble,a), x, convert(Cdouble,p), dy, dx) :
    error("$T not supported")
    gpusync()
end

# Utility fn, used by rgen
axpb!(x; a=1, p=1, b=0)=(axpbforw(a,x,p,b,x); x)


text:
using PyThermo
using Test
using PyThermo.ShockTube: shockcalc
using Unitful
using Conda

@testset "PyThermo.jl" begin
    @test Conda.version("thermo") <= v"0.1.40"

    @testset "Species" begin
        SF6 = Species("SF6")
        @test isapprox(ustrip(density(SF6)), 6.0383, rtol=2e-3)
        SF6.calculate(T = 500)
        @test isapprox(ustrip(density(SF6)), 3.5657, rtol=2e-3)
    end
    @testset "Mixture" begin
        HeAce = Mixture(["Helium" => 0.95, "Acetone" => 0.05])
        ρ_HeAce = HeAce.rho
        @test !isnothing(ρ_HeAce)
        @test isapprox(ρ_HeAce, 0.2747138795604815, rtol=2e-3)
    end
end

@testset "ShockTube.jl" begin
    st = shockcalc(Species("N2"), Species("Ar"), 1.8)
    @test isapprox(ustrip(pressure(st.driver)), 1.426894182009389e6, rtol=2e-3)
    @test isapprox(ustrip(st.u2), 300.0825207636464, rtol=2e-3)
    @test isapprox(ustrip(density(st.shocked)), 3.3867027987446874, rtol=2e-3)
    @test isapprox(ustrip(soundspeed(st.reflected)), 544.5317089952606, rtol=5e-3)
end

ext: jl


text:
include("StaticAnalyzer.jl")

struct ExecutionException <: Exception
  msg::String
  line::Int
end

Base.showerror(io::IO, e::LexicalException) = print(io,
  "A lexical error occurred: $(e.msg)\n")
Base.showerror(io::IO, e::SyntaxException) = print(io,
  "A syntax error occurred: $(e.msg)\n")
Base.showerror(io::IO, e::StaticAnalysisException) = print(io,
  "Static analysis produced an error: $(e.msg)\n")
Base.showerror(io::IO, e::ExecutionException) = print(io,
  "Static analysis produced an error on line $(e.line): $(e.msg)\n")

function run(program::String)
  AST = parseInput(scanInput(program))
  staticAnalysis(AST)
  executeProgram(AST)
end

function evaluate(l::Literal)
  if l.token.class == int_literal
    value = parse(Int, l.token.lexeme)
    return SValue(MInt, value)
  elseif l.token.class == string_literal
    return SValue(MString, l.token.lexeme)
  elseif l.token.class == kw_true
    return SValue(MBool, true)
  elseif l.token.class == kw_false
    return SValue(MBool, false)
  end
end

evaluate(l::Literal, st) = evaluate(l)

function evaluate(v::Var, st::SymbolTable)
  return getValue(st, v.variable.lexeme)
end

function evaluate(node::UnaryOperation, st::SymbolTable)
  operation = operator_to_function[node.operator.class]
  return operation(evaluate(node.operand, st))
end

function evaluate(node::BinaryOperation, st::SymbolTable)
  operation = operator_to_function[node.operator.class]
  left = evaluate(node.leftOperand, st)
  right = evaluate(node.rightOperand, st)
  if (operation == (÷) && right.value == 0)
    throw(ExecutionException("Divide by zero", node.line))
  end
  return operation(left, right)
end

function executeStatements(statements::Array{Statement,1}, st::SymbolTable)
  for s in statements
    executeStatement(s, st)
  end
end

executeStatements(s::Statements, st::SymbolTable) =
  executeStatements(s.statements, st)

executeStatements(s::Array{Statement,1}) =
  executeStatements(s, SymbolTable())

executeProgram(p::Statements) = executeStatements(p.statements)

function executeStatement(d::Declaration, st::SymbolTable)
  addOrUpdateSymbol(st, d.variable.lexeme, SValue(d.type.class))
end

function executeStatement(d::DecAssignment, st::SymbolTable)
  addOrUpdateSymbol(st, d.variable.lexeme, evaluate(d.value, st))
end

function executeStatement(a::Assignment, st::SymbolTable)
  addOrUpdateSymbol(st, a.variable.lexeme, evaluate(a.value, st))
end

function executeStatement(f::For, st::SymbolTable)
  varName = f.variable.lexeme
  from = evaluate(f.from, st).value
  to = evaluate(f.to, st).value
  for i in from:to
    setIterationVariable(st, varName, i)
    executeStatements(f.body, st)
  end
  releaseIterationVariable(st, varName)
end

function executeStatement(p::Print, st::SymbolTable)
    print(evaluate(p.argument, st).value)
end

function executeStatement(r::Read, st::SymbolTable)
  varName = r.variable.lexeme
  type = getValue(st, varName).type
  println("\n<MiniPL-Julia is waiting for input>")
  rawInput = split(readline())[1]
  if type == MInt
    addOrUpdateSymbol(st, varName, SValue(MInt, parse(Int, rawInput)))
  elseif type == MString
    addOrUpdateSymbol(st, varName, SValue(MString, string(rawInput)))
  elseif type == MBool
    addOrUpdateSymbol(st, varName, SValue(MBool, parse(Bool, rawInput)))
  end
end

function executeStatement(a::Assert, st::SymbolTable)
  if !(evaluate(a.argument, st).value)
    println("Assertion on line $(a.line) failed.")
  end
end

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
"""
    abs_val(num)

Program to find the absolute value of a number

# Example
```julia
abs_val(-100) # returns 100
abs_val(0) # returns 0
abs(123.1) # returns 123.1
-1000 == abs_val(-1000) #returns false
1000 == abs_val(1000) #returns true
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee) and [Rratic](https://github.com/Rratic)
"""
function abs_val(num::N)where N<:Real
    return num < 0 ? -num : num
end

"""
    abs_max(x)

Program to find the max absolute value in a vector

# Example
```julia
abs_max([1,3,4]) # returns 4
abs_max([-3,1,2]) # returns -3
abs_max([-7,-3,6]) #returns -7
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee) and [Rratic](https://github.com/Rratic)
"""
function abs_max(x::Vector{N})where N<:Real
    _max=x[1]
    for i in x
        if abs(i)>abs(_max) _max=i end
    end
    return _max
end

"""
    abs_min(num)

Program to find the min absolute value in a vector

# Example
```julia
abs_min([1,3,4]) # returns 1
abs_min([-3,1,2]) # returns 1
abs_min([-7,-3,6]) #returns -3
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee) and [Rratic](https://github.com/Rratic)
"""
function abs_min(x::Vector{N})where N<:Real
    _min=x[1]
    for i in x
        if abs(i)<abs(_min) _min=i end
    end
    return _min
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
module Rainflow
import Base.show

export sort_peaks, find_boundary_vals, count_cycles, sum_cycles

""" This function sort out points where the slope is changing sign."""
function sort_peaks(signal::AbstractArray{Float64,1}, dt=[0.:length(signal)-1.])
    slope = diff(signal)
    # Determines if the point is local extremum
    is_extremum = [true, (slope[1:end-1].*slope[2:end]).<=0., true]
    return signal[is_extremum] , dt[is_extremum]
end

immutable Cycle  # This is the information stored for each cycle found
    count::Float64
    range::Float64
    mean::Float64  # value
    v_s::Float64   # value start
    t_s::Float64   # time start
    v_e::Float64   # value end
    t_e::Float64   # time end
end

show(io::IO,x::Cycle) = print(io, "Cycle: count=", x.count, ", range=",x.range, ", mean=",x.mean)

function cycle(count::Float64, v_s::Float64, t_s::Float64, v_e::Float64, t_e::Float64)
    Cycle(count, abs(v_s-v_e), (v_s+v_e)/2, v_s, t_s, v_e, t_e)
end

""" Count the cycles from the data. """
function count_cycles(ext_in::Array{Float64,1},t::Array{Float64,1})
    ext = copy(ext_in) # Makes a copy because there is going to be sorted in the vectors
    time = copy(t)
    i = 1
    j = 2
    cycles = Cycle[]
    sizehint(cycles, length(ext)) # This reduces the memory consumption a bit
    @inbounds begin
    while length(ext)>(i+1)
        Y = abs(ext[i+1]-ext[i])
        X = abs(ext[j+1]-ext[j])
        if X>=Y
            if i == 1 # This case counts a half and cycle deletes the poit that is counted
                #println("Half cycle $(ext[i]), $(ext[i+1])")
                push!(cycles,cycle(0.5 ,ext[i], time[i], ext[i+1],time[i+1]))
                shift!(ext) # Removes the first entrance in ext and time
                shift!(time)
            else # This case counts one cycle and deletes the poit that is counted
                #println("One cycle $(ext[i]), $(ext[i+1])")
                push!(cycles,cycle(1. ,ext[i], time[i], ext[i+1],time[i+1]))
                splice!(ext,i+1)  # Removes the i and i+1 entrance in ext and time
                splice!(ext,i)
                splice!(time,i+1)
                splice!(time,i)
            end
            i = 1
            j = 2
        else
            i += 1
            j += 1
        end
    end
    for i=1:length(ext)-1 # This counts the rest of the points that have not been counted as a half cycle
        #println("Half cycle $(ext[i]), $(ext[i+1])")
        push!(cycles,cycle(0.5 ,ext[i], time[i], ext[i+1],time[i+1]))
    end
    end
    return cycles
end

type Cycles_bounds #
    min_mean::Float64
    max_mean::Float64
    max_range::Float64
end

show(io::IO,x::Cycles_bounds) = print(io, "Cycles_bounds : min mean value=", x.min_mean, ", max mean value=", x.max_mean, ", max range=",x.max_range)

""" Find the minimum and maximum mean value and maximum range from a vector of cycles. """
function find_boundary_vals(cycles::Array{Cycle,1})
    bounds = Cycles_bounds(Inf, -Inf, -Inf)
    for cycle in cycles
        cycle.mean > bounds.max_mean && setfield!(bounds, :max_mean, cycle.mean)
        cycle.mean < bounds.min_mean && setfield!(bounds, :min_mean, cycle.mean)
        cycle.range > bounds.max_range && setfield!(bounds, :max_range, cycle.range)
    end
    return bounds
end

""" Returns the range index the value is belonging in """
function find_range{T<:Real}(interval::Array{T,1},value)
    issorted(interval) || error("The array needs to be sorted in raising order")
    for i=1:length(interval)-1
        if interval[i] <= value <= interval[i+1]
            return i
        end
    end
    error("The value where not in range")
end

if v"0.4.0-dev+4986" <  VERSION
    typealias Interval{T} Union{Array{T,1}, LinSpace{T}}
    """ Returns the range index the value is belonging in """
    function find_range{T<:Real}(interval::LinSpace{T}, value)
        issorted(interval) || error("The array needs to be sorted in raising order")
        start = interval.start
        stop = interval.stop
        (start < value < stop) || error("The value where not in range, see if the vectors in calc_sum(cycles::Array{Cycle,1}, range_intervals::Array{T,1}, mean_intervals::Array{T,1}) are continious increasing in value, or adjust the nr_digits parameter")
        inc = (interval.stop - start) / interval.divisor
        i = int(fld(value - start, inc) ) + 1
    end
else
    typealias Interval{T} Array{T,1}
end

""" Sums the cycle count given intervals of range_intervals and mean_intervals. The range_intervals and mean_intervals is given in fraction of range size"""
function sum_cycles{T<:Real}(cycles::Array{Cycle,1}, range_intervals::Interval{T}, mean_intervals::Interval{T})
    bounds = find_boundary_vals(cycles)
    bins = zeros(length(range_intervals)-1, length(mean_intervals)-1)
    range_intervals *= bounds.max_range/100
    #show(range_intervals)
    mean_intervals *= (bounds.max_mean-bounds.min_mean)/100
    mean_intervals += bounds.min_mean
    nr_digits = 14  # The rounding is performed due to numerical noise in the floats when comparing
    if v"0.4.0-dev+4986" >  VERSION || isa(mean_intervals,LinSpace)
        mean_intervals = round(mean_intervals, nr_digits)
    elseif v"0.4.0-dev+4986" >  VERSION || isa(range_intervals,LinSpace)
        range_intervals = round(range_intervals, nr_digits)
    end
    #show(mean_intervals)
    for cycle in cycles
        i = find_range(range_intervals,round(cycle.range, nr_digits))
        j = find_range(mean_intervals,round(cycle.mean, nr_digits))
        bins[i,j] += cycle.count
    end
    return bins
end

function sum_cycles(cycles::Array{Cycle,1}, nr_ranges::Int=10, nr_means::Int=1)
    range_intervals = linspace(0,100,nr_ranges+1)
    mean_intervals = linspace(0,100,nr_means+1)
    sum_cycles(cycles, range_intervals, mean_intervals)
end

try
    include("plot.jl")
catch
    println("""For added plotting features do: Pkg.add("PyPlot")""")
end

end # module


text:
using Random, LinearAlgebra

# Standardised tests for the eigen decomposition of a square kronecker product
function eigen_tests(rng, C::AbstractKroneckerProduct)

    # Check approximate correctness of decomposition
    λ, Γ = eigen(C)
    @test Γ * Diagonal(λ) * Γ' ≈ C

    # Check approximate correctness of shifted decomposition
    σ² = abs2(randn(rng))
    λ′, Γ′ = (eigen(C) + σ² * I)
    @test Γ′ * Diagonal(λ′) * Γ′' ≈ Matrix(C + σ² * I)
    @test σ² * I + C == C + σ² * I

    # Test with decomposition
    v = randn(rng, length(λ))
    @test eigen(C) \ v ≈ Float64.(Matrix(C)) \ v

    # Test various linear algebra operations with decomposition
    C_dense = Float64.(Matrix(C))
    @test det(eigen(C)) ≈ det(C_dense)
    @test logdet(eigen(C)) ≈ logdet(C_dense)
    @test Matrix(inv(eigen(C))) ≈ inv(C_dense)
end

@testset "eigen" begin
    rng, P, Q, R = MersenneTwister(123456), 3, 5, 7

    # Generate some positive definite matrices so that logdet can be tested.
    A_, B_, C_ = randn(rng, P, P), randn(rng, Q, Q), randn(rng, R, R)
    A, B, C = Symmetric(A_ * A_' + I), Symmetric(B_ * B_' + I), Symmetric(C_ * C_' + I)

    D = kronecker(A, B)
    eigen_tests(rng, D)
    eigen_tests(rng, kronecker(D, C))
    eigen_tests(rng, kronecker(C, D))
    eigen_tests(rng, kronecker(kronecker(A, B), kronecker(B, A)))
    eigen_tests(rng, kronecker(A, 4))
end


text:
import Base: eltype, size, length

"""
    Observation{D,T}

Types inheriting from this struct provide all information about the observation
of the stochastic process being made. `D` denotes the dimension of the
observation whereas `T` its eltype.
"""
abstract type Observation{D,T} end

"""
    ParametrizedObservation{D,T} <: Observation{D,T}

Same as observations but additionally, types inheriting from this struct depend
on some parameters (that may in practice be unknown).
"""
abstract type ParametrizedObservation{D,T} <: Observation{D,T} end

"""
    eltype(::Type{K}) where {K<:Observation{D,T}}

Type of each entry of the collection holding an observation.
"""
eltype(::Type{K}) where {K<:Observation{D,T}} where {D,T} = T

eltype(::K) where K<:Observation = eltype(K)

"""
    size(::Type{K}) where {K<:Observation{D}}

Size of the observation
"""
size(::Type{K}) where {K<:Observation{D}} where D = D

size(::K) where K<:Observation = size(K)

"""
    length(::Type{K}) where {K<:Observation{D}}

Length of the observation (equal to a number of entries in a vector holding the
data)
"""
length(::Type{K}) where {K<:Observation{D}} where D = _size_to_length(D)

length(::K) where K<:Observation = length(K)

_size_to_length(s::Number) = s
_size_to_length(s::Tuple) = prod(s)

#TODO check if can be removed, doesn't seem to be used
"""
    parameter_names

Return names of the parameters of a given object
"""
function parameter_names end

# TODO check if can be safely removed
function clone end

function set_parameters!(obs::Observation, η, entries)
    for (i_η, i_obs) in entries
        obs.θ[i_obs] = η[i_η]
    end
end


text:
using TensorOperations
using Printf

"""
replace this with table lookup
"""
function get_nchk(n,k)
    return calc_nchk(n,k)
end

function get_fiedler(A)

    dim = size(A)[1]
    dim > 1 || return [1.0]
    size(A)[1] == size(A)[2] || error(" a not square")
    D = zeros(size(A)[1])
    for i in 1:dim
        D[i] = sum(A[i,:])
    end
    L = Diagonal(D) - A
    F = eigen(L)

    @printf(" eigenvalues of laplacian\n")
    for (ii,i) in enumerate(F.values)
        @printf(" %4i = %12.8f\n", ii, i)
    end
    return F.vectors[:,2]
end

function fiedler_sort(C,K)
    fvec = get_fiedler(abs.(C'*K*C))
    perm, = bubble_sort(fvec)
    return C[:,perm]
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
abstract type GridFile end

mutable struct CurvilinearSphericalGridFile <: GridFile

    R    :: Float64
    Ω    :: Float64

    Nx    :: Integer
    Ny    :: Integer
    lsize :: Integer

    xc    :: Array{Float64, 2}
    yc    :: Array{Float64, 2}

    xv    :: Array{Float64, 3}
    yv    :: Array{Float64, 3}

    mask  :: Array{Float64, 2}
    area  :: Array{Float64, 2}
    frac  :: Array{Float64, 2}

    missing_value :: Float64

    function CurvilinearSphericalGridFile(
        filename::String;
        R       :: Float64,
        Ω       :: Float64,
        missing_value::Float64 = 1e20
    )

        ds = Dataset(filename, "r")
        _mask = ds["mask"][:, :]
        _area = ds["area"][:, :]
        _frac = ds["frac"][:, :]
        _xc  = ds["xc"][:, :]
        _yc  = ds["yc"][:, :]
        _xv  = ds["xv"][:]
        _yv  = ds["yv"][:]

        _Nx  = ds.dim["ni"]
        _Ny  = ds.dim["nj"]
        close(ds)

        return new(
            R, Ω,
            _Nx, _Ny, _Nx * _Ny,
            _xc, _yc, _xv, _yv,
            _mask, _area, _frac,
            missing_value,
        )

    end

end

mutable struct CylindricalGridFile <: GridFile

    R     :: Float64
    Ω     :: Float64

    Nx    :: Integer
    Ny    :: Integer

    Ly    :: Float64
    lat0  :: Float64
    β     :: Float64

    mask  :: Array{Float64, 2}

    function CylindricalGridFile(;
        R :: Float64,
        Ω :: Float64,
        Nx :: Integer,
        Ny :: Integer,
        Ly :: Float64,
        lat0 :: Float64,
        β    :: Float64,
        mask :: Union{AbstractArray{Float64, 2}, Nothing} = nothing,
    )

        if mask == nothing
            mask = ones(Float64, Nx, Ny)
        end

        if size(mask) != (Nx, Ny)
            throw(ErrorException("Size of mask does not match Nx and Ny"))
        end

        return new(
            R, Ω,
            Nx, Ny, Ly,
            lat0, β,
            mask,
        )

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Dates

struct Transaction
  sender::String
  recipient::String
  amount::Int
end

struct Block
  index::Int
  timestamp::Dates.DateTime
  transaction::Array{Transaction}
  proof::Int
  previous_hash::String
end

mutable struct blockchain
  chain::Array{Block}
  current_transaction::Array{Transaction}
end

max_stars_repo_name: timo-cmd/Bisquit


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    @assert precompile(Tuple{typeof(watch_manifest), String})
    @assert precompile(Tuple{typeof(watch_file), String, Int})
    @assert precompile(Tuple{Rescheduler{typeof(watch_manifest), String}})
    @assert precompile(Tuple{Rescheduler{typeof(revise_dir_queued),Tuple{String}}})
    @assert precompile(Tuple{typeof(revise)})
    @assert precompile(Tuple{typeof(setindex!), ExprsSigs, Nothing, RelocatableExpr})
    @assert precompile(Tuple{typeof(setindex!), ExprsSigs, Vector{Any}, RelocatableExpr})
    @assert precompile(Tuple{typeof(setindex!), ModuleExprsSigs, ExprsSigs, Module})
    @assert precompile(Tuple{typeof(setindex!), Dict{PkgId,PkgData}, PkgData, PkgId})
    @assert precompile(Tuple{typeof(setindex!), Dict{String,WatchList}, WatchList, String})
    # @assert precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)), NamedTuple{(:define,),Tuple{Bool}}, typeof(Revise.methods_by_execution!), Function, CodeTrackingMethodInfo, Dict{Module,Vector{Expr}}, Frame})
    # m = bodymethod(which(methods_by_execution!, (Function, CodeTrackingMethodInfo, Dict{Module,Vector{Expr}}, Module, Expr)))
    # @assert precompile(getfield(Revise, m.name), (Bool, typeof(methods_by_execution!), Function, CodeTrackingMethodInfo, Dict{Module,Vector{Expr}}, Module, Expr))
    @assert precompile(Tuple{typeof(get_def), Method})
    # precompile(Tuple{typeof(run_backend), REPL.REPLBackend})
    # precompile(Tuple{typeof(Revise._watch_package), Base.PkgId})
    # precompile(Tuple{typeof(Revise.watch_package), Base.PkgId})
    # precompile(Tuple{typeof(Revise.sig_type_exprs), Module, Expr})

    # # Here are other methods that require >10ms for inference but which do
    # # not successfully precompile
    # for dct in (watched_files, pkgdatas)
    #     precompile(Tuple{typeof(setindex!), typeof(dct), valtype(dct), keytype(dct)})
    # end
    # precompile(Tuple{typeof(setindex!), ModuleExprsSigs, FMMaps, Module})
    # precompile(Tuple{typeof(setindex!), Dict{String,FileInfo}, FileInfo, String})
    # precompile(Tuple{typeof(watch_file), String, Int})
    # precompile(Tuple{typeof(empty!), Dict{Tuple{PkgData,String},Nothing}})

    return nothing
end


text:
abstract type GeometricStatement end

struct Hypothesis <: GeometricStatement
    points::OrderedDict{Symbol, PointStatus}
    constraints::Vector{GeometricConstraint}
end

function show(io::IO, hp::Hypothesis)
    println(io, "POINTS:")
    println(io, "------------")
    for (p, s) in hp.points
        println(io, p, " : ", s)
    end
    println(io, "")
    println(io, "HYPOTHESIS:")
    println(io, "------------")
    for (i, c) in enumerate(hp.constraints)
        println(io, "(", i, ") ", c)
    end
end

struct Thesis <: GeometricStatement
    constraints::Vector{GeometricConstraint}
end

function show(io::IO, th::Thesis)
    println(io, "THESIS:")
    println(io, "------------")
    for c in th.constraints
        println(io, c)
    end
end

# macros
"""
    @hp(block)

macro to construct the hypothesis of the theorem.

### Input

block -- An expression containing the hypothesis of the theorem, can be a single statement
         or a sequence of statements between `begin...end`.

### Output

An object of type `Hypothesis`.

### Examples

```jldoctest
julia> @hp O := Circle(A, B, C)
POINTS:
------------
A : free
B : free
C : free
O : dependent by (1)

HYPOTHESIS:
------------
(1) OA ≅ OB ≅ OC
```
"""
macro hp(block)
    _parse_hp(block)
end

"""
    @th(block)

macro to construct the thesis of the theorem.

### Input

block -- An expression containing the thesis of the theorem, can be a single statement
         or a sequence of statements between `begin...end`.

### Output

An object of type `Thesis`.

### Examples

```jldoctest
julia> @th Segment(A, O) ≅ Segment(O, C)
THESIS:
------------
AO ≅ OC
```
"""
macro th(block)
    _parse_th(block)
end

#####################
## Prover interface #
#####################

abstract type AbstractGeometricProver end
abstract type AbstractGeometricProof end

struct NoProof <: AbstractGeometricProof end
show(io::IO, ::NoProof) = println(io, "PROOF:\n------\nProof missing")
isproved(::NoProof) = false

"""
    prove(hp::Hypothesis, th::Thesis[, method=RittWuMethod])
    prove(theorem::Theorem[, method=RittWuMethod])

Proves a thereom with given hypothesis `hp` and thesis `th`.

### Input

- `hp::Hypothesis`                  -- hypothesis of the theorem
- `th::Thesis`                      -- thesis of the theorem
- `method<:AbstractGeometricProver` -- method to prove the theorem, default `RittWuMethod`.

Alternatively, instead of passing hypothesis and thesis as two distinct parameters, it is
possible to pass them together as a Theorem object.

### Output

A proof of the the theorem. The type of the output depends on the chosen method

### Algorithms

Currently the following provers are supported

- RittWuMethod

### Examples

```jldoctest
julia> hp = @hp D := Parallelogram(A, B, C, D)
POINTS:
------------
A : free
B : free
C : free
D : dependent by (1)

HYPOTHESIS:
------------
(1) ABCD parallelogram

julia> th = @th Segment(A, B) ≅ Segment(C, D)
THESIS:
------------
AB ≅ CD

julia> prove(hp, th)
Assigned coordinates:
---------------------
A = (0, 0)
B = (u₁, 0)
C = (u₂, u₃)
D = (x₁, x₂)

Goal 1: success

Nondegeneracy conditions:
-------------------------
```
"""
prove(hp::Hypothesis, th::Thesis) = prove(hp, th, RittWuMethod)

#####################
# Theorem interface #
#####################
struct Theorem{T<:AbstractGeometricProof}
    hp::Hypothesis
    th::Thesis
    p::T
end

function show(io::IO, thm::Theorem)
    show(io, thm.hp)
    println(io)
    show(io, thm.th)
    println(io)
    show(io, thm.p)
end

Theorem(hp::Hypothesis, th::Thesis) = Theorem(hp, th, NoProof())

prove(t::Theorem, prover=RittWuMethod) = Theorem(t.hp, t.th, prove(t.hp, t.th, prover))

isproved(t::Theorem) = isproved(t.p)


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

mutable struct IkeEncryption <: SwaggerModel

    function IkeEncryption(;)
        o = new()
        o
    end
end # type IkeEncryption

const _property_map_IkeEncryption = Dict{Symbol,Symbol}()
const _property_types_IkeEncryption = Dict{Symbol,String}()
Base.propertynames(::Type{ IkeEncryption }) = collect(keys(_property_map_IkeEncryption))
Swagger.property_type(::Type{ IkeEncryption }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IkeEncryption[name]))}
Swagger.field_name(::Type{ IkeEncryption }, property_name::Symbol) =  _property_map_IkeEncryption[property_name]

function check_required(o::IkeEncryption)
    true
end

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


text:
"""
    add_loopless_constraints(;
        max_flux_bound = _constants.default_reaction_bound,
        strict_inequality_tolerance = _constants.loopless_strict_inequality_tolerance,
    )

Add quasi-thermodynamic constraints to the model to ensure that no thermodynamically
infeasible internal cycles can occur. Adds the following constraints to the problem:
```
-max_flux_bound × (1 - yᵢ) ≤ xᵢ ≤ max_flux_bound × yᵢ
-max_flux_bound × yᵢ + strict_inequality_tolerance × (1 - yᵢ) ≤ Gᵢ
Gᵢ ≤ -strict_inequality_tolerance × yᵢ + max_flux_bound × (1 - yᵢ)
Nᵢₙₜ' × G = 0
yᵢ ∈ {0, 1}
Gᵢ ∈ ℝ
i ∈ internal reactions
Nᵢₙₜ is the nullspace of the internal stoichiometric matrix
```
Note, this modification introduces binary variables, so an optimization solver capable of
handing mixed integer problems needs to be used. The arguments `max_flux_bound` and
`strict_inequality_tolerance` implement the "big-M" method of indicator constraints.

For more details about the algorithm, see `Schellenberger, Lewis, and, Palsson. "Elimination
of thermodynamically infeasible loops in steady-state metabolic models.", Biophysical
journal, 2011`.
"""
add_loopless_constraints(;
    max_flux_bound = _constants.default_reaction_bound, # needs to be an order of magnitude bigger, big M method heuristic
    strict_inequality_tolerance = _constants.loopless_strict_inequality_tolerance,
) =
    (model, opt_model) -> begin

        internal_rxn_idxs = [
            ridx for (ridx, rid) in enumerate(reactions(model)) if
            !is_boundary(reaction_stoichiometry(model, rid))
        ]

        N_int = nullspace(Array(stoichiometry(model)[:, internal_rxn_idxs])) # no sparse nullspace function

        y = @variable(opt_model, y[1:length(internal_rxn_idxs)], Bin)
        G = @variable(opt_model, G[1:length(internal_rxn_idxs)]) # approx ΔG for internal reactions

        x = opt_model[:x]
        for (cidx, ridx) in enumerate(internal_rxn_idxs)
            @constraint(opt_model, -max_flux_bound * (1 - y[cidx]) <= x[ridx])
            @constraint(opt_model, x[ridx] <= max_flux_bound * y[cidx])

            @constraint(
                opt_model,
                -max_flux_bound * y[cidx] + strict_inequality_tolerance * (1 - y[cidx]) <= G[cidx]
            )
            @constraint(
                opt_model,
                G[cidx] <=
                -strict_inequality_tolerance * y[cidx] + max_flux_bound * (1 - y[cidx])
            )
        end

        @constraint(opt_model, N_int' * G .== 0)
    end


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a native test compiler, and generate reflection methods for it

function native_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = NativeCompilerTarget()
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function native_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function native_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function native_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function native_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    native_method = Symbol("native_$(method)")
    @eval begin
        $native_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $native_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function native_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end


text:
using Test
using SimpleTropical

x = Tropical(3.5)
y = Tropical(4)
inf = TropicalInf

@testset "Comparisons" begin
    @test x == x
    @test x != y

    @testset "Infinity" begin
        @test x != inf
        @test inf != y
        @test inf == inf
        @test inf != Tropical(0)
        @test Tropical(0) != inf

        @test !isequal(x, inf)
        @test !isequal(inf, y)
        @test isequal(inf, inf)
        @test !isequal(inf, Tropical(0))
        @test !isequal(Tropical(0), inf)
    end

    @testset "Not-a-number" begin
        nan = Tropical(NaN)
        @test nan != nan
        @test isequal(nan, nan)
    end
end

@testset "Sum" begin
    @test x + y == Tropical(3.5)
    @test x + inf == x
    @test inf + y == y
end

@testset "Product" begin
    z = Tropical(0)
    @test x * y == Tropical(7.5)
    @test x * z == x
    @test y * z == y
    @test inf * x == inf
    @test y * inf == inf
end

@testset "Inv and power" begin
    @test inv(x) == Tropical(-3.5)
    @test y^-1 == inv(y)
    @test x^10 == Tropical(35.0)
    @test y^0 == Tropical(0)
    @test inf^5 == inf
    @test_throws AssertionError inf^-3
    @test_throws AssertionError inv(inf)
end

@testset "Conversions" begin
    @test convert(Tropical, 5) === Tropical(5)
    @test convert(Tropical{Float64}, 5) === Tropical(5.0)
    @test convert(Tropical, x) === x
    @test convert(Tropical{Int}, y) === y
    @test convert(Tropical{Float32}, x) === Tropical(convert(Float32, x.val))
    @test convert(Tropical{Float64}, inf) === Tropical{Float64}(0.0, true)
end

@testset "Identity elements" begin
    a = Tropical(5)
    @test a + zero(a) == a
    @test a * one(a) == a

    x = ones(Tropical{Int}, 5)
    y = zeros(Tropical, 5)
    @test x + y == x
    @test x .* y == y
end

@testset "⊕ and ⊗ notation" begin
    @test 3 ⊕ 5 == 3
    @test 3 ⊗ 5 == 8

    @test Inf ⊕ 5 == 5
    @test Inf ⊗ 5 == Tropical(Inf)
    @test Inf ⊗ 5 == Inf
end

nothing


text:
module RandomExtensionsTests

using ReTest
include("tests.jl")

end


text:
import Random
function example_design_matrix(;num_rows = 50_000)
    num_categorical_features = 100
    num_continuous_features = 10
    X = zeros(Float64, num_rows, num_categorical_features + num_continuous_features)
    for j = 1:num_categorical_features # binary categorical features
        X[:, j] = rand(num_rows) .> 0.5
    end
    for j = (num_categorical_features+1):(num_categorical_features+num_continuous_features) # continuous features that we have normalized
        X[:, j] = randn(num_rows)
    end
    return X
end


text:
using LazyBandedMatrices, BlockBandedMatrices, BandedMatrices, LazyArrays,
            ArrayLayouts, MatrixFactorizations, LinearAlgebra, Random, Test
import LazyArrays: Applied, resizedata!, FillLayout, MulAddStyle
import LazyBandedMatrices: MulBandedLayout, VcatBandedMatrix, BroadcastBandedLayout
import BandedMatrices: BandedStyle, _BandedMatrix, AbstractBandedMatrix

Random.seed!(0)

struct PseudoBandedMatrix{T} <: AbstractMatrix{T}
    data::Array{T}
    l::Int
    u::Int
end

Base.size(A::PseudoBandedMatrix) = size(A.data)
function Base.getindex(A::PseudoBandedMatrix, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k]
    else
        zero(eltype(A.data))
    end
end
function Base.setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k] = v
    else
        error("out of band.")
    end
end

struct PseudoBandedLayout <: AbstractBandedLayout end
Base.BroadcastStyle(::Type{<:PseudoBandedMatrix}) = BandedStyle()
BandedMatrices.MemoryLayout(::Type{<:PseudoBandedMatrix}) = PseudoBandedLayout()
BandedMatrices.isbanded(::PseudoBandedMatrix) = true
BandedMatrices.bandwidths(A::PseudoBandedMatrix) = (A.l , A.u)
BandedMatrices.inbands_getindex(A::PseudoBandedMatrix, j::Int, k::Int) = A.data[j, k]
BandedMatrices.inbands_setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int) = setindex!(A.data, v, j, k)
LinearAlgebra.fill!(A::PseudoBandedMatrix, v) = fill!(A.data,v)

@testset "Mul" begin
    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = PseudoBandedMatrix(rand(4, 4), 2, 3)
    C = PseudoBandedMatrix(zeros(5, 4), 3, 4)
    D = zeros(5, 4)

    @test (C .= Mul(A, B)) ≈ (D .= Mul(A, B)) ≈ A*B
end

@testset "Vcat Zeros special case" begin
    A = _BandedMatrix((1:10)', 10, -1,1)
    x = Vcat(1:3, Zeros(10-3))
    @test A*x isa Vcat{Float64,1,<:Tuple{<:Vector,<:Zeros}}
    @test length((A*x).args[1]) == length(x.args[1]) + bandwidth(A,1) == 2
    @test A*x == A*Vector(x)

    A = _BandedMatrix(randn(3,10), 10, 1,1)
    x = Vcat(randn(10), Zeros(0))
    @test A*x isa Vcat{Float64,1,<:Tuple{<:Vector,<:Zeros}}
    @test length((A*x).args[1]) == 10
    @test A*x ≈ A*Vector(x)
end

@testset "MulMatrix" begin
    A = brand(6,5,0,1)
    B = brand(5,5,1,0)
    M = ApplyArray(*,A,B)

    @test isbanded(M) && isbanded(Applied(M))
    @test bandwidths(M) == bandwidths(Applied(M))
    @test BandedMatrix(M) == A*B
    @test MemoryLayout(typeof(M)) isa MulBandedLayout
    @test colsupport(M,1) == colsupport(Applied(M),1) == 1:2
    @test rowsupport(M,1) == rowsupport(Applied(M),1) == 1:2

    @test Base.BroadcastStyle(typeof(M)) isa BandedStyle
    @test M .+ A isa BandedMatrix

    V = view(M,1:4,1:4)
    @test bandwidths(V) == (1,1)
    @test MemoryLayout(typeof(V)) == MemoryLayout(typeof(M))
    @test M[1:4,1:4] isa BandedMatrix

    A = brand(5,5,0,1)
    B = brand(6,5,1,0)
    @test_throws DimensionMismatch ApplyArray(*,A,B)

    A = brand(6,5,0,1)
    B = brand(5,5,1,0)
    C = brand(5,6,2,2)
    M = Mul(A,B,C)
    @test @inferred(eltype(M)) == Float64
    @test bandwidths(M) == (3,3)
    @test M[1,1] ≈ (A*B*C)[1,1]

    M = @inferred(ApplyArray(*,A,B,C))
    @test @inferred(eltype(M)) == Float64
    @test bandwidths(M) == (3,3)
    @test BandedMatrix(M) ≈ A*B*C

    M = ApplyArray(*, A, Zeros(5))
    @test colsupport(M,1) == colsupport(Applied(M),1)
    @test_skip colsupport(M,1) == 1:0
end

@testset "Cat" begin
    A = brand(6,5,2,1)
    H = Hcat(A,A)
    @test H[1,1] == applied(hcat,A,A)[1,1] == A[1,1]
    @test isbanded(H)
    @test bandwidths(H) == (2,6)
    @test BandedMatrix(H) == BandedMatrix(H,(2,6)) == hcat(A,A) == hcat(A,Matrix(A)) ==
            hcat(Matrix(A),A) == hcat(Matrix(A),Matrix(A))
    @test hcat(A,A) isa BandedMatrix
    @test hcat(A,Matrix(A)) isa Matrix
    @test hcat(Matrix(A),A) isa Matrix
    @test isone.(H) isa BandedMatrix
    @test bandwidths(isone.(H)) == (2,6)
    @test Base.replace_in_print_matrix(H,4,1,"0") == "⋅"

    H = Hcat(A,A,A)
    @test isbanded(H)
    @test bandwidths(H) == (2,11)
    @test BandedMatrix(H) == hcat(A,A,A) == hcat(A,Matrix(A),A) == hcat(Matrix(A),A,A) ==
            hcat(Matrix(A),Matrix(A),A) == hcat(Matrix(A),Matrix(A),Matrix(A))
    @test hcat(A,A,A) isa BandedMatrix
    @test isone.(H) isa BandedMatrix
    @test bandwidths(isone.(H)) == (2,11)

    V = Vcat(A,A)
    @test V isa VcatBandedMatrix
    @test isbanded(V)
    @test bandwidths(V) == (8,1)
    @test BandedMatrix(V) == vcat(A,A) == vcat(A,Matrix(A)) == vcat(Matrix(A),A) == vcat(Matrix(A),Matrix(A))
    @test vcat(A,A) isa BandedMatrix
    @test vcat(A,Matrix(A)) isa Matrix
    @test vcat(Matrix(A),A) isa Matrix
    @test isone.(V) isa BandedMatrix
    @test bandwidths(isone.(V)) == (8,1)
    @test Base.replace_in_print_matrix(V,1,3,"0") == "⋅"

    V = Vcat(A,A,A)
    @test bandwidths(V) == (14,1)
    @test BandedMatrix(V) == vcat(A,A,A) == vcat(A,Matrix(A),A) == vcat(Matrix(A),A,A) ==
            vcat(Matrix(A),Matrix(A),A) == vcat(Matrix(A),Matrix(A),Matrix(A))
    @test vcat(A,A,A) isa BandedMatrix
    @test isone.(V) isa BandedMatrix
    @test bandwidths(isone.(V)) == (14,1)
end

@testset "BroadcastMatrix" begin
    A = BroadcastMatrix(*, brand(5,5,1,2), brand(5,5,2,1))
    @test eltype(A) == Float64
    @test bandwidths(A) == (1,1)
    @test LazyArrays.colsupport(A, 1) == 1:2
    @test A == broadcast(*, A.args...)

    B = BroadcastMatrix(+, brand(5,5,1,2), 2)
    @test B == broadcast(+, B.args...)

    C = BroadcastMatrix(+, brand(5,5,1,2), brand(5,5,3,1))
    @test bandwidths(C) == (3,2)
    @test MemoryLayout(typeof(C)) == BroadcastBandedLayout{typeof(+)}()
    @test isbanded(C) == true
    @test BandedMatrix(C) == C
end

@testset "Cache" begin
    A = _BandedMatrix(Fill(1,3,10_000), 10_000, 1, 1)
    C = cache(A);
    @test C.data isa BandedMatrix{Int,Matrix{Int},Base.OneTo{Int}}
    @test colsupport(C,1) == rowsupport(C,1) == 1:2
    @test bandwidths(C) == bandwidths(A) == (1,1)
    @test isbanded(C)
    resizedata!(C,1,1);
    @test C[1:10,1:10] == A[1:10,1:10]
    @test C[1:10,1:10] isa BandedMatrix
end

@testset "NaN Bug" begin
    C = BandedMatrix{Float64}(undef, (1,2), (0,2)); C.data .= NaN;
    A = brand(1,1,0,1)
    B = brand(1,2,0,2)
    C .= Mul(A,B)
    @test C == A*B

    C.data .= NaN
    C .= @~ 1.0 * A*B + 0.0 * C
    @test C == A*B
end

@testset "Applied" begin
    A = brand(5,5,1,2)
    @test applied(*,Symmetric(A),A) isa Applied{MulAddStyle}
    B = apply(*,A,A,A)
    @test B isa BandedMatrix
    @test all(B .=== (A*A)*A)
    @test bandwidths(B) == (3,4)
end

@testset "QL tests" begin
    for T in (Float64,ComplexF64,Float32,ComplexF32)
        A=brand(T,10,10,3,2)
        Q,L=ql(A)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        @test Matrix(Q)*Matrix(L) ≈ A
        b=rand(T,10)
        @test mul!(similar(b),Q,mul!(similar(b),Q',b)) ≈ b
        for j=1:size(A,2)
            @test Q' * A[:,j] ≈ L[:,j]
        end

        A=brand(T,14,10,3,2)
        Q,L=ql(A)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        @test_broken Matrix(Q)*Matrix(L) ≈ A

        for k=1:size(A,1),j=1:size(A,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        A=brand(T,10,14,3,2)
        Q,L=ql(A)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        @test Matrix(Q)*Matrix(L) ≈ A

        for k=1:size(Q,1),j=1:size(Q,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        A=brand(T,10,14,3,6)
        Q,L=ql(A)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        @test Matrix(Q)*Matrix(L) ≈ A

        for k=1:size(Q,1),j=1:size(Q,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        A=brand(T,100,100,3,4)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        b=rand(T,100)
        @test ql(A)\b ≈ Matrix(A)\b
        b=rand(T,100,2)
        @test ql(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch ql(A) \ randn(3)
        @test_throws DimensionMismatch ql(A).Q'randn(3)

        A=brand(T,102,100,3,4)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        b=rand(T,102)
        @test_broken ql(A)\b ≈ Matrix(A)\b
        b=rand(T,102,2)
        @test_broken ql(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch ql(A) \ randn(3)
        @test_throws DimensionMismatch ql(A).Q'randn(3)

        A=brand(T,100,102,3,4)
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        b=rand(T,100)
        @test_broken ql(A)\b ≈ Matrix(A)\b

        A = Tridiagonal(randn(T,99), randn(T,100), randn(T,99))
        @test ql(A).factors ≈ ql!(Matrix(A)).factors
        @test ql(A).τ ≈ ql!(Matrix(A)).τ
        b=rand(T,100)
        @test ql(A)\b ≈ Matrix(A)\b
        b=rand(T,100,2)
        @test ql(A)\b ≈ Matrix(A)\b
        @test_throws DimensionMismatch ql(A) \ randn(3)
        @test_throws DimensionMismatch ql(A).Q'randn(3)
    end

    @testset "lmul!/rmul!" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            A = brand(T,100,100,3,4)
            Q,R = qr(A)
            x = randn(T,100)
            b = randn(T,100,2)
            @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
            @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
            @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
            @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
            c = randn(T,2,100)
            @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
            @test rmul!(copy(c), Q') ≈ c*Matrix(Q')

            A = brand(T,100,100,3,4)
            Q,L = ql(A)
            x = randn(T,100)
            b = randn(T,100,2)
            @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
            @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
            @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
            @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
            c = randn(T,2,100)
            @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
            @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
        end
    end

    @testset "Mixed types" begin
        A=brand(10,10,3,2)
        b=rand(ComplexF64,10)
        Q,L=ql(A)
        @test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b
        @test Q*L ≈ A

        A=brand(ComplexF64,10,10,3,2)
        b=rand(10)
        Q,L=ql(A)
        @test Q*L ≈ A
        @test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b

        A = BandedMatrix{Int}(undef, (2,1), (4,4))
        A.data .= 1:length(A.data)
        Q, L = ql(A)
        @test_broken Q*L ≈ A
    end
end

@testset "kron" begin
    A = brand(5,5,2,2)
    B = brand(2,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (5,4)
    @test Matrix(K) == kron(Matrix(A), Matrix(B))

    A = brand(3,4,1,1)
    B = brand(3,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (7,2)
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B))
    K = kron(B,A)
    @test Matrix(K) ≈ kron(Matrix(B), Matrix(A))

    K = kron(A, B')
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B'))
    K = kron(A', B)
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A'), Matrix(B))
    K = kron(A', B')
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A'), Matrix(B'))

    A = brand(5,6,2,2)
    B = brand(3,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (12,4)
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B))

    n = 10; h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))
    D_xx = kron(D², Eye(n))
    @test D_xx isa BandedMatrix
    @test bandwidths(D_xx) == (10,10)
    D_yy = kron(Eye(n), D²)
    @test D_yy isa BandedMatrix
    @test bandwidths(D_yy) == (1,1)
    Δ = D_xx + D_yy
    @test Δ isa BandedMatrix
    @test bandwidths(Δ) == (10,10)

    @testset "#87" begin
        @test kron(Diagonal([1,2,3]), Eye(3)) isa Diagonal{Float64,Vector{Float64}}
    end
end

struct FiniteDifference{T} <: AbstractBandedMatrix{T}
    n::Int
end

FiniteDifference(n) = FiniteDifference{Float64}(n)

Base.getindex(F::FiniteDifference{T}, k::Int, j::Int) where T =
    if k == j
        -2*one(T)*F.n^2
    elseif abs(k-j) == 1
        one(T)*F.n^2
    else
        zero(T)
    end

BandedMatrices.bandwidths(F::FiniteDifference) = (1,1)
Base.size(F::FiniteDifference) = (F.n,F.n)

@testset "Misc" begin
    @testset "Block banded Kron" begin
        n = 10
        h = 1/n
        D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2

        @time D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))
        @time D_yy = BandedBlockBandedMatrix(Kron(Eye(n),D²))
        @time Δ = D_xx + D_yy

        @test Δ isa BandedBlockBandedMatrix
        @test blockbandwidths(Δ) == subblockbandwidths(Δ) == (1,1)
        @test Δ == kron(Matrix(D²), Matrix(I,n,n)) + kron(Matrix(I,n,n), Matrix(D²))

        n = 10
        D² = FiniteDifference(n)
        D̃_xx = Kron(D², Eye(n))
        @test blockbandwidths(D̃_xx) == (1,1)
        @test subblockbandwidths(D̃_xx) == (0,0)

        V = view(D̃_xx, Block(1,1))
        @test bandwidths(V) == (0,0)

        @test BandedBlockBandedMatrix(D̃_xx) ≈ D_xx

        D̃_yy = Kron(Eye(n), D²)
        @test blockbandwidths(D̃_yy) == (0,0)
        @test subblockbandwidths(D̃_yy) == (1,1)

        V = view(D̃_yy, Block(1,1))
        @test bandwidths(V) == (1,1)

        @test BandedBlockBandedMatrix(D̃_yy) ≈ D_yy
    end

    @testset "Diagonal interface" begin
        n = 10
        h = 1/n
        D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
        D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))

        D = Diagonal(randn(n^2))
        @test D_xx + D isa BandedBlockBandedMatrix
        @test blockbandwidths(D_xx + D) == blockbandwidths(D_xx)
        @test subblockbandwidths(D_xx + D) == subblockbandwidths(D_xx)
        @test D_xx + D == Matrix(D_xx) + D

        @test D_xx - D isa BandedBlockBandedMatrix
        @test blockbandwidths(D_xx - D) == blockbandwidths(D_xx)
        @test subblockbandwidths(D_xx - D) == subblockbandwidths(D_xx)
        @test D_xx - D == Matrix(D_xx) - D

        @test D_xx*D == Matrix(D_xx)*D
        @test D_xx*D isa BandedBlockBandedMatrix
        @test blockbandwidths(D_xx*D) == blockbandwidths(D_xx)
        @test subblockbandwidths(D_xx*D) == subblockbandwidths(D_xx)

        @test D*D_xx == D*Matrix(D_xx)
        @test D*D_xx isa BandedBlockBandedMatrix
        @test blockbandwidths(D*D_xx) == blockbandwidths(D_xx)
        @test subblockbandwidths(D*D_xx) == subblockbandwidths(D_xx)
    end
end


text:
include("../src/EasyStocks.jl")
using Test

@testset "EasyStocks.jl" begin
    # Write your tests here.
    @test 1+1==2
end

@testset "Make sure all functions run" begin
    @test try
        mp = ModPar()
        np = NumPar(mp)
        MS = ModelSolution(mp,np)
        SolveModel!(MS)
        true
    catch
        false
    end
end

@testset "Sales Costs" begin
    @test SalesCost(1,2,0.0)==0.0 # Test that cost = 0  when sizes are different, but cost is 0
    @test SalesCost(1,1,0.1)==0.0 # Test that cost = 0  when sizes are the same
    @test SalesCost(1,2,0.1)>0.0 # Test that cost > 0  when sizes are not the same
end


text:
using BinaryProvider
using Libdl: dlext

function compile(libname, tarball_url, hash; prefix=BinaryProvider.global_prefix, verbose=false)
    # download to tarball_path
    tarball_path = joinpath(prefix, "downloads", "src.tar.gz")
    download_verify(tarball_url, hash, tarball_path; force=true, verbose=verbose)

    # unpack into source_path
    tarball_dir = joinpath(prefix, "downloads", split(first(list_tarball_files(tarball_path)), '/')[1]) # e.g. "zeromq-4.2.5"
    source_path = joinpath(prefix, "downloads", "src")
    verbose && @info("Unpacking $tarball_path into $source_path")
    rm(tarball_dir, force=true, recursive=true)
    rm(source_path, force=true, recursive=true)
    unpack(tarball_path, dirname(tarball_dir); verbose=verbose)
    mv(tarball_dir, source_path)

    install_dir = joinpath(source_path, "julia_install")
    verbose && @info("Compiling in $source_path...")
    cd(source_path) do
        run(`./configure --prefix=$install_dir --without-docs --disable-libunwind --disable-perf --disable-eventfd --without-gcov --disable-curve-keygen`)
        run(`make`)
        run(`make install`)
        mkpath(libdir(prefix))
        cp("$install_dir/lib/libzmq.$dlext", joinpath(libdir(prefix), libname*"."*dlext),
           force=true, follow_symlinks=true)
    end
end


text:
module GridModule

using UserExceptionModule: UserError
using CoordsModule: Coords,DirectionIndicator,LatLonCoords,get_next_cell_coords
using CoordsModule: LatLonSectionCoords,Generic1DCoords
using InteractiveUtils

abstract type Grid end

function for_all(function_on_point::Function,
                 grid::Grid)
  throw(UserError())
end

function for_all_fine_cells_in_coarse_cell(function_on_point::Function,
                                           fine_grid::Grid,coarse_grid::Grid,
                                           coarse_cell_coords::Coords)
  throw(UserError())
end

function find_coarse_cell_containing_fine_cell(fine_grid::Grid,coarse_grid::Grid,
                                               fine_cell_coords::Coords)
  throw(UserError())
end

function find_downstream_coords(grid::Grid,
                                flow_direction::DirectionIndicator,
                                coords::Coords)
  throw(UserError())
end

function coords_in_grid(grid::Grid,coords::Coords)
  throw(UserError())
end

function wrap_coords!(grid::Grid,coords::Coords)
  throw(UserError())
end

function get_number_of_cells(grid::Grid)
  throw(UserError)
end

get_number_of_dimensions(obj::T) where {T <: Grid} =
  obj.number_of_dimensions::Int64

struct LatLonGrid <: Grid
  nlat::Int64
  nlon::Int64
  wrap_east_west::Bool
  number_of_dimensions::Int64
  function LatLonGrid(nlat::Int64,
                      nlon::Int64,
                      wrap_east_west::Bool)
    return new(nlat,nlon,wrap_east_west,2)
  end
end

struct UnstructuredGrid <: Grid
  ncells::Int64
  number_of_dimensions::Int64
  clat::Array{Float64,1}
  clon::Array{Float64,1}
  clat_bounds::Array{Float64,2}
  clon_bounds::Array{Float64,2}
  mapping_to_coarse_grid::Array{Int64,1}
  function UnstructuredGrid(ncells::Int64,clat::Array{Float64,1},clon::Array{Float64,1},
                            clat_bounds::Array{Float64,2},clon_bounds::Array{Float64,2},
                            mapping_to_coarse_grid::Array{Int64,1})
    return new(ncells,1,clat,clon,clat_bounds,clon_bounds,mapping_to_coarse_grid)
  end
end

UnstructuredGrid(ncells::Int64,clat::Array{Float64,1},clon::Array{Float64,1},
                 clat_bounds::Array{Float64,2},clon_bounds::Array{Float64,2}) =
  UnstructuredGrid(ncells::Int64,clat::Array{Float64,1},clon::Array{Float64,1},
                   clat_bounds::Array{Float64,2},clon_bounds::Array{Float64,2},
                   [i for i = 1:ncells]::Array{Int64,1})

UnstructuredGrid(ncells::Int64) = UnstructuredGrid(ncells,Array{Float64,1}(undef,ncells),Array{Float64,1}(undef,ncells),
                                                   Array{Float64,2}(undef,3,ncells),Array{Float64,2}(undef,3,ncells))

UnstructuredGrid(ncells::Int64) = UnstructuredGrid(ncells,Array{Float64,1}(undef,ncells),Array{Float64,1}(undef,ncells),
                                                   Array{Float64,2}(undef,3,ncells),Array{Float64,2}(undef,3,ncells))

UnstructuredGrid(ncells::Int64,
                 mapping_to_coarse_grid::Array{Int64,1}) =
  UnstructuredGrid(ncells,Array{Float64,1}(undef,ncells),Array{Float64,1}(undef,ncells),
                   Array{Float64,2}(undef,3,ncells),Array{Float64,2}(undef,3,ncells),
                   mapping_to_coarse_grid)

LatLonGridOrUnstructuredGrid = Union{LatLonGrid,UnstructuredGrid}

function for_all(function_on_point::Function,
                 grid::LatLonGrid)
  for j = 1:grid.nlon
    for i = 1:grid.nlat
        function_on_point(LatLonCoords(i,j))
      end
  end
end

function for_all(function_on_point::Function,
                 grid::UnstructuredGrid)
  for i = 1:grid.ncells
    function_on_point(Generic1DCoords(i))
  end
end

function for_all_with_line_breaks(function_on_point::Function,
                 grid::LatLonGrid)
  for i = 1:grid.nlat
    for j = 1:grid.nlon
      function_on_point(LatLonCoords(i,j))
    end
    println("")
  end
end

function for_all_with_line_breaks(function_on_point::Function,
                                  grid::UnstructuredGrid)
  if grid.ncells >= 16
    block_size::Int64 = round(grid.ncells/16)
    print("                              ")
    for i = 1:block_size
      function_on_point(Generic1DCoords(i))
    end
    println("")
    print("            ")
    for i = (block_size+1):(4*block_size)
      function_on_point(Generic1DCoords(i))
    end
    println("")
    for i = (4*block_size+1):(8*block_size)
      function_on_point(Generic1DCoords(i))
    end
    println("")
    for i = (8*block_size+1):(12*block_size)
      function_on_point(Generic1DCoords(i))
    end
    println("")
    print("            ")
    for i = (12*block_size+1):(15*block_size)
      function_on_point(Generic1DCoords(i))
    end
    println("")
    print("                              ")
    for i = (15*block_size+1):(16*block_size)
      function_on_point(Generic1DCoords(i))
    end
  else
    for i = 1:grid.ncells
      function_on_point(Generic1DCoords(i))
    end
  end
end

function for_section_with_line_breaks(function_on_point::Function,
                 grid::LatLonGrid,section_coords::LatLonSectionCoords)
  for i = section_coords.min_lat:section_coords.max_lat
    for j = section_coords.min_lon:section_coords.max_lon
      function_on_point(LatLonCoords(i,j))
    end
    println("")
  end
end

function for_all_fine_cells_in_coarse_cell(function_on_point::Function,
                                           fine_grid::LatLonGrid,
                                           coarse_grid::LatLonGrid,
                                           coarse_cell_coords::LatLonCoords)
  nlat_scale_factor = fine_grid.nlat/coarse_grid.nlat
  nlon_scale_factor = fine_grid.nlon/coarse_grid.nlon
  for j = 1+(coarse_cell_coords.lon - 1)*nlon_scale_factor:coarse_cell_coords.lon*nlon_scale_factor
    for i = 1+(coarse_cell_coords.lat - 1)*nlat_scale_factor:coarse_cell_coords.lat*nlat_scale_factor
      function_on_point(LatLonCoords(i,j))
    end
  end
end

function for_all_fine_cells_in_coarse_cell(function_on_point::Function,
                                           fine_grid::UnstructuredGrid,
                                           coarse_grid::UnstructuredGrid,
                                           coarse_cell_coords::Generic1DCoords)
  for i = 1:fine_grid.ncells
    if (fine_grid.mapping_to_coarse_grid[i] == coarse_cell_coords.index)
      function_on_point(Generic1DCoords(i))
    end
  end
end

function find_coarse_cell_containing_fine_cell(fine_grid::LatLonGrid,coarse_grid::LatLonGrid,
                                               fine_cell_coords::LatLonCoords)
  fine_cells_per_coarse_cell_lat::Int64 = fine_grid.nlat/coarse_grid.nlat
  fine_cells_per_coarse_cell_lon::Int64 = fine_grid.nlon/coarse_grid.nlon
  coarse_lat::Int64 = ceil(fine_cell_coords.lat/fine_cells_per_coarse_cell_lat);
  coarse_lon::Int64 = ceil(fine_cell_coords.lon/fine_cells_per_coarse_cell_lon);
  return LatLonCoords(coarse_lat,coarse_lon);
end

function find_coarse_cell_containing_fine_cell(fine_grid::UnstructuredGrid,
                                               coarse_grid::UnstructuredGrid,
                                               fine_cell_coords::Generic1DCoords)
  return Generic1DCoords(fine_grid.mapping_to_coarse_grid[fine_cell_coords.index])
end

function find_downstream_coords(grid::T,
                                flow_direction::DirectionIndicator,
                                coords::Coords) where {T<:LatLonGridOrUnstructuredGrid}
  return get_next_cell_coords(flow_direction,coords)
end

function coords_in_grid(grid::LatLonGrid,coords::LatLonCoords)
  return ((1 <= coords.lat <= grid.nlat) || wrap_east_west) &&
          (1 <= coords.lon <= grid.nlon)
end

function wrap_coords(grid::LatLonGrid,coords::LatLonCoords)
  wrapped_lon::Int64 = coords.lon
  if coords.lon < 1
    wrapped_lon += grid.nlon
  elseif coords.lon > grid.nlon
    wrapped_lon -= grid.nlon
  end
  return LatLonCoords(coords.lat,wrapped_lon)
end

function get_number_of_cells(grid::LatLonGrid)
  return grid.nlat*grid.nlon
end

function get_number_of_cells(grid::UnstructuredGrid)
  return grid.ncells
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
# Some general utilities

using Base.Threads

export thread_map, spawn_for_each

function thread_map(f::Function, x)
  buffer = Array{Any}(undef, size(x))

  input = collect(x)

  @threads for ii in 1:length(x)
    buffer[ii] = f(input[ii])
  end

  # Infer the return type using the union of element types
  Array{Union{typeof.(buffer)...}}(buffer)
end

function thread_map(f::Function, x, output_type::Type{T}) where T
  buffer = Array{T}(undef, size(x))

  input = collect(x)

  @threads for ii in 1:length(x)
    buffer[ii] = f(input[ii])
  end

  buffer
end

# Runs a simple worker on a number of tasks
# Returns *nothing*
function spawn_for_each(f::Function, X; threaded=true)
  next_job = Atomic{Int64}(1)

  # Spawn simple workers to run the tasks
  if threaded
    @threads for _ in 1:nthreads()
      while (index = atomic_add!(next_job, 1)) <= length(X)
        f(X[index])
      end
    end
  else
    for x in X
      f(x)
    end
  end

  nothing
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using NReco
using ATools
using Test
using DataFrames
using Distributions
using Statistics
using Logging
using DataFrames
#using StatsModels
# Lower function verbosity
logger = global_logger(SimpleLogger(stdout, Logging.Warn))

fname = "testdata/n3-window-1m-LXe-20mm-1-20.h5"
pdf   = NReco.read_abc(fname)
dconf = NReco.DetConf(0.3f0, 0.05f0, 2.0f0, 100.0f0, 5000.0f0, 7)
sxyz  = pdf.sensor_xyz
wfm   = pdf.waveform
vdf   = pdf.vertices

@testset "util" begin
    @test NReco.select_by_index(sxyz,
    "sensor_id", 0) == NReco.select_by_column_value(sxyz, "sensor_id", 0)

    wevt = NReco.select_event(wfm, 4995000)
    @test mean(wevt.time) ≈ 42.39652f0
    @test NReco.select_by_column_value(wevt, "sensor_id", 1783).sensor_id[1] == 1783
    @test mean(NReco.select_by_column_value_lt(wevt, "time", 5.0).time) < 5.0
    @test mean(NReco.select_by_column_value_gt(wevt, "time", 5.0).time) > 5.0
    @test mean(NReco.select_by_column_value_interval(wevt, "time", 5.0, 10.0).time) >5.0
    @test mean(NReco.select_by_column_value_interval(wevt, "time", 5.0, 10.0).time) <10.0

    _, imx, _ = NReco.find_max_xy(wevt, "sensor_id", "time")
    m, i      = findmax(wevt.time)
    @test wevt.sensor_id[i] == imx
end

@testset "recof" begin
    @test all(NReco.primary_in_lxe(vdf).parent_id .== 0)

    df          = DataFrame(:q1=>Float32[100.0, 120.5, 132.6],
                            :q2=>Float32[123.6, 122.9,  99.9])
    filtered_df = NReco.filter_energies(df, 100.0f0, 150.0f0)
    @test nrow(filtered_df) == 1
    @test all(in(Array(filtered_df[1, :])).([120.5f0, 122.9f0]))

    df[!, :zstd1] = Float32[2.5, 3.2, 4.9]
    df[!, :zstd2] = Float32[7.9, 4.5, 2.3]
    NReco.calculate_interaction_radius!(df, x -> 2*x, "zstd")
    @test all(in(propertynames(df)).([:r1x, :r2x]))
    @test all(isapprox.(2 * df[!, :zstd1], df[!, :r1x]))
    @test all(isapprox.(2 * df[!, :zstd1], df[!, :r1x]))
end

@testset "nemareco" begin
    exp_keys = [:event_id, :phot1, :phot2, :nsipm1, :nsipm2, :q1, :q2,
	            :E1, :E2, :r1,  :r2, :r1x, :r2x,
                :phistd1, :zstd1, :widz1, :widphi1, :corrzphi1,
                :phistd2, :zstd2, :widz2, :widphi2, :corrzphi2,
			    :xs, :ys, :zs, :ux, :uy, :uz, :xt1, :yt1, :zt1,
                :t1, :xt2, :yt2, :zt2, :t2, :x1, :y1, :z1,
                :x2, :y2, :z2, :xr1, :yr1, :zr1, :tr1,
                :xr2, :yr2, :zr2, :tr2, :xb1, :yb1, :zb1, :ta1,
			    :xb2, :yb2, :zb2, :ta2]
    _, result = NReco.nemareco([fname], dconf)
    result_fields = fieldnames(typeof(result[1]))
    @test length(result_fields) == length(exp_keys)
    @test all(in(exp_keys).(result_fields))
    filter_func = ismissing, isnothing, isnan
    corrzphi1   = filter(c -> !any(f -> f(c), filter_func), getfield.(result, :corrzphi1))
    corrzphi2   = filter(c -> !any(f -> f(c), filter_func), getfield.(result, :corrzphi2))
    @test all((corrzphi1 .<= 1.0) .& (corrzphi1 .>= -1.0))
    @test all((corrzphi2 .<= 1.0) .& (corrzphi2 .>= -1.0))
end

@testset "classify_events" begin
    evt_keys   = [:total, :empty, :single, :prompt, :single_prompt, :good_prompt]
    evt_counts = NReco.event_classifier([fname], dconf)
    evt_fields = fieldnames(typeof(evt_counts))
    @test length(evt_fields) == length(evt_keys)
    @test all(in(evt_keys).(evt_fields))

    @test evt_counts.total         == 20
    @test evt_counts.empty         ==  7
    @test evt_counts.single        == 11
    @test evt_counts.prompt        ==  2
    @test evt_counts.single_prompt == 11
    @test evt_counts.good_prompt   ==  2
end

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
void main(){
    String x;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using DataFrames, CSV, Tables
using StanSample
using Test

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "WaffleDivorce.csv"), DataFrame);

stan5_1_t = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] D; // Outcome
 vector[N] A; // Predictor
}

parameters {
 real a; // Intercept
 real bA; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

transformed parameters {
    vector[N] mu;
    mu = a + + bA * A;
}

model {
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ student_t( 2, mu , sigma );
}
generated quantities{
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = student_t_lpdf(D[i] | 2, mu[i], sigma);
}
";

begin
    data = (N=size(df, 1), D=df.Divorce, A=df.MedianAgeMarriage,
        M=df.Marriage)
    m5_1s_t = SampleModel("m5.1s_t", stan5_1_t)
    rc5_1s_t = stan_sample(m5_1s_t; data)

    if success(rc5_1s_t)
        post5_1s_t_df = read_samples(m5_1s_t, :dataframe)
    end
end

if success(rc5_1s_t)
    nt5_1s_t = read_samples(m5_1s_t, :namedtuple)
    df5_1s_t = read_samples(m5_1s_t, :dataframe)
    log_lik_1_t = nt5_1s_t.log_lik'
    a5_1s_t, cnames = read_samples(m5_1s_t, :array; return_parameters=true);
end

st5_1_t = convert_a3d(a5_1s_t, cnames, Val(:table));
@test names(st5_1_t)[end] == Symbol("log_lik.50")
@test size(DataFrame(st5_1_t)) == (4000, 103)

mu = matrix(st5_1_t, "mu")
@test size(mu) == (4000, 50)


text:
using Bridge
using StaticArrays
import Bridge: b, σ, B, β, a, constdiff
const ℝ = SVector{N,T} where {N,T}

"""
    RadialOU{T} <: ContinuousTimeProcess{ℝ{1,T}}

Struct defining Radial Ornstein Uhlenbeck processs
"""
struct RadialOU{T} <: ContinuousTimeProcess{ℝ{1,T}}
    η::T
    σ::T
    RadialOU(η::T, σ::T) where T = new{T}(η, σ)
end

b(t, x, P::RadialOU) = ℝ{1}(-P.η*x[1] + 0.5*P.σ^2/x[1])
σ(t, x, P::RadialOU) = ℝ{1}(P.σ^2)

domain(::RadialOU{T}) where T = LowerBoundedDomain((zero(T),), (1,))
constdiff(::RadialOU) = true
clone(P::RadialOU, θ) = RadialOU(θ...)
params(P::RadialOU) = [P.η, P.σ]

"""
    RadialOUAux{S} <: ContinuousTimeProcess{ℝ{1,S}}
"""
struct RadialOUAux{S} <: ContinuousTimeProcess{ℝ{1,S}}
    trgtDomain::LowerBoundedDomain{S,1}
    η::S
    σ::S
    t::Float64
    u::S
    T::Float64
    v::S

    function RadialOUAux(η::S, σ::S, t::Float64, u::S, T::Float64, v::S) where S
        domain = LowerBoundedDomain((zero(S),), (1,))
        new{S}(domain, η, σ, t, u, T, v)
    end

    function RadialOUAux(domain::LowerBoundedDomain{S,1}, η::S, σ::S,
                         t::Float64, u::S, T::Float64, v::S) where S
        new{S}(domain, η, σ, t, u, T, v)
    end
end

B(t, P::RadialOUAux) = @SMatrix[0.0]
β(t, P::RadialOUAux) = ℝ{1}(0.0)
σ(t, P::RadialOUAux) = ℝ{1}(P.σ^2)
dependsOnParams(::RadialOUAux) = (2,)
constdiff(::RadialOUAux) = true
b(t, x, P::RadialOUAux) = B(t,P)*x + β(t,P)
a(t, P::RadialOUAux) = σ(t,P)*σ(t,P)'
clone(P::RadialOUAux, θ) = RadialOUAux(P.trgtDomain, θ..., P.t, P.u, P.T, P.v)
params(P::RadialOUAux) = [P.η, P.σ]


text:
module test_misc

using Test, XtalsPyTools

@testset "pydep functions" begin
    @test isnothing(XtalsPyTools.load_pydep("bogus_python_package"))
    @test ismissing(XtalsPyTools.check_pydep(:bogus_python_package => "bogus_python_package"))
end

end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
@testset "838.push-dominoes.jl" begin
    @test push_dominoes("RR.L") == "RR.L"
    @test push_dominoes(".L.R...LR..L..") == "LL.RR.LLRRLL.."
    @test push_dominoes(".L.R.") == "LL.RR"
end


text:
# DAE kernel, defined for each time ∈ t.
function f!_simple(resid,du,u,p,t)
    @unpack Ξ₁, L_1, L_2, z, r, μ, υ, π, T, P, x, ω = p

    v_t = u[1:P]
    g_t = u[P+1]
    A = (r(t) - μ - υ^2/2)*I - (μ + υ^2 - g_t)*L_1 - υ^2/2 * L_2 # (19)

    resid[1:P] .= A * v_t - π.(t, z) # (28)
    resid[1:P] .-= du[1:P] # discretized system of ODE for v (28)
    resid[P+1] = Ξ₁*v_t[1] - dot(ω, v_t) + x(t) # value matching condition (29)
end

# DAE constructor
function simpleDAE(parameters, settings)
    @unpack μ, υ, θ, r, x, π = parameters
    @unpack z_ex, T = settings

    @assert z_ex[1] == 0.0 && issorted(z_ex) # validate grid
    z = z_ex[2:end-1]
    P = length(z)

    ω = ω_weights(z_ex, θ, 1)

    bc = (Mixed(ξ = 1), Mixed(ξ = 1)) # boundary conditions for differential operators
    L_1_minus = L₁₋bc(z_ex, bc) # use backward difference as the drift is negative
    L_2 = L₂bc(z_ex, bc)

    stationary = stationary_numerical_simple(parameters, settings)
    g_T = stationary.g
    v_T = stationary.v

    Ξ₁ = 1/(1 - 1*(z[1] - 0.0)) # (24)
    p = (Ξ₁ = Ξ₁, L_1 = L_1_minus, L_2 = L_2, z = z, r = r, υ = υ, π = π, T = T, μ = μ, g_T = g_T, P = P, x = x, ω = ω)

    u_T = [v_T; g_T]
    resid_M1 = zeros(P+1)
    return DAEProblem(f!_simple, resid_M1, u_T, (T, 0.0), differential_vars = [fill(true, P); false], p)
end

function solve_simple_transition(parameters, settings)
    ode_prob = simpleDAE(parameters, settings)
    @unpack ω, x, Ξ₁, P = ode_prob.p
    @unpack ts, ode_solve_algorithm = settings
    sol = Sundials.solve(ode_prob, ode_solve_algorithm, tstops = ts)

    # Calculate the residual at each time point
    residuals = zeros(length(ts))
    v_ts = zeros(P, length(ts))
    g_ts = zeros(length(ts))
    for (i, t) in enumerate(ts)
        v_t = sol(t)[1:P] # i.e., the value function at the point-in-time.
        residuals[i] = Ξ₁*v_t[1] + x(t) - dot(ω, v_t) # value matching condition (29)
        v_ts[:,i] = v_t
        g_ts[i] = sol(t)[end]
    end
    return (residuals = residuals, v_ts = v_ts, g_ts = g_ts)
end


text:
"""
ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля (присутствуют перегородки)
РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке, следующий - за исключением двух последних клеток на Востоке, и т.д.
"""
function stepper(r::Robot) #Маркировка "лесенкой"
    steps=0
    arr_of_steps=[]
    while ((isborder(r,Sud))&&(isborder(r,West)))==false
        push!(arr_of_steps,moves!(r,West))
        push!(arr_of_steps,moves!(r,Sud))
        steps+=2
    end
    putmarker!(r)
    while isborder(r,Ost)==false
        move!(r,Ost)
        putmarker!(r)
    end
    counter=moves!(r,West)
    while isborder(r,Nord)==false && counter>0
        point=counter
        move!(r,Nord)
        while point>0
            putmarker!(r)
            point-=find_decision(r,Ost)
        end
        user_decision(r,West)
        counter-=1
    end
    moves!(r,Sud)
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for a ∈ 1:arr_of_steps[steps]
            move!(r,side)
        end
        steps-=1
    end
end

function find_decision(r::Robot, side::HorizonSide) #Обход локальных преград
    num_steps=0
    while isborder(r,side)==true&&isborder(r,next(side))==false
        move!(r,next(side))
        num_steps+=1
    end
    x=0 #Счётчик
    if isborder(r,side)==false
        move!(r,side)
        x+=1
    end
    if num_steps !=0
        while isborder(r,inverse(next(side)))==true
            move!(r,side)
            x+=1
        end
        for a ∈ 1:num_steps
            move!(r,inverse(next(side)))
        end
    end
    return x #Возвращение значения счётчика
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function user_decision(r::Robot, side::HorizonSide)
    x=1
    while x!=0
        x=find_decision(r,side)
    end
end

function next(side::HorizonSide)
    (HorizonSide(mod(Int(side)+1,4)))
end
function inverse(side::HorizonSide)
    (HorizonSide(mod(Int(side)+2,4)))
end

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
using Test,
      RandomQuantum,
      QuantumInfo,
      Cliffords,
      Qlab

import SchattenNorms.trnorm

function do_unitary2choi()
    choi = unitary2choi(str2unitary("1QX90p"))
    ideal = [0.25+0.0im    0.0+0.25im   0.0+0.25im  0.25+0.0im;
             0.0-0.25im  0.25+0.0im   0.25+0.0im    0.0-0.25im;
             0.0-0.25im  0.25+0.0im   0.25+0.0im    0.0-0.25im;
             0.25+0.0im    0.0+0.25im   0.0+0.25im  0.25+0.0im]
    @test trnorm(choi-ideal) < 3e-2
end

function do_choi2pauliMap()
    choi = unitary2choi(str2unitary("1QX90p"))
    p_map = choi2pauliMap(choi)
    ideal = [1.0 0.0 0.0 0.0;
             0.0 1.0 0.0 0.0;
             0.0 0.0 0.0 1.0;
             0.0 0.0 -1.0 0.0]
    @test isapprox(choi, ideal, atol=1e-14)
end

function do_choi2chi()
    choi = unitary2choi(str2unitary("1QX90p"))
    chi = choi2chi(choi)
    ideal = [0.5 0.0-0.5im 0.0 0.0;
             0.0-0.5im 0.5 0.0 0.0;
             0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0]
    @test isapprox(chi, ideal, atol=1e-14)
end

function test_processFidelity()
    gates = Dict("CNOT12" => Matrix{ComplexF64}([[1,0,0,0] [0,1,0,0] [0,0,0,1] [0,0,1,0]]),
                "InvCNOT12" => Matrix{ComplexF64}([[0,1,0,0] [1,0,0,0] [0,0,1,0] [0,0,0,1]]),
                "CNOT21" => Matrix{ComplexF64}([[1,0,0,0] [0,0,0,1] [0,0,1,0] [0,1,0,0]]),
                "1QId" => I,
                "1QX90p" => exp(-1im*pi*X/4),
                "1QX90m" => exp(1im*pi*X/4),
                "1QY90p" => exp(-1im*pi*Y/4),
                "1QY90m" => exp(1im*pi*Y/4),
                "1QXp" => exp(-1im*pi*X/2),
                "1QYp" => exp(-1im*pi*Y/2),
                "1QX45p" => exp(-1im*pi*X/8),
                "1QX22p" => exp(-1im*pi*X/16),
                "1QHad" => exp(-1im*(pi/2)*(1/sqrt(2))*(X+Z)),
                "1QZ90" => exp(-1im*(pi/4)*Z),
                "1QT" => exp(-1im*(pi/8)*Z),
                "II" => kron(I,I),
                "XI" => exp(-1im*kron(X,I)*pi/2),
                "IX" => exp(-1im*kron(I,X)*pi/2),
                "YI" => exp(-1im*kron(Y,I)*pi/2),
                "IY" => exp(-1im*kron(I,Y)*pi/2),
                "X_2I" => exp(-1im*kron(X,I)*pi/4),
                "IX_2" => exp(-1im*kron(I,X)*pi/4),
                "Y_2I" => exp(-1im*kron(Y,I)*pi/4),
                "IY_2" => exp(-1im*kron(I,Y)*pi/4),
                "XZ90" => exp(-1im*kron(X,Z)*pi/4),
                "ZX90" => exp(-1im*kron(Z,X)*pi/4),
                "XZm90" => exp(1im*kron(X,Z)*pi/4),
                "ZXm90" => exp(1im*kron(Z,X)*pi/4),
                "ZXm90_2" => exp(1im*kron(Z,I)*pi/2)*exp(1im*kron(I,X)*pi/4),
                "X_8I" => exp(-1im*kron(X,I)*pi/16),
                "X_4I" => exp(-1im*kron(X,I)*pi/8),
                "IX_8" => exp(-1im*kron(I,X)*pi/16),
                "IX_4" => exp(-1im*kron(I,X)*pi/8),
                "CNOT12_y" => Matrix{ComplexF64}([[1,0,0,0] [0,1,0,0] [0,0,-1im,0] [0,0,0,1im]]),
                "X90ZX90" => Matrix{ComplexF64}([[1,0,0,0] [0,0,0,-1im] [0,0,1,0] [0,-1im,0,0]]),
                "XX" => exp(-1im*kron(X,I)*pi/2)*exp(-1im*kron(I,X)*pi/2),
                "X90X90" => exp(-1im*kron(X,I)*pi/4)*exp(-1im*kron(I,X)*pi/4),
                "Y90Y90" => exp(-1im*kron(Y,I)*pi/4)*exp(-1im*kron(I,Y)*pi/4)
            );
    for gate in keys(gates)
        proces_fid, gate_fid = getProcessFidelity(str2unitary(gate), gate)
        @test isapprox(proces_fid, 0.0, atol=1e-14)
end

do_unitary2choi()
do_choi2pauliMap()
do_choi2chi()
test_processFidelity()


text:
"""
    get_lattice_timers() :: Nothing

Function to test current performance of lattice implementation for L = 6.
"""
function get_lattice_timers() :: Nothing

    # init timer
    to = TimerOutput()

    # init list of lattices
    lattices = String["square",
                      "honeycomb",
                      "kagome",
                      "triangular",
                      "cubic",
                      "fcc",
                      "bcc",
                      "hyperhoneycomb",
                      "hyperkagome",
                      "pyrochlore",
                      "diamond"]

    # time lattice building
    for name in lattices
        @timeit to "=> " * name begin
            for reps in 1 : 10
                @timeit to "-> build"  l = get_lattice(name, 6, verbose = false)
                @timeit to "-> reduce" r = get_reduced_lattice("heisenberg", [[0.0]], l, verbose = false)
            end
        end
    end

    show(to)

    return nothing
end


text:
"""
    virusfilter()

Returns a `namefinder` limited to the viral division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments. Note that phage are covered by `phagefinder`.
"""
virusfilter() = namefilter(:VRL)

"""
    bacteriafilter()

Returns a `namefinder` limited to the bacterial division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments.
"""
bacteriafilter() = namefilter(:BCT)

"""
    plantfilter()

Returns a `namefinder` limited to the "plant and fungi" division of the NCBI
taxonomy. See the documentation for `namefinder` and `taxid` for more
information about arguments.
"""
plantfilter() = namefilter(:PLN)

"""
    primatefilter()

Returns a `namefinder` limited to the primate division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments.
"""
primatefilter() = namefilter(:PRI)

"""
    rodentfilter()

Returns a `namefinder` limited to the rodent division of the NCBI
taxonomy. See the documentation for `namefinder` and `taxid` for more
information about arguments.
"""
rodentfilter() = namefilter(:ROD)

"""
    mammalfilter(;inclusive::Bool=true)

Returns a `namefinder` limited to the mammal division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments.

If the keyword argument `inclusive` is set to `false`, this will *not* search
for organisms assigned to a lower division, in this case rodents (covered by
`rodentfinder`) and primates (covered by `primatefinder`). The default behavior
is to include these groups.
"""
mammalfilter(inclusive::Bool=true) = inclusive ? namefilter([:MAM, :ROD, :PRI]) : namefilter(:MAM)

"""
    vertebratefilter(;inclusive::Bool=true)

Returns a `namefinder` limited to the vertebrate division of the NCBI taxonomy.
See the documentation for `namefinder` and `taxid` for more information about
arguments.

If the keyword argument `inclusive` is set to `false`, this will *not* search
for organisms assigned to a lower division, in this case mammals (covered by
`mammalfinder`). The default behavior is to include these groups, which also
include the groups covered by `mammalfinder` itself.
"""
vertebratefilter(inclusive::Bool=true) = inclusive ? namefilter([:VRT, :MAM, :ROD, :PRI]) : namefilter(:VRT)

"""
    invertebratefilter()

Returns a `namefinder` limited to the invertebrate division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments.

Note that this is limited organisms not covered by `plantfinder`,
`bacteriafinder`, and `virusfinder`.
"""
invertebratefilter() = namefilter(:INV)

"""
    phagefilter()

Returns a `namefinder` limited to the phage division of the NCBI taxonomy. See
the documentation for `namefinder` and `taxid` for more information about
arguments.
"""
phagefilter() = namefilter(:PHG)

"""
    environmentalsamplesfilter()

Returns a `namefinder` limited to the environmental samples division of the NCBI
taxonomy. See the documentation for `namefinder` and `taxid` for more
information about arguments.
"""
environmentalsamplesfilter() = namefilter(:ENV)


text:
#!/usr/bin/julia

# Trizen
# 14 October 2016
# https://github.com/trizen

# Algorithm from:
#   http://oeis.org/wiki/User:Peter_Luschny/ComputationAndAsymptoticsOfBernoulliNumbers#Seidel

using Printf
const BONE = big"1"

function bernoulli_seidel(n)

    n == 0   && return 1
    n == 1   && return 1//2
    isodd(n) && return 0

    D = append!([0, BONE], [0 for i in 1:(n/2-1)])

    h = 2
    b = false
    for i in 1:n
        if b
            for k in 2:h-1
                D[k] += D[k-1]
            end
        else
            w = h-1
            h += 1
            while w > 1
                D[w] += D[w+1]
                w -= 1
            end
        end
        b = !b
    end

    D[h-1] // ((BONE << (n+1)) - 2) * ((n % 4) == 0 ? -1 : 1)
end

for i in 0:50
    @printf("B%-3d = %s\n", 2*i, bernoulli_seidel(2*i))
end

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
# This file was generated, do not modify it. # hide
ctry_selec = r"BEL|FRA|DEU"
tech_selec = r"Solar"

cap_sum_plot = cap_sum[occursin.(ctry_selec, cap_sum.country) .& occursin.(tech_selec, cap_sum.primary_fuel), :]


text:
using Distributed
using exactQmodule
using Arpack
using PauliStrings
using MPOmodule
using finiteMPS
using FFTW
#@everywhere using ProgressMeter
@everywhere using DelimitedFiles
@everywhere include("finiteTDVP_SSE.jl")

#
#	get frequency bins
#
function fftfreq(n,d=1)
	if n%2 == 0
		bins = [0:1:Int(n/2-1)...,(Int(-n/2):1:-1)...]
	else
		bins = [0:1:Int((n-1)/2)...,(Int(-(n-1)/2):1:-1)...]
	end

	return bins/(d*n)
end

function get_time_domain_func(f, freqvec; sorted = false)
	df = freqvec[2]-freqvec[1]
	N = size(freqvec)[1]

	f_vec = f.(freqvec)

	fft_val = fft(f_vec)
	tvec = fftfreq(N, df)

	fft_val = fft_val.*[(-1)^i for i in 0:N-1]

	dt = tvec[2]-tvec[1]

	fft_val = fft_val*df

	if sorted == true

		idx = sortperm(tvec)
		fft_val = fft_val[idx]
		tvec = tvec[idx]
	end

	return fft_val, tvec
end

function constructMPO(L::Vector{PauliString}, Basis)

	N = L[1].N
	bondDim = length(L)

	mpo_ID_start = MPO((1,bondDim), [localOp(Array{Complex{Float64}}(I,2,2)) for i in 1:bondDim], [i for i in 1:bondDim])
	mpo_ID = MPO(bondDim, [localOp(Array{Complex{Float64}}(I,2,2)) for i in 1:bondDim], [1+bondDim*i+i for i in 0:bondDim-1])
	mpo_ID_end = MPO((bondDim, 1), [localOp(Array{Complex{Float64}}(I,2,2)) for i in 1:bondDim], [i for i in 1:bondDim])

	MPOvec = [mpo_ID_start, [deepcopy(mpo_ID) for i in 2:N-1]..., mpo_ID_end]

	for (i, pString) in enumerate(L)

		l_op = localOp.(Basis[pString.baseIdx.+1]) # +1 since Basis contains ID at index 1 and pString.baseIdx is 0 for id
		l_op[1] = coef(pString)*l_op[1]

		mpos = MPOvec[pString.sites]
		for (j, op) in enumerate(l_op)
			mpo = mpos[j]
			mpo.Operator[i] = op

		end

	end

	MPOvec = depara!(MPOvec)
	return MPOvec
end

function getJump_from_Comutator(a::Vector{PauliString}, H, G_time, time_vec, order)
	@show a[1].sites
	pauli_string = a
	flip_sign(x) = x != 0 ? sign(x) : 1
	L = Vector{PauliString}([])
	dt = time_vec[2]-time_vec[1]
	@show dt
	for n in 0:order
		#  ((-1i)ⁿ/n!) * ∫g(t)⋅tⁿ dt
		c =(2*pi)^(n+1)*sum(flip_sign.(time_vec).*G_time.*(time_vec.^n).*dt)*(((-1im)^n)/factorial(n))

		@show c #just for sanity check to see if real
		push!(L, [PauliString(x.N, x.sites, x.baseIdx, coef(x)*c)
			  for x in pauli_string]...)
		pauli_string = simplify_stringList(vcat([adH(H,p,1) for p in pauli_string]...))
		#pauli_string = vcat([adH(H,p,1) for p in pauli_string]...)
	end
	L = simplify_stringList(L)

	return L
end

#
#
#
include("MPOcollection.jl")

#
#	sys size
N = 10
Np = 5
"""****************************** initialize ham  ************************************"""
V = 1.5
U = 0.5
J = 1.0 #intra
J′ = 0.5 #inter
μ = 0.0
Temp = 2.0

"""**************************** get Jump Operator  *********************************"""
#
#	get Ham as pauli string
#
HRM = Vector{PauliString}([])
for site in 1:N-1 # just has to be build up to order

	push!(HRM, PauliString(N,[site,site+1],[1,1],[-J′,-J][Int(isodd(site))+1]*0.5))
	push!(HRM, PauliString(N,[site,site+1],[2,2],[-J′,-J][Int(isodd(site))+1]*0.5))

	push!(HRM, PauliString(N,[site,site+1],[2,1],1im*[-J′,-J][Int(isodd(site))+1]*0.25))
	push!(HRM, PauliString(N,[site,site+1],[1,2],-1im*[-J′,-J][Int(isodd(site))+1]*0.25))

	push!(HRM, PauliString(N,[site,site+1],[1,2], 1im*[-J′,-J][Int(isodd(site))+1]*0.25))
	push!(HRM, PauliString(N,[site,site+1],[2,1], -1im*[-J′,-J][Int(isodd(site))+1]*0.25))

	push!(HRM, PauliString(N,[site],[3],0.5*[-V, V][Int(isodd(site))+1]))
	push!(HRM, PauliString(N,[site],[0],0.5*[-V, V][Int(isodd(site))+1]))

	# unit cell repulsive interaction  U*nᵢnᵢ₊1 = 0.25*U*(sᶻᵢsᶻᵢ₊₁ + sᶻᵢIᵢ₊₁ + Iᵢsᶻᵢ₊₁ + IᵢIᵢ₊₁)
	if isodd(site)
		push!(HRM, PauliString(N, [site, site+1], [3,3], 0.25*U))
		push!(HRM, PauliString(N, [site, site+1], [3,0], 0.25*U))
		push!(HRM, PauliString(N, [site, site+1], [0,3], 0.25*U))
		push!(HRM, PauliString(N, [site, site+1], [0,0], 0.25*U))  # is irrelevant for expansion
	end

end
push!(HRM, PauliString(N,[N],[3],0.5*[-V, V][Int(isodd(N))+1]))
push!(HRM, PauliString(N,[N],[0],0.5*[-V, V][Int(isodd(N))+1]))

#for Temp in 0.4:0.05:6.0
@show Temp
GHz = 1.0
function J_CutOff(ω, Λ, T)
	res = ω == 0 ?  T/GHz : (ω/(1*GHz))*exp(-ω^2/(2*Λ^2)) / (1-exp(-ω/T))
	return res
end

Λ = 40
γ = 0.02
function S(x) return  J_CutOff(x, Λ, Temp) end
g(x) = sqrt(S(x)/(2*pi))
freqvec = -20*Λ:0.01:20*Λ
G_time, time_vec = get_time_domain_func(g, freqvec, sorted = true) # get g in time

order = 4
tol = nothing #10e-10
jump_pauli = []
jump_op = []
#for i in 1:2:N
for i in 1:2:Int(N)
#for i in Int(N/2)

	@show i
	pauli_string = [PauliString(N, [i], [3], 1.0)]

	#
	#	from commutator
	#
	L = getJump_from_Comutator(pauli_string, HRM, G_time, time_vec, order)
	L .*= sqrt(γ)
	L = simplify_stringList(Vector{PauliString}(L), tol = tol)
	L = simplify_stringList(Vector{PauliString}(L), tol = tol)

	push!(jump_pauli, Vector{PauliString}(L))

end

save("pauli_string_T$(Temp)_γ$(γ)_order$(order).jld", "string", jump_pauli)
#end

sz = Array{ComplexF64}([-1.0 0.0; 0.0 1.0])
sy = Array{ComplexF64}([0.0 -1.0im; 1.0im 0.0])
sx = Array{ComplexF64}([0.0 1.0; 1.0 0.0])
Id(n) = Array{ComplexF64}(I,n,n)
Basis = [Id(2), sx, sy, sz]
constructMPO(x) = constructMPO(x, Basis)
Jump_MPO_List = []
for x in jump_pauli
	@show length(x)

	#testMPS = MPSvec
	#leftCanMPS(testMPS)
	#rightCanMPS(testMPS)
	x = constructMPO(x)
	#Lenv = [Array{ComplexF64,2}(I,1,1)]
	#@show applyTM_MPO(testMPS[1:end], x, Lenv; left = false)

	push!(Jump_MPO_List, x)
end
#construct_MPO(x) =  MPOTensor_to_MPO(Op_to_MPO(N, x; tol = 1e-4))
#Jump_MPO_List = construct_MPO.(jump_op)

# check if particle number conserving
for x in Jump_MPO_List
	@show [y.bDim for y in x]
	#testMPS, __ = finiteMPS(N, 2, 10)
	#testMPS = MPSvec
	#leftCanMPS(testMPS)
	#rightCanMPS(testMPS)
	#Lenv = [Array{ComplexF64,2}(I,1,1)]
	#Nbefore = applyTM_MPO(testMPS[1:end], MPOvec_TotalN[1:end], Lenv; left = false)
	#testMPS = new_apply_MPO(testMPS, x)
	#testMPS = iter_applyMPO(testMPS, x, 10)
	#leftCanMPS(testMPS)
	#rightCanMPS(testMPS)
	#Lenv = [Array{ComplexF64,2}(I,1,1)]
	#Nafter = applyTM_MPO(testMPS[1:end], MPOvec_TotalN[1:end], Lenv; left = false)
	#@show Nbefore, Nafter
end

"""**************************** get Decay Operator  *********************************"""
Decay_MPO_List = []
Decay_pauli_List = []
for i in 1:size(jump_pauli)[1]

	#
	#	construct by commutator
	#

	L_decay = []
	for x in Base.product(conj.(jump_pauli[i]), jump_pauli[i])
		push!(L_decay, *(x...))

	end
	L_decay = simplify_stringList(Vector{PauliString}(L_decay))
	L_decay = simplify_stringList(Vector{PauliString}(L_decay))
	filter!(x->abs.(coef(x))>1e-4,L_decay)
	push!(Decay_pauli_List, L_decay)
	#decay_mpo = constructMPO(L_decay)

	decay_mpo = constructMPO(L_decay)

	#
	#	which one should be used
	#
	#push!(Decay_MPO_List, Decay_1_MPO) #exact
	push!(Decay_MPO_List, decay_mpo) # commutator
	println("size of Decay operators: ", [size(m.Op_index) for m in Decay_MPO_List[end]])

	# check if particle number conserving
	#testMPS, __ = finiteMPS(N, 2, 10)
	#testMPS = MPSvec
	#rightCanMPS(testMPS)
	#Lenv = [Array{ComplexF64,2}(I,1,1)]
	#Nbefore = applyTM_MPO(testMPS[1:end], MPOvec_TotalN[1:end], Lenv; left = false)
	#testMPS = iter_applyMPO(testMPS, Decay_MPO_List[end], 100)
	#rightCanMPS(testMPS)
	#Lenv = [Array{ComplexF64,2}(I,1,1)]
	#Nafter = applyTM_MPO(testMPS[1:end], MPOvec_TotalN[1:end], Lenv; left = false)
	#@show Nbefore, Nafter
end

test = vcat([-0.5im*vcat(Decay_pauli_List...)...]...)
test = simplify_stringList(Vector{PauliString}(test))
test = simplify_stringList(Vector{PauliString}(test))
heff_mpo = constructMPO(test)
#heff_mpo = sum(Decay_MPO_List)
heff_mpo = depara!(heff_mpo)
#heff_mpo = H_MPOvec
@show [size(m.Op_index) for m in heff_mpo]


text:
module Gabi

#Included Files
include(joinpath(Pkg.dir(),"Gabi", "src","models","svm","LibSvm.jl"))
include(joinpath(Pkg.dir(),"Gabi","src","models","svm","SvmModel.jl"))
include(joinpath(Pkg.dir(),"Gabi","src","models","knn","knn.jl"))
# Exports
export Svm
export predict
export train!
export SvmType
export SvmKernel
export Knn
export normalize
end # module


text:
using DataStructures
import LightGraphs: SimpleGraph, add_edge!, neighbors

sorted(i::Int, j::Int) = (min(i, j), max(i, j))

function fill_if_scalar(x::Union{Vector{T}, T}, dims) where T
    (x isa Vector) ? x : fill(x, dims)
end

function invcount(ptr::Vector{Int}, brks::Vector{Int})
    m = length(ptr) - 1
    cnt = counter([sorted(ptr[i], ptr[i + 1]) for i in 1:m  if i + 1 ∉ brks])
    edgecnts = [float(cnt[sorted(ptr[i], ptr[i + 1])]) for i in 1:m]
    edgewt = [x <= 1. ? x : 1 ./ x for x in edgecnts]
    push!(edgewt, 1.)
    edgewt
end

function plateaus(
        β::Vector{Float64},
        g::SimpleGraph,
        tol::Float64)
    n = length(β)
    # plateaus = Vector{Vector{Int}}

    # checked if we have already walked through node, to_check has missing nodes
    nplateaus = 0.
    checked = falses(n)
    to_check = collect(1:n)

    while !isempty(to_check)
        # find next unchecked index
        idx = pop!(to_check)
        while !isempty(to_check) && checked[idx]
            idx = pop!(to_check)
        end

        # edge case, exit loop if there's no unchecked element
        checked[idx] && break

        # current plateau and member conditions
        cur_unchecked = [idx]

        # check every possible boundary of the plateau
        while !isempty(cur_unchecked)
            idx = pop!(cur_unchecked)

            # check the index of unchecked neighbours
            for local_idx in neighbors(g, idx)
                if !checked[local_idx] && abs(β[idx] - β[local_idx]) < tol
                    # check neighbour
                    checked[local_idx] = true

                    # add it to current plateau
                    push!(cur_unchecked, local_idx)
                end
            end
        end
        nplateaus += 1.
    end

    nplateaus
end

# ----------------
function bincount(
    y::Vector{Float64},
    node::Vector{Int},
    splits::Vector{Float64};
    priorcount::Float64 = 0.,
    numnodes::Int = maximum(node)
)
    n, m = length(y), length(splits) - 1
    countmat = fill(priorcount, numnodes, m)
    for j in 1:m
        @simd for i in 1:n
            if splits[j] <= y[i] <  splits[j + 1]
                countmat[node[i], j] += 1.
            end
        end
    end
    countmat
end

function generate_bins2(levels::Int)
    bins = Tuple{Int, Int}[]
    for l = 1:levels
        currlevel_size = 2^(l - 1)
        currlevel_width = 2^(levels - l + 1)
        for k = 1: currlevel_size
            left = currlevel_width * (k - 1) + 1
            right = currlevel_width * k + 1
            push!(bins, (left, right))
        end
    end
    bins
end

function generate_bins(levels::Int)
    bins = Tuple{Int, Int}[]
    for l = 1:levels
        currlevel_size = 2^(l - 1)
        currlevel_width = 2^(levels - l + 1)
        for k = 1: currlevel_size
            left = currlevel_width * (k - 1) + 1
            right = currlevel_width * k
            push!(bins, (left, right))
        end
    end
    bins
end

function pred(k::Int)
    @assert k >= 1 "k must be positive"
    parent, direction = k ÷ 2, k % 2 == 1
end
function allpred(k::Int)
    if k == 1
        Tuple{Int,Int}[]
    else
        p, d = pred(k)
        [allpred(p); (p, d)]
    end
end
function level_numnodes(l::Int)
    [2^(l - 1) + k for k in 0:2^(l - 1) - 1]
end
# -------------------