text:
using OrdinaryDiffEq, DiffEqSensitivity, Zygote, ForwardDiff, Test

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

# Initial condition
u0 = [1.0, 1.0]

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

# LV equation parameter. p = [α, β, δ, γ]
p = [1.5, 1.0, 3.0, 1.0]

# Setup the ODE problem, then solve
prob = ODEProblem(lotka_volterra!, u0, tspan, p)

function loss(p)
  sol = solve(prob, Tsit5(), p=p, save_idxs=[2], saveat = tsteps, abstol=1e-14, reltol=1e-14)
  loss = sum(abs2, sol.-1)
  return loss
end

grad1 = Zygote.gradient(loss,p)[1]
grad2 = ForwardDiff.gradient(loss,p)
@test grad1 ≈ grad2


text:
#Verigen.jl - stuff pertaining to verilog generation and processing.

type Verigen
  module_name::Symbol
  inputs::Vector{Pair{Symbol,VerilogRange}}
  wires ::Dict{Symbol,Tuple{VerilogRange, Tuple}}
  assignments::Array{String}
  modulecalls::Array{String}
  dependencies::Set{Tuple}
  last_assignments::Vector{Pair{Symbol,VerilogRange}}
  Verigen(s::Symbol) = new(s, Tuple{Symbol,VerilogRange}[], Dict{Symbol,Tuple{VerilogRange, Tuple}}(),String[], String[], Set{Tuple}(),[])
end

#formatting for input declarations only.
function v_decl_fmt(r::VerilogRange)
  r == 0:0 && return ""
  return "[$(r.stop):$(r.start)] "
end

function v_fmt(r::VerilogRange)
  r.stop == r.start && return "[$(r.stop)] "
  return "[$(r.stop):$(r.start)] "
end

function v_fmt(i::Integer)
  return "[$i] "
end

type ModuleObject
  moduleparams::Tuple
  modulename::Symbol
  inputlist::Vector{String}
  outputlist::Vector{Pair{Symbol, VerilogRange}}
end

#type for module cache members
type ModuleCache
  txt::String
  module_name::Symbol
  inputs::Vector{Symbol}
  outputlist::Vector{Pair{Symbol, VerilogRange}}
end

const __global_definition_cache = Dict{Tuple,ModuleCache}()
#caches the generated code.  Keys are Tuple{Symbol, Pairs...} where ... are
#pairs of parameter_symbol => parameter.  Values are

const __global_dependency_cache = Dict{Tuple, Set{Tuple}}()
#caches a list of dependencies.  Keys are the same tuple structure as above,
#where pairs are a parameter_symbol =>  parameter.

function wire_shape(w::Tuple{VerilogRange, Tuple})
  v_decl_fmt(w[1])
end
function array_shape(w::Tuple{VerilogRange, Tuple})
  join(["[$(l-1):0]" for l in w[2]])
end

function describe(v::Verigen)
  outputs = v.last_assignments

  inpnames = Set{Symbol}([inp[1] for inp in v.inputs])

  io_declarations = string("\n  ",
    join([string("input ", v_decl_fmt(input[2]), input[1]) for input in v.inputs], ",\n  "),
    ",\n  ",
    join([string("output ", v_decl_fmt(output[2]), output[1]) for output in outputs], ",\n  "))

  wirestrings = [string("  wire ", wire_shape(v.wires[wire]), wire, array_shape(v.wires[wire]), ";") for wire in keys(v.wires) if
    (!(wire in inpnames)) && !(wire in [output[1] for output in outputs])]

  wire_declarations = string(join(wirestrings, "\n"), length(wirestrings) > 0 ? "\n\n" : "")

  assignments = join(v.assignments, "\n")

  modulecalls = string(join(v.modulecalls, "\n"), length(v.modulecalls) > 0 ? "\n" : "")
"""
module $(v.module_name)($io_declarations);

$wire_declarations$modulecalls$assignments
endmodule
"""
end

function modulecall(moduleparams, params)
  plist = [p.lexical_representation for p in params]
  gdef = __global_definition_cache[moduleparams]
  slist = [".$(gdef.inputs[idx]) ($(plist[idx]))" for idx = 1:length(plist)]
  ModuleObject(moduleparams, gdef.module_name, slist, gdef.outputlist)
end

macro verimode(s, p...)
  if length(p) == 1
    assign_call_parameters = :(__call_parameters = [$(p[1])])
  elseif length(p) > 1
    assign_call_parameters = :(__call_parameters = [])
    for idx = 1:length(p)
      q = p[idx]
      assign_call_parameters = :($assign_call_parameters; push!(__call_parameters, $(q)))
    end
  else
    assign_call_parameters = :()
  end

  m = esc(quote
    __synth_mode = $s
    if __synth_mode == :modulecall
      $assign_call_parameters
    end
  end)
  return m
end

macro verigen(s, p...)
  module_name = QuoteNode(s)

  if length(p) > 0
    listgen = :(ppairs = [])
    for parameter in p
      psym = QuoteNode(parameter)
      listgen = :($listgen; push!(ppairs, $psym => $parameter))
    end
    paramgen = :($listgen; __module_params = ($module_name, ppairs...))
  else
    paramgen = :(__module_params = ($module_name,))
  end

  esc(quote
    $paramgen

    if haskey(Verilog.__global_definition_cache, __module_params)
      if __synth_mode == :verilog
        return Verilog.__global_definition_cache[__module_params].txt
      elseif __synth_mode == :modulecall
        return Verilog.modulecall(__module_params, __call_parameters)
      end
    end

    __verilog_state = Verilog.Verigen($module_name)
  end)
end

macro verifin()
  esc(quote
    txt = Verilog.describe(__verilog_state)
    #set the global definition cache.
    Verilog.__global_definition_cache[__module_params] = Verilog.ModuleCache(txt,
      __verilog_state.module_name,
      [vi[1] for vi in __verilog_state.inputs],
      __verilog_state.last_assignments)

    #populate the global dependency cache.
    Verilog.__global_dependency_cache[__module_params] = __verilog_state.dependencies

    #return the correct output type depending on what we've set out.
    if __synth_mode == :verilog
      return txt
    elseif __synth_mode == :modulecall
      return Verilog.modulecall(__module_params, __call_parameters)
    end
  end)
end

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

doc"""
  `@input identifier rangedescriptor`

  binds an input to a given range descriptor.
"""
macro input(identifier, rangedescriptor)
  ident_symbol = QuoteNode(identifier)
  esc(quote
    if (__synth_mode == :verilog || __synth_mode == :modulecall)

      push!(__verilog_state.inputs, ($ident_symbol => $rangedescriptor))
      #also put it in the wires object.
      __verilog_state.wires[$ident_symbol] = ($rangedescriptor, ())
      $identifier = Verilog.WireObject{$rangedescriptor}(string($ident_symbol))
    else
      #in the general case
      if isa($identifier, Integer)
        $identifier = Wire{$rangedescriptor}($identifier)
      elseif isa($identifier, Wire)
        $rangedescriptor == range($identifier) || throw(Verilog.SizeMismatchError())
      end
    end
  end)
end

export @input

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

type AssignError <: Exception; s::String; end

