text:
using FISH_MDS

filename = "data_fish/combined_Sen_merged_chr4_arm1.csv"
data = load_data(filename, scale=1, exponent=1/3)
#mds = MDSProblem(data, 1.0)#, lads, nads, 0.1)
mds = new_mds(data, 1.0, "data_fish/chr4_sen_probes.txt")
mds = remove_infs(mds)
println(filename)
#println(mds)
ipopt_problem = make_ipopt_problem(mds, radius_constraint=true,
    fish_constraint=true)
solveProblem(ipopt_problem)
coords = ipopt_problem.x
coords = reshape(coords, 3, int(length(coords)/3))'
interp_coords = interp_3d(coords, 10)
output_txt(coords, "test_out/chr4_sen_arm1_with_constraint.txt")
output_txt(interp_coords, "test_out/interp_chr4_sen_arm1_with_constraint.txt")

# tests: no constraint, with constraint, no interp,
# with interp, etc.

filename = "data_fish/combined_Sen_merged_chr4_arm1.csv"

println("no interp, with constraint")

run_mds(filename, interp=false,
    output_name = "test_out/sen_chr4_constrained_interp.txt",
    constraint="data_fish/chr4_sen_probes.txt")

println("no interp, with constraint, with autoscale")

run_mds(filename, interp=false, auto_scale=true,
    output_name = "test_out/sen_chr4_constrained_interp_autoscale.txt",
    constraint="data_fish/chr4_sen_probes.txt")

println("interp, with constraint")

run_mds(filename, interp=true,
    output_name = "test_out/sen_chr4_constrained_interp.txt",
    constraint="data_fish/chr4_sen_probes.txt")

println("no interp, no constraint")

run_mds(filename, interp=false,
    output_name = "test_out/sen_chr4_default.txt")

println("no interp, no constraint, fixed start")

run_mds(filename, interp=false,
    output_name = "test_out/sen_chr4_default_2.txt",
    starting_points_file="test_out/sen_chr4_default.txt")
# new file
filename = "data_fish/combined_Young_Qui_merged_chr18.csv"

run_mds(filename, output_name = "test_out/qui_chr18_default.txt")


text:
# This file is part of InvertedFiles.jl

include("invfiles.jl")
include("knr.jl")

ext: jl