macro assign(ident, expr)

  ##############################################################################
  ## AN ASSIGNMENT THAT IS AN ASSIGNMENT TO A SINGLE IDENTIFIER
  if isa(ident, Symbol)
    ident_symbol = QuoteNode(ident)

    esc(quote
      #check to see if our assignment is passing a wire.
      assign_temp = $expr
      if isa(assign_temp, Verilog.WireObject)
        #check to see if the identifier in the wire table.
        if $ident_symbol in keys(__verilog_state.wires)
          push!(__verilog_state.assignments, string("  assign ", $ident_symbol, " = ", assign_temp.lexical_representation, ";"))
        else
          __verilog_state.wires[$ident_symbol] = (range(assign_temp), ())
          push!(__verilog_state.assignments, string("  assign ", $ident_symbol, " = ", assign_temp.lexical_representation, ";"))
          $ident = Verilog.WireObject{range(assign_temp)}(string($ident_symbol))
        end
        #remember the last assignment
        __verilog_state.last_assignments = [$ident_symbol => range(assign_temp)]
      elseif isa(assign_temp, Wire)
        #if we're passing it a wire object, then it must be either a direct
        #wire declaration or some sort of wire constant.
        if Verilog.assigned(assign_temp)
          __verilog_state.wires[$ident_symbol] = (range(assign_temp), ())
          push!(__verilog_state.assignments, string("  assign ", $ident_symbol, " = ", Verilog.wo_concat(assign_temp), ";"))
          $ident = Verilog.WireObject{range(assign_temp)}(string($ident_symbol))
        else
          __verilog_state.wires[$ident_symbol] = (range(assign_temp), ())
          $ident = Verilog.WireObject{range(assign_temp)}(string($ident_symbol))
        end
      elseif isa(assign_temp, Verilog.ModuleObject)
        mname = assign_temp.modulename
        idsym = $ident_symbol
        mcaller = string(assign_temp.modulename, "_", idsym)
        iplist = assign_temp.inputlist
        mout = assign_temp.outputlist[1].first
        push!(__verilog_state.modulecalls, string("  $mname $mcaller(\n    ", join(iplist, ",\n    "), ",\n    .$mout ($idsym));\n"))
        #add this to the list of dependencies.
        push!(__verilog_state.dependencies, assign_temp.moduleparams)
        #create the wire associated with this module call.
        assign_range = assign_temp.outputlist[1].second
        __verilog_state.wires[$ident_symbol] = (assign_range, ())
        #this could be the last assignment.
        __verilog_state.last_assignments = [$ident_symbol => assign_range]
        #and also instantiate a new variable with this parameter.
        $ident = Verilog.WireObject{assign_range}(string($ident_symbol))
      elseif isa(assign_temp, Array) && (typeof(assign_temp).parameters[1] <: Wire)
        wtype = typeof(assign_temp).parameters[1]
        wdim = typeof(assign_temp).parameters[2]
        wrange = wtype.parameters[1]
        #make sure this array object has a defined type (and not just a blank wire)
        @assert isa(wrange, Verilog.VerilogRange)
        wa_size = size(assign_temp)

        __verilog_state.wires[$ident_symbol] = (wrange, wa_size)
        #go ahead and overwrite all of these with empty wires of the correct type.
        #reassign assign_temp as an array of wireobjects with the same dimensions
        #as the array of wires.
        assign_temp = Array{Verilog.WireObject{wrange}, wdim}(wa_size...)
        for idx = 1:length(assign_temp)
          #fill the wireobject with the appropriate dereferencing indexes.
          name = string($ident_symbol, join(map((n) -> "[$(n-1)]", ind2sub(wa_size, idx))))
          assign_temp[idx] = Verilog.WireObject{wrange}(name)
        end
        $ident = assign_temp
      else
        #just pass the value to ident, without touching it.
        $ident = assign_temp
      end
    end)

  ##############################################################################
  ## ASSIGNING TO AN ARRAY REFERENT.  THIS COULD EITHER BE A WIRE SUBRANGE OR IT
  ## COULD BE AN ARRAY OF WIRES.

  elseif ident.head == :ref
    ident_base = ident.args[1]
    ident_symbol = QuoteNode(ident.args[1])
    ident_reference = ident.args[2]
    esc(quote
      ##########################################################################
      # array-of-wires case.
      if isa($ident_base, Array) && (typeof($ident_base).parameters[1] <: Verilog.WireObject)
        @assert isa($ident_reference, Integer)
        if isa($expr, Wire)
          push!(__verilog_state.assignments, string("  assign ", ($ident).lexical_representation, " = ", Verilog.wo_concat($expr), ";"))
        elseif isa($expr, Verilog.WireObject)
          push!(__verilog_state.assignments, string("  assign ", ($ident).lexical_representation, " = ", ($expr).lexical_representation, ";"))
        elseif isa($expr, Verilog.ModuleObject)
          mname = $expr.modulename
          #create the output "look"
          idsym = string($ident_symbol, "[", ($ident_reference - 1), "]")
          #create verilog's caller assignment name, which can't have brackets,
          #use an underscore to make it a subscript.
          mcaller = string($expr.modulename, "_", $ident_symbol, "_", ($ident_reference - 1))
          iplist = $expr.inputlist
          mout = $expr.outputlist[1].first
          push!(__verilog_state.modulecalls, string("  $mname $mcaller(\n    ", join(iplist, ",\n    "), ",\n    .$mout ($idsym));\n"))
          #add this to the list of dependencies.
          push!(__verilog_state.dependencies, $expr.moduleparams)
          #asserting that the sizes match.
          if length(($expr).outputlist[1].second) != length($ident)
            print("length of ", $expr.outputlist[1], " mismatches ")
            print("length of ", $ident, ".")
            throw(Verilog.SizeMismatchError())
          end
        end
      else
      ##########################################################################
      # wire subrange case
        assign_temp = $expr
        if isa($ident_reference, Verilog.RelativeRange)
          parsed_reference = Verilog.parse_msb($ident_reference, (__verilog_state.wires[$ident_symbol])[1])
        elseif isa($ident_reference, Verilog.msb)
          parsed_reference = (__verilog_state.wires[$ident_symbol])[1].stop - ($ident_reference).value
        elseif isa($ident_reference, Type{Verilog.msb})
          parsed_reference = (__verilog_state.wires[$ident_symbol])[1].stop
        else
          parsed_reference = $ident_reference
        end

        if isa(assign_temp, Verilog.WireObject)
          if $ident_symbol in keys(__verilog_state.wires)
            push!(__verilog_state.assignments, string("  assign ", $ident_symbol, Verilog.v_fmt(parsed_reference), "= ", assign_temp.lexical_representation, ";"))
          else
            throw(AssignError("can't make a partial assignment to a nonexistent wire."))
          end
        elseif isa(assign_temp, Verilog.ModuleObject)
          if $ident_symbol in keys(__verilog_state.wires)
            mname = assign_temp.modulename
            mcaller = string(assign_temp.modulename, "_", $ident_symbol, "_", (parsed_reference).stop, "_", (parsed_reference).start)
            mout = assign_temp.outputlist[1].first
            iplist = assign_temp.inputlist
            idsym = string($ident_symbol, Verilog.v_fmt(parsed_reference))
            push!(__verilog_state.modulecalls, string("  $mname $mcaller(\n    ", join(iplist, ",\n    "), ",\n    .$mout ($idsym));\n"))

            #this could be the last assignment.
            __verilog_state.last_assignment = [$ident_symbol => range(assign_temp)]
          else
            throw(AssignError("can't make a partial assignment to a nonexistent wire."))
          end
        elseif isa(assign_temp, Verilog.Wire)
          #assume that naked wire objects that are tried to be assigned must be
          #constant values.
          if Verilog.assigned(assign_temp)
            push!(__verilog_state.assignments, string("  assign ", $ident_symbol, Verilog.v_fmt($ident_reference), "= ", Verilog.wo_concat(assign_temp), ";"))
          else
            throw(UnassignedError())
          end
        else
          $ident = $expr
        end
      end
    end)

  ##############################################################################
  ## ASSIGNING FROM A MODULE CALL TO A TUPLE.  WIRE ARRAY ASSIGNMENTS NOT
  ## SUPPORTED (FOR NOW)

  elseif ident.head == :tuple
    ident_list = ident.args
    assign_temp = expr

    #unroll the process of assigning the arguments.  This gets created no matter what,
    #it might not get used in the case that some other set of variables is using
    #a tuple assignment.
    assignment_code = :()
    for idx = 1:length(ident_list)
      ident = ident_list[idx]
      ident_symbol = QuoteNode(ident_list[idx])
      if isa(ident, Symbol)
        #unroll this!
        assignment_code = quote
          $assignment_code
          #update the full list of identifiers.
          push!(idlist, string($ident_symbol))
          push!(oplist, string(".", output_list[$idx].first, " (", $ident_symbol, ")"))
          #create the wire associated with this module call.
          __verilog_state.wires[$ident_symbol] = (output_list[$idx].second, ())
          #instantiate a new variable with this parameter.
          $ident = Verilog.WireObject{output_list[$idx].second}(string($ident_symbol))
        end

      elseif isa(ident, Expr) && (ident.head == :ref)

        #the case where it's a reference, it could either be a wire subreference
        #or it could be a wire array.

        ident_ref = ident.args[1]
        ident_symbol = QuoteNode(ident.args[1])
        ident_range = ident.args[2]

        assignment_code = quote
          $assignment_code

          #check if it's a wire object.
          if isa($ident_ref, Array) && (typeof($ident_ref).parameters[1] <: Verilog.WireObject)
            #assert that the index is not something strange.
            @assert isa($ident_range, Integer)
            #push this into the identifier list
            push!(idlist, string($ident_symbol, "_", $ident_range - 1))
            #push this onto the output list.
            push!(oplist, string(".", output_list[$idx].first, " (", $ident_symbol,"[", $ident_range - 1,"])"))
          elseif isa($ident_ref, Verilog.WireObject)
            wirerange = output_list[$idx].second
            dest_range = $ident_range
            if isa(dest_range, Integer)
              @assert (length(wirerange) == 1)
              push!(idlist, string($ident_symbol, "_", dest_range))
            else
              @assert (length(wirerange) == length(dest_range))
              push!(idlist, string($ident_symbol, "_", dest_range.start == dest_range.stop ? dest_range.start : "$(wirerange.stop)_$(wirerange.start)"))
            end
            push!(oplist, string(".", output_list[$idx].first, " (", $ident_symbol, strip(Verilog.v_fmt(dest_range)),")"))
          end
        end

      else #do nothing.
      end
    end

    #we're trying to make a tuple assignment.
    esc(quote
      if isa($assign_temp, Verilog.ModuleObject)

        #assign the module name.
        mname = $assign_temp.modulename
        idlist = String[]
        input_list = $assign_temp.inputlist
        output_list = $assign_temp.outputlist
        ident_list = $ident_list
        iplist = $assign_temp.inputlist
        oplist = String[]

        #do the unrolled assignment.
        $assignment_code

        mcaller = string($assign_temp.modulename, "_", join(idlist, "_"))

        #add this to the list of module calls.
        push!(__verilog_state.modulecalls, string("  $mname $mcaller(\n    ", join(iplist, ",\n    "), ",\n    ",
                                                                              join(oplist, ",\n    "), ");\n"))
        #add this to the list of dependencies.
        push!(__verilog_state.dependencies, $assign_temp.moduleparams)

      else
        #just do the boring assignment
        $ident = $assign_temp
      end
    end)
  else
    #transparently pass on the assigment without intercepting it.
    esc(:($ident = $expr))
  end
end

macro suffix(stringvalue)
  esc(quote
    __verilog_state.module_name = string(__verilog_state.module_name, "_", $stringvalue)
  end)
end

macro final(identifiers...)
  if length(identifiers) == 1
    ident_symbol = QuoteNode(identifiers[1])
    esc(:(__verilog_state.last_assignments = [$ident_symbol => (__verilog_state.wires[$ident_symbol])[1]]))
  else
    #println(:($identifiers))
    esc(quote
      pairlist = []
      for sym in $identifiers
        push!(pairlist, sym => (__verilog_state.wires[sym])[1])
      end
      __verilog_state.last_assignments = pairlist
    end)
  end
end

export @suffix

ext: jl


text:
# # Duffing Oscillator
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/DuffingOscillator.ipynb)
#
#md # !!! note "Overview"
#md #     System type: nonlinear\
#md #     State dimension: 2\
#md #     Application domain: Nonlinear physics
#
# ## Implementation

using ReachabilityAnalysis, Plots

@taylorize function duffing!(du, u, p, t)
    local α = -1.0
    local β = 1.0
    local δ = 0.3
    local γ = 0.37

    x, v = u
    f = γ * cos(ω * t)

    du[1] = u[2]
    du[2] = - α*x - δ*v - β*x^3 + f
end

ω = 1.2
T = 2*pi / ω
X0 = Singleton([1.0, 0.0]) ⊕ BallInf(zeros(2), 0.1)
prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim=2);

# ## Results

sol = solve(prob, tspan=(0.0, 20*T), alg=TMJets21a());

#-

plot(sol, vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue)

#-

plot(sol, vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue)

#-

plot(sol, vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red)


text:
using JLD

abstract type AbstractSolver end

struct RK4 <: AbstractSolver
    param
end

struct Problem

    param   :: NamedTuple
    solver  :: AbstractSolver

end

struct ProblemSave

    param   :: Dict
    solver  :: Symbol

    ProblemSave( param, solver) = new( param, solver)

end

import Base.convert

function convert(::Type{ProblemSave}, p :: Problem )

    param   = Dict(pairs(p.param))
    solver  = Symbol(typeof(p.solver))

    ProblemSave(param, solver)

end

function convert(::Type{Problem}, p :: ProblemSave)

    @show param = (;p.param...)
    @show p.solver

    if p.solver == :RK4
        solver = RK4(param)
    end

    Problem(param, solver)

end

import JLD.save

function save(p::Problem, name::String)

    filename = string(name,".jld")
    JLD.save(filename, name, convert(ProblemSave, p))

end

import JLD.load

function load(name::String)

    convert(Problem, JLD.load(string(name,".jld"), name ))

end

param = ( ϵ  = 1/2,
          N  = 2^10,
          L  = 10,
          T  = 5,
          dt = 0.01,
          theta = 1)

problem = Problem( param, RK4(param) )

p_dump  = convert(ProblemSave, problem )
p_load  = convert(Problem, p_dump )

save(problem, "problem")

p_load = load("problem")


text:
module HPFEM

using Jacobi
using Base.LinAlg
using Polynomials
using JuBLAS

import Base.length
import Base.call
import Base.size

# package code goes here

include("mesh1d.jl")
include("localnum.jl")
include("basis_fun1d.jl")
include("quadrature1d.jl")
include("basis1d.jl")
include("element1d.jl")
include("operator.jl")
include("dofmap.jl")
include("dirichilet_lift.jl")
include("bbmatrix.jl")
include("static_cond_solver.jl")
include("plotmatrix.jl")

export isbndryint, ModalC01d, Legendre1d, Lagrange1d
export QuadType
export nbndry, ninterior, bndry_idx, interior_idx
export seq2bi!, seq2bi, bi2seq!, bi2seq, seq2b!, seq2b, seq2i!, seq2i
export Basis1d, qnodes, nmodes, nquad, qweights, basis, qbasis, dqbasis, basis1d, basis1d!
export SpecElem1d, project, project!
export mass_matrix, stiff_matrix, plot_matrix
end # module


text:
## Before we got started
# I added the CSV package via
# typing "]" then "add CSV"
# alternatively, we could use
# using Pkg; Pkg.add("CSV")

## Let's start by reading our house prices
using CSV
data = CSV.read("1-julia-intro/houses.csv")

## Let's get the first house info
data[1,:]'

## See the column names
names(data)

##
scatter(data[:latitude], data[:longitude])

## Let's add house prices
scatter(data[:latitude], data[:longitude], marker_z=Int.(data[:price]))

## Let's see zipcodes
# how many are there
@show unique(Int.(data[:zip]))

## TODO - show a true categorical vector plot of zipcodes.
## Plot by zip
scatter(data[:latitude], data[:longitude], marker_z=Int.(data[:zip]))

## Let's plot house-prices by sqft.
scatter(data[:latitude], data[:longitude],
    marker_z=Int.(data[:price])./Int.(data[:sq__ft]))

##
datar = data[data[:type] .== "Residential",:]
datar = datar[datar[:sq__ft] .> 0, :]
##
scatter(datar[:latitude], datar[:longitude],
    marker_z=Int.(datar[:price])./Int.(datar[:sq__ft]))

## Get the average price per sqft based on zipcodes.
# Let's use Julia's data frames first
using Statistics
using DataFrames
datar.price_per_sqft = Int.(datar[:price])./Int.(datar[:sq__ft])
by(datar, :zip, :price_per_sqft => mean)

##
"""
This takes in a set of data in terms of zipcodes
and price_per_sqrt data for a house sold in that zipcode.
    So the two vectors zipcodes, price_per_sqft should have
    the same length. Then we will return a Dictionary
    which is the average price per sqft in each zipcode.
"""
function avgprice_by_zip(zipcodes, price_per_sqft)
    zipprices = Dict{Int,Float64}()
    zipcounts = Dict{Int,Int}() # this counts the number of vals
    @assert length(zipcodes) == length(price_per_sqft)
    for i in 1:length(zipcodes)
        if !(haskey(zipprices, zipcodes[i]))
            # initialize
            zipprices[zipcodes[i]] = 0
            zipcounts[zipcodes[i]] = 0
        end
        zipprices[zipcodes[i]] += price_per_sqft[i]
        zipcounts[zipcodes[i]] += 1 # add one to our count
    end
    for k in keys(zipprices)
        zipprices[k] /= zipcounts[k]
    end
    return zipprices
end
avgs = avgprice_by_zip(Int.(datar[:zip]), datar[:price_per_sqft])

##
bar(collect(keys(avgs)), collect(values(avgs)))

## Our goal -- is there really a relationship between
## prices and zip codes?
## LEt's pick an arbitrary zipcode
using Random
testzip = 95690
# Here the avg prices
testavg = avgs[testzip]

function permtest_on_zip(N::Int)
    # note that this is going to use
    # global variables from outside,
    # you can just use them!
    zipcodes = Int.(datar[:zip])
    prices = datar[:price_per_sqft]
    rand_prices = zeros(N)
    for i=1:N
        p = randperm(length(zipcodes))
        zipcodes = zipcodes[p]
        randavgs = avgprice_by_zip(zipcodes, prices)
        rand_prices[i] = randavgs[testzip]
    end
    return rand_prices
end
rand_prices = permtest_on_zip(1000)
@show sum(rand_prices .>= testavg)
histogram(rand_prices)
vline!([testavg])

## Let's look at a max-vs-min price by zip
diffs = maximum(values(avgs)) - minimum(values(avgs))

## Differs under perm test
function permtest_diff_on_zip(N::Int)
    # note that this is going to use
    # global variables from outside,
    # you can just use them!
    zipcodes = Int.(datar[:zip])
    prices = datar[:price_per_sqft]
    rand_diffs = zeros(N)
    for i=1:N
        p = randperm(length(zipcodes))
        zipcodes = zipcodes[p]
        randavgs = avgprice_by_zip(zipcodes, prices)
        rand_diffs[i] = maximum(values(randavgs)) - minimum(values(randavgs))
    end
    return rand_diffs
end
rand_diffs = permtest_diff_on_zip(1000)
@show sum(rand_diffs .>= diffs)
histogram(rand_prices)
vline!([diffs])


text:
import Codex
import Codex.TransformExport
import CSV

using DataFrames
using Test

data_dir = joinpath(Codex.PROJECT_ROOT, "test", "data")
export_dir = joinpath(data_dir, "2020-08")