text:
# Autogenerated wrapper script for ORC_jll for x86_64-w64-mingw32
export liborc, liborc_test, orc_bugreport, orcc

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Relative path to `liborc`
const liborc_splitpath = ["bin", "liborc-0.4-0.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
liborc_path = ""

# liborc-specific global declaration
# This will be filled out by __init__()
liborc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const liborc = "liborc-0.4-0.dll"

# Relative path to `liborc_test`
const liborc_test_splitpath = ["bin", "liborc-test-0.4-0.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
liborc_test_path = ""

# liborc_test-specific global declaration
# This will be filled out by __init__()
liborc_test_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const liborc_test = "liborc-test-0.4-0.dll"

# Relative path to `orc_bugreport`
const orc_bugreport_splitpath = ["bin", "orc-bugreport.exe"]

# This will be filled out by __init__() for all products, as it must be done at runtime
orc_bugreport_path = ""

# orc_bugreport-specific global declaration
function orc_bugreport(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ';', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ';', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(orc_bugreport_path)
    end
end

# Relative path to `orcc`
const orcc_splitpath = ["bin", "orcc.exe"]

# This will be filled out by __init__() for all products, as it must be done at runtime
orcc_path = ""

# orcc-specific global declaration
function orcc(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ';', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ';', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(orcc_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"ORC")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global liborc_path = normpath(joinpath(artifact_dir, liborc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global liborc_handle = dlopen(liborc_path)
    push!(LIBPATH_list, dirname(liborc_path))

    global liborc_test_path = normpath(joinpath(artifact_dir, liborc_test_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global liborc_test_handle = dlopen(liborc_test_path)
    push!(LIBPATH_list, dirname(liborc_test_path))

    global orc_bugreport_path = normpath(joinpath(artifact_dir, orc_bugreport_splitpath...))

    push!(PATH_list, dirname(orc_bugreport_path))
    global orcc_path = normpath(joinpath(artifact_dir, orcc_splitpath...))

    push!(PATH_list, dirname(orcc_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
include("../../src/IpyJL.jl")
using Main.IpyJL
using Test
    # Doesn't properly test now, ideally this will
    #    read the file after and test output, but I am lazy.
ipynbjl("ipynbtestbook.ipynb", "example.jl")


text:
export GumbelSoftmaxExplorer

using Random
using Flux: logsoftmax

struct GumbelSoftmaxExplorer <: AbstractExplorer
    rng::AbstractRNG
end

GumbelSoftmaxExplorer(; rng = Random.GLOBAL_RNG) = GumbelSoftmaxExplorer(rng)

function (p::GumbelSoftmaxExplorer)(v::AbstractVector{T}) where {T}
    logits = logsoftmax(v)
    u = rand(p.rng, T, length(logits))
    argmax(logits .- log.(-log.(u)))
end

function (p::GumbelSoftmaxExplorer)(
    v::AbstractVector{T},
    mask::AbstractVector{Bool},
) where {T}
    v[.!mask] .= typemin(T)
    p(v)
end


text:
# Run the test from HypoTest folder
using HypoTest
using Test
import Distributions
import Random

# List of test files:
tests = ["momentestim_test.jl",
		 "histdata_test.jl",
		 "disttest_test.jl",
		 "othertest_test.jl"]

# Run all tests in the list
for i in tests
	include(i)
end

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
module MaxHelpingHandCoreModeSpikes

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CoreModeSpikesUp" CoreModeSpikesUp(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, hotType::String="MaxHelpingHand/heatspike", coldType::String="cliffside")
@mapdef Entity "MaxHelpingHand/CoreModeSpikesDown" CoreModeSpikesDown(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, hotType::String="MaxHelpingHand/heatspike", coldType::String="cliffside")
@mapdef Entity "MaxHelpingHand/CoreModeSpikesLeft" CoreModeSpikesLeft(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, hotType::String="MaxHelpingHand/heatspike", coldType::String="cliffside")
@mapdef Entity "MaxHelpingHand/CoreModeSpikesRight" CoreModeSpikesRight(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, hotType::String="MaxHelpingHand/heatspike", coldType::String="cliffside")

const spikeTypes = String[
    "default",
    "outline",
    "cliffside",
    "reflection",
    "MaxHelpingHand/heatspike"
]

const entities = Dict{String, Type}(
    "up" => CoreModeSpikesUp,
    "down" => CoreModeSpikesDown,
    "left" => CoreModeSpikesLeft,
    "right" => CoreModeSpikesRight
)

const spikesUnion = Union{CoreModeSpikesUp, CoreModeSpikesDown, CoreModeSpikesLeft, CoreModeSpikesRight}

const placements = Ahorn.PlacementDict()
for (dir, entity) in entities
    key = "Core Mode Spikes ($(uppercasefirst(dir))) (max480's Helping Hand)"
    placements[key] = Ahorn.EntityPlacement(
        entity,
        "rectangle"
    )
end

Ahorn.editingOptions(entity::spikesUnion) = Dict{String, Any}(
    "hotType" => spikeTypes,
    "coldType" => spikeTypes
)

const directions = Dict{String, String}(
    "MaxHelpingHand/CoreModeSpikesUp" => "up",
    "MaxHelpingHand/CoreModeSpikesDown" => "down",
    "MaxHelpingHand/CoreModeSpikesLeft" => "left",
    "MaxHelpingHand/CoreModeSpikesRight" => "right"
)

const offsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (4, -4),
    "down" => (4, 4),
    "left" => (-4, 4),
    "right" => (4, 4)
)

const rotations = Dict{String, Number}(
    "up" => 0,
    "right" => pi / 2,
    "down" => pi,
    "left" => pi * 3 / 2
)

const resizeDirections = Dict{String, Tuple{Bool, Bool}}(
    "up" => (true, false),
    "down" => (true, false),
    "left" => (false, true),
    "right" => (false, true),
)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::spikesUnion)
    direction = get(directions, entity.name, "up")
    theta = rotations[direction] - pi / 2

    width = Int(get(entity.data, "width", 0))
    height = Int(get(entity.data, "height", 0))

    x, y = Ahorn.position(entity)
    cx, cy = x + floor(Int, width / 2) - 8 * (direction == "left"), y + floor(Int, height / 2) - 8 * (direction == "up")

    Ahorn.drawArrow(ctx, cx, cy, cx + cos(theta) * 24, cy + sin(theta) * 24, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::spikesUnion)
    if haskey(directions, entity.name)
        x, y = Ahorn.position(entity)

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        direction = get(directions, entity.name, "up")
        variant = get(entity.data, "hotType", "default")

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        ox, oy = offsets[direction]

        return Ahorn.Rectangle(x + ox - 4, y + oy - 4, width, height)
    end
end

Ahorn.minimumSize(entity::spikesUnion) = (8, 8)

function Ahorn.resizable(entity::spikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "hotType", "default")
        direction = get(directions, entity.name, "up")

        return resizeDirections[direction]
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::spikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "hotType", "default")
        direction = get(directions, entity.name, "up")

        width = get(entity.data, "width", 8)
        height = get(entity.data, "height", 8)

        for ox in 0:8:width - 8, oy in 0:8:height - 8
            drawX = ox + offsets[direction][1]
            drawY = oy + offsets[direction][2]

            Ahorn.drawSprite(ctx, "danger/spikes/$(variant)_$(direction)00", drawX, drawY)
        end
    end
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module Tabulars

import Base: @pure, @propagate_inbounds, @_pure_meta, @_inline_meta,
    @_propagate_inbounds_meta, get, getindex, setindex!, view, indices, size, length,
    summary, show, length, start, next, done, transpose, permutedims, ndims, show, print,
    endof

import Base: ==

export IndexError
export AbstractTabular, AbstractTable, AbstractSeries
export Tabular, Table, Series
export PermutedDimsTabular, PermutedDimsTable, PermutedDimsSeries
export SubTabular, SubTable, SubSeries
export tabulate

export Label, @l_str

include("util.jl")
include("Label.jl")

include("Tabular.jl")
include("indexing.jl")
include("PermutedDimsTabular.jl")
include("SubTabular.jl")
include("iteration.jl")
include("generic.jl")
include("show.jl")
include("tabulate.jl")

end # module

# endof (endof(series), endof(table, i))
# didn't throw for non-unique indices of tuples
# order of indices is weird
# fancy getindex matrix table

# TODO LIST
#
# Basic indexing interface:
# * slice/fancy setindex!
# * setindex! works when inner type is immutable
# * setindex for immutable types, and related trait?
#
# New containers
# * sort-based Associative or Tabular (instead of hash-based Dict everywhere)
# * TupleTabular (move current TupleTabular -> TypedTabular)
# * GroupedTabular - nested indexing (non-Cartesian)
# * SparseTabular - list of indices => values, with default value?
# * Support Generator
#
# Features / methods
# * rows(table), cols(table)
# * map (maybe map(series), map(rows(table)), map(cols(table))
# * reduce / reducedim / reducerows / reducecols
# * filter / filterrows / filtercols
# * groupby
# * indexable values (accelerated find/search)
#
# Performance issues
# * Long TupleTabulars (generated functions?)
# * Fast constructors (elide sameindices, maybe use @boundscheck?)
# * Avoid wrapping and reallocating data upon construction, where possible


text:
module FileUtils

using DeepDiffs
using Pkg

export replace_in_file, replace_in_files, map_filenames, replace_filenames, max_version

function _replace_in_files(filenames, replacement; showdiffs)
  for filename in filenames
    !isfile(filename) && continue
    replace_in_file(filename, replacement; showdiffs)
  end
  return nothing
end

function replace_in_file(filename, replacement; showdiffs=true)
  txt = read(filename, String)
  open(filename, "w") do f
    write(f, replace(txt, replacement))
  end
  if showdiffs
    txt_final = read(filename, String)
    diff = deepdiff(txt, txt_final)
    println(filename)
    removed_lines = removed(diff)
    added_lines = added(diff)
    txt_lines = split(txt, "\n")
    txt_final_lines = split(txt_final, "\n")
    @assert removed_lines == added_lines
    for n in 1:length(removed_lines)
      line_diff = deepdiff(txt_lines[removed_lines[n]], txt_final_lines[added_lines[n]])
      println("Line $(removed_lines[n]): ", line_diff)
    end
    println()
  end
  return nothing
end

function _replace_in_files_recursive(path, replacement; ignore_dirs, showdiffs)
  for (root, dirs, files) in walkdir(path)
    if !any(ignore_dir -> contains(root, ignore_dir), ignore_dirs)
      _replace_in_files(joinpath.(root, files), replacement; showdiffs=showdiffs)
    end
  end
  return nothing
end

"""
    replace_in_files(path::String, replacement::Pair{String,String}; ignore_dirs=[], recursive=false, showdiffs=true)

# Examples

```julia
replace_in_files(".", "@show" => "@test"; recursive=true, ignore_dirs=[".git"])
```
"""
function replace_in_files(path, replacement; ignore_dirs=[], recursive=false, showdiffs=true)
  if recursive
    _replace_in_files_recursive(path, replacement; ignore_dirs=ignore_dirs, showdiffs=showdiffs)
  else
    files = readdir(path)
    _replace_in_files(joinpath.(path, files), replacement; showdiffs=showdiffs)
  end
  return nothing
end

"""
    map_filenames(f::Function, path=pwd(); filter_filenames=Returns(true), force=false)

# Examples

```julia
map_filenames(filename -> replace(filename, ".txt" => ".txt.backup"); filter_filenames=startswith("example_"))
```
"""
function map_filenames(f, path=pwd(); filter_filenames=Returns(true), force=false)
  files = Base.filter(filter_filenames, readdir(path))
  for (i, file) in enumerate(files)
    old_filename = joinpath(path, file)
    new_filename = joinpath(path, f(file))
    println(old_filename, " => ", new_filename)
    if old_filename ≠ new_filename
      mv(old_filename, new_filename; force=force)
    end
  end
  return nothing
end

"""
    replace_filenames(path=pwd(), replacement::Pair{String,String}; filter=Returns(true), force=false)

# Examples

```julia
replace_filenames(pwd(), ".txt" => ".txt.backup")
```
"""
function replace_filenames(replacement, args...; kwargs...)
  return map_filenames(file -> replace(file, replacement), args...; kwargs...)
end

# TODO: move to a package called PkgUtils?
function max_version(depot_path::AbstractString, pkgname::AbstractString)
  path = joinpath(depot_path, "registries", "General",
                  first(pkgname, 1), pkgname, "Versions.toml")
  if isfile(path)
    return maximum(VersionNumber.(keys(Pkg.TOML.parse(read(path, String)))))
  end
  return v"0.0.0"
end

"""
    max_version(pkgname)

# Examples

```julia
max_version("ITensors")
```
"""
max_version(pkgname::AbstractString) = maximum(max_version.(DEPOT_PATH, pkgname))

end


text:
"""

## stability

Function spabac! performs Cholesky forward and back-substitution on
a symmetric skyline global matrix. The loads vector is updated.

###Arguments
```julia
stability(gsm, ggm, tol, limit)
```

### Arguments
```julia
* gsm::SparseMatrixCSC{Float64,Int}   : Factored global stiffness matrix
* ggm::SparseMatrixCSC{Float64,Int}   : Factored geometric matrix
* tol::Float64                          : Convergence tolerance
* limit::Int                          : Iteration limit
```
"""
function stability(gsm::SparseMatrixCSC{Float64,Int},
  ggm::SparseMatrixCSC{Float64,Int}, tol::Float64, limit::Int)
  #
  # This subroutine computes the smallest eigenvalue in a beam
  # stability analysis. This is the Julia version. Notice no !.
  #
  # Updateds iters, evec, ival
  #
  converged = true
  neq=size(gsm, 1)
  x0 = zeros(neq+1)
  x1 = zeros(neq+1)
  big = 0.0

  cfgsm = cholesky(gsm)

  iters = 0
  x0[2] = 1.0

  while true
   iters = iters+1
   x1[2:end] = ggm * x0[2:end]
   x1[2:end] = cfgsm \ x1[2:end]
   big = maximum(x1[2:end])
   if abs(minimum(x1[2:end])) > big
     big = minimum(x1[2:end])
   end
   x1 = x1 / big
   converged = (maximum(abs.(x1[2:end]-x0[2:end]))/maximum(abs.(x1[2:end])) < tol)
   x0 = copy(x1)
   if converged || iters==limit
     break
   end
  end
  x1[2:end]=x1[2:end] / sqrt(sum(x1[2:end] .^ 2))
  ival = 1.0 / big
  (iters, x1, ival)
end


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

@doc raw"""PersistentVolumeClaimSpec describes the common attributes of storage devices and allows a Source for provider-specific attributes

    IoK8sApiCoreV1PersistentVolumeClaimSpec(;
        accessModes=nothing,
        dataSource=nothing,
        resources=nothing,
        selector=nothing,
        storageClassName=nothing,
        volumeMode=nothing,
        volumeName=nothing,
    )

    - accessModes::Vector{String} : AccessModes contains the desired access modes the volume should have. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#access-modes-1
    - dataSource::IoK8sApiCoreV1TypedLocalObjectReference : This field requires the VolumeSnapshotDataSource alpha feature gate to be enabled and currently VolumeSnapshot is the only supported data source. If the provisioner can support VolumeSnapshot data source, it will create a new volume and data will be restored to the volume at the same time. If the provisioner does not support VolumeSnapshot data source, volume will not be created and the failure will be reported as an event. In the future, we plan to support more data source types and the behavior of the provisioner may change.
    - resources::IoK8sApiCoreV1ResourceRequirements : Resources represents the minimum resources the volume should have. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#resources
    - selector::IoK8sApimachineryPkgApisMetaV1LabelSelector : A label query over volumes to consider for binding.
    - storageClassName::String : Name of the StorageClass required by the claim. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#class-1
    - volumeMode::String : volumeMode defines what type of volume is required by the claim. Value of Filesystem is implied when not included in claim spec. This is a beta feature.
    - volumeName::String : VolumeName is the binding reference to the PersistentVolume backing this claim.
"""
mutable struct IoK8sApiCoreV1PersistentVolumeClaimSpec <: SwaggerModel
    accessModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: accessModes
    dataSource::Any # spec type: Union{ Nothing, IoK8sApiCoreV1TypedLocalObjectReference } # spec name: dataSource
    resources::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ResourceRequirements } # spec name: resources
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    storageClassName::Any # spec type: Union{ Nothing, String } # spec name: storageClassName
    volumeMode::Any # spec type: Union{ Nothing, String } # spec name: volumeMode
    volumeName::Any # spec type: Union{ Nothing, String } # spec name: volumeName

    function IoK8sApiCoreV1PersistentVolumeClaimSpec(;accessModes=nothing, dataSource=nothing, resources=nothing, selector=nothing, storageClassName=nothing, volumeMode=nothing, volumeName=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("accessModes"), accessModes)
        setfield!(o, Symbol("accessModes"), accessModes)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("dataSource"), dataSource)
        setfield!(o, Symbol("dataSource"), dataSource)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("storageClassName"), storageClassName)
        setfield!(o, Symbol("storageClassName"), storageClassName)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("volumeMode"), volumeMode)
        setfield!(o, Symbol("volumeMode"), volumeMode)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimSpec, Symbol("volumeName"), volumeName)
        setfield!(o, Symbol("volumeName"), volumeName)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaimSpec

const _property_map_IoK8sApiCoreV1PersistentVolumeClaimSpec = Dict{Symbol,Symbol}(Symbol("accessModes")=>Symbol("accessModes"), Symbol("dataSource")=>Symbol("dataSource"), Symbol("resources")=>Symbol("resources"), Symbol("selector")=>Symbol("selector"), Symbol("storageClassName")=>Symbol("storageClassName"), Symbol("volumeMode")=>Symbol("volumeMode"), Symbol("volumeName")=>Symbol("volumeName"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaimSpec = Dict{Symbol,String}(Symbol("accessModes")=>"Vector{String}", Symbol("dataSource")=>"IoK8sApiCoreV1TypedLocalObjectReference", Symbol("resources")=>"IoK8sApiCoreV1ResourceRequirements", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("storageClassName")=>"String", Symbol("volumeMode")=>"String", Symbol("volumeName")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaimSpec }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaimSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaimSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaimSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaimSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaimSpec[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaimSpec)
    true
end

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


text:
using RegressionAndOtherStories
using Test

include(joinpath(@__DIR__, "test_hibbs.jl"))
include(joinpath(@__DIR__, "test_model_summary.jl"))


text:
"""

# cmdline

Recursively parse the model to construct command line.

### Method
```julia
cmdline(m)
```

### Required arguments
```julia
* `m::SampleModel`                : CmdStanSampleModel
```
"""
function cmdline(m::Union{DiagnoseModel, Diagnose, Gradient}, id)

  #=
  `./bernoulli diagnose test=gradient epsilon=1.0e-6 error=1.0e-6
  random seed=-1 id=1 data file=bernoulli_1.data.R
  output file=bernoulli_diagnose_1.csv refresh=100`
  =#

  cmd = ``
  # parse top level
  if isa(m, DiagnoseModel)
    # Handle the model name field for unix and windows
    cmd = `$(m.exec_path)`

    # Sample() specific portion of the model, might be recursive
    cmd = `$cmd $(cmdline(getfield(m, :method), id))`

    # Common to all models, not recursive
    cmd = `$cmd random seed=$(getfield(m, :seed).seed)`

    # Init file required?
    if length(m.init_file) > 0 && isfile(m.init_file[id])
      cmd = `$cmd init=$(m.init_file[id])`
    else
      cmd = `$cmd init=$(m.init.bound)`
    end

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file[id])
      cmd = `$cmd id=$(id) data file=$(m.data_file[id])`
    end

    # Output files
    cmd = `$cmd output`
    if length(m.sample_file) > 0
      cmd = `$cmd file=$(m.sample_file[id])`
    end
    if length(m.diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(m.diagnostic_file[id])`
    end

    cmd = `$cmd refresh=$(string(getfield(m, :output).refresh))`

  else # The 'recursive' part, currently only Gradient <: Diagnostics

    if isa(m, Diagnostics)
      # Inset 'test=gradient' into cmdline
      cmd = `$cmd test=$(split(lowercase(string(typeof(m))), '.')[end])`
    else
      # Insert initial `diagnose` into cmdline
      cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    end
    for name in fieldnames(typeof(m))
      if length(fieldnames(typeof(getfield(m, name)))) == 0
        cmd = `$cmd $(name)=$(getfield(m, name))`
      else
        # Composite (Gradient) object, handle all fields
        # (by recursively calling cmdline)
        cmd = `$cmd $(cmdline(getfield(m, name), id))`
      end
    end
  end

  cmd

end


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

using LazyArtifacts
using Test

mktempdir() do tempdir
    LazyArtifacts.Artifacts.with_artifacts_directory(tempdir) do
        socrates_dir = artifact"socrates"
        @test isdir(socrates_dir)
        ex = @test_throws ErrorException artifact"c_simple"
        @test startswith(ex.value.msg, "Artifact \"c_simple\" was not installed correctly. ")
    end
end

# Need to set depwarn flag before testing deprecations
@test success(run(setenv(`$(Base.julia_cmd()) --depwarn=no --startup-file=no -e '
    using Artifacts, Pkg
    using Test
    mktempdir() do tempdir
        Artifacts.with_artifacts_directory(tempdir) do
            socrates_dir = @test_logs(
                    (:warn, "using Pkg instead of using LazyArtifacts is deprecated"),
                    artifact"socrates")
            @test isdir(socrates_dir)
        end
    end'`,
    dir=@__DIR__)))

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
# Adapted from https://ccrma.stanford.edu/~jos/filters/Parallel_SOS_Transfer_Function.html

function psos2tf(sos,g=1)
    nsecs = size(sos,1)
    if nsecs<1, B=[]; A=[]; return; end

    Bs = sos[:,1:3];
    As = sos[:,4:6];
    B = Bs[1,:];
    A = As[1,:];

    for i=2:nsecs
        B = conv(B,As[i,:]) + conv(A,Bs[i,:]);
        A = conv(A,As[i,:]);
    end
    return [B,A]
end # function


text:
module IsConstFoldable

using InteractiveUtils
using Rematch

export @is_const_foldable

"""
    @is_const_foldable foo(1, 2)

A macro which will return true if we are *sure* that an expression const-folds, and false
otherwise. Note that there may be false negatives (because LLVM does its own optimization
passes), but there should not be false positives.

Note that it tests this by checking the @code_typed of a helper function whose body is your
expression, allowing the optimizer the chance to fold away the logic with awareness of the
context of being called with Const inputs.

You can use this within a Unit Test, e.g.:

```julia
@test @is_const_foldable zero(Int)
@test @is_const_foldable zero(Rational)  # Fails
```
"""
macro is_const_foldable(expr)
    return quote
        # We embed the expr in a wrapper() function, which ensures that during compliation
        # and optimization, the optimizer has access to all static constants provided.
        function wrapper() $(esc(expr)) end
        function _test_wrapper()
            # If the code in wrapper() is reduced to only a single return statement, then
            # the logic in `expr` has been const-folded away.
            # This single return statment would be something like this: `:(return 0)`
            codeinfo = ($InteractiveUtils.@code_typed wrapper())[1]
            if length(codeinfo.code) == 1
                # Ensure that the statement is a return line of some kind
                # (This also works for empty functions)
                codeinfo.code[end].head == :return && return true
            end
            # Otherwise, the optimizer was not able to eliminate the entirety of the operation.
            return false
        end
        _test_wrapper()
    end
end

end


text:
module MultivariateCalculusTest

using VarianceComponentModels, Base.Test

srand(123)

# test vech
@testset "vech" begin
  a = reshape(1:9, 3, 3)
  vecha = vech(a)
  @test vecha[1] == a[1, 1]
  @test vecha[2] == a[2, 1]
  @test vecha[3] == a[3, 1]
  @test vecha[4] == a[2, 2]
  @test vecha[5] == a[3, 2]
  @test vecha[6] == a[3, 3]
  a = 1.0
  @test vech(a) == a
end

# test trilind
@testset "trilind" begin
  n = 3
  A = randn(n, n)
  @test vecnorm(vech(A) - A[trilind(A)]) ≈ 0.0
end

# test triuind
@testset "triuind" begin
  n = 2
  idx = triuind(n, n)
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(n)
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(randn(n, n))
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(randn(n, n), 1)
  @test length(idx) == binomial(n, 2)
  @test idx[1] == 3
end

# test commutation
@testset "commutation" begin
  m, n = 3, 2
  A = randn(m, n)
  @test vecnorm(commutation(m, n) * vec(A) - vec(A')) ≈ 0.0
  @test vecnorm(commutation(A) * vec(A) - vec(A')) ≈ 0.0
  @test vecnorm(spcommutation(m, n) * vec(A) - vec(A')) ≈ 0.0
  @test vecnorm(spcommutation(A) * vec(A) - vec(A')) ≈ 0.0
  B = randn(m, m)
  @test vecnorm(commutation(m) * vec(B) - vec(B')) ≈ 0.0
  @test vecnorm(commutation(eltype(B), m) * vec(B) - vec(B')) ≈ 0.0
  @test vecnorm(spcommutation(eltype(B), m) * vec(B) - vec(B')) ≈ 0.0
end

# test duplication
@testset "duplication" begin
  n = 3
  A = randn(n, n)
  A = 0.5(A + A') # symmetrize
  @test vecnorm(duplication(n) * vech(A) - vec(A)) ≈ 0.0
  @test vecnorm(duplication(A) * vech(A) - vec(A)) ≈ 0.0
  @test vecnorm(spduplication(n) * vech(A) - vec(A)) ≈ 0.0
  @test vecnorm(spduplication(A) * vech(A) - vec(A)) ≈ 0.0
end

# test chol_gradient
@testset "chol_gradient" begin
  n = 3
  A = randn(n, n)
  A = 0.5(A + A') # symmetrize
  L = tril(randn(n, n))
  # calculate gradient wrt L using chol_gradient
  dL1 = chol_gradient(vec(A), L)
  # alternative way to calculate gradient wrt L
  dL2 = 2.0vech((A * L) + L' * A' - diagm(diag(A * L)))
  @test vecnorm(dL1 - dL2) ≈ 0.0
end

# test kron_gradient
@testset "kron_gradient" begin
  n, q, p, r = 2, 3, 4, 5
  X = randn(n, q)
  Y = randn(p, r)
  M = kron(X, Y)
  # calculate gradient wrt X using kron_gradient
  dX1 = kron_gradient(vec(M), Y, n, q)
  # alternative way to calculate gradient wrt X
  dX2 = zeros(X)
  for j = 1:q
    for i = 1:n
      dX2[i, j] = vecdot(Y, M[(i-1)*p+1:i*p, (j-1)*r+1:j*r])
    end
  end
  @test vecnorm(dX1 - vec(dX2)) < 1e-8
end

# test duplication
@testset "kronaxpy" begin
  m, n, p, q = 3, 4, 5, 6
  A = randn(m, n)
  X = randn(p, q)
  Y = zeros(m * p, n * q)
  kronaxpy!(A, X, Y)
  @test vecnorm(Y - kron(A, X)) ≈ 0.0
end

# test bump_diagonal
@testset "bump_diag" begin
  A = randn(3, 3)
  B = copy(A)
  bump_diagonal!(B, 1e-3)
  @test B[1, 1] - A[1, 1] ≈ 1e-3
  @test B[2, 2] - A[2, 2] ≈ 1e-3
  @test B[3, 3] - A[3, 3] ≈ 1e-3
end

# test clamp_diagonal
@testset "clamp_diag" begin
  n = 3
  A = randn(n, n)
  B = copy(A)
  clamp_diagonal!(B, 0.0, 1.0)
  @test all(0.0 .≤ diag(B) .≤ 1.0)
end

end # module MultivariateCalculusTest


text:
using Cairo

function common_header()
    ## header to provide surface and context
    c = CairoRGBSurface(256,256);
    cr = CairoContext(c);

    save(cr);
    set_source_rgb(cr,0.8,0.8,0.8);    # light gray
    rectangle(cr,0.0,0.0,256.0,256.0); # background
    fill(cr);
    restore(cr);

    return c,cr
end

function write_out_picture(filename::AbstractString,c::CairoSurface)
    ## mark picture with current date
    move_to(cr,0.0,12.0);
    set_source_rgb(cr, 0,0,0);
    show_text(cr,Libc.strftime(time()));
    ##
    write_to_png(c,filename);
    nothing
end

function example_current_point(cr)
    save(cr);

    x=25.6;  y=128.0;
    x1=102.4; y1=30.4;
    x2=123.6; y2=45.6;
    x3=130.4; y3=128.0;

    move_to(cr, x, y);
    curve_to(cr, x1, y1, x2, y2, x3, y3);

    if has_current_point(cr)
        x,y = get_current_point(cr);
        save(cr)
        move_to(cr,x,y)
        set_source_rgb(cr, 0,0,1.0);
        show_text(cr,"current point")
        restore(cr)
    end

set_line_width(cr, 10.0);
stroke(cr);

    restore(cr);
    save(cr);

    nothing
end

c,cr = common_header();
save(cr);
example_current_point(cr);
restore(cr);
write_out_picture("sample_current_point.png",c);


text:
module SimilarityMetrics
	using Distance

	export cosine, expl1, expl2, jaccard, matching, overlap, pearson

	typealias AbstractSet Union(Set, IntSet)

	include("cosine.jl")
	include("expl1.jl")
	include("expl2.jl")
	include("jaccard.jl")
	include("matching.jl")
	include("overlap.jl")
	include("pearson.jl")
end


text:
#https://www.mathworks.com/help/matlab/ref/matlab.graphics.chart.primitive.histogram.html
function histogram(X;options=Dict())
    MATLAB.@mput X
    optionString = namevaluepairs(options)
    MATLAB.eval_string("histogram(X "*optionString*")")
end


text:
"""
makes some bitwise operations easier
"""

module BitWiseUtils

export llvmPassOnes,@bitDilatate,bitPassOnesUINt,bitDilatateUINt,@setBitTo1, @setBitTo , @setBitTo1UINt32, isBit1AtPos, bitDilatate,@bitPassOnes,bitPassOnesUINt
"""
sets given bit of supplied number to 1
"""
macro setBitTo1(numb,pos)
    return esc(quote
        $numb|= 1 << ($pos-1)
    end)

end#setBitTo1

macro setBitTo1UINt32(numb,pos)
    return esc(quote
        $numb|= UInt32(1) << ($pos-1)
    end)

end#setBitTo1

"""
set bit of number numb in position pos to value val
"""
macro setBitTo(numb,pos,val)
    return esc(quote
    # Suppose you want to change bit N of x, where N=0 means the least-significant bit. You don’t care what the old value was, but the new value is B (either 0 or 1). This will do the trick:
     ($numb) =(($numb) & ~(1<<($pos-1))) | ($(val)<<($pos-1))
    end)

end#setBitTo1

"""
returning boolean indicating whether we have 1 at given spot in the integer
"""
function isBit1AtPos(numb,pos)::Bool
  return  (numb & (1 << (pos - 1)))>0
end#isBit1AtPos

"""
needed for bit dilatation in the kernel
    Given 32 bit integer source and target i woul like to do sth like one dimensional mathematical morphology dilatation so given bit at position x is 1 in source i would like to set bits at position x+1 and x-1 at target to 1 Other bits of target should not be changed so for example

Given 1 0 0 0 1 0 0 0 ... Set target integer into 1 1 0 1 1 1 0 0 ...

https://stackoverflow.com/questions/70134566/bit-wise-dilatation
"""
function bitDilatate(x::UInt32,locArr::UInt32)::UInt32
        locArr = x
        if(isBit1AtPos(locArr,1))
            @setBitTo1UINt32(x,2)
        end
        for i in 2:31
            if(isBit1AtPos(locArr,i))
                @setBitTo1UINt32(x,i+1)
                @setBitTo1UINt32(x,i-1)
            end
        end
        if(isBit1AtPos(locArr,32))
            @setBitTo1UINt32(x,31)
        end
        return x
end
macro bitDilatate(x)
    return esc(quote
    reinterpret(UInt32,( Int32( (($x )>> 1) | ($x) | (($x) << 1))))
    end)#quote
end

function bitDilatateUINt(x::UInt32)::UInt32
    return ((x )>> 1) | (x) | ((x) << 1)
end

"""
Given 32 bit integers x  and y  i would like to set bits of x to 1 if in corresponding position of y there is 1  without modyfing other bits of y
    So for example  if x is 1 0 0 0 0 1 ...And y is 0 0 0 1 0 0  ...
    I would like a result that would be  1 0 0 1 0 1
"""
macro bitPassOnes(source,target)
    return esc(quote
    reinterpret(UInt32,( Int32( (($target)|($source)) )))
    end)#quote
end

function bitPassOnesUINt(source::UInt32,target::UInt32) ::UInt32
    for i in 15:20
        if(isBit1AtPos(source,i))
            @setBitTo1UINt32(target,i)
        end
    end
    return reinterpret(UInt32,( UInt32(target )))
    # return ((target)|(source))
end

"""
increments given UINT16 given both boolGold and boolSegm are true
"""
function llvmPassOnes(source::UInt32,target::UInt32)::UInt32
    Base.llvmcall("""
    %3 = or i32 %0, %1
    ret i32 %3""", UInt32, Tuple{UInt32,UInt32}, source, target)
end

end#BitWiseUtils

# for i in 1:16
#     div,remm = divrem(i-1,4)
#     # aFrag[fld(div+1,2)+1,rem+1 ] = dataShmem[rem+1,fld(div+1,2)+1 ]

#     a =  ((i-1) & (3))+1
#     b = ~((i+3)>>2 & 1)+3
#     c = ((i-1)>>2 )+1
#     aFrag[a,c] = dataShmem[b,a]*( ((i>4 && i<13)*-2)+1 )
# @info "i$(i)   a $(a)   b $(b)   c $(c)"

# end

# numb = UInt32(0)
# #settingcorrectly
# numb |= 1 << 1
# numb |= 1 << 2
# numb |= 1 << 0
# numb |= 1 << 0
# numb |= 1 << 5
# numb

#reading...
# numb>>1 & UInt32(1)
# numb>>2 & UInt32(1)
# numb>>3 & UInt32(1)
# numb>>4 & UInt32(1)
# numb>>5 & UInt32(1)


text:
# TODO overload the correct AbstractPlotting functions with the below js fragments

function mouse_position(scene::Scene, window::ThreeScreen)
    on_mousedown = @js function (e, context)
        $painting[] = true
        @var el = context.dom.querySelector("#surface")
        @var context = el.getContext("2d")
        @var rect = el.getBoundingClientRect();
        @var x = e.clientX - rect.left;
        @var y = e.clientY - rect.top;
        window.addclick(x, y, false);
        window.redraw(context, $paintbrush_ob[], rect, false);
    end
    on_mouseup = @js function (e, context)
        $painting[] = false
end
    disconnect!(event); disconnect!(window, mouse_position)
    event[] = correct_mouse(window, GLFW.GetCursorPos(window)...)
    GLFW.SetCursorPosCallback(window, cursorposition)
end

function disconnect!(window::GLFW.Window, ::typeof(mouse_position))
    GLFW.SetCursorPosCallback(window, nothing)
end
document.addEventListener( 'mousemove', onDocumentMouseMove, false );
window.addEventListener( 'resize', onWindowResize, false );
function onWindowResize() {
    camera.aspect = window.innerWidth / window.innerHeight;
    camera.updateProjectionMatrix();
    renderer.setSize( window.innerWidth, window.innerHeight );
}
function onDocumentMouseMove( event ) {
    event.preventDefault();
    mouse.x = ( event.clientX / window.innerWidth ) * 2 - 1;
    mouse.y = - ( event.clientY / window.innerHeight ) * 2 + 1;
}

camera.updateMatrixWorld();
// find intersections
raycaster.setFromCamera( mouse, camera );
var intersects = raycaster.intersectObjects( parentTransform.children, true);
if ( intersects.length > 0 ) {
	if ( currentIntersected !== undefined ) {
		currentIntersected.material.linewidth = 1;
	}
	currentIntersected = intersects[ 0 ].object;
	currentIntersected.material.linewidth = 5;
	sphereInter.visible = true;
	sphereInter.position.copy( intersects[ 0 ].point );
} else {
	if ( currentIntersected !== undefined ) {
		currentIntersected.material.linewidth = 1;
	}
	currentIntersected = undefined;
	sphereInter.visible = false;
}
renderer.render( scene, camera );
}


text:
module adapters_sqlite_structured_test

using Test # @test
using Octo.Adapters.SQLite # Schema from SELECT FROM WHERE

struct User
end

u = from(User)
@test_throws Schema.TableNameError to_sql([FROM u])

Schema.model(User, table_name="users")

@test to_sql([FROM u]) == "FROM users"
@test to_sql([SELECT * FROM u]) == "SELECT * FROM users"
@test to_sql([SELECT u.id FROM u]) == "SELECT id FROM users"
@test to_sql([SELECT (u.id, u.name) FROM u]) == "SELECT id, name FROM users"
@test to_sql([WHERE u.id == 2]) == "WHERE id = 2"

u = from(User, :u)
@test to_sql([FROM u]) == "FROM users u" #
@test to_sql([SELECT (u.id, u.name) FROM u]) == "SELECT u.id, u.name FROM users u" #
@test to_sql([WHERE u.id == 2]) == "WHERE u.id = 2"
@test to_sql([SELECT (u.id, u.name) FROM u WHERE u.id == 2]) == "SELECT u.id, u.name FROM users u WHERE u.id = 2" #

# using Octo.Adapters.SQLite # INSERT INTO VALUES UPDATE SET
using Octo.Adapters.SQLite: Enclosed, PlaceHolder

struct Temp
end
Schema.model(Temp, table_name="temp")

using Octo.Adapters.SQLite: placeholders

temp = from(Temp)
@test to_sql([INSERT INTO temp VALUES placeholders(3)]) == "INSERT INTO temp VALUES (?, ?, ?)"
@test to_sql([UPDATE temp SET (title="Texas",) WHERE temp.AlbumId == 6]) == "UPDATE temp SET title = 'Texas' WHERE AlbumId = 6"

end # module adapters_sqlite_structured_test


text:
struct A end
struct B end
struct C end
struct D end
struct E end
struct F end

abstract type T1 end
abstract type T2 <: T1 end
abstract type T3 <: T2 end
abstract type T4 <: T3 end
abstract type T5 <: T4 end
abstract type T6 <: T5 end
abstract type T7 <: T6 end
abstract type T8 <: T7 end

f(::Any) = Union{}
f(::Type{Any}) = Union{}
f(::Type{T}) where T = f(supertype(T))
f(::Type{T1}) = Union{A, f(supertype(T1))}
f(::Type{T2}) = Union{B, f(supertype(T2))}

trait(::Type{Float64}, ::Type{T1}) = A
trait(::Type{Float64}, ::Type{T2}) = B
trait(::Type{T}, ::Type{U}) where {T, U} =
  trait(T, supertype(U))

struct Trait{T} end

# Can we get the lowest type, at the type level?

traits(::Type{T}) = Union{trait(T, Union{}), traits()}


text:
# The integrand. p encapsulates the passed parameters.
"The idea here is to have a function which steps in dtau, not ds."
function f_tilde(du, u, args, t)

    # tilde_alpha is a function
    zstart, v0, pars, interp, gpre, tilde_alpha_func, tilde_tau_inv = args

    # t is now tilde_tau, not s
    zprime = tilde_tau_inv(t)
    println("tilde tau ", t, " tau ", u[1], " zprime ", zprime/AU)

    # Write in coordinate conversions for querying alpha and sfunc
    rcyl, z, vlos = get_coords(zprime, gpre)

    # Calculate Delta v at this position
    Deltav = v0 * 1.e5 - vlos # get_vlos(xprime, yprime, zprime, pars)

    # Look up RT quantities from nearest-neighbor interp.
    DV2, sfunc, Upsilon = interp(rcyl, z)

    # Calculate alpha
    alpha = Upsilon * exp(-Deltav^2/DV2)
    tilde_alpha = tilde_alpha_func(zprime)
    println("alpha / tilde_alpha, ", alpha/tilde_alpha)

    # Update dtau/ds and dI/ds
    du[1] = alpha/tilde_alpha
    du[2] = exp(-u[1]) * sfunc * alpha/tilde_alpha
end

"Use estimates of the disk structure to take better steps in tau."
function trace_pixel_tilde(xprime, yprime, v0, mol::Molecule, pars::AbstractParameters, interp)
    rmax = 700.0 * AU

    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    z0s = geometry.get_zps(xprime, yprime, pars, v0)

    # create a tilde_alpha_func
    # Get the minimum rcyl
    rcyl_min = xprime

    # Occurs at what value of zprime?
    zprime_rcyl_min = cosd(pars.incl)/sind(pars.incl) * yprime

    # Actually, we want the zprime that corresponds to the midplane crossing
    zprime_midplane = -tand(pars.incl) * yprime
    # Calculate the rcyl_mid here

    # Get amplitude at midplane crossing
    # Corresponds to r_cyl_mid, 0.0
    r_cyl_mid = sqrt(xprime^2 + (cosd(pars.incl) * yprime + sind(pars.incl) * zprime_midplane)^2)

    # Get scale height at rcyl_min (will be an underestimate)
    H_rcyl_min = model.Hp(rcyl_min, pars)
    H_rcyl_mid = model.Hp(r_cyl_mid, pars)

    sigma_Upsilon = H_rcyl_mid / cosd(pars.incl)
    a_Upsilon = Upsilon_nu(r_cyl_mid, 0.0, mol.nu_0, pars, mol)
    mu_Upsilon = zprime_midplane

    DeltaV = sqrt(2 * kB * model.temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)

    mu_upsilon1, mu_upsilon2 = z0s
    sigma_upsilon = DeltaV * sqrt(2)/3 * (xprime * sqrt(G * pars.M_star * M_sun) * sind(pars.incl))^(4./3) / (mu_upsilon1 * (v0 * 1.e5)^(7/3))
    # a_upsilon = sqrt(2pi) * sigma_upsilon

    sigma_alpha2 = 1/(sigma_Upsilon^(-2) + sigma_upsilon^(-2)) # actually sigma^2, hence the 2 in the name
    mu_alpha1 = sigma_alpha2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon1/sigma_upsilon^2)
    mu_alpha2 = sigma_alpha2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon2/sigma_upsilon^2)

    a_alpha1 = a_Upsilon * exp(- (mu_Upsilon - mu_upsilon1)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))
    a_alpha2 = a_Upsilon * exp(- (mu_Upsilon - mu_upsilon2)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))

    function tilde_alpha_func(zprime)
        if zprime > 0
            return a_alpha1 * exp(-0.5 * (zprime - mu_alpha1)^2/sigma_alpha2)
        else
            return a_alpha2 * exp(-0.5 * (zprime - mu_alpha2)^2/sigma_alpha2)
        end

    end

    "Given tilde_tau, invert to find zprime"
    function tilde_tau_inv(tilde_tau)
        zprime = mu_alpha1 + sqrt(2 * sigma_alpha2) * erfcinv(tilde_tau * sqrt(2)/(a_alpha1 * sqrt(sigma_alpha2 * pi)))
    end

    # Calculate the bounding positions for emission
    # If it's two, just trace it.
    # If it's four, break it up into two integrals.

    bounds = get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(bounds) == 4
        z1start, z1end, z2start, z2end = bounds
    elseif length(bounds) == 2
        z1start, z1end = bounds
    end

    # Start the ode to integrate the first path.
    u0 = Float64[0.0; 0.0]

    gpre = precalc_geo(xprime, yprime, pars)

    args = (z1start, v0, pars, interp, gpre, tilde_alpha_func, tilde_tau_inv)

    println("z1start ", z1start/AU)
    tilde_tau_start = tilde_alpha_func(z1start)
    println("Starting at tilde_tau ", tilde_tau_start)

    # stop when we get to dtau of a large value (will trigger completion before)
    tspan = (tilde_tau_start, 100.)

    prob = ODEProblem(f_tilde, u0, tspan, args)

    # sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dtmax=5*AU, dt=dt)
    sol = solve(prob, alg, dtmax=0.5, dt=0.2, callback=cb)
    # println(sol)

    tau_final, int_final = sol.u[end]
    return tau_final, int_final

    # if (tau_final < tau_thresh) & (length(bounds) == 4)
    #     #We didn't reach the threshold in the first bounding region, so continue in the second.
    #     tspan = (0, z2start - z2end)
    #
    #     # Start the ode to integrate the first path from the end of the second.
    #     u0 = Float64[tau_final; int_final]
    #     args = (z2start, v0, pars, interp, gpre)
    #     prob = ODEProblem(f, u0, tspan, args)
    #     # sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dense=false, save_everystep=false, dtmax=3*AU, dt=dt)
    #     sol = solve(prob, alg, reltol=rt, abstol=at, dtmax=3*AU, dt=dt)
    #     tau_final, int_final = sol.u[end]
    #
    #     zps = z2start - sol.t
    #
    #     return zps, sol
    # else
    #     return zps, sol
    # end

end

"v0 is the central velocity of the channel, relative to the disk systemic velocity."
function trace_pixel_debug(xprime, yprime, v0, mol::Molecule, pars::AbstractParameters, interp)

    rmax = 700.0 * AU

    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # Calculate the bounding positions for emission
    # If it's two, just trace it.
    # If it's four, break it up into two integrals.

    bounds = get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(bounds) == 4
        z1start, z1end, z2start, z2end = bounds
    elseif length(bounds) == 2
        z1start, z1end = bounds
    end

    # Start the ode to integrate the first path.
    u0 = Float64[0.0; 0.0]

    gpre = precalc_geo(xprime, yprime, pars)

    args = (z1start, v0, pars, interp, gpre)

    tspan = (0.0, z1start - z1end)

    prob = ODEProblem(f, u0, tspan, args)

    # sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dtmax=5*AU, dt=dt)
    sol = solve(prob, alg, reltol=rt, abstol=at, dt=dt)
    # println(sol)

    zps = z1start - sol.t

    tau_final, int_final = sol.u[end]

    if (tau_final < tau_thresh) & (length(bounds) == 4)
        #We didn't reach the threshold in the first bounding region, so continue in the second.
        tspan = (0, z2start - z2end)

        # Start the ode to integrate the first path from the end of the second.
        u0 = Float64[tau_final; int_final]
        args = (z2start, v0, pars, interp, gpre)
        prob = ODEProblem(f, u0, tspan, args)
        # sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dense=false, save_everystep=false, dtmax=3*AU, dt=dt)
        sol = solve(prob, alg, reltol=rt, abstol=at, dtmax=3*AU, dt=dt)
        tau_final, int_final = sol.u[end]

        zps = z2start - sol.t

        return zps, sol
    else
        return zps, sol
    end

end


text:
# ---
# title: 786. K-th Smallest Prime Fraction
# id: problem786
# author: Indigo
# date: 2021-01-28
# difficulty: Hard
# categories: Binary Search, Heap
# link: <https://leetcode.com/problems/k-th-smallest-prime-fraction/description/>
# hidden: true
# ---
#
# A sorted list `A` contains 1, plus some number of primes.  Then, for every p <
# q in the list, we consider the fraction p/q.
#
# What is the `K`-th smallest fraction considered?  Return your answer as an
# array of ints, where `answer[0] = p` and `answer[1] = q`.
#
#
#
#     **Examples:**
#     Input: A = [1, 2, 3, 5], K = 3
#     Output: [2, 5]
#     Explanation:
#     The fractions to be considered in sorted order are:
#     1/5, 1/3, 2/5, 1/2, 3/5, 2/3.
#     The third fraction is 2/5.
#
#     Input: A = [1, 7], K = 1
#     Output: [1, 7]
#
#
# **Note:**
#
#   * `A` will have length between `2` and `2000`.
#   * Each `A[i]` will be between `1` and `30000`.
#   * `K` will be between `1` and `A.length * (A.length - 1) / 2`.
#
#
## @lc code=start
using LeetCode

function kth_smallest_prime_fraction(arr::Vector{Int}, k::Int)
    hp = BinaryMinHeap(1 .// arr[2:end])
    next_prime = Dict((arr[i], arr[i+1]) for i in 1:length(arr)-1)
    while (k -= 1) > 0
        tp = pop!(hp)
        if next_prime[tp.num] < tp.den
            push!(hp, next_prime[tp.num] // tp.den)
        end
    end
    return first(hp)
end
## @lc code=end


text:
# Autogenerated wrapper script for HiGHS_jll for x86_64-w64-mingw32-cxx11
export highs, libhighs

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("HiGHS")
JLLWrappers.@declare_executable_product(highs)
JLLWrappers.@declare_library_product(libhighs, "libhighs.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        highs,
        "bin\\highs.exe",
    )

    JLLWrappers.@init_library_product(
        libhighs,
        "bin\\libhighs.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Vector3D # this is simpler and faster than the ImmutableArrays

export Vec, +, *, -, /, dot, norm, unitize, myconvert

immutable Vec # a three dimensional vector
    x::Float64
    y::Float64
    z::Float64
end
Vec(xyz::Vector{Float64}) = Vec(xyz...)
# various functions for manipulatin the vector type:
+(a::Vec, b::Vec) = Vec(a.x+b.x, a.y+b.y, a.z+b.z)
*(a::Float64, b::Vec) = Vec(a*b.x, a*b.y, a*b.z)
-(a::Vec, b::Vec) = Vec(a.x-b.x, a.y-b.y, a.z-b.z)
-(a::Vec) = *(-1.,a)
/(a::Vec, b::Float64) = Vec(a.x/b, a.y/b, a.z/b)
# like in linear algebra:
dot(a::Vec, b::Vec) = (a.x*b.x + a.y*b.y + a.z*b.z)
norm(a::Vec) = sqrt(dot(a, a))
unitize(a::Vec) = (a/norm(a))
myconvert(a::Vec) = [a.x,a.y,a.z] # just a helper function to convert this to a matrix
end

max_stars_repo_name: JuliaTagBot/borat

max_issues_repo_name: JuliaTagBot/borat

max_forks_repo_name: JuliaTagBot/borat


text:
include("../src/LarSurf.jl")
using Plasm
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
# using LarSurf

block_size = [5, 5, 5]
margin_size = 0

## Artifical data
data3d = ones(Int16, 19, 15, 31)
data3d[2:14,8:13,9:11] .= 10
# data3d

voxelsize_mm = [0.5, 0.9, 0.8]
threshold=0

blocks_number, blocks_number_axis = LarSurf.number_of_blocks_per_axis(
    size(data3d), block_size)

block1 = LarSurf.get_block(
    1, data3d, block_size, margin_size, blocks_number_axis, false
)

verts, trifaces = LarSurf.import_data3d(block1, voxelsize_mm)

println(verts)
println(trifaces)

println("======= To LAR ========")

V, EV, FE = LarSurf.to_lar(verts, trifaces)

println(V)
println(EV)
println(FE)

println("V: ", typeof(V), size(V))
println("EV: ", typeof(EV))
print("FV: ", typeof(FE))

@show EV;


text:
julia> join(reverse(collect(graphemes("as⃝df̅"))))
"f̅ds⃝a"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Documenter, MessyTimeSeriesOptim;
makedocs(
    sitename="MessyTimeSeriesOptim.jl",
    pages = [
        "index.md",
        "Estimation"  => ["man/estimation.md"],
        "Validation"  => ["man/validation.md"],
    ]
);

deploydocs(
    repo = "github.com/fipelle/MessyTimeSeriesOptim.jl.git",
)


text:
#= Currency info lookup functions =#

# Macro to define flexible lookup functions that accept type or Monetary object.
macro flexible(assignment)
    @assert assignment.head == :(=)
    symb = assignment.args[1].args[1]
    quote
        $(esc(assignment))
        $(esc(symb))(::Type{U}) where U <: Monetary{S} where S = $(esc(symb))(S)
        $(esc(symb))(m::Monetary) = $(esc(symb))(typeof(m))
    end
end

"""
    currency(m::Monetary) → Symbol

Return a symbol corresponding to the ISO 4217 currency code of the currency that
the given monetary amount is representing. For example, `currency(80USD)` will
return `:USD`. If the given monetary value is of a non-ISO 4217 currency, then
the returned symbol should contain only lowercase letters.

Prefer `iso4217alpha` to this function if a string is desired.
"""
currency(::Monetary{T}) where {T} = T

"""
    decimals(m::Monetary) → Int
    decimals(s::Symbol)   → Int
    decimals(d::DataType) → Int

Get the precision, in terms of the number of decimal places after the major
currency unit, of the given `Monetary` value or type. Alternatively, if given a
symbol, gets the default exponent (the number of decimal places to represent the
minor currency unit) for that symbol. Return `-1` if there is no sane minor
unit, such as for several kinds of precious metal.
"""
decimals(c::Symbol) = ISO4217[c][1]
decimals(::Monetary{T,U,V}) where {T,U,V} = V
decimals(::Type{Monetary{T,U,V}}) where {T,U,V} = V
decimals(::Type{T}) where {T<:Monetary} = decimals(filltype(T))

"""
Get a brief human-readable English-language description of the currency. The
description should begin with the common name of the currency, which should
describe it unambiguously (up to variations on the same currency). Optionally,
parentheses following the main description may include additional information
(such as the unit of a major currency unit).

This function may be called with either a symbol, a `Monetary` type, or a
`Monetary` object.
"""
function currencyinfo end
@flexible currencyinfo(s::Symbol) = ISO4217[s][2]

"""
    iso4217num(s::Symbol)   → Int
    iso4217num(m::Monetary) → Int
    iso4217num(t::DataType) → Int

Get the ISO 4217 numeric code for a currency. For custom currencies, a value of
`0` will be returned. This function may be called with either a symbol, a
`Monetary` type, or a `Monetary` object. Note that most applications should
zero-pad this code to three digits.
"""
function iso4217num end
@flexible iso4217num(s::Symbol) = ISO4217[s][3]

"""
    iso4217alpha(s::Symbol)   → String
    iso4217alpha(m::Monetary) → String
    iso4217alpha(t::DataType) → String

Get the ISO 4217 alphabetic code for a currency. For custom currencies, a
lowercase string will be returned, and this should not be interpreted as an ISO
4217 code. Otherwise, a three-letter uppercase string will be returned. This
function may be called with either a symbol, a `Monetary` type, or a `Monetary`
object.
"""
function iso4217alpha end
@flexible iso4217alpha(s::Symbol) = string(s) |> String

"""
    shortsymbol(s::Symbol)   → String
    shortsymbol(m::Monetary) → String
    shortsymbol(t::DataType) → String

Get a short, possibly ambiguous, commonly-used symbol for a currency. This
function may be called with either a symbol, a `Monetary` type, or a `Monetary`
object.
"""
function shortsymbol end
@flexible shortsymbol(s::Symbol) = get(SHORT_SYMBOL, s, iso4217alpha(s))

"""
    longsymbol(s::Symbol)   → String
    longsymbol(m::Monetary) → String
    longsymbol(t::DataType) → String

Get a commonly-used currency symbol for a currency, with at least enough
disambiguation to be non-ambiguous. This function may be called with either a
symbol, a `Monetary` type, or a `Monetary` object.
"""
function longsymbol end
@flexible longsymbol(s::Symbol) = get(LONG_SYMBOL, s, iso4217alpha(s))

"""
    majorunit(s::Symbol)   → Monetary{s}
    majorunit(m::Monetary) → typeof(m)
    majorunit(t::DataType) → t

Get the major unit of the currency. This function may be called with either a
symbol, a `Monetary` type, or a `Monetary` object.
"""
function majorunit end
majorunit(s::Symbol) = majorunit(Monetary{s})
majorunit(::Type{Monetary{T,U,V}}) where {T,U,V} = Monetary{T,U,V}(convert(U, 10)^V)
majorunit(::Type{T}) where {T<:Monetary} = majorunit(filltype(T))
majorunit(x::Monetary) = majorunit(typeof(x))


text:
"""
Abstract type that allows users to pass in their preferred Algorithm
"""
abstract AbstractFlowAlgorithm

"""
Forces the maximum_flow function to use the Edmonds–Karp algorithm.
"""
type EdmondsKarpAlgorithm <: AbstractFlowAlgorithm
end

"""
Forces the maximum_flow function to use Dinic\'s algorithm.
"""
type DinicAlgorithm <: AbstractFlowAlgorithm
end

"""
Forces the maximum_flow function to use the Boykov-Kolmogorov algorithm.
"""
type BoykovKolmogorovAlgorithm <: AbstractFlowAlgorithm
end

"""
Forces the maximum_flow function to use the Push-Relabel algorithm.
"""
type PushRelabelAlgorithm <: AbstractFlowAlgorithm
end

"""
Type that returns 1 if a forward edge exists, and 0 otherwise
"""

type DefaultCapacity <: AbstractArray{Int, 2}
    flow_graph::DiGraph
    nv::Int
    DefaultCapacity(flow_graph::DiGraph) = new(flow_graph, nv(flow_graph))
end

getindex(d::DefaultCapacity, s::Int, t::Int) = if has_edge(d.flow_graph, s , t) 1 else 0 end
size(d::DefaultCapacity) = (d.nv, d.nv)
transpose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))
ctranspose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))

"""
Constructs a residual graph for the input flow graph. Creates a new graph instead
of modifying the input flow graph.

The residual graph comprises of the same Vertex list, but ensures that for each
edge (u,v), (v,u) also exists in the graph. (to allow flow in the reverse direction).

If only the forward edge exists, a reverse edge is created with capacity 0. If both
forward and reverse edges exist, their capacities are left unchanged. Since the capacities
in DefaultDistance cannot be changed, an array of ones is created. Returns the
residual graph and the modified capacity_matrix (when DefaultDistance is used.)

Requires arguments:

- flow_graph::DiGraph,                    # the input graph
- capacity_matrix::AbstractArray{T,2}     # input capacity matrix
"""

residual(flow_graph::DiGraph) = DiGraph(Graph(flow_graph))

# Method for Edmonds–Karp algorithm

function maximum_flow{T<:Number}(
    flow_graph::DiGraph,                   # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    algorithm::EdmondsKarpAlgorithm        # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return edmonds_karp_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Dinic's algorithm

function maximum_flow{T<:Number}(
    flow_graph::DiGraph,                   # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    algorithm::DinicAlgorithm              # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return dinic_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Boykov-Kolmogorov algorithm

function maximum_flow{T<:Number}(
    flow_graph::DiGraph,                   # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    algorithm::BoykovKolmogorovAlgorithm   # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return boykov_kolmogorov_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Push-relabel algorithm

function maximum_flow{T<:Number}(
    flow_graph::DiGraph,                   # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2},   # edge flow capacities
    algorithm::PushRelabelAlgorithm        # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return push_relabel(residual_graph, source, target, capacity_matrix)
end

"""
Generic maximum_flow function. Requires arguments:

- flow_graph::DiGraph                   # the input graph
- source::Int                           # the source vertex
- target::Int                           # the target vertex
- capacity_matrix::AbstractArray{T,2}   # edge flow capacities
- algorithm::AbstractFlowAlgorithm      # keyword argument for algorithm
- restriction::T                        # keyword argument for a restriction

The function defaults to the Push-relabel algorithm. Alternatively, the algorithm
to be used can also be specified through a keyword argument. A default capacity of 1
is assumed for each link if no capacity matrix is provided.
If the restriction is bigger than 0, it is applied to capacity_matrix.

All algorithms return a tuple with 1) the maximum flow and 2) the flow matrix.
For the Boykov-Kolmogorov algorithm, the associated mincut is returned as a third output.

### Usage Example:

```julia

# Create a flow-graph and a capacity matrix
flow_graph = DiGraph(8)
flow_edges = [
    (1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
    (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
    (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
]
capacity_matrix = zeros(Int, 8, 8)
for e in flow_edges
    u, v, f = e
    add_edge!(flow_graph, u, v)
    capacity_matrix[u,v] = f
end

# Run default maximum_flow without the capacity_matrix
f, F = maximum_flow(flow_graph, 1, 8)

# Run default maximum_flow with the capacity_matrix
f, F = maximum_flow(flow_graph, 1, 8)

# Run Endmonds-Karp algorithm
f, F = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=EdmondsKarpAlgorithm())

# Run Dinic's algorithm
f, F = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=DinicAlgorithm())

# Run Boykov-Kolmogorov algorithm
f, F, labels = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=BoykovKolmogorovAlgorithm())

```
"""

function maximum_flow{T<:Number}(
    flow_graph::DiGraph,                   # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2} =  # edge flow capacities
        DefaultCapacity(flow_graph);
    algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
        PushRelabelAlgorithm(),
    restriction::T = zero(T)               # keyword argument for restriction max-flow
    )
    if restriction > zero(T)
      return maximum_flow(flow_graph, source, target, min(restriction, capacity_matrix), algorithm)
    end
    return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm)
end


text:
"""
Code for creating, loading and maniuplating line lists and masks.

Author: Eric Ford
Created: August 2020
"""

using DataFrames, CSV

"""
    Read line list in ESPRESSO csv format.
ESPRESSO format: lambda and weight.
Warning: ESPRESSO masks don't provide line depth and sometimes include one entry for a blend of lines.
"""
function read_linelist_espresso(fn::String) #; calcΔ::CCWT = default_calc_chunk_width) where { CCWT<:AbstractCalcChunkWidth }
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","weight"],delim=' ',ignorerepeated=true)
    @assert hasproperty(df, :lambda)
    @assert hasproperty(df, :weight)
    df[!,:lambda] .= λ_air_to_vac.(df[!,:lambda])
    return df
end

""" Read line list in VALD csv format.
   VALD format: lambda_lo, lambdaa_hi and depth.
"""
function read_linelist_vald(fn::String) # ; calcΔ::CCWT = default_calc_chunk_width) where { CCWT<:AbstractCalcChunkWidth }
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda_lo","lambda_hi","depth"])
    @assert hasproperty(df, :lambda_lo)
    @assert hasproperty(df, :lambda_hi)
    @assert hasproperty(df, :depth)
    df[!,:lambda_lo] .= λ_air_to_vac.(df[!,:lambda_lo])
    df[!,:lambda_hi] .= λ_air_to_vac.(df[!,:lambda_hi])
    df[!,:lambda] = sqrt.(df[!,:lambda_lo].*df[!,:lambda_hi])
    df[!,:weight] = df[!,:depth] # TODO: Decide out what we want to do about tracking depths and weights sepoarately
    return df
end

"""  Convert vacuum wavelength (in Å) to air wavelength
Ref: Donald Morton (2000, ApJ. Suppl., 130, 403) via
     https://www.astro.uu.se/valdwiki/Air-to-vacuum%20conversion
"""
function λ_vac_to_air(λ_vac::Real)
    @assert 3500 < λ_vac < 13000  # Making sure in Å for optical/NIR spectra.
    local s = 10000/λ_vac
    local n = 1 + 0.0000834254 + 0.02406147 / (130 - s^2) + 0.00015998 / (38.9 - s^2)
    return λ_vac/n
end

""" Convert air wavelength (in Å) to vacuum wavelength
Ref: https://www.astro.uu.se/valdwiki/Air-to-vacuum%20conversion
     VALD3 tools use the following solution derived by N. Piskunov
"""
function λ_air_to_vac(λ_air::Real)
    @assert 3500 < λ_air < 13000  # Making sure in Å for optical/NIR spectra.
    local s = 10000/λ_air
    local n = 1 + 0.00008336624212083 + 0.02408926869968 / (130.1065924522 - s^2) + 0.0001599740894897 / (38.92568793293 - s^2)
    λ_air*n
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
include("../../../good_plots/line_and_render.jl")

angle_vars = [15.0, 40.0, 85.0]

mass_fixed = 10
spin_fixed = 0.998
angle_fixed = 40.0

for angle in angle_vars
    plt = combined_plot(;
                        mass=mass_fixed,
                        spin=spin_fixed,
                        obs_angle=angle,
                        tolerance=1e-12,
                        dtmax=0.5,
                        size_multiplier=10,
                        resolution=3000,
                        nbins=200,
                        fov=10
                        )
    png(plt, "iron_line_combined angle=$angle.png")
end


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

using JuMP, Ipopt

m = Model()
@variable(m, x[1:3])

eps = 1e-4
scaler = 1.1
@objective(m, Min, sum(x[i] for i = 1:3) )
@NLconstraint(m, scaler * x[1] <= 0)
@NLconstraint(m, x[1] >= eps)
@NLconstraint(m, scaler^2 * x[2] >= 0)
@NLconstraint(m, x[2] <= eps)
@NLconstraint(m, scaler^3 * x[3] >= 0)
@NLconstraint(m, x[3] <= eps)

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
it = one_phase_solve(nlp_raw,my_pars);

IpoptSolve(nlp_raw)

#it.point.x


text:
using Hyperopt, Nonconvex, NLopt
using ChainRulesCore, ForwardDiff
using Plots
using CUTEst

function ChainRulesCore.rrule(::typeof(CUTEst.obj), nlp::CUTEstModel, x::AbstractVector)
    val = CUTEst.obj(nlp, x)
    grad = CUTEst.grad(nlp, x)
    val, Δ -> (NoTangent(), NoTangent(), Δ * grad)
end

function multistart(options)
    start_time = time()
    epoch_time, epoch_best = [], []
    function callback(results)
        t = time()-start_time
        push!(epoch_best, minimum(results))
        push!(epoch_time, t)
    end
    function f(nlp)
        if options == BOHB_options
            options = HyperoptOptions(sub_options = max_iter -> IpoptOptions(first_order = true, max_iter = max_iter),
                        sampler = Nonconvex.Hyperband(R=100, η=3, inner=BOHB(dims=[Continuous() for _ in 1:length(nlp.meta.x0)])))
        end
        model = Model(x -> CUTEst.obj(nlp, x))
        addvar!(model, nlp.meta.lvar, nlp.meta.uvar)
        alg = HyperoptAlg(IpoptAlg())
        result = Nonconvex.optimize(model, alg, nlp.meta.x0, options = options, callback=callback)
        result, epoch_time, epoch_best
    end
    return f
end

function get_fn(dims)
    function counting_ones(x::AbstractVector)
        target = 0.
        for (i, _x, dim) in zip(1:length(dims), x, dims)
            if dim isa Continuous
                target += _x
            else
                target += _x == 2 ? 1 : 0
            end
        end
        target
    end
    counting_ones
end

function multistart_counting_ones()
    dims = Array{Hyperopt.DimensionType}([Continuous() for _ in 1:50])
    append!(dims, [Categorical(2) for _ in 1:25])
    append!(dims, [UnorderedCategorical(2) for _ in 1:25])
    fn = get_fn(dims)
    options = HyperoptOptions(sub_options=max_iter->IpoptOptions(first_order=true, max_iter=max_iter),
                                sampler=Nonconvex.Hyperband(R=100, η=3, inner=BOHB(dims=dims)))
    model = Model(fn)
    addvar!(model, [0 for _ in 1:100], [1 for _ in 1:100])
    alg = HyperoptAlg(IpoptAlg())
    # define callback
    start_time = time()
    epoch_time, epoch_best = [], []
    function callback(results)
        t = time()-start_time
        push!(epoch_best, minimum(results))
        push!(epoch_time, t)
    end
    # optim
    result = Nonconvex.optimize(model, alg, [0 for _ in 1:100], options=options, callback=callback)
    result, epoch_time, epoch_best
end

function time_vs_regret(algorithm, problem)
    println("problem: $(problem)")
    if problem == "counting_ones"
        result, epoch_time, epoch_best = multistart_counting_ones()
    else
        nlp = CUTEstModel(problem)
        result, epoch_time, epoch_best = multistart(algorithm)(nlp)
        finalize(nlp)
    end
    return (result, epoch_time, epoch_best)
end

problems =
["counting_ones",
"DALLASL",
"HIER13",
"HIER16",
"READING3",
"READING1",
"SOSQP2",
"SREADIN3",
"SOSQP1",
"ZAMB2",]

random_options = HyperoptOptions(sub_options = IpoptOptions(), sampler = Hyperopt.RandomSampler())
LH_options = HyperoptOptions(sub_options = IpoptOptions(), sampler = Nonconvex.LHSampler())
hyperband_options = HyperoptOptions(sub_options = max_iter -> IpoptOptions(first_order = true, max_iter = max_iter),
                            sampler = Hyperopt.Hyperband(R=100, η=3, inner=Hyperopt.RandomSampler()))
BOHB_options = HyperoptOptions(sub_options = max_iter -> IpoptOptions(first_order = true, max_iter = max_iter), sampler = Hyperopt.Hyperband(R=100, η=3, inner=Hyperopt.RandomSampler()))

algorithms = Dict( "hyperband" => hyperband_options, "BOHB" => BOHB_options, "random" => random_options, "LH" => LH_options)

results = Dict()
epoch_times = Dict()
epoch_bests = Dict()
for problem in problems
    plot(title=problem, xaxis=:log, xlabel="time elapsed", ylabel="regret")
    for algorithm in keys(algorithms)
        println("algorithm: $(algorithm), problem $(problem): ")
        linewidth = (problem ∈ ["random", "LH"] ? 10 : 2)
        results[algorithm, problem], epoch_times[algorithm, problem], epoch_bests[algorithm, problem] = time_vs_regret(algorithms[algorithm], problem)
        plot!(epoch_times[algorithm, problem], epoch_bests[algorithm, problem], label=algorithm, lw=linewidth)
    end
    savefig("$(problem).svg")
end

max_stars_repo_name: noilreed/notebooks

max_issues_repo_name: noilreed/notebooks

max_forks_repo_name: noilreed/notebooks


text:
"""
Clamped cylinder, transient vibration.

The cylinder is banged at the initial time (i.e. it is given an initial velocity).
"""
module clamp_cyl_transient_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using SymRCM
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot; #@gp "clear"
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection
using ThreadedSparseCSR
using UnicodePlots
using InteractiveUtils
using BenchmarkTools

const E = 200e9;
const nu = 0.3;
const rho = 7800.0
const R = 0.1
const L = 0.8
const ksi = 0.0
const omegad = 2*pi*1000
const color = "red"
const omegaf = 2*pi*10^5
const tend = 25 * (2*pi/omegaf)
const visualize = true
const visualizeclear = true
const visualizevtk = !true
const distributedforce = !true

function _pwr(K, M, maxit = 30, rtol = 1/10000)
    invM = fill(0.0, size(M, 1))
    invM .= 1.0 ./ (vec(diag(M)))
    v = rand(size(M, 1))
    w = fill(0.0, size(M, 1))
    everyn = Int(round(maxit / 50)) + 1
    lambda = lambdap = 0.0
    for i in 1:maxit
        ThreadedSparseCSR.bmul!(w, K, v)
        wn = norm(w)
        w .*= (1.0/wn)
        v .= invM .* w
        vn = norm(v)
        v .*= (1.0/vn)
        if i % everyn  == 0
            lambda = sqrt((v' * (K * v)) / (v' * M * v))
            # @show i, abs(lambda - lambdap) / lambda
            if abs(lambda - lambdap) / lambda  < rtol
                break
            end
            lambdap = lambda
        end
    end
    return lambda
end

function _cd_loop!(M, K, ksi, U0, V0, tend, dt, force!, peek)
    # Central difference integration loop, for mass-proportional Rayleigh damping.
    U = similar(U0)
    V = similar(U0)
    A = similar(U0)
    F = similar(U0)
    E = similar(U0)
    C = similar(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = similar(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end
    dt2_2 = ((dt^2)/2)
    dt_2 = (dt/2)

    t = 0.0
    @. U = U0; @. V = V0; # Initial Conditions
    A .= invMC .* force!(F, t); # Compute initial acceleration
    peek(0, U, V, t)
    for step in 1:nsteps
        t = t + dt
        @. U += dt*V + dt2_2*A;   # Displacement update
        force!(F, t);  # External forces are computed
        ThreadedSparseCSR.bmul!(E, K, U)  # Evaluate elastic restoring forces
        @. F -= E + C * (V + dt_2 * A)  # Calculate total force
        @. V += dt_2 * A; # Update the velocity, part one: add old acceleration
        @. A = invMC * F # Compute the acceleration
        @. V += dt_2 * A; # Update the velocity, part two: add new acceleration
        peek(step, U, V, t)
    end
end

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _set_up_model(n = 64, thickness = 0.01, drilling_stiffness_scale = 1.0)
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = drilling_stiffness_scale

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    return fens, fes, geom0, dchi, K, M
end

function _execute(n = 64, thickness = 0.01)
    fens, fes, geom0, dchi, K, M = _set_up_model(n, thickness)

    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Compute the maximum eigenvalue to determine the stable time step
    omega_max = _pwr(K, M, 300)

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        if t < 2 * (2*pi/omegaf)
            F .= sin(omegaf*t) .* Fmag
        else
            F .= 0.0
        end
        # @show norm(F)
        return F
    end

    # Suddenly applied constant force at a node
    # Fmag = fill(0.0, dchi.nfreedofs)
    # Fmag[qpointdof] = -1.0
    # function force!(F, t)
    #     F .= Fmag
    #     return F
    # end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, V, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    _cd_loop!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        if visualizevtk
            @info "Dumping visualization"
            times = Float64[]
            vectors = []
            for i in 1:length(displacements)
                scattersysvec!(dchi, displacements[i])
                push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
                push!(times, i*dt*nbtw)
            end
            vtkwritecollection("clamp_cyl_forc_damp_expl_$n", fens, fes, times; vectors = vectors)
        end
    end
end

function _test_stable_time_step(n = 64, thickness = 0.01, drilling_stiffness_scale = 1.0)
    fens, fes, geom0, dchi, K, M = _set_up_model(n, thickness, drilling_stiffness_scale)

    omega_max = _pwr(K, M, 300)

    # evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)

    # omega_max = sqrt(evals[1])

    dt = Float64(0.9* 2/omega_max)
    return dt
end

function test(ns = [4*64])
    @info "Clamped cylinder: transient vibration"
    for n in ns
        _execute(n, 0.01)
    end
    return true
end

function test_stable_time_step()
    @info "Clamped cylinder: test stable time step"
    results = []
    for t in [0.01, 0.002, 0.001, 0.0001, 0.00001]
        dts = []
        for d in [0.01, 0.1, 1.0, 1.5, 2.0, 2.5, 5.0, 10.0]
            dt = _test_stable_time_step(2*64, t, d)
            push!(dts, dt)
        end
        push!(results, (t, dts))
    end
    return results
end

function allrun(ns = [2*64])
     println("#####################################################")
    println("# test ")
    test(ns)
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
function fixedpointiteration(g::Function, xi::Number, ea::Number, maxIter::Int64)
  # Estiamte roots for an equation using the fixed-point iteration method
  for i = 1:maxIter
    xiold = xi;
    xi = g(xi);

    if (abs(xi - xiold) / xi) < ea
      break;
    end;
  end

  return xi;
end

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_name: alexjohnj/numerical-methods


text:
#neuralNet type: Defines the structure of the neural net
#and allocates arrays that will be used over and over again
#during training
mutable struct neuralNet
    nNeurons::Array{Int64} #number of neurons by level (including input and output levels)
    nLevels::Int64         #number of levels

    x::Array{Float64} #activations (one per neuron, including inputs and outputs)
    w::Array{Float64} #weights
    b::Array{Float64} #biases (one per neuron, excluding inputs)

    xIdx::Array{UnitRange{Int64},1} #indices of x values by level (includes input level)
    wIdx::Array{UnitRange{Int64},1} #indices of weights by level (excludes input level)
    bIdx::Array{UnitRange{Int64},1} #indices of biases by level (excludes input level)

    nWeights::Int64               #total number of weights
    nWeightsByLevel::Array{Int64} #number of weights by level

    da::Array{Float64}   #derivative of the activation function at each non-input neuron
    dCdb::Array{Float64} #derivative of the cost function wrt biases
    dCdw::Array{Float64} #derivative of the cost function wrt weights

    neuralNet() = new()

    function neuralNet(nNeurons::Array{Int64})
        nn = neuralNet();
        nn.nNeurons = nNeurons;

        nn.nLevels = length(nNeurons);

        #x value indices
        xIdx    = [0; cumsum(nNeurons)[1:end-1]].+1;
        nn.xIdx = [ xIdx[l]:(xIdx[l]+nn.nNeurons[l]-1) for l=1:nn.nLevels ];

        #weight indices
        nn.nWeightsByLevel = [ nNeurons[l]*nNeurons[l-1] for l=2:nn.nLevels ];
        nn.nWeights = sum( nn.nWeightsByLevel );
        wIdx = [0; cumsum(nn.nWeightsByLevel)[1:end-1]].+1;
        nn.wIdx = [ wIdx[l]:(wIdx[l]+nn.nWeightsByLevel[l]-1) for l=1:nn.nLevels-1 ];

        #biase indices
        nBiases = sum(nNeurons[2:end]);
        bIdx    = [0; cumsum(nNeurons[2:end-1])].+1;
        nn.bIdx = [ bIdx[l]:(bIdx[l]+nn.nNeurons[l+1]-1) for l=1:nn.nLevels-1 ];

        #allocate
        nn.x    = zeros(sum(nNeurons)); #activations & inputs
        nn.w    = zeros(nn.nWeights);   #weights
        nn.b    = zeros(nBiases);       #biases

        #allocate (use for back propagation)
        nn.da   = zeros(nBiases);      #derivative of activation function
        nn.dCdb = zeros(nBiases);      #derivative of cost wrt biases (delta)
        nn.dCdw = zeros(nn.nWeights);  #derivative of cost wrt weights

        return nn;
    end
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
function rand_datum()
  x = rand()
  return x > 0.5 ? "WGS-84 " : x > 0.25 ? "ETRS89 " : x > 0.1 ? "GRS 80 " : "JGD2011"
end

# Things that work for both regularly and irregularly sampled data
function randResp(n::Int64=0)
  if n > 0
    T = Float64
    i = max(1, div(n,2))
  else
    T = rand() < 0.5 ? Float32 : Float64
    i = rand(2:2:8)
  end
  zstub = zeros(T, 2*i)
  pstub = 10 .*rand(T, i)
  if T == Float32
    resp = PZResp(a0 = 1.0f0, f0 = 1.0f0, p = vcat(pstub .+ pstub.*im, pstub .- pstub*im), z = complex.(zstub))    # resp
  else
    resp = PZResp64(a0 = 1.0, f0 = 1.0, p = vcat(pstub .+ pstub.*im, pstub .- pstub*im), z = complex.(zstub))    # resp
  end
  return resp
end

function rand_t(Lx::Int64, n::Int64=0)
  ts = time()-86400+randn()
  L = n == 0 ? rand(0:9) : n
  t = zeros(2+L, 2)

  # first row is always start time
  t[1,:] = [1 round(Int64, ts/μs)]

  # rest are random time gaps
  if L > 0
    gaps = unique(rand(2:Lx, L, 1))

    # try to create exactly L gaps
    j = 0
    while j < 5
      j = j + 1
      gaps_tmp = unique(rand(2:Lx, L, 1))
      if length(gaps_tmp) > length(gaps)
        gaps = gaps_tmp
      end
      (length(gaps) == L) && break
    end
    L = length(gaps)

    t[2:L+1,1] = gaps
    for i = 2:L+1
      t[i,2] = round(Int64, (rand(1:100)+rand())*sμ)
    end

    # control for gap in last sample
    if any(t[:,1].==Lx) == true
      t = t[1:L+1,:]
    else
      t[L+2,:] = [Lx 0]
    end
  else
    t[2,:] = [Lx 0]
  end
  return sortslices(t, dims=1)
end

function randLoc(randtype::Bool=true)
  y = randtype == true ? rand() : 1.0
  datum = rand_datum()
  loc = [ 180*(rand()-0.5),
          360*(rand()-0.5),
          1000*rand(),
          100*rand(),
          360*(rand()-0.5),
          180*(rand()-0.5) ]
  if y > 0.5
    return GeoLoc(datum, loc...)
  else
    g = GenLoc(loc)
    g.datum = datum
    return g
  end
end

function pop_chan_tail!(Ch::SeisChannel)
  setfield!(Ch, :name, randstring(rand(12:64)))       # name
  setfield!(Ch, :gain,  rand()*10^rand(0:10))         # gain
  setfield!(Ch, :loc, randLoc())                      # loc
  pop_rand_dict!(Ch.misc, rand(4:24))                 # misc
  note!(Ch, "Created by function populate_chan!.")    # notes
end

# Populate a channel with irregularly-sampled (campaign-style) data
function populate_irr!(Ch::SeisChannel; nx::Int64=0)
  chan = "OY"*randstring('A':'Z',1)
  net = ur2()
  sta = uppercase(randstring('A':'Z', rand(1:5)))
  loc = ur2()

  ts = round(Int, sμ*(time()-86400+randn()))
  if nx == 0
    L = 2^rand(6:12)
  else
    L = nx
  end
  Ls = rand(1200:7200)

  Ch.id     = join([net,sta,loc,chan],'.')
  Ch.fs     = 0
  Ch.src    = string("randSeisChannel(c=true, nx=",  nx, ")")
  Ch.units  = rand(irregular_units)
  Ch.t      = hcat(collect(1:1:L), ts.+sort(rand(UnitRange{Int64}(1:Ls), L)))
  Ch.x      = (rand(L) .- (rand(Bool) == true ? 0.5 :  0.0)).*(10 .^ (rand(1:10, L)))
  pop_chan_tail!(Ch)
  return nothing
end

# Populate a channel with regularly-sampled (time-series) data
function populate_chan!(Ch::SeisChannel; s::Bool=false, nx::Int64=0)

  bcodes  = Char['V', 'L', 'M', 'M', 'B', 'S', 'S', 'S', 'S', 'S', 'S', 'H', 'S', 'E', 'E', 'C']

  # Ch.fs
  fs = rand(fs_vals)
  fc = rand(fc_vals[fc_vals .< fs/2])

  # An empty ID generates codes and units to match values real data might have
  bc = getbandcode(fs)
  (ic, cc, units) = getyp2codes(bc, s)
  sta = randstring('A':'Z', rand(3:5))
  loc = rand() < 0.5 ? "" : ur2()
  cha = string(bc, ic, cc)

  Ch.id     = join([ur2(), sta, loc, cha],'.')                        # id
  Ch.fs     = fs                                                      # fs
  Ch.units  = units                                                   # units
  Ch.resp   = randResp()                                              # resp
  Ch.src    = string("randSeisChannel(c=false, nx=",  nx, ")")        # src

  # random noise for data, with random time gaps
  if nx == 0
    Ls = rand(1200:7200)
    Lx = ceil(Int, Ls*Ch.fs)
  else
    Lx = nx
  end
  Ch.x = randn(rand() < 0.5 ? Float32 : Float64, Lx)                  # x
  Ch.t = rand_t(Lx)                                                   # t

  pop_chan_tail!(Ch)
  return nothing
end

"""
    randSeisChannel()

Generate a random channel of geophysical time-series data as a SeisChannel.

    randSeisChannel(c=true)

Generate a random channel of irregularly-sampled data.

    randSeisChannel(s=true)

Generate a random channel of regularly-sampled seismic data.
"""
function randSeisChannel(; c::Bool=false, s::Bool=false, nx::Int64=0)
  Ch = SeisChannel()
  if c == true
    populate_irr!(Ch, nx=nx)
  else
    populate_chan!(Ch, s=s, nx=nx)
  end
  return Ch
end


text:
# this is for presynaptic , covers also Poisson inputs generators, etc
abstract type PSSpikingType{NT} <: PopulationState{NT} end

# this is for postsynaptic that integrate inputs over currents
# LIF, EIF, QIF, etc
abstract type PSGeneralCurrentIFType{NT} <: PSSpikingType{NT} end

# connection for IF neurons
# must keep track of internal trace (current/conductance)
struct ConnGeneralIF{N,TP<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  post_trace::Vector{Float64}
  plasticities::TP
end
function ConnGeneralIF(weights::SparseMatrixCSC)
  npost=size(weights,1)
  return ConnGeneralIF(weights,zeros(Float64,npost),())
end
function ConnGeneralIF(weights::SparseMatrixCSC,
    (plasticities::PlasticityRule)...)
  npost=size(weights,1)
  return ConnGeneralIF(weights,zeros(Float64,npost),plasticities)
end
function reset!(conn::ConnGeneralIF)
  reset!.(conn.plasticities)
  fill!(conn.post_trace,0.0)
  return nothing
end

# or two internal traces (plus and minus)
struct ConnGeneralIF2{N,TP<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  post_trace1::Vector{Float64}
  post_trace2::Vector{Float64}
  plasticities::TP
end
function ConnGeneralIF2(weights::SparseMatrixCSC)
  npost=size(weights,1)
  return ConnGeneralIF2(weights,zeros(Float64,npost),zeros(Float64,npost),())
end
function ConnGeneralIF2(weights::SparseMatrixCSC,
    (plasticities::PlasticityRule)...)
  npost=size(weights,1)
  return ConnGeneralIF2(weights,zeros(Float64,npost),zeros(Float64,npost),
    plasticities)
end
function reset!(conn::ConnGeneralIF2)
  reset!.(conn.plasticities)
  fill!(conn.post_trace1,0.0)
  fill!(conn.post_trace2,0.0)
  return nothing
end

# This is an "infinitely strong" voltage connection
# if pre spikes, and *any* weight is present
# then  post will also spike (no delay for now)
# it cannot be plastic, but one can still modify the connectivity matrix
struct ConnSpikeTransfer{N,TP<:NTuple{N,PlasticityRule}} <: Connection{N}
  weights::SparseMatrixCSC{Float64,Int64}
  plasticities::TP
  function ConnSpikeTransfer(weights::SparseMatrixCSC)
    return  new{0,Tuple{}}(weights,())
  end
end

# see ConnSpikeTransfer comments
function forward_signal!(t_now::Real,dt::Real,
      pspost::PSGeneralCurrentIFType,
      conn::ConnSpikeTransfer,pspre::PSSpikingType)
	post_idxs = rowvals(conn.weights) # postsynaptic neurons
	for _pre in findall(pspre.isfiring)
		_posts_nz = nzrange(conn.weights,_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
      post_idx = post_idxs[_pnz]
			pspost.input[post_idx] = Inf
		end
	end
  return nothing
end

function forward_signal!(t_now::Real,dt::Real,
      pspost::PSGeneralCurrentIFType,
      conn::ConnGeneralIF,pspre::PSSpikingType)
	post_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
	τ_decay = pspre.neurontype.τ_output_decay
	for _pre in findall(pspre.isfiring)
		_posts_nz = nzrange(conn.weights,_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
			post_idx = post_idxs[_pnz]
			# update the post currents even when refractory
			conn.post_trace[post_idx] +=
						weightsnz[_pnz] / τ_decay
		end
	end
	# add the currents to postsynaptic input
	# ONLY non-refractory neurons
	post_refr = findall(pspost.isrefractory) # refractory ones
	@inbounds @simd for i in eachindex(pspost.input)
		if !(i in post_refr)
			pspost.input[i] += conn.post_trace[i]
		end
	end
  # finally, postsynaptic currents decay in time
	@inbounds @simd for i in eachindex(conn.post_trace)
		conn.post_trace[i] -= dt*conn.post_trace[i] / τ_decay
	end
  return nothing
end

# examines state_now, if above threshold, isfiring and refractory turn true
# and potential is reset. Finally, removes expired refractoriness
function _spiking_state_update!(state_now::Vector{R},
    isfiring::BitArray{1},isrefractory::BitArray{1},
    last_fired::Vector{R},
    t_now::R, t_refractory::R,
    v_threshold::R, v_reset::R) where R<:Real
  reset_spikes!(isfiring)
	@inbounds @simd for i in eachindex(state_now)
		if state_now[i] > v_threshold
			state_now[i] =  v_reset
			isfiring[i] = true
			last_fired[i] = t_now
			isrefractory[i] = true
		# check only when refractory
		elseif isrefractory[i] &&
				( (t_now-last_fired[i]) >= t_refractory)
			isrefractory[i] = false
		end
	end
  return nothing
end