@testset "TransformExport" begin
    export_dir = joinpath(data_dir, "2020-08")
    dfs = TransformExport.responses(export_dir)
    @test typeof(dfs) == Dict{String,DataFrame}
    @test "first" in keys(dfs)
    @test "second" in keys(dfs)

    @test TransformExport.rm_timing!(DataFrame(:v1_timing => 1)) == DataFrame()
    @test TransformExport.rm_boring_timestamps!(DataFrame(completed_at = 0, opened_at = 1, open_from = 2)) == DataFrame(completed_at = 0)
    @test TransformExport.rm_boring_foreign_keys!(DataFrame(t_for_id = 1, t_by_id = 2)) == DataFrame(t_by_id = 2)
    @test TransformExport.rm_empty_rows!(DataFrame(filled_out_by_id = [1, missing])) == DataFrame(filled_out_by_id = 1)
    @test TransformExport.rename_id_col!(DataFrame(filled_out_by_id = 1)) == DataFrame(id = 1)

    df = TransformExport.read_csv(joinpath(export_dir, "responses", "first.csv"); delim=';')
    without_descr = TransformExport.rm_descriptions(df)
    @test TransformExport.simplify(without_descr) ==
        DataFrame(id = "aaaa", completed_at = "22-08-2020 14:49:47", v2 = 2, v3 = 3)
    @test string(TransformExport.parsedatetime(df[1, :completed_at])) == "2020-08-22T14:49:47"

    @test size(df) == (3, 17)
    simple = TransformExport.simplify(df)
    @test simple == DataFrame(id = "aaaa", completed_at = "22-08-2020 14:49:47", v2 = "2
 (lorem)", v3 = "3 (heel erg)")

    @test TransformExport._rm_description("1") == 1
    @test TransformExport._rm_description("23 (lorem)") == 23
    @test TransformExport._rm_description("45
 (lorem)") == 45
    @test TransformExport._rm_descriptions(["6"]) == [6]
    @test TransformExport._contains_description(["7 (ipsum)"]) == true
    @test TransformExport._contains_description(["2018"]) == false

    people = DataFrame(:person_id => ["aaaa", "aaab"], :first_name => ["0001", "0002"])
    with_names = DataFrame(id = "0001", completed_at = "22-08-2020 14:49:47", v2 = "2
 (lorem)", v3 = "3 (heel erg)")
    @test TransformExport.substitute_names(simple, people) == with_names

    id_username = TransformExport.read_csv(joinpath(data_dir, "id-username.csv"))
    @test TransformExport.names2usernames(with_names, id_username)[!, :id] == ["jackson"]

    third = TransformExport.read_csv(joinpath(export_dir, "responses", "third.csv"))
    @test Set(first(TransformExport.split_datetime(third, :completed_at))) == Set(["kim", Date("2020-01-01"), "10:00:00", 2])
end


text:
module ConcurrentCollectionsTests

include("utils.jl")
include("test_bench_dict_histogram.jl")
include("test_bench_smoke.jl")
include("test_crq.jl")
include("test_dict.jl")
include("test_dlcrq.jl")
include("test_doctest.jl")
include("test_lcrq.jl")
include("test_mpcrq.jl")
include("test_msqueue.jl")
include("test_ssqueue.jl")
include("test_tsstack.jl")
include("test_work_stealing_deque.jl")
include("test_aqua.jl")

function __init__()
    @info "Starting test" Threads.nthreads()
end

end  # module ConcurrentCollectionsTests


text:
import Dates: Date, year, month, day

cd(dirname(@__FILE__))

function jl_am_put_crr(S,K,r,c,σ,dt1,dt2,N)
    t = Int(dt2 - dt1)/365
    Δt = t/N
    u = exp(σ*√Δt)
    d = exp(-σ*√Δt)
    B = exp((r-c)*Δt)
    p = (B-d)/(u-d)
    q = (u-B)/(u-d)
    iR = exp(-Δt*r)

    # final time
    Z = map(0:N) do i
        x = K - S*exp((2*i-N)σ*√Δt)
        max(x,0)
    end

    for n = N-1:-1:0
        for i = 0:n
            x = K - S*exp((2*i-n)σ*√Δt)
            y = iR*(q*Z[i+1] + p*Z[i+2])
            Z[i+1] = max(x, y)
        end
    end
    return Z[1]
end

function ql_am_put_crr(S,K,r,c,σ,dt1::Date,dt2::Date, n)
    ccall((:ql_am_put_crr, :quantlib), Cdouble,
          (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
           Cint, Cint, Cint,
           Cint, Cint, Cint,
           Cint),
          S, K, r, c, σ,
          year(dt1), month(dt1), day(dt1),
          year(dt2), month(dt2), day(dt2),
          n)
end

using BenchmarkTools
using Miletus

println("CRR Model")
println("=========")
println()

println("Miletus")
println(@benchmark value(CRRModel(Date(2010,01,01), Date(2011,01,01), 801, 36.0, 0.06, 0.0, 0.2),
                 AmericanPut(Date(2011,01,01), SingleStock(), 40.0)))

println("Miletus Currency")
println(@benchmark value(CRRModel(Date(2010,01,01), Date(2011,01,01), 801, 36.0USD, 0.06, 0.0, 0.2),
                 AmericanPut(Date(2011,01,01), SingleStock(), 40.0USD)))

if isfile("quantlib")
    println()
    println("QuantLib")
    println(@benchmark ql_am_put_crr(40.0, 36.0, 0.06, 0.0, 0.2,
                                     Date(2010,01,01), Date(2011,01,01), 801))
end

println()
println("Pure Julia")
println(@benchmark jl_am_put_crr(40.0, 36.0, 0.06, 0.0, 0.2,
           Date(2010,01,01), Date(2011,01,01), 801))
println()
println()

println("Monte Carlo Model")
println("=========")
println()

d1 = Dates.today()
d2 = d1 + Dates.Day(120)
m = GeomBMModel(d1, 100.00, 0.05, 0.0, 0.3)
o = EuropeanCall(d2, SingleStock(), 100.00)

println("Miletus")
println(@benchmark value($m, $o, $MonteCarloModel, 10_000))
println()

m = GeomBMModel(d1, 100.00USD, 0.05, 0.0, 0.3)
o = EuropeanCall(d2, SingleStock(), 100.00USD)

println("Miletus Currency")
println(@benchmark value($m, $o, $MonteCarloModel, 10_000))
println()

println("Least Squares Monte Carlo Model")
println("=========")
println()

d1 = Dates.today()
d2 = d1 + Dates.Day(120)
m = GeomBMModel(d1, 100.00, 0.05, 0.0, 0.3)
o = AmericanPut(d2, SingleStock(), 100.00)

println("Miletus")
println(@benchmark value($m, $o, $LeastSquaresMonteCarlo, 10_000, 3))
println()


text:
function search_msg(
    search::String,
    field::Field,
    sort_by::SortBy,
    sort_order::SortOrder,
    max_results::Integer,
)
    println("\narXiv.jl: processing request...")
    println("Searching \"$(Symbol(field))\" for \"$(search)\"...")
    println("$(max_results) results are sorted by $(sort_by) with $(sort_order) order.")
end


text:
abstract type AbstractBand end

"""
Energy band from DFT calculation.
"""
mutable struct Band <: AbstractBand
    k_points_cart  :: Vector{Vec3}
    k_points_cryst :: Vector{Vec3{Float64}}
    eigvals        :: Vector{Float64}
    extra          :: Dict{Symbol,Any}
end
function Band(k_points_cart, k_points_cryst, eigvals)
    return Band(k_points_cart, k_points_cryst, eigvals, Dict{Symbol,Any}())
end

function kpoints(band::Band, kind = :cryst)
    return kind == :cart ? band.k_points_cart : band.k_points_cryst
end
eigvals(band::Band) = band.eigvals
StructTypes.StructType(::Type{<:AbstractBand}) = StructTypes.Mutable()

"""
    bandgap(bands::AbstractVector{Band}, fermi=0.0)

Calculates the bandgap (possibly indirect) around the fermi level.
"""
function bandgap(bands::Union{Iterators.Flatten,AbstractVector{<:AbstractBand}},
                 fermi = 0.0)
    max_valence = -Inf
    min_conduction = Inf
    for b in bands
        max = maximum(eigvals(b) .- fermi)
        min = minimum(eigvals(b) .- fermi)
        if max_valence <= max <= 0.0
            max_valence = max
        end
        if 0.0 <= min <= min_conduction
            min_conduction = min
        end
    end
    return min_conduction - max_valence
end
function bandgap(u_d_bands::Union{NamedTuple,Tuple}, args...)
    return bandgap(Iterators.flatten(u_d_bands), args...)
end

mutable struct TimingData
    name::String
    cpu::Dates.CompoundPeriod
    wall::Dates.CompoundPeriod
    calls::Int
    children::Vector{TimingData}
end
StructTypes.StructType(::Type{TimingData}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Dates.CompoundPeriod}) = StructTypes.Struct()
StructTypes.StructType(::Type{Dates.Millisecond}) = StructTypes.Struct()
StructTypes.StructType(::Type{Dates.Second}) = StructTypes.Struct()


text:
import Base.*

#=

* : matrix multiply between two Assoc.

If there are Strings in Assoc Value mapping, the Assoc will be reduced to logical.

=#
function *(A::Assoc,B::Assoc)
    #Check A,B, if string => Logical
    At = A
    Bt = B
    if(!isa(A.val[1],Number))
        At = logical(A)
    end
    if(!isa(B.val[1],Number))
        Bt = logical(B)
    end
    ## A*B operation
    #Improve intersect with unique sorted sets
    ABintersect = sortedintersect(At.col,Bt.row)

#    if (size(ABintersect,1) == 0)
#        return Assoc([1],[1],0,(+))
#    end

#    AintMap = searchsortedmapping(ABintersect,At.col)
#    BintMap = searchsortedmapping(ABintersect,Bt.row)

    Aref = @spawn searchsortedmapping(ABintersect,At.col)
    Bref = @spawn searchsortedmapping(ABintersect,Bt.row)
    AintMap = fetch(Aref)
    BintMap = fetch(Bref)
#    AintMap, BintMap = sortedintersectmapping(At.col,Bt.row)

#    AintMap,BintMap = map(x -> searchsortedmapping(x[1],x[2]),[(ABintersect,At.col) (ABintersect, Bt.row)])
    Aref = @spawn At.A[:,AintMap]
    Bref = @spawn Bt.A[BintMap,:]
    AA = fetch(Aref)
    BB = fetch(Bref)

    ABA = AA*BB

    AB = Assoc(At.row,Bt.col,Array{Union{AbstractString,Number}}([1.0]),ABA)

    return condense(AB)
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
function IodLambert(arg0::jdouble)
    return IodLambert((jdouble,), arg0)
end

function estimate(obj::IodLambert, arg0::Frame, arg1::jboolean, arg2::jint, arg3::Vector3D, arg4::AbsoluteDate, arg5::Vector3D, arg6::AbsoluteDate)
    return jcall(obj, "estimate", KeplerianOrbit, (Frame, jboolean, jint, Vector3D, AbsoluteDate, Vector3D, AbsoluteDate), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end


text:
"""
    coast(cmd0::String=""; kwargs...)

Plot continents, shorelines, rivers, and borders on maps.
Plots grayshaded, colored, or textured land-masses [or water-masses] on
maps and [optionally] draws coastlines, rivers, and political
boundaries. A map projection must be supplied.

Full option list at [`coast`]($(GMTdoc)coast.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- **A** | **area** :: [Type => Str or Number]

    Features with an area smaller than min_area in km^2 or of
    hierarchical level that is lower than min_level or higher than
    max_level will not be plotted.
    ($(GMTdoc)coast.html#a)
- $(GMT.opt_B)
- **C** | **river_fill** :: [Type => Str]

    Set the shade, color, or pattern for lakes and river-lakes.
    ($(GMTdoc)coast.html#c)
- **D** | **res** | **resolution** :: [Type => Str]		``Arg = c|l|i|h|f|a``

    Selects the resolution of the data set to use ((f)ull, (h)igh, (i)ntermediate, (l)ow, (c)rude), or (a)uto).
    ($(GMTdoc)coast.html#d)
- **E** | **DCW** :: [Type => Str]

    Select painting or dumping country polygons from the Digital Chart of the World.
    ($(GMTdoc)coast.html#e)
    + Tuple("code", Str); Tuple("code" [,"fill"], (pen)); Tuple((...),(...),...)
    + ex: ("PT",(0.5,"red","--")); (("PT","gblue",(0.5,"red"),("ES",(0.5,"yellow")))
- **F** | **box** :: [Type => Str]

    Draws a rectangular border around the map scale or rose.
    ($(GMTdoc)coast.html#f)
- **G** | **land** :: [Type => Str]

    Select filling or clipping of “dry” areas.
    ($(GMTdoc)coast.html#g)
- **I** | **rivers** :: [Type => Str]

    Draw rivers. Specify the type of rivers and [optionally] append pen attributes.
    ($(GMTdoc)coast.html#i)
- **L** | **map_scale** :: [Type => Str]

    Draw a map scale.
    ($(GMTdoc)coast.html#l)
- **M** | **dump** :: [Type => Str]

    Dumps a single multisegment ASCII output. No plotting occurs.
    ($(GMTdoc)coast.html#m)
- **N** | **borders** :: [Type => Str]

    Draw political boundaries. Specify the type of boundary and [optionally] append pen attributes
    ($(GMTdoc)coast.html#n)
- $(GMT.opt_P)
- **clip** :: [Type => Str]		``Arg = land|water|end``

    To clip land do *clip=:land*, *clip=:water* clips water. Use *end* to mark end of existing clip path.
    No projection information is needed.
    ($(GMTdoc)coast.html#q)
- **S** | **water** | **ocean** :: [Type => Str]

    Select filling or clipping of “wet” areas.
    ($(GMTdoc)coast.html#s)
- **Td** | **rose`** :: [Type => Str]

    Draws a map directional rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)coast.html#t)
- **Tm** | **compass** :: [Type => Str]

    Draws a map magnetic rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)coast.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **shore** :: [Type => Str]
    Draw shorelines [Default is no shorelines]. Append pen attributes.
    ($(GMTdoc)coast.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function coast(cmd0::String=""; clip=nothing, first=true, kwargs...)

	length(kwargs) == 0 && clip === nothing && return monolitic("pscoast", cmd0)

	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	maybe_more = false				# If latter set to true, search for lc & lc pen settings
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12cd/0")
	cmd = parse_common_opts(d, cmd, [:F :JZ :UVXY :bo :c :p :t :params], first)
	cmd = auto_JZ(cmd)		# Add -JZ if perspective for the case -R.../z_min/z_max
	cmd = parse_these_opts(cmd, d, [[:A :area], [:C :river_fill], [:D :res :resolution], [:M :dump]])
	#cmd = parse_TdTmL(d, cmd)
    cmd = parse_type_anchor(d, cmd, [[:Td :rose], [:Tm :compass], [:L :map_scale]])
	cmd = add_opt_fill(cmd, d, [:G :land], 'G')
	cmd = add_opt_fill(cmd, d, [:S :water :ocean], 'S')

	if (clip !== nothing)
		clip = string(clip)
		if     (clip == "land")    cmd *= " -Gc"
		elseif (clip == "water" || clip == "ocean") cmd *= " -Sc"
		elseif (clip == "end")     cmd *= " -Q"
		else
			@warn("The 'clip' argument can only be a string with 'land', 'water' or 'end'. Ignoring it.")
		end
	end

	# Parse these three options that can be made to respond to same code`
	symbs = [[:I :rivers], [:N :borders], [:W :shore]];	flags ="INW"
	for k = 1:3
		if ((val = find_in_dict(d, symbs[k])[1]) !== nothing)
			if (isa(val, NamedTuple) || (isa(val, Tuple) && isa(val[1], NamedTuple)))
				cmd = add_opt(cmd, flags[k], d, symbs[k], (type="/#", level="/#", pen=("", add_opt_pen)))
			elseif (isa(val, Tuple))  cmd *= " -" * flags[k] * parse_pen(val)
			else                      cmd *= " -" * flags[k] * arg2str(val)	# Includes Str, Number or Symb
			end
		end
	end

	if ((val = find_in_dict(d, [:E :DCW])[1]) !== nothing)
		if (isa(val, String) || isa(val, Symbol))
			cmd = string(cmd, " -E", val)			# Simple case, ex E="PT,+gblue"
		elseif (isa(val, NamedTuple))
			cmd = add_opt(cmd, "E", d, [:DCW :E], (country="", name="", continent="=",
			                                       pen=("+p", add_opt_pen), fill=("+g", add_opt_fill)))
		elseif (isa(val, Tuple))
			cmd = parse_dcw(cmd, val)
		end
	end

	if (!occursin("-C",cmd) && !occursin("-E",cmd) && !occursin("-G",cmd) && !occursin("-I",cmd) &&
		!occursin("-M",cmd) && !occursin("-N",cmd) && !occursin("-Q",cmd) && !occursin("-S",cmd) && !occursin("-W",cmd))
		cmd *= " -W0.5p"
	end
	if (!occursin("-D",cmd))  cmd *= " -Da"  end		# Then pick automatic
	finish = !occursin("-M ",cmd) ? true : false		# Otherwise the dump would be redirected to GMTjl_tmp.ps

	return finish_PS_module(d, "pscoast " * cmd, "", K, O, finish)
end

# ---------------------------------------------------------------------------------------------------
function parse_dcw(cmd::String, val::Tuple)
	for k = 1:length(val)
		if (isa(val[k], NamedTuple))
			cmd *= add_opt("", "E", Dict(:DCW => val[k]), [:DCW],
				(country="", name="", continent="=", pen=("+p", add_opt_pen), fill=("+g", add_opt_fill)))
		elseif (isa(val[k], Tuple))
			cmd *= parse_dcw(val[k])
		else
			cmd *= parse_dcw(val)
			break
		end
	end
	return cmd
end

function parse_dcw(val::Tuple)
    t = string("", " -E", val[1])
	if (length(val) > 1)
		if (isa(val[2], Tuple))  t *= "+p" * parse_pen(val[2])
		else                     t *= string(val[2])
		end
		if (length(val) > 2)
			if (isa(val[3], Tuple))  t *= add_opt_fill("+g", Dict(fill => val[3]), [:fill])
			else                     t *= string(val[3])
			end
		end
	end
	return t
end

# ---------------------------------------------------------------------------------------------------
coast!(cmd0::String=""; clip=nothing, first=false, kw...) = coast(cmd0; clip=clip, first=first, kw...)

const pscoast  = coast			# Alias for GMT5
const pscoast! = coast!			# Alias for GMT5


text:
using PowerModels
using Ipopt
using DataFrames
using CSV

filename = ARGS[1]

network_data = PowerModels.parse_file(filename)
result = run_ac_opf(network_data,IpoptSolver())

gen_nums = [parse(Int64,i) for (i,gen) in result["solution"]["gen"]]
pg = [gen["pg"] for (i,gen) in result["solution"]["gen"]]
qg = [gen["qg"] for (i,gen) in result["solution"]["gen"]]

df = DataFrame()
df[:gen_num] = gen_nums
df[:active_power] = pg
df[:reactive_power] = qg

CSV.write("output.csv",df)

max_stars_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_issues_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_forks_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    pix_read_gif(
        filename::AbstractString
    )::Union{Pix, Nothing}

Read a GIF image from the specified file.  Returns `nothing` on error.

__Arguments:__

| T | Name     | Default | Description
|---| :------- | ------- | -----------
| R | filename |         | The name of the GIF file to load.
"""
function pix_read_gif(
            filename::AbstractString
        )::Union{Pix, Nothing}
    local retval = nothing

    try
        local data = read(filename)
        retval = pix_read_gif(data)
    catch ex
        if isa(ex, SystemError) == false
            rethrow(ex)
        end
        retval = nothing
    end

    return retval
end

# =================================================================================================
"""
    pix_read_gif(
        stream::IO
    )::Union{Pix, Nothing}

Read a GIF image from the specified stream.  Returns `nothing` on error.

__Arguments:__

| T | Name   | Default  | Description
|:--| :----- | :------- | :----------
| R | stream |          | The IO stream to read the GIF file from.

__Details:__

This implementation mirrors the API provided by Leptonica when you pass in a FILE pointer. This
assumes that the remainder of the stream contains a GIF image.
"""
function pix_read_gif(
            stream::IO
        )::Union{Pix, Nothing}
    local data = read(stream)
    return pix_read_gif(data)
end

# =================================================================================================
"""
    pix_read_gif(
        data::AbstractArray{UInt8}
    )::Union{Pix, Nothing}

Read a GIF image from the byte array.  Returns `nothing` on error.

__Arguments:__

| T | Name | Default  | Description
|:--| :--- | :------- | :----------
| R | data |          | The byte array to read the GIF image from.
"""
function pix_read_gif(
            data::AbstractArray{UInt8}
        )::Union{Pix, Nothing}
    local ptr = ccall(
        (:pixReadMemGif, LEPTONICA),
        Ptr{Cvoid},
        (Ptr{Cvoid}, Csize_t),
        data,
        length(data)
    )

    if ptr == C_NULL
        return nothing
    end

    return Pix(ptr)
end

# =================================================================================================
"""
    pix_write_gif_i(
        pix::Pix
    )::Tuple{Ptr{UInt8}, Csize_t}

Write an image to memory in the GIF image format.  If there is an error `(C_NULL, 0)` is returned.

__Arguments:__

| T | Name     | Default  | Description
|:--| :------- | :------- | :----------
| R | pix      |          | The image to write to memory.

__Details:__

**This method is meant for internal use only.**  It returns the memory pointer and number of bytes
allocated by the Leptonica library to write out the image.  The caller is responsible for freeing
this memory with a call to `lept_free()`.
"""
function pix_write_gif_i(
            pix::Pix
        )::Tuple{Ptr{UInt8}, Csize_t}
    local temp = Ref(Ptr{UInt8}(C_NULL))
    local size = Ref(Csize_t(0))

    if is_valid(pix) == false
        @error "Pix has been freed."
        return (C_NULL, 0)
    end

    local retval = ccall(
        (:pixWriteMemGif, LEPTONICA),
        Cint,
        (Ptr{Ptr{UInt8}}, Ptr{Csize_t}, Ptr{Cvoid}),
        temp,
        size,
        pix.ptr
    )

    if retval == 1
        return (C_NULL, 0)
    end

    return (temp[], size[])
end

# =================================================================================================
"""
    pix_write_gif(
        filename::AbstractString,
        pix::Pix
    )::Bool

Write an image to a file in the GIF image format.  If there is an error `false` is returned.

__Arguments:__

| T | Name     | Default  | Description
|:--| :------- | :------- | :----------
| R | filename |          | The name of the file to write to.
| R | pix      |          | The image to write to the file.

__Details:__

If the file exists it will be overwritten.
"""
function pix_write_gif(
            filename::AbstractString,
            pix::Pix
        )::Bool
    local result = false
    local file = nothing
    local data, size = pix_write_gif_i(pix)

    if data != C_NULL
        try
            if size > 0
                file = open(filename, create=true, write=true)
                unsafe_write(file, data, size)
                result = true
            end
        catch ex
            if isa(ex, SystemError) == false
                rethrow(ex)
            end
        finally
            lept_free(data)

            if file != nothing
                close(file)
            end
        end
    end
    return result
end

# =================================================================================================
"""
    pix_write_gif(
        stream::IO,
        pix::Pix
    )::Bool

Write an image to an IO stream in the GIF image format.  If there is an error `false` is returned.

__Arguments:__

| T | Name   | Default  | Description
|:--| :----- | :------- | :----------
| R | stream |          | The IO stream to write to.
| R | pix    |          | The image to write to the stream.
"""
function pix_write_gif(
            stream::IO,
            pix::Pix
        )::Bool
    local result = false
    local data, size = pix_write_gif_i(pix)

    if data != C_NULL
        try
            if size > 0
                unsafe_write(stream, data, size)
                result = true
            end
        finally
            lept_free(data)
        end
    end

    return result
end

# =================================================================================================
"""
    pix_write_gif(
        pix::Pix
    )::Union{Vector{UInt8}, Nothing}

Write an image to a byte array in the GIF image format.  If there is an error `nothing` is returned.

__Arguments:__

| T | Name | Default  | Description
|:--| :--- | :------- | :----------
| R | pix  |          | The image to write to the byte array.
"""
function pix_write_gif(
            pix::Pix
        )::Union{Vector{UInt8}, Nothing}
    local output     = nothing
    local data, size = pix_write_gif_i(pix)

    if data != C_NULL
        try
            if size > 0
                output = Vector{UInt8}(undef, size)
                unsafe_copyto!(pointer(output), data, size)
            end
        finally
            lept_free(data)
        end
    end

    return output
end


text:
@json struct Project
    key::String
    links::NamedTuple
    name::String
    uuid::UUID
end
@json_struct Project

@json struct Branch
    name::String
end
@json_struct Branch

@json struct Workspace
    links::NamedTuple
    name::String
    slug::String
    uuid::UUID
end
@json_struct Workspace

@json mutable struct Repo
    created_on::DateTime
    description::String
    fork_policy::String
    full_name::String
    has_issues::Bool
    has_wiki::Bool
    is_private::Bool
    language::String
    links::NamedTuple
    mainbranch::Branch
    name::String
    owner::User
    project::Project
    size::Int
    slug::String
    updated_on::DateTime
    uuid::UUID
    website::String
    workspace::Workspace
end
@json_struct Repo

@not_implemented(::BitbucketAPI, ::typeof(get_user_repos))
@not_implemented(::BitbucketAPI, ::typeof(get_user_repos), ::String)
@not_implemented(::BitbucketAPI, ::typeof(get_user_repos), ::UUID)
into(::BitbucketAPI, ::typeof(get_user_repos)) = Vector{Repo}

endpoint(api::BitbucketAPI, ::typeof(get_repo), workspace::AStr, repo::AStr) =
    Endpoint(:GET, "/repositories/$workspace/$repo")
@not_implemented(::BitbucketAPI, ::typeof(get_repo), ::String, ::String, ::String)
into(::BitbucketAPI, ::typeof(get_repo)) = Repo

endpoint(::BitbucketAPI, ::typeof(create_repo), workspace::AStr, repo::AStr) =
    Endpoint(:POST, "/repositories/$workspace/$repo")
@not_implemented(::BitbucketAPI, ::typeof(create_repo), org::AStr)
@not_implemented(::BitbucketAPI, ::typeof(create_repo), ::UUID)
into(::BitbucketAPI, ::typeof(create_repo)) = Repo

function is_bitbucket_collaborator(resp::HTTP.Response)
    if HTTP.status(resp) == 200
        d = JSON2.read(String(HTTP.body(resp)), Dict)
        return haskey(d, "values") && length(d["values"]) > 0 &&
            d["values"][1]["permissions"] in ["write" "admin"]
    end
    false
end

endpoint(::BitbucketAPI, ::typeof(is_collaborator), wrksp::AStr, repo::AStr) =
    Endpoint(:GET, "/user/permissions/repositories";
             allow_404=true,
             query=Dict(
                 :q=> "repository.full_name=\"$wrksp/$repo\""
             ),
             )
endpoint(::BitbucketAPI, ::typeof(is_collaborator), wrksp::AStr, repo::AStr, user::AStr) =
    Endpoint(:GET, "/workspaces/$wrksp/permissions/repositories/$repo";
             allow_404=true,
             query=Dict(
                 :q=> "user.nickname=\"$user\""
             )
             )
endpoint(::BitbucketAPI, ::typeof(is_collaborator), wrksp::String, repo::String, user::UUID) =
    Endpoint(:GET, "/workspaces/$wrksp/permissions/repositories/$repo";
             allow_404=true,
             query=Dict(
                 :q=> "user.uuid=\"{$user}\""
             )
             )
postprocessor(::BitbucketAPI, ::typeof(is_collaborator)) = DoSomething(is_bitbucket_collaborator)
into(::BitbucketAPI, ::typeof(is_collaborator)) = Bool

endpoint(api::BitbucketAPI, ::typeof(get_file_contents), wrksp::AStr, repo::AStr, path::AStr) =
    Endpoint(:GET, "/repositories/$wrksp/$repo/src/$path")
@not_implemented(::BitbucketAPI, ::typeof(get_file_contents), repo::UUID, path::AStr)
postprocessor(::BitbucketAPI, ::typeof(get_file_contents)) = DoSomething() do req
    String(req.body)
end
into(::BitbucketAPI, ::typeof(get_file_contents)) = String

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

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

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

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

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

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


text:
"""
A graph to test greedy multi path routing
"""

greedy_test = QNetwork()
A = BasicNode("A")
B = BasicNode("B")
M1 = BasicNode("M1")
M2 = BasicNode("M2")
M3 = BasicNode("M3")

AM1 = BasicChannel(A, M1, exp_cost=false)
AM2 = BasicChannel(A, M2, exp_cost=false)
AM3 = BasicChannel(A, M3, exp_cost=false)
M1B = BasicChannel(M1, B, exp_cost=false)
M2B = BasicChannel(M2, B, exp_cost=false)
M3B = BasicChannel(M3, B, exp_cost=false)

AM1.costs = Dict("loss"=>0.5, "Z"=>0.5)
AM2.costs = Dict("loss"=>1.0, "Z"=>1.0)
AM3.costs = Dict("loss"=>1.5, "Z"=>1.5)
M1B.costs = Dict("loss"=>0.5, "Z"=>0.5)
M2B.costs = Dict("loss"=>1.0, "Z"=>1.0)
M3B.costs = Dict("loss"=>1.5, "Z"=>1.5)

for i in [A, B, M1, M2, M3, AM1, AM2, AM3, M1B, M2B, M3B]
    add(greedy_test, i)
end

refresh_graph!(greedy_test)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
function sparse_to_csr_test()
    i = [1;2;3]
    j = [3;4;4]
    v = [8;9;10]
    (rp,ci,ai,m) = sparse_to_csr(i,j,v)

    # more tests added here
    for t = 1:100
        A = sprand(100,80,0.01)
        (rp,ci,ai,m)=sparse_to_csr(A)
        i = zeros(Int64,length(ai))
        j = zeros(Int64,length(ai))
        a = zeros(Float64,length(ai))
        n = length(rp)-1
        nz = 0
        for cr = 1:n
            for ri = rp[cr]:rp[cr+1]-1
                nz=nz+1
                i[nz]=cr
                j[nz]=ci[ri]
                a[nz]=ai[ri]
            end
        end
        A2 = sparse(i,j,a,n,m)
        if ~isequal(A,A2)
            error("sparse_to_csr_test failed")
        end
    end
    return true
end


text:
export TwoStreamInstability

struct TwoStreamInstability
    grid :: OneDGrid
end

function sample!(pg::ParticleGroup{1,1}, d::TwoStreamInstability)

    np = pg.n_particles
    xmin = d.grid.xmin
    xmax = d.grid.xmax
    dx = d.grid.dx

    wp = +1.0  # plasma frequency
    qm = -1.0  # charge/mass ratio
    qp = wp^2 / (qm * np / (xmax - xmin))
    xp = collect(LinRange(xmin, xmax, np + 1))[1:end-1]  # Particle positions
    v0 = 0.9       # Stream velocity
    vt = 0.0000001 # Thermal speed

    # Particle momentum, initially Maxwellian
    vp = vt .* (1 .- vt .^ 2) .^ (-0.5) .* randn(rng, Float64, np)
    pm = collect(0:np-1)
    pm = 1 .- 2 * mod.(pm .+ 1, 2)
    vp .+= pm .* (v0 * (1 - v0^2)^(-0.5)) # Momentum + stream velocity

    # Add electron perturbation to excite the desired mode perturbation
    xp1 = 1.0
    mode = 1
    xp .+= xp1 * dx * sin.(2π .* xp ./ (xmax - xmin) .* mode)
    xp .= mod.(xp, xmax - xmin)

    pg.array[:,1] .= xp
    pg.array[:,2] .= vp

end


text:
function parseGtsLine( s::AbstractString, C, T=eltype(C) )
    firstInd = findfirst( isequal(' '), s )
    secondInd = findnext( isequal(' '), s, firstInd+1 )
    firstNum = parse( T, s[1:firstInd] )
    if secondInd != nothing
        secondNum = parse( T, s[firstInd:secondInd] )
        thirdNum = parse( T, s[secondInd:end] )
        return C([firstNum, secondNum, thirdNum])
    else
        secondNum = parse( T, s[firstInd:end] )
        return C([firstNum, secondNum])
    end
end

function load( st::Stream{format"GTS"}, MeshType=GLNormalMesh )
    io = stream(st)
    head = readline( io )
    FT = facetype(MeshType)
    VT = vertextype(MeshType)

    nVertices, nEdges, nFacets = parseGtsLine( head, Tuple{Int,Int,Int} )
    iV = iE = iF = 1
    vertices = Vector{VT}(undef, nVertices)
    edges = Vector{Vector{Int}}(undef, nEdges)
    facets = Vector{Vector{Int}}(undef, nFacets)
    for full_line::String in eachline(io)
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(full_line))
        !isascii(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !all(iscntrl, line) #ignore comments
            if iV <= nVertices
                vertices[iV] = parseGtsLine( line, VT )
                iV += 1
            elseif iV > nVertices && iE <= nEdges
                edges[iE] = parseGtsLine( line, Array{Int} )
                iE += 1
            elseif iE > nEdges && iF <= nFacets
                facets[iF] = parseGtsLine( line, Array{Int} )
                iF += 1
            end # if
        end # if
    end # for
    faces = [ FT( union( edges[facets[i][1]], edges[facets[i][2]], edges[facets[i][3]] ) ) for i in 1:length(facets) ]  # orientation not guaranteed
    return MeshType( vertices, faces )
end

function save( st::Stream{format"GTS"}, mesh::AbstractMesh )
    # convert faces to edges and facets
    edges = [[ 0, 0 ]] # TODO
    facets = [[ 0, 0 ]] # TODO
    # write to file
    io = stream( st )
    println( io, length(mesh.vertices), legth(edges), length(mesh.faces) )
    for v in mesh.vertices
        println( io, v[1], v[2], v[3] )
    end
    for e in edges
        println( io, e[1], e[2] )
    end
    for f in facets
        println( io, f[1], f[2], f[3] )
    end
end

max_stars_repo_name: JuliaIO/MeshIO

max_issues_repo_name: JuliaIO/MeshIO

max_forks_repo_name: JuliaIO/MeshIO


text:
v06 = v"0.6.2"
v07 = v"0.7.0-alpha.0"

supported() = VERSION >= v07 ||
             (VERSION >= v06 && VERSION < v"0.7.0-DEV")

compat_warn() = warn("""
    HTTP.jl has not been tested with Julia version $VERSION
    Supported versions are $v06 and $v07 and above.
    """)

__init__() = supported() || compat_warn()

@static if VERSION >= v07

    using Base64
    import Dates

    const bytesavailable = Base.bytesavailable
    const compat_findfirst = Base.findfirst
    const compat_replace = Base.replace
    const compat_occursin = Base.occursin
    const compat_parse = Base.parse
    const compat_string = Base.string

    compat_stdout() = stdout

    compat_search(s::AbstractString, c::Char) = Base.findfirst(isequal(c), s)
    using Sockets

    sprintcompact(x) = sprint(show, x; context=:compact => true)

else

    supported() || compat_warn()

    eval(:(module Base64 end))
    const Dates = Base.Dates

    const bytesavailable = Base.nb_available
    compat_findfirst(a...) = (r = Base.findfirst(a...); r == 0 ? nothing : r)
    compat_replace(s, p) = Base.replace(s, p.first, p.second)
    compat_occursin(needle, haystack) = contains(haystack, needle)
    compat_parse(s, T; base::Int=10) = Base.parse(s, T, base)
    compat_string(s; base::Int=16, pad::Int=0) = hex(s, pad)

    compat_stdout() = STDOUT

    compat_search(s::AbstractString, c::Char) = compat_findfirst(s, c)

    const Nothing = Void
    const isnumeric = isnumber
    const isletter = isalpha

    Base.SubString(s) = SubString(s, 1)
    Base.String(x::SubArray{UInt8,1}) = String(Vector{UInt8}(x))

    macro debug(s) DEBUG_LEVEL > 0 ? :(("D- ", $(esc(s)))) : :() end
    macro info(s)  DEBUG_LEVEL >= 0 ? :(println("I- ", $(esc(s)))) : :() end
    macro warn(s)  DEBUG_LEVEL >= 0 ? :(println("W- ", $(esc(s)))) : :() end
    macro error(m, args...)
        args = [:(print("|  "); @show $a) for a in args]
        esc(:(println("E- ", $m); $(args...); nothing))
    end
    if !isdefined(Main, :Compat)
        Base.fetch(t::Task) = wait(t)
    end

    eval(:(module Sockets
        export TCPSocket
        import Base: TCPSocket, TCPServer, IPAddr, @ip_str, DNSError,
            getsockname, getaddrinfo, connect, listen, DNSError, InetAddr
        const localhost = ip"127.0.0.1"
        end))
    using .Sockets

    sprintcompact(x) = sprint(showcompact, x)
end

#https://github.com/JuliaWeb/MbedTLS.jl/issues/122
@static if !applicable(bytesavailable, MbedTLS.SSLContext())
    Base.bytesavailable(ssl::MbedTLS.SSLContext) = nb_available(ssl)
end

@static if !applicable(Sockets.getsockname, MbedTLS.SSLContext())
    Sockets.getsockname(ctx::SSLContext) = Sockets.getsockname(ctx.bio)
end


text:
# Autogenerated wrapper script for ModularForms_jll for i686-linux-gnu
export libmodularforms_jll

using Arb_jll
using FLINT_jll
JLLWrappers.@generate_wrapper_header("ModularForms")
JLLWrappers.@declare_library_product(libmodularforms_jll, "libmodularforms_jll.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Arb_jll, FLINT_jll)
    JLLWrappers.@init_library_product(
        libmodularforms_jll,
        "lib/libmodularforms_jll.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StatsBase                 # just for sample()

# testing a-priori algorithm
transactions = [sample(1:10, 5, replace = false) for x in 1:100_000]
fk = frequent(transactions, 0.1)

groceries = ["milk", "bread", "eggs", "apples", "oranges", "beer"]
transactions = [sample(groceries, 4, replace = false) for x in 1:1000]
fk = frequent(transactions, 0.1)

rules1 = apriori(transactions, 0.1, 0.4, false);        # false for only single-item consequents
display(rules1)

# reading data from .csv
adult_data = dataset("adult")
adult_trans = make_transactions(adult_data[1:1000, :])  # take only sub-set of data for convenience
@time rules2 = apriori(adult_trans, 0.2, 0.8);          # takes about 26 sec.
display(rules2)

# testing SPADE algorithm
zaki_data = dataset("zaki_data")
seqs = make_sequences(zaki_data, sid_col = 2, eid_col = 3, item_col = 1)
@time res = spade(seqs, 0.2, 20);

@assert length(res[1]) == 8
@assert length(res[2]) == 57
@assert length(res[3]) == 191
@assert length(res[4]) == 444
@assert length(res[5]) == 743
@assert length(res[6]) == 897
@assert length(res[7]) == 784
@assert length(res[8]) == 494
@assert length(res[9]) == 220
@assert length(res[10]) == 66
@assert length(res[11]) == 12
@assert length(res[12]) == 1

# testing SPADE algorithm with subset of Zaki data
zaki_data = dataset("zaki_subset_data")
seqs2 = make_sequences(zaki_data, sid_col = 2, eid_col = 3, item_col = 1)
@time res2 = spade(seqs2, 0.2, 6);

# test sequential-rule building algorithm
seqrules = sequential_rules(res, 0.01, 12)
@assert length(seqrules) == 3886

seqrules2 = sequential_rules(res2, 0.01, 5)
@assert length(seqrules2) == 28

# testing SPADE algorithm with subset of Zaki data
zaki_data = readcsv("../data/zaki_data_extra.csv", skipstart = 1)
seqs3 = make_sequences(zaki_data, sid_col = 2, eid_col = 3, item_col = 1)
@time res3 = spade(seqs3, 0.2, 5);

@code_warntype spade!(res3[3], Array{IDList,1}[], 20, 0.2)

# test sequential-rule building algorithm
@time seqrules3 = sequential_rules(res3, 0.01, 5);
length(seqrules3)

letters = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N"]
n = 1000
nseq = 20
ntime = 50
zaki_rep_data = hcat(rand(letters, n), rand(collect(1:nseq), n), rand(collect(1:ntime), n))

seqs4 = make_sequences(zaki_rep_data, sid_col = 2, eid_col = 3, item_col = 1)
@time res4 = spade(seqs4, 0.2, 5);

function benchmark1(n)
    nseq = 10
    ntime = 10
    maxdepth = 10
    letters = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T"]

    dat = hcat(rand(letters, n), rand(collect(1:nseq), n), rand(collect(1:ntime), n))
    seqs = make_sequences(dat, sid_col = 2, eid_col = 3, item_col = 1);
    @time x = spade(seqs, 0.2, maxdepth);
    return nothing
end

benchmark1(100)


text:
function gmv(mus::DataFrame, covMatr::DataFrame)
    ## global minimum variance portfolio
    nAss = size(mus, 2)

    oneArr = ones(nAss)
    invCov = inv(array(covMatr))
    wgts = 1/(oneArr'*invCov*oneArr).*invCov*oneArr

    return wgts'
end

function gmv(mus::Array{Float64, 2}, covMatr::Array{Float64, 2})
    ## global minimum variance portfolio
    nAss = size(mus, 2)

    oneArr = ones(nAss)
    invCov = inv(covMatr)
    wgts = 1./(oneArr'*invCov*oneArr).*invCov*oneArr

    return wgts'
end

function maxSharpeRatio(mus::Array{Float64, 2},
                        covMatr::Array{Float64, 2},
                        r::Float64)
    ## calculate portfolio weights that maximize Sharpe ratio
    ##
    ## Inputs:
    ## 	mus			1 x nAss Array{Float64, 2} of asset mus
    ## 	covMatr		nAss x nAss Array{Float64, 2} of asset
    ## 					covariances
    ## 	r				Float64 for annualized interest rate
    ##
    ## Outputs:
    ## 	1 x nAss Array{Float64, 2} of optimal portfolio weights

    ## transform to daily interest rate
    intRate = (1 + r)^(1/250) - 1

    ## number of assets
    nAss = size(mus, 2)

    ## Array{Float64, 2} expected daily excess returns
    muExcess = (mus .- intRate)[:]

    oneArr = ones(nAss)
    invCov = inv(covMatr)
    wgts = 1./(oneArr'*invCov*muExcess).*invCov*muExcess

    return wgts'
end

## global minimum variance, no short-selling
function gmvNoSS(mus::Array{Float64, 2}, covMatr::Array{Float64, 2})
    ## numerical calculation of gmv without shortselling
    nAss = length(mus)
    env = Gurobi.Env()
    setparams!(env; IterationLimit=100, Method=1)
    minVarModel = gurobi_model(env;
                               name = "minimumVariance",
                               H = covMatr,
                               f = zeros(nAss),
                               Aeq = ones(1, nAss),
                               beq = [1.],
                               lb = zeros(nAss),
                               ub = ones(nAss))
    optimize(minVarModel)
    wgts = get_solution(minVarModel)
    return wgts'
end


text:
using Agents

@agent Particle ContinuousAgent{3} begin
    radius::Float64
end

function initialize_model(;
    N = 100, space_extents = (10.0, 10.0, 10.0),
    )
    space = ContinuousSpace(space_extents, 1.0; periodic = true)
    model = ABM(Particle, space)
    for i in 1:N
        particle = Particle(i, random_position(model), (0.0, 0.0, 0.0), rand())
        add_agent_pos!(particle, model)
    end
    return model
end

model = initialize_model()

using InteractiveDynamics, GLMakie

as(agent) = 0.1agent.radius + 0.1
fig, step = abm_plot(model; as)
fig


text:
using DataTables
using Test

@testset "DataTables" begin

    dt1 = DataTable(a=[1,2,3], b=[4.,5.,6.], c=["John", "Sally", "Jim"])

    @testset "Core" begin

        @test Base.IndexStyle(DataTable) isa Base.IndexLinear
        @test Base.IndexStyle(typeof(dt1)) isa Base.IndexLinear

        @test length(dt1) == 3
        @test dt1.a == [1,2,3]
        @test dt1.b == [4.,5.,6.]
        @test dt1.c == ["John", "Sally", "Jim"]
        @test dt1[1] == (a = 1, b = 4., c = "John")
        @test dt1[2] == (a = 2, b = 5., c = "Sally")
        @test dt1[3] == (a = 3, b = 6., c = "Jim")

        @test collect(dt1) == [(a = 1, b = 4., c = "John"), (a = 2, b = 5., c = "Sally"), (a = 3, b = 6., c = "Jim")]

        dt2 = DataTable([(a = 1, b = 4., c = "John"), (a = 2, b = 5., c = "Sally"), (a = 3, b = 6., c = "Jim")])

        @test dt1 == dt2

    end

    @testset "show" begin

        @test sprint(show, dt1) ==
    "3x3 DataTable\na │ b   │ c    \n──┼─────┼──────\n1 │ 4.0 │ John \n2 │ 5.0 │ Sally\n3 │ 6.0 │ Jim  "

        @test sprint((stream, data) -> show(stream, "text/plain", data), dt1) ==
    "3x3 DataTable\na │ b   │ c    \n──┼─────┼──────\n1 │ 4.0 │ John \n2 │ 5.0 │ Sally\n3 │ 6.0 │ Jim  "

        @test sprint((stream, data) -> show(stream, "text/html", data), dt1) ==
    "

a	b	c

1	4.0	"John"

2	5.0	"Sally"

3	6.0	"Jim"

"

        @test sprint((stream, data) -> show(stream, "application/vnd.dataresource+json", data), dt1) ==
    "{\"schema\":{\"fields\":[{\"name\":\"a\",\"type\":\"integer\"},{\"name\":\"b\",\"type\":\"number\"},{\"name\":\"c\",\"type\":\"string\"}]},\"data\":[{\"a\":1,\"b\":4.0,\"c\":\"John\"},{\"a\":2,\"b\":5.0,\"c\":\"Sally\"},{\"a\":3,\"b\":6.0,\"c\":\"Jim\"}]}"

        @test showable("text/html", dt1) == true
        @test showable("application/vnd.dataresource+json", dt1) == true

    end

end


text:
const generated_functions = []

function load_generated_functions()
    for function_defn in generated_functions
        Core.eval(Main, function_defn)
    end
end

# built-in generator types
include("vector_trace.jl")
include("lightweight/lightweight.jl")
include("plate/plate.jl")
include("plate_of_dists/plate_of_dists.jl")
include("markov/markov.jl")
include("basic/basic.jl")
include("at_dynamic/at_dynamic.jl")

export load_generated_functions

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
function find_low_risk_path(chitons)
    cavern = CartesianIndices(chitons)
    moves = CartesianIndex.(((1, 0), (-1, 0), (0, 1), (0, -1)))
    start, finish = extrema(cavern)
    prev_position = Dict{CartesianIndex{2}, CartesianIndex{2}}()
    tentative_risk = fill(typemax(Int), size(chitons))
    tentative_risk[start] = 0
    queue = [(start, 0)]
    while !isempty(queue)
        position, risk = popat!(queue, last(findmin(last, queue)))
        if position == finish
            path = [finish]
            while first(path) != start
                pushfirst!(path, prev_position[first(path)])
            end
            return risk, path
        end
        if risk <= tentative_risk[position]
            tentative_risk[position] = risk
            for neighbor in [position + move for move in moves if position + move in cavern]
                neighbor_risk = risk + chitons[neighbor]
                if neighbor_risk < tentative_risk[neighbor]
                    tentative_risk[neighbor] = neighbor_risk
                    prev_position[neighbor] = position
                    push!(queue, (neighbor, neighbor_risk))
                end
            end
        end
    end
    nothing  # finish is not reachable
end

function expand_chitons(chitons)
    chitons = [
        chitons    chitons.+1 chitons.+2 chitons.+3 chitons.+4
        chitons.+1 chitons.+2 chitons.+3 chitons.+4 chitons.+5
        chitons.+2 chitons.+3 chitons.+4 chitons.+5 chitons.+6
        chitons.+3 chitons.+4 chitons.+5 chitons.+6 chitons.+7
        chitons.+4 chitons.+5 chitons.+6 chitons.+7 chitons.+8
    ]
    chitons .= mod1.(chitons, 9)
    chitons
end

chitons = collect.(readlines("day15/chitons.txt"))
chitons = parse.(Int, permutedims(hcat(chitons...)))

p1, _ = find_low_risk_path(chitons)

@show p1

p2, _ = find_low_risk_path(expand_chitons(chitons))

@show p2


text:
module NeutronTransport

using UnPack
using Setfield
using Reexport
using StaticArrays
@reexport using RayTracing
using Gridap: get_face_labeling
using Gridap.Geometry: get_face_tag, get_tag_from_name

import Base: show, eltype

abstract type TransportFormulation end

struct MethodOfCharacteristics <: TransportFormulation end
struct CollisionProbability <: TransportFormulation end
struct DiscreteOrdinates <: TransportFormulation end
struct Diffusion <: TransportFormulation end

const MoC = MethodOfCharacteristics
const CP = CollisionProbability
const SN = DiscreteOrdinates

include("cross_section.jl")
include("polar_quad.jl")
include("quadrature.jl")

abstract type TransportProblem{Dim,NRegions,NGroups} end

dimension(::TransportProblem{Dim}) where {Dim} = Dim
nregions(::TransportProblem{Dim,NRegions}) where {Dim,NRegions} = NRegions
ngroups(::TransportProblem{Dim,NRegions,NGroups}) where {Dim,NRegions,NGroups} = NGroups

include("mocproblem.jl")

abstract type TransportSolution end

include("mocsolver.jl")

export TabuchiYamamoto, GaussLegendre
export CrossSections
export MoCProblem
export solve

end


text:
using OnlineAI, FactCheck
import OnlineStats: @LOG

A = SigmoidMapping()
# A = TanhMapping()
C = CrossentropyCostModel()
# C = L2CostModel()

nin, nout = 3, 1
hidden = [2]
net = buildClassificationNet(nin, nout, hidden;
                              finalMapping = A,
                              params = NetParams(mloss = C)
                            )

x = rand(nin) - 0.5
y = [1.0]

ε = 1e-4
numgrads = []
ls = net.layers

for l in ls
  w = l.w
  numgrad = similar(w)  # numerical gradient := (C(wij + ε) - C(wij - ε)) / (2ε)
  for i in 1:nrows(w)
    for j in 1:ncols(w)
      wij = w[i,j]

      w[i,j] = wij + ε
      cp = cost(net, x, y)

      w[i,j] = wij - ε
      cm = cost(net, x, y)

      numgrad[i,j] = (cp - cm) / (2ε)

      # put wij back
      w[i,j] = wij
    end
  end
  push!(numgrads, numgrad)
end

yhat = forward!(net, x, true)
multiplyDerivative = OnlineAI.costMultiplier!(net.params.mloss, net.costmult, y, yhat)
OnlineAI.updateSensitivities!(ls[end], net.costmult, multiplyDerivative)
for i in length(ls)-1:-1:1
  OnlineAI.updateSensitivities!(ls[i:i+1]...)
end

δs = [l.δ for l in ls]
xs = [l.x for l in ls]
Σs = [l.Σ for l in ls]
ws = [l.w for l in ls]
bs = [l.b for l in ls]

@show x
@show y
@show yhat
@show net.costmult
@show multiplyDerivative

facts("gradient") do
  for i in 1:length(ls)
    println()
    @show i
    @show xs[i]
    @show δs[i]
    @show ws[i]
    @show bs[i]
    @show Σs[i]
    @show numgrads[i]
    @show δs[i] * xs[i]'

    @fact δs[i] * xs[i]' --> roughly(numgrads[i], rtol=1e-3)
  end
end

# now update the weights
for l in ls
  OnlineAI.updateWeights!(l, net.params.updater)
end

ws = [l.w for l in ls]
bs = [l.b for l in ls]

for i in 1:length(ls)
  println()
  @show i
  @show ws[i]
  @show bs[i]
end


text:
using NeuronBuilder, ModelingToolkit, OrdinaryDiffEq, Plots

const Cm = 10.0 # specific capacitance cₘ is a biological constant (around) 10 nF/mm^2
Liu_conv = Cm

channels = [Liu.Na(700.0 * Liu_conv), Liu.CaS(4.0 * Liu_conv), Liu.CaT(2.0 * Liu_conv), Liu.Ka(50.0 * Liu_conv), Liu.KCa(40.0 * Liu_conv),
    Liu.Kdr(70.0 * Liu_conv), Liu.H(0.03 * Liu_conv), Liu.leak(0.01 * Liu_conv)]

τmRNAs = [5e4 5e4 5e4 5e4 5e4 5e4 5e4]
τg = 1e4
Ca_tgt = 5.0

#average calcium at the start is 10, Ca_tgt is half as much so conductances have to decrease
reg_chs = [OLeary_reg(channel,taumrna,τg, Ca_tgt) for (channel, taumrna) in zip(channels[1:end-1],τmRNAs)]
Oreg_channels = [reg_chs... , channels[end]]

τCa = 20.0
Ca∞ = 0.05
fxarea = 14.96 * 0.0628

defaults = Dict(Voltage => BasicVoltageDynamics(),
    Calcium => Liu.CalciumDynamics(τCa, Ca∞, fxarea),
    Reversal{Calcium} => Liu.CaReversalDynamics())

somatic_parameters = Dict(
    Reversal{Sodium} => 50.0,
    Reversal{Potassium} => -80.0,
    Reversal{Leak} => -50.0,
    Reversal{Proton} => -20.0,
    Voltage => -60.0,
    Calcium => 0.05,
    Reversal{Calcium} => 0.0)

b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, Oreg_channels, :regul_Liu)

neur = b()

prob = ODEProblem(neur, [], (0.0, 15000.0), [])

sol = @time solve(prob, AutoTsit5(Rosenbrock23()))

p1 = plot(sol, vars=[neur.V, neur.Ca])
p2 = plot(sol, vars=[neur.Na_regul₊mRNA, neur.CaS_regul₊mRNA, neur.CaT_regul₊mRNA, neur.Ka_regul₊mRNA, neur.KCa_regul₊mRNA, neur.Kdr_regul₊mRNA, neur.H_regul₊mRNA], title="mRNAs", titlefontsize=9)
p3 = plot(sol, vars=[neur.Na_regul₊gNa, neur.CaS_regul₊gCa, neur.CaT_regul₊gCa, neur.Ka_regul₊gK, neur.KCa_regul₊gK, neur.Kdr_regul₊gK, neur.H_regul₊gH], title="conductances", titlefontsize=9)
plot(p1,p2,p3,layout=(3,1), legend=false)


text:
# -------- Loss -------------------------------------------------------------- #

"""
Loss type that determines how the output of a NeuralNetwork model is compared
to policy, value, and feature labels.
"""
struct Loss

  value          :: NamedTuple{(:loss, :weight), Tuple{Function, Float32}}
  policy         :: NamedTuple{(:loss, :weight), Tuple{Function, Float32}}
  reg            :: NamedTuple{(:loss, :weight), Tuple{Function, Float32}}

  features :: Vector{Feature}
  fweights :: Vector{Float32}

end

# Default value loss
default_loss(x, y :: Union{Vector, Knet.KnetVector}) = sum(abs2, x .- y)

# Default policy loss
default_loss(x, y :: Union{Matrix, Knet.KnetMatrix}) = -sum(y .* log.(x.+1f-10))

# Default regularization loss for parameters
default_loss(param) = sum(abs2, param)

# Convert Constructor arguments
conv_loss_arg(w :: Real) = (loss = default_loss, weight = Float32(w))
conv_loss_arg(f :: Function) = (loss = f, weight = 1f0)
conv_loss_arg(t :: Tuple{Function, Real}) = (loss = t[1], weight = Float32(t[2]))
conv_loss_arg(nt) = nt

"""
    Loss([; value, policy, reg, features, feature_weights])

Construct a loss.

The arguments `value`, `policy`, and `reg` take either a float (the
respective weight value) or tuples `(f, w)` of a loss function and
a regularization value. The argument `features` is a list of supported features,
and `feature_weights` are the corresponding weights.
"""
function Loss(
             ; value = 1f0
             , policy = 1f0
             , reg = 0f0
             , features = Feature[]
             , feature_weights = ones(Float32, length(features)) )

  Loss( conv_loss_arg(value)
      , conv_loss_arg(policy)
      , conv_loss_arg(reg)
      , Feature[f for f in features]
      , feature_weights )

end

function caption(l :: Loss)
  base = [:value, :policy, :reg]
  features = feature_name.(l.features)
  vcat(base, features)
end

features(l :: Loss) = l.features

# -------- Loss Calculation -------------------------------------------------- #

function loss( l :: Loss
             , model :: NeuralModel{G, GPU}
             , cache :: DataCache{G, GPU}
             ) where {G, GPU}

  n = length(cache)
  use_features = !isnothing(cache.flabel)

  # Get the model output for value/policy labels and feature labels
  v, p, f = model(cache.data, use_features)

  # Calculate the value and policy loss
  vloss = l.value.weight * l.value.loss(v, cache.vlabel) / n
  ploss = l.policy.weight * l.policy.loss(p, cache.plabel) / n

  # Calculate the regularization loss
  rloss = l.reg.weight * sum(Knet.params(model)) do param
    s = size(param)
    maximum(s) < prod(s) ? l.reg.loss(param) : 0f0
  end

  # Calculate the feature losses
  if use_features

    feats = features(l)
    indices = feature_indices(feats, G)

    flosses = map(feats, indices) do feat, sel
      feature_loss(feat, f[sel,:], cache.flabel[sel,:]) / n
    end

  else

    flosses = zeros(Float32, length(features(l)))

  end

  # Combine all losses in one array
  vcat(vloss, ploss, rloss, flosses...)

end

"""
    loss(l, model, dataset [; maxbatch = 1024])

Calculate the loss determined by `l` for `model` on `dataset` while evaluating
model with at most `maxbatch` game states at once.
"""
function loss( l :: Loss
             , model :: NeuralModel{G, GPU}
             , dataset :: DataSet{G}
             ; maxbatch = 1024
             ) where {G, GPU}

  # Check if features can be used
  use_features = check_features(l, model, dataset)

  # Cut the dataset in batches if it is too large
  batches = Batches(dataset, maxbatch, gpu = GPU, use_features = use_features)

  # Get the total loss
  sum(batches) do cache
    loss(l, model, cache) .* length(cache)
  end ./ length(dataset)

end

"""
    loss(l, model, game, label)

Calculate the loss determined by `l` between `label` and `model(game)`
"""
function loss(l :: Loss, model :: Model, game, label)

  loss(l, model, DataSet([game], [label], [zeros(Float32, 0)]))

end

max_stars_repo_name: roSievers/jtac


text:
mutable struct Optimizer <: MOI.AbstractOptimizer
    form::DMIPFormulation
    obj_sense::MOI.OptimizationSense
    primal_bound::Float64
    config::AlgorithmConfig
    result::Union{Nothing,Result}

    function Optimizer(config::AlgorithmConfig = AlgorithmConfig())
        return new(
            DMIPFormulation(),
            MOI.FEASIBILITY_SENSE,
            NaN,
            config,
            nothing,
        )
    end
end

_is_max_sense(opt::Optimizer) = opt.obj_sense == MOI.MAX_SENSE

function _get_primal_bound(opt::Optimizer)
    return isnan(opt.primal_bound) ? Inf : opt.primal_bound
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Cerberus"
MOI.get(opt::Optimizer, ::MOI.RawSolver) = opt.form

# TODO: Figure out what to say here
function MOI.is_empty(opt::Optimizer)
    return isempty(opt.form) &&
           opt.obj_sense == MOI.FEASIBILITY_SENSE &&
           isnan(opt.primal_bound) &&
           opt.result === nothing
end

function MOI.empty!(opt::Optimizer)
    opt.form = DMIPFormulation()
    opt.obj_sense = MOI.FEASIBILITY_SENSE
    opt.primal_bound = NaN
    opt.result = nothing
    return nothing
end

function MOI.optimize!(opt::Optimizer)
    orig_obj = opt.form.obj
    if _is_max_sense(opt)
        opt.form.obj =
            CSAF(-orig_obj.coeffs, copy(orig_obj.indices), -orig_obj.constant)
    end
    opt.result = optimize!(opt.form, opt.config, _get_primal_bound(opt))
    if _is_max_sense(opt)
        opt.form.obj = orig_obj
    end
    return nothing
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(opt::Optimizer, ::MOI.Silent, value::Bool)
    opt.config.silent = value
    return nothing
end
MOI.get(opt::Optimizer, ::MOI.Silent) = opt.config.silent

MOIU.supports_default_copy_to(opt::Optimizer, copy_names::Bool) = !copy_names
function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(dest, src; kws...)
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(opt::Optimizer, ::MOI.TimeLimitSec, value::Real)
    opt.config.time_limit_sec = value
    return nothing
end
MOI.get(opt::Optimizer, ::MOI.TimeLimitSec) = opt.config.time_limit_sec


text:
@inline sqnorm(x) = sum(x.^2)

"""
pinvli{n,m,T}(a::Mat{n,m,T})

Compute pseudo inverse of matrix with linear independent columns.
"""
function pinvli(a::Mat{n,m,T}) where {n,m,T}
    @assert n >= m
    return inv(a'*a)*a'
end


text:
using NanoclusterModeler
using Test

@testset "NanoclusterModeler.jl" begin
    # Write your own tests here.
end


text:
@testset "linfree_00001.tif" begin
filename = joinpath(testbase, "linfree_00001.tif")
    # $ ScanImageTiffReader image shape linfree_00001.tif
    # Shape: 512 x 512 x 10 @ i16
    tsize = ScanImageTiffReader.open(filename) do io
        size(io)
    end
    @test tsize == (512, 512, 10)
    ttype = ScanImageTiffReader.open(filename) do io
        pxtype(io)
    end
    @test ttype == Int16

    # $ ScanImageTiffReader image bytes linfree_00001.tif
    # 5 MB
    dat = ScanImageTiffReader.open(filename) do io
        data(io)
    end
    @test abs(sizeof(dat)/2^20 - 5.0) < 1e-5
    # @test ScanImageTiffReader.open(filename, length) == size(dat, 3) # TODO: check this assumption

    # $ ScanImageTiffReader descriptions --frame 0 linfree_00001.tif
    # "frameNumbers = 1"
    # "acquisitionNumbers = 1"
    # "frameNumberAcquisition = 1"
    # "frameTimestamps_sec = 0.000000"
    # "acqTriggerTimestamps_sec = "
    # "nextFileMarkerTimestamps_sec = "
    # "endOfAcquisition =  0"
    # "endOfAcquisitionMode = 0"
    # "dcOverVoltage = 0"
    # "epoch = [2016 6 4 13 52 6.7667]"
    desc = ScanImageTiffReader.open(filename) do io
        split(description(io, 1), "\n")
    end
    @test desc[1] == "frameNumbers = 1"
    @test desc[2] == "acquisitionNumbers = 1"
    @test desc[3] == "frameNumberAcquisition = 1"
    @test desc[4] == "frameTimestamps_sec = 0.000000"
    @test desc[5] == "acqTriggerTimestamps_sec = "
    @test desc[6] == "nextFileMarkerTimestamps_sec = "
    @test desc[7] == "endOfAcquisition =  0"
    @test desc[8] == "endOfAcquisitionMode = 0"
    @test desc[9] == "dcOverVoltage = 0"
    @test desc[10] == "epoch = [2016 6 4 13 52 6.7667]"

    # $ ScanImageTiffReader metadata linfree_00001.tif
    # SI.TIFF_FORMAT_VERSION = 3
    # SI.VERSION_MAJOR = '2015'
    # SI.VERSION_MINOR = '4'
    # SI.acqState = 'grab'
    # SI.acqsPerLoop = 1
    # SI.extTrigEnable = false
    # SI.hBeams.beamCalibratedStatus = [false false false]
    # SI.hBeams.beamStatus = [0 0 0]
    # SI.hBeams.directMode = [false false false]
    # SI.hBeams.enablePowerBox = false
    # [...]
    # didn't parse for one reason or another...
    md = ScanImageTiffReader.open(filename) do io
        split(metadata(io), "\n")
    end
    @test md[1] == "SI.TIFF_FORMAT_VERSION = 3"
    @test md[2] == "SI.VERSION_MAJOR = '2015'"
    @test md[3] == "SI.VERSION_MINOR = '4'"
    @test md[4] == "SI.acqState = 'grab'"
    @test md[5] == "SI.acqsPerLoop = 1"
    @test md[6] == "SI.extTrigEnable = false"
    @test md[7] == "SI.hBeams.beamCalibratedStatus = [false false false]"
    @test md[8] == "SI.hBeams.beamStatus = [0 0 0]"
    @test md[9] == "SI.hBeams.directMode = [false false false]"
    @test md[10] == "SI.hBeams.enablePowerBox = false"
end


text:
using Documenter
using ConcurrentCollections

makedocs(
    sitename = "ConcurrentCollections",
    format = Documenter.HTML(),
    modules = [ConcurrentCollections],
)

for (root, dirs, files) in walkdir(joinpath(@__DIR__, "build"))
    for name in files
        path = joinpath(root, name)
        if endswith(name, ".html") || name == "search_index.js"
            html = replace(
                read(path, String),
                "ConcurrentCollections.Implementations" => "ConcurrentCollections",
            )
            write(path, html)
        end
    end
end

deploydocs(
    repo = "github.com/JuliaConcurrent/ConcurrentCollections.jl",
    push_preview = true,
)