text:
##
## Phenotype file (not standardized, with header and family-ID): phenotypes.reordered.csv
## Final phenotype file (standardized, without header/family-ID): phenotypes.reordered.standardized.csv
## Final covariates file (with family ID, not standardized): covariates.reordered.csv
## Final covariates file (without family ID, standardized, extra intercept column): covariates.reordered.standardized.csv
## Final genotype file: ukb.merged.metabolic.subset.european.400K.QC.bed
##

#
# mIHT test run (K = 5000, no cross validation)
#
using DelimitedFiles, MendelIHT, SnpArrays, LinearAlgebra
BLAS.set_num_threads(1)
y = readdlm("phenotypes.reordered.standardized.csv", ',')
z = readdlm("covariates.reordered.standardized.csv", ',')
standardize!(@view(z[:, 2:end]))
x = SnpArray("ukb.merged.metabolic.subset.european.400K.QC.bed")
xla = SnpLinAlg{Float64}(x, center=true, scale=true, impute=true)
Yt = Matrix(Transpose(y))
Zt = Matrix(Transpose(z))
fit_iht(Yt, Transpose(xla), Zt; verbose=true,
    k=5000, d = MvNormal(), max_iter=10, max_step=3, init_beta=false)

#
# running mIHT
#
using MendelIHT, Random, LinearAlgebra
BLAS.set_num_threads(1)
Random.seed!(2022)
plinkfile = "ukb.merged.metabolic.subset.european.400K.QC"
phenotypes = "phenotypes.reordered.standardized.csv"
covariates = "covariates.reordered.standardized.csv"

path = 1000:1000:10000
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, min_iter=10,
    cv_summaryfile="cviht.summary.roughpath1.txt")

k_rough_guess = path[argmin(mses)]
path = (k_rough_guess - 900):100:(k_rough_guess + 900)
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, min_iter=10,
    cv_summaryfile="cviht.summary.roughpath2.txt")

k_rough_guess = path[argmin(mses)]
path = (k_rough_guess - 90):10:(k_rough_guess + 90)
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, min_iter=10,
    cv_summaryfile="cviht.summary.roughpath3.txt")

k_rough_guess = path[argmin(mses)]
path = (k_rough_guess - 9):(k_rough_guess + 9)
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, min_iter=10,
    cv_summaryfile="cviht.summary.final.txt")

K = path[argmin(mses)]
@time iht_result = iht(plinkfile, K, MvNormal,
    summaryfile = "iht.final.summary.txt",
    betafile = "iht.final.beta.txt",
    covariancefile = "iht.final.cov.txt",
    covariates=covariates, phenotypes=phenotypes, max_iter=2000)


text:
using RunMetrics, Statistics, Plots, LinearAlgebra, GLM
old_dir = pwd()
example_dir = joinpath(@__DIR__, "..", "sample_data", "multi")

files = readdir(example_dir)

cd(example_dir)

tcx_files = files[findall(x -> endswith(x, ".tcx") && isrun(x), files)]
deleteat!(tcx_files, 6)
deleteat!(tcx_files, 9)

run_sums = [parse_run(tcx) for tcx in tcx_files]

rs = run_sums[1]
# window = 90s seems best (mse = 29.6)
# α = 0.006 seems best (mse = 37.0)
y, X = RunMetrics.design_mat(run_sums, 90.0,
                             decay_function = RunMetrics.sliding_scale,
                             standardize = true)

keep_inds = findall(x -> x > 60.0, X[:, 1])

# X[:, 1] .= log.(X[:, 1])

y = y[keep_inds]
X = X[keep_inds, :]

X = [log.(X[:, 1]) X]

first_run = 5
n, p = size(X)

function logitize(x::Float64, bounds::Tuple{Float64, Float64})
    return (x - bounds[1]) / (bounds[2] - bounds[1])
end

y = logitize.(y, Ref((50.0, 200.0)))

df = DataFrame([y X], [["hr", "log time", "time", "h_speed", "v_speed"]; ["run$i" for i = first_run:p]])
nms = Symbol.(names(df))
form = FormulaTerm(Term(nms[1]), Tuple(Term(nms[i]) for i = 2:p))

ft = glm(form, df, Binomial(), LogitLink())
y_hat = predict(ft)

# X_cont = @view X[:, 1:(first_run - 1)]
# X_cont .= (X_cont - ones(size(X, 1)) * mean(X_cont, dims=1)) * inv(Diagonal(vec(std(X_cont, dims=1))))

# β = inv(X' * X) * X' * y

# y_hat = X * β
mse = var(y_hat - y)
r2 = 1.0 - mse / var(y)
@show mse
@show r2
cd(@__DIR__)

cd(old_dir)

scatter(y_hat, y)

run_inds = [findall(isequal(1.0), @view(X[:, i])) for i = first_run:size(X, 2)]
errs = y - y_hat
p = plot()
for i = 1:length(run_inds)
    p = plot!(run_inds[i], errs[run_inds[i]], legend=false)
end
p = plot!(x -> 0)

ext: jl


text:
const _FESPACE_ALLOWED_KWARGS = [
  :reffe,
  :conformity,
  :order,
  :labels,
  :valuetype,
  :model, :triangulation,
  :dirichlet_tags, :dirichlet_masks,
  :dof_space,
  :constraint, :zeromean_trian, :zeromean_quad,
]

supported_kwargs(::Type{FESpace}) = _FESPACE_ALLOWED_KWARGS

function check_all_kwargs_supported(f, given_kw)
  allowed_kw = supported_kwargs(f)
  for key in keys(given_kw)
    if !(key in allowed_kw)
      msg = """$f got unsupported keyword $(sprint(print, key)). Allowed keywords are:
      $allowed_kw
      """
      throw(ArgumentError(msg))
    end
  end
end

"""
    $(SIGNATURES)

Construct a `FESpace`. Supported keywords are:
$(supported_kwargs(FESpace))
"""
function FESpace(;kwargs...)
  check_all_kwargs_supported(FESpace, kwargs)

  constraint = _get_kwarg(:constraint,kwargs,nothing)
  reffe = _get_kwarg(:reffe,kwargs)

  fespace = nothing

  if isa(reffe,ReferenceFE)

    fespace = _setup_generic_space(kwargs)

  elseif reffe in [:Lagrangian,:QLagrangian,:PLagrangian,:SLagrangian]

    fespace = _setup_lagrange_spaces(kwargs)

  elseif reffe == :RaviartThomas

    fespace = _setup_hdiv_space(kwargs)

  elseif reffe == :Nedelec

    fespace = _setup_hcurl_space(kwargs)

  else
    @notimplemented "Unsupported reffe $reffe."
  end

  @assert fespace != nothing

  restricted_at = _get_restricted_triangulation(kwargs)
  if restricted_at == nothing
    _fespace = fespace
  else
    @assert isa(restricted_at,RestrictedTriangulation)
    _fespace = ExtendedFESpace(fespace,restricted_at)
  end

  if constraint == nothing
    return _fespace

  elseif constraint == :zeromean
    order = _get_kwarg(:order,kwargs)
    model = _get_kwarg(:model,kwargs,nothing)
    zeromean_trian = _get_kwarg(:zeromean_trian,kwargs,nothing)
    zeromean_quad = _get_kwarg(:zeromean_quad,kwargs,nothing)
    if zeromean_trian !== nothing
      trian = zeromean_trian
    elseif model === nothing
      trian = _get_kwarg(:triangulation,kwargs)
    else
      trian = get_triangulation(model)
    end
    if zeromean_quad === nothing
      quad = CellQuadrature(trian,order)
    else
      quad = zeromean_quad
    end
    return ZeroMeanFESpace(_fespace,trian,quad)

  else
    @unreachable "Unknown constraint value $constraint"
  end

end

function _get_restricted_triangulation(kwargs)
  model  = _get_kwarg(:model,kwargs,nothing)
  _trian = _get_kwarg(:triangulation,kwargs,nothing)
  if _trian == nothing
    if model == nothing
      @unreachable "either a model or a triangulation has to be provided for building a FESpace"
    end
    trian = get_triangulation(model)
  else
    if model != nothing
      @unreachable "either a model or a triangulation BUT NOT BOTH has to be provided for building a FESpace"
    end
    trian = _trian
  end
  _get_restricted_triangulation_or_nothing(trian)
end

function  _get_restricted_triangulation_or_nothing(trian)
  nothing
end

function  _get_restricted_triangulation_or_nothing(trian::RestrictedTriangulation)
  trian
end

"""
"""
function TestFESpace(;kwargs...)
  FESpace(;kwargs...)
end

function _setup_generic_space(kwargs)

  reffe = _get_kwarg(:reffe,kwargs)
  model = _get_kwarg(:model,kwargs)
  labels = _get_kwarg(:labels,kwargs,get_face_labeling(model))
  conformity = _get_kwarg(:conformity,kwargs,true)
  diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
  order = _get_kwarg(:order,kwargs,nothing)
  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
  ( dofspace == :reference ? true : false )

  is_ref = (dofspace==:reference)

  Tf = _get_kwarg(:valuetype,kwargs,VectorValue{1,Float64})
  T = eltype(Tf)

  polytopes = get_polytopes(model)
  if length(polytopes) > 2 || polytopes[1] != get_polytope(reffe)
    @unreachable "RefFE and geometrical model are not compatible in generic FESpace constructor"
  end

  reffes = [ reffe ]

  if conformity in [true, :default]
    conf = get_default_conformity(reffe)
    V =  ConformingFESpace(reffes,conf,model,labels,diritags,nothing,is_ref)
  elseif isa(conformity,Conformity)
    V =  ConformingFESpace(reffes,conformity,model,labels,diritags,nothing,is_ref)
  else
    s = "Unsuported conformity = $conformity"
    @unreachable s
  end

  V
end

function _setup_hdiv_space(kwargs)

  reffe = _get_kwarg(:reffe,kwargs)
  model = _get_kwarg(:model,kwargs)
  labels = _get_kwarg(:labels,kwargs,get_face_labeling(model))
  conformity = _get_kwarg(:conformity,kwargs,true)
  diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
  order = _get_kwarg(:order,kwargs,nothing)
  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
  ( dofspace == :reference ? true : false )

  is_ref = (dofspace==:reference)

  Tf = _get_kwarg(:valuetype,kwargs,VectorValue{1,Float64})
  T = eltype(Tf)

  if order == nothing
    @unreachable "order is a mandatory keyword argument in FESpace constructor for RaviartThomas reference FEs"
  end

  polytopes = get_polytopes(model)
  reffes = [RaviartThomasRefFE(T,p,order) for p in polytopes]

  if conformity in [true, :default, :HDiv, :Hdiv]
    V =  ConformingFESpace(reffes,DivConformity(),model,labels,diritags,nothing,is_ref)
  else
    s = "Conformity $conformity not implemented for $reffe reference FE on polytopes $(polytopes...)"
    @unreachable s
  end

  V
end

function _setup_hcurl_space(kwargs)

  reffe = _get_kwarg(:reffe,kwargs)
  model = _get_kwarg(:model,kwargs)
  labels = _get_kwarg(:labels,kwargs,get_face_labeling(model))
  conformity = _get_kwarg(:conformity,kwargs,true)
  diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
  order = _get_kwarg(:order,kwargs,nothing)
  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
  Tf = _get_kwarg(:valuetype,kwargs,VectorValue{1,Float64})
  T = eltype(Tf)

  is_ref = dofspace==:reference

  if order == nothing
    @unreachable "order is a mandatory keyword argument in FESpace constructor for Nedelec reference FEs"
  end

  polytopes = get_polytopes(model)
  reffes = [NedelecRefFE(T,p,order) for p in polytopes]

  if conformity in [true, :default, :HCurl, :Hcurl]
    V =  ConformingFESpace(reffes,CurlConformity(),model,labels,diritags,nothing,is_ref)
  else
    s = "Conformity $conformity not implemented for $reffe reference FE on polytopes $(polytopes...)"
    @unreachable s
  end

  V
end

function _setup_lagrange_spaces(kwargs)

  conformity = _get_kwarg(:conformity,kwargs,true)
  reffe = _get_kwarg(:reffe,kwargs)
  order = _get_kwarg(:order,kwargs)
  T = _get_kwarg(:valuetype,kwargs,nothing)
  diritags = _get_kwarg(:dirichlet_tags,kwargs,Int[])
  dirimasks = _get_kwarg(:dirichlet_masks,kwargs,nothing)
  # dirimasks = nothing
  dofspace = _get_kwarg(:dof_space,kwargs,:reference)
  labels = _get_kwarg(:labels,kwargs,nothing)
  # labels = nothing
  model = _get_kwarg(:model,kwargs,nothing)

  is_ref = (dofspace==:reference)

  if T == nothing
    @unreachable "valuetype is a mandatory keyword argument in FESpace constructor for Lagrangian reference FEs"
  end

  if conformity in [false, :L2]

    s = "Strong dirichlet conditions cannot be imposed in discontinuous spaces for the moment"
    @notimplementedif diritags != Int[] s
    @notimplementedif dirimasks != nothing s

    _trian = _get_kwarg(:triangulation,kwargs,nothing)
    if _trian == nothing
      if model == nothing
        @unreachable "either a model or a triangulation has to be provided for building a discontinuous Lagrangian space"
      end
      trian = get_triangulation(model)
    else
      if model != nothing
        @unreachable "either a model or a triangulation BUT NOT BOTH has to be provided for building a discontinuous Lagrangian space"
      end
      trian = _trian
    end

    polytopes = [get_polytope(r) for r in get_reffes(trian)]

    if _is_reffe_lagrangian_compatible_with_polytopes(reffe,polytopes)
      if reffe == :SLagrangian
        _reffes = [SerendipityRefFE(T,p,order) for p in polytopes]
      else
        _reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
      end
    else
      if reffe == :PLagrangian
        _reffes = [LagrangianRefFE(T,p,order;space=:P) for p in polytopes]
      else
        @unreachable "Not possible to use a $reffe reffe on polytopes $(polytopes...)"
      end
    end

    return  DiscontinuousFESpace(_reffes,trian,is_ref)

  elseif conformity in [true, :default, :H1, :C0]

    if model == nothing
      @unreachable "model is a mandatory keyword argument in FESpace constructor for conforming Lagrangian reference FEs"
    end
    polytopes = get_polytopes(model)
    trian = get_triangulation(model)
    if ! _is_reffe_lagrangian_compatible_with_polytopes(reffe,polytopes)
      s = "Conformity $conformity not possible for $reffe reference FE on polytopes $(polytopes...)"
      @unreachable s
    end
    if reffe == :SLagrangian
      _reffes = [SerendipityRefFE(T,p,order) for p in polytopes]
    else
      _reffes = [LagrangianRefFE(T,p,order) for p in polytopes]
    end
    if labels == nothing
      return GradConformingFESpace(_reffes,model,diritags,dirimasks,is_ref)
    else
      return ConformingFESpace(_reffes,GradConformity(),model,labels,diritags,dirimasks,is_ref)
    end

  else
    s = "Conformity $conformity not implemented for lagrangian reference FEs"
    @unreachable s
  end

end

function _get_kwarg(kwarg,kwargs)
  try
    return kwargs[kwarg]
  catch
    s = "The key-word argument $(kwarg) is mandatory in the FESpace constructor"
    error(s)
  end
end

function _get_kwarg(kwarg,kwargs,value)
  try
    return kwargs[kwarg]
  catch
    return value
  end
end

function _is_reffe_lagrangian_compatible_with_polytopes(reffe,polytopes)
  a = true
  for p in polytopes
    a = a && _is_reffe_lagrangian_compatible_with_polytope(reffe,p)
  end
  a
end

function _is_reffe_lagrangian_compatible_with_polytope(reffe,polytope)
  v = (reffe == :Lagrangian)
  v = v || (reffe == :QLagrangian && is_n_cube(polytope))
  v = v || (reffe == :SLagrangian && is_n_cube(polytope))
  v = v || (reffe == :PLagrangian && is_simplex(polytope))
  v
end


text:
@testset "Arc-Standard" begin

    tree = deptree(x -> deptoken(x[1], x[2]), [
        ("the", 2),
        ("cat", 3),
        ("slept", 0),
        ("on", 3),
        ("the", 6),
        ("bed", 4)
    ])

    @testset "Static" begin
        oracle = Oracle(ArcStandard(), static_oracle)
        for (cfg, t) in oracle(tree)
            @test t in possible_transitions(cfg)
        end
        @test initconfig(oracle, tree) == initconfig(oracle, tree)
    end
end


text:
# This file provides functions for reading files into an hdf5a db
using Compat
using GZip
using Libz
using Lumberjack
using Docile
using HDF5,JLD
include( Pkg.dir("GenomicStore","src","crud.jl") )
include( Pkg.dir("GenomicStore","src","util.jl") )

function read_gzip_file(filename)
   gzio=gzopen(filename,"r")
   line_num =0
   lines=ASCIIString[]
   for line in eachline(gzio)
        line=chomp(line)
        push!(lines,line)
        line_num +=1
        if line_num % 1000000 == 0
            Lumberjack.info("Procesed $line_num")
        end
    end
    return lines
end

function memory_read_file(filename)
    io = open(filename)
    Lumberjack.info("reading all")
    file=readall(io)
    Lumberjack.info("finished reading all")
    Lumberjack.info("split line")
    lines=split(file,'\n')
    Lumberjack.info("done split")
    close(io)
    return lines
end

function isgzip(filepath)
    (path,ext) = splitext(filepath)
    gzip = (ext==".gz" || ext==".gzip" ) ? true : false
    return gzip
end

"""
  read_and_parse_point_file
  This function expects to have the first three columns to contain
  seq id,start position, and strand.
  e.g.
  chr1    3000826 +
  chr1    3000827 -

  Returns a tuple of vectors
  1. sequence_ids
  2. start_positions
"""
function read_and_parse_point_file(filename;gzip=false,strand_filter_str="-")
    if gzip
        # reset on end must be done for bgzip files
        stream=ZlibInflateInputStream(open(filename),reset_on_end=true)
    else
        stream= open(filename)
    end
    seq_ids=UTF8String[]
    starts=Int64[]
    count=0
    pass_filter_count=-0
    for line in eachline(stream)
        (seq_id,start,strand)=split(line,'\t')
        strand=chomp(strand)
        count+=1
        if count % 1000000 == 0
           Lumberjack.info("Lines read: $count")
        end
        if strand == strand_filter_str
            continue
        end
        pass_filter_count +=1
        push!(starts,parse(Int64,start) )
        push!(seq_ids,seq_id)
    end
    Lumberjack.info("Finished parsing and assigning to array.\n
    Number of lines in file: $count\n
    Number of lines passing filter: $pass_filter_count\n
    First seq_id: $(seq_ids[1])\n
    Last seq_id: $(seq_ids[end])")
    return (seq_ids,starts)
end

"""
  read_and_parse_methpipe_cg_bed_file

 Reads into memory a file and parses a file in the format:

  chr1    3000826 +       CG      1       17
  chr1    3001006 +       CG      0.8     10

 extracting out either the 5th or 6th column (the coverage column)

 returns tuple of 4 arrays

 1. seq_ids
 2, starts
 3. stops
 4. levels/coverage

"""
function read_and_parse_methpipe_cg_bed_file(filename;gzip=false,value_column=5)
    if gzip
        # reset on end must be done for bgzip files
        stream=ZlibInflateInputStream(open(filename),reset_on_end=true)
    else
        stream= open(filename)
    end
    seq_ids=UTF8String[]
    starts=Int64[]
    stops=Int64[]
    if  value_column == 5
        scores=Float32[]
        T=Float32
    elseif value_column == 6
        scores=Int32[]
        T=Int32
    else
        Lumberjack.error("invalid column number (valid: 5 or 6)")
    end
    count=0
    for line in eachline(stream)
        # This splits the line into the following:
        #(seq_id,start,strand,context,score,coverage) = split(line,'\t')
        split_line=split(line,'\t')
        count+=1
        if count % 1000000 == 0
           Lumberjack.info("Lines read: $count")
        end
        push!(seq_ids,split_line[1])
        start_int=parse(Int64,split_line[2])
        push!(starts,start_int)
        push!(stops,start_int+1)
        push!(scores,parse(T,split_line[value_column]))
    end
    Lumberjack.info("Finished parsing and assigning to array.\n
    Number of lines in file: $count\n
    First seq_id: $(seq_ids[1])\n
    Last seq_id: $(seq_ids[end])")
    return (seq_ids,starts,stops,scores)
end

"""
memory_read_and_parse_interval_file

 This function expects a file to have the first three columns to contain
 a sequence id, start position, and stop/end position

 Initially used to parse files like below:

 chr1    3670331 3672761 HYPO0   184     +
 chr1    3993551 3993737 HYPO1   5       +

 Returns a tuple of vectors

 1. sequence_ids
 2. start_positions
 3. stop_positions
"""

function  read_and_parse_interval_file(filename; gzip=false, only_start=false,strand_filter_char=nothing)
    if gzip
        # reset on end must be done for bgzip files
        stream=ZlibInflateInputStream(open(filename),reset_on_end=true)
    else
        stream= open(filename)
    end
    seq_ids=UTF8String[]
    starts=Int64[]
    stops=Int64[]
    count=0
    for line in eachline(stream)
        # This splits the line into the following:
        #(seq_id,start,stop,id,cpg_count,strand) = split(line,'\t')
        split_line=split(line,'\t')
        count+=1
        if count % 1000000 == 0
           Lumberjack.info("Lines read: $count")
        end
        push!(seq_ids,split_line[1])
        push!(starts,parse(Int64,split_line[2]))
        push!(stops, parse(Int64,split_line[3]))
    end
    Lumberjack.info("Finished parsing and assigning to array.\n
    Number of lines in file: $count\n
    First seq_id: $(seq_ids[1])\n
    Last seq_id: $(seq_ids[end])")
    if only_start
        return (seq_ids,starts)
    end
    return  (seq_ids,starts,stops)
end

"""
  read_and_parse_bedgraph

  This function expects to have the first three columns to contain
  seq id,start position, and stop/end postion.

  e.g.

  chr1 3670331 3670333 0.05

  Returns a tuple of vectors

  1. sequence_ids
  2. start_positions
  3. stop_positions
  4. score
"""
function read_and_parse_bedgraph(filename;gzip=false)
    if gzip
        # reset on end must be done for bgzip files
        stream=ZlibInflateInputStream(open(filename),reset_on_end=true)
    else
        stream= open(filename)
    end
    seq_ids=UTF8String[]
    starts=Int64[]
    stops=Int64[]
    scores=Float32[]
    count=0
    for line in eachline(stream)
           (seq_id,start,stop,score)=split(line,'\t')
           push!(starts,parse(Int64,start) )
           push!(scores,parse(Float32,score))
           push!(stops,parse(Int64,stop))
           push!(seq_ids,seq_id)
           count+=1
           if count % 1000000 == 0
             Lumberjack.info("Read $count")
           end
       end
    Lumberjack.info("finished parsing and assigning to array")
    return (seq_ids,starts,stops,scores)
end

function track_exists(genomic_store_path, track_id)
   db=getdb(genomic_store_path)
   return _istrack(db,track_id)
end

#=

save_bed_track

 This function saves "bed" like files to hdf5 storage. The name should more accurately
 be "save_bedlike_track"

 It handles saving values specified in

 a) bedgraph format
 b) methpipe_cpg_bed
     i) extracting the methylation levels value
     ii) extracting the coverage levels

 Current requirement that input file is pre-sorted
=#

function save_bed_track(genomic_store_path,input_file,track_id,chr_sizes_path;
                             start_coord_shift=0,
                             OUT_OF_RANGE_VAL=0.0,
                             gzip=false,
                             val_type="float32",
                             bedtype="bedgraph")
    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)
    #--- sort in temporary dir
    #--- load
    if bedtype == "bedgraph_cpg"
        (seq_ids,starts,stops,scores)=read_and_parse_bedgraph(input_file,gzip=gzip)
    elseif bedtype == "methpipe_bed_levels"
        (seq_ids,starts,stops,scores)=read_and_parse_methpipe_cg_bed_file(input_file,gzip=gzip,value_column=5)
    elseif bedtype == "methpipe_bed_coverage"
        (seq_ids,starts,stops,scores)=read_and_parse_methpipe_cg_bed_file(input_file,gzip=gzip,value_column=6)
    else
         Lumberjack.error("Invalid bedtype specified $bedtype")
    end
    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")

    db=getdb(genomic_store_path) # from crud.jl

    # initialise for the first sequence_id
    seq_id=seq_ids[1]
    seq_len  = chr_sizes_dict[seq_id]
    seq_vals = fill(OUT_OF_RANGE_VAL,seq_len)

    # initialise strucutres for checking sorted data
    seen_seq_ids = Set()
    last_start_pos = 0
    println("sequence id: $seq_id")
    for i=1:bedgraph_file_length
        # if there is a change in sequence_id
        if seq_id != seq_ids[i]
            println("change: $seq_id, $(seq_ids[i])")
            if( in(seq_id,seen_seq_ids ) )
                Lumberjack.error("Unordered sequence identifier found at $i ($seq_id)")
            end

            write_track(db,track_id,seq_id,seq_vals)

            # -- indicate we have processed the current seq_id
            push!(seen_seq_ids,seq_id)

            # -- reset last start position
            last_start_pos=0

            # -- update seq_id
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --
            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(OUT_OF_RANGE_VAL,seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end

        #--- check that the start positions are sorted
        if last_start_pos > starts[i]
            Lumberjack.error("Unsorted position found at $i ($last_start_pos > $(starts[i])")
        end
        # update last_start_pos
        last_start_pos = starts[i]

        #--- check that the stop position is not less than the start position
        if starts[i] > stops[i]
            Lumberjack.error("Start position greater than stop position line $i $(starts[i]) > $(stops[i])")
        end

        if val_type == "float32"
            try
                seq_vals[ starts[i] + start_coord_shift ] = Float32( scores[i] )
            catch y
                println("processing this line:\n")
                curr_row = join([ seq_ids[i], starts[i],stops[i], scores[i] ],"\t" )
                println(curr_row)
                println("shift: $start_coord_shift")
                Lumberjack.error(y)
            end
        end
        if val_type == "int32"
            seq_vals[ starts[i] + start_coord_shift ] = Int32( scores[i] )
        end
        if val_type == "int16"
            seq_vals[ starts[i] + start_coord_shift ] = Int16( scores[i] )
        end
        if val_type == "int8"
            seq_vals[ starts[i] + start_coord_shift ] = Int8( scores[i] )
        end
    end
    # write the final track
    write_track(db,track_id,seq_id,seq_vals)
    Lumberjack.info("Written track $seq_id")
end

#=
 save_point_start_track

 This function saves a "bed" like file containing just a start position and strand to hdf5 storage.
 The tab delimited format describing CpG positions looks like this:
 chr1    3000826 +
 chr1    3000827 -

 Current requirement that input file is pre-sorted

=#

function save_start_point_track(genomic_store_path,point_file,track_id,chr_sizes_path;start_coord_shift=0,strand_filter_str=nothing)
    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)

    (seq_ids,starts)=read_and_parse_point_file(point_file, strand_filter_str=strand_filter_str)
    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")

    db=getdb(genomic_store_path)

    # initialise for the first sequence_id
    seq_id=seq_ids[1]
    seq_len  = chr_sizes_dict[seq_id]
    seq_vals = fill(Int8(0),seq_len)

    # initialise structures for checking sorted data
    seen_seq_ids = Set()
    last_start_pos = 0

    for i=1:bedgraph_file_length
        if seq_id != seq_ids[i]

            if( in(seq_id,seen_seq_ids ) )
                Lumberjack.error("Unordered sequence identifier found at $i ($seq_id)")
            end

            write_track(db,track_id,seq_id,seq_vals)

            # -- indicate we have processed the current seq_id
            push!(seen_seq_ids,seq_id)

            # -- reset last start position
            last_start_pos=0

            # -- update seq_id
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --
            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            # -- check we have sorted data by chr
            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(Int8(0),seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end

        #--- check that the start positions are sorted
        if last_start_pos > starts[i]
            Lumberjack.error("Unsorted position found at $i ($last_start_pos > $(starts[i])")
        end
        # update last_start_pos
        last_start_pos = starts[i]

        seq_vals[ starts[i] + start_coord_shift ] = 1
    end
    # write the final track
    write_track(db,track_id,seq_id,seq_vals)
end

#=

 save_interval_track

# this assumes each interval needs to be recorded
# as existing or not existing - use 1 or 0
# this function just records the start position at the moment.
# interval_label_type may be either "binary" uses 1 or 0 to indicate an interval covers a basepair
# or multiclass_int8 which allows an 256 possible integer values (-128 .. 127)
# this last option is not supported anymore.
 Current requirement that input file is pre-sorted
=#
function save_interval_track(genomic_store_path,interval_file,track_id,chr_sizes_path;
                             start_coord_shift=0,strand_filter_char=nothing, interval_label_type="binary")

    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)

    if interval_label_type == "binary"
        (seq_ids,starts,stops)=read_and_parse_interval_file(interval_file, strand_filter_char=strand_filter_char)
    elseif interval_label_type == "multiclass_int8"
        #(seq_ids,starts,stops,scores)=read_and_integer_parse_bedgraph(interval_file)
       Lumberjack.error("multiclass_int8 file type not supported")
    else
        Lumberjack.error("Interval file type not supported")
    end

    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")

    db=getdb(genomic_store_path)

    # initialise for the first sequence_id
    seq_id=seq_ids[1]
    seq_len  = chr_sizes_dict[seq_id]
    seq_vals = fill(Int8(0),seq_len)

    # initialise structures for checking sorted data
    seen_seq_ids = Set()
    last_start_pos = 0

    for i=1:bedgraph_file_length
        if seq_id != seq_ids[i]

            if( in(seq_id,seen_seq_ids ) )
                Lumberjack.error("Unordered sequence identifier found at $i ($seq_id)")
            end

            write_track(db,track_id,seq_id,seq_vals)

            # -- indicate we have processed the current seq_id
            push!(seen_seq_ids,seq_id)

            # -- reset last start position
            last_start_pos=0

            # -- update seq_id
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --
            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            # -- check we have sorted data by chr
            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(Int8(0),seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end

        #--- check that the start positions are sorted
        if last_start_pos > starts[i]
            Lumberjack.error("Unsorted position found at $i ($last_start_pos > $(starts[i])")
        end
        # update last_start_pos
        last_start_pos = starts[i]

        #--- check that the stop position is not less than the start position
        if starts[i] > stops[i]
            Lumberjack.error("Start position greater than stop position line $i $(starts[i]) > $(stops[i])")
        end

        starting_pos=starts[i] + start_coord_shift
        stopping_pos=stops[i]  + start_coord_shift
        if interval_label_type == "binary"
            score=1
        elseif interval_label_type == "multiclass_int8"
            score = scores[i]
        end
        for pos=starting_pos:stopping_pos
            seq_vals[ pos ] = score
        end
    end
    # write the final track
    write_track(db,track_id,seq_id,seq_vals)
end

# -- Exposed methods

#=
 store_methpipe_bed_points
   exported method
=#

function store_methpipe_bed_points(filepath::AbstractString,
                                   genomic_store_path::AbstractString,
                                   track_id::AbstractString,
                                   chr_sizes_path::AbstractString;
                                   na_val=0,
                                   coord_shift=0,
                                   measurement="levels")

    gzip = isgzip(filepath) ? true : false

    if measurement == "levels"
        val_type = "float32"
        bedtype = "methpipe_bed_levels"
        na_val=Float32(na_val)
    elseif measurement == "coverage"
        val_type = "int32"
        bedtype = "methpipe_bed_coverage"
        na_val=Int32(na_val)
    end

    save_bed_track(genomic_store_path,filepath,track_id,chr_sizes_path,
                                start_coord_shift=coord_shift, OUT_OF_RANGE_VAL=na_val, gzip=gzip,
                                val_type=val_type,
                                bedtype=bedtype
                           )
end

#=
 store_cpg_points
  default is to filter out lines on negative strand
=#
function store_cpg_points(filepath::AbstractString,
                                   genomic_store_path::AbstractString,
                                   track_id::AbstractString,
                                   chr_sizes_path::AbstractString;
                                   strand_filter_str="-",
                                   coord_shift=0)

  save_start_point_track(genomic_store_path, filepath,track_id, chr_sizes_path,
                         start_coord_shift=coord_shift,
                         strand_filter_str=strand_filter_str
                         )

end

#=
 store_methpipe_bed_intervals
=#
function store_methpipe_bed_intervals(filepath::AbstractString,
                                   genomic_store_path::AbstractString,
                                   track_id::AbstractString,
                                   chr_sizes_path::AbstractString;
                                   coord_shift=0)

   save_interval_track(genomic_store_path,filepath,track_id,chr_sizes_path,
                       start_coord_shift=0, interval_label_type="binary")

end

#=
 store_bedgraph
=#

function store_bedgraph_cpg( filepath::AbstractString,
                         genomic_store_path::AbstractString,
                         track_id::AbstractString,
                         chr_sizes_path::AbstractString;
                         na_val=0,
                         coord_shift=0,
                         measurement="levels")

    gzip = isgzip(filepath) ? true : false
    val_type = "float32"
    na_val=Float32(na_val)
    save_bed_track(genomic_store_path,filepath,track_id,chr_sizes_path,
                   start_coord_shift=coord_shift, OUT_OF_RANGE_VAL=na_val,
                   gzip=gzip, val_type=val_type, bedtype="bedgraph_cpg" )

end


text:
insts = OrderedDict()

insts["test"] = [
    ((3, 0.7, 0.5, true),),
    ((3, 0.7, 0.5, false),),
    ((6, 0.3, 0.7, true),),
    ((6, 0.3, 0.7, false),),
    ((10, 0.7, 0.2, true),),
    ((10, 0.7, 0.2, false),),
    ((15, 0.2, 0.3, true),),
    ((15, 0.2, 0.3, false),),
]

function completablepsd_insts(use_nat::Bool, d_max::Int)
    return [((d, inv(sqrt(d)), 0.7, use_nat),) for d in vcat(10, 10:5:d_max)]
end

insts["nat"] = completablepsd_insts(true, 65)
insts["ext"] = completablepsd_insts(false, 55)

return (CompletablePSD, insts)


text:
using GeoPhyInv
using Test

model = Medium(:acou_homo2);
update!(model, [:vp,:rho], randn_perc=1)

model0 = Medium(:acou_homo2);
update!(model0, [:vp,:rho], randn_perc=1)

ageom=AGeom(model.mgrid, :surf)
wav, tgrid=ricker(model, 3, 1.0)
srcwav=SrcWav(tgrid,ageom,[:P])
update!(srcwav,[:P], wav)

parameterization=[:χvp, :χrho, :null]

mgrid=model.mgrid

@testset "test parallel implementation during gradient" begin
	for attrib_mod in [Fdtd(), FdtdBorn()]
		global pa=SeisInvExpt(attrib_mod, Migr(), srcwav=srcwav, ageom=ageom, tgrid=tgrid, modm=model0,
				     modm_obs=model,
				     modm0=model0,
				     igrid_interp_scheme=:B2,
				     igrid=broadcast(x->range(x[1],stop=x[end],step=300.),mgrid),
				     parameterization=parameterization,   verbose=false,
				     nworker=1)

		global pa_parallel=SeisInvExpt(attrib_mod, Migr(), srcwav=srcwav, ageom=ageom, tgrid=tgrid, modm=model0,
				     modm_obs=model,
				     modm0=model0,
				     igrid_interp_scheme=:B2,
				     igrid=broadcast(x->range(x[1],stop=x[end],step=300.),mgrid),
				     parameterization=parameterization,   verbose=false,
				     nworker=nothing)

		result=update!(pa)

		result_parallel=update!(pa_parallel)

		@test result[2] ≈ result_parallel[2]
	end
end

@testset "Testing Born Modeling and its gradient" begin
	global expt=x->SeisInvExpt(FdtdBorn(), x,
			     srcwav=srcwav, ageom=ageom, tgrid=tgrid, modm=model0,
	     		     modm0=model0,
			     modm_obs=model,
			     igrid_interp_scheme=:B2,
			     igrid=broadcast(x->range(x[1],stop=x[end],step=350.),mgrid),
			     parameterization=parameterization,   verbose=false)

	global pa=expt(LS())

	update!(pa, bounded_flag=true, solver=:ipopt,
			ipopt_options=[["max_iter", 0],["derivative_test", "first-order"]])
	pa=expt(Migr())

	result=update!(pa)

	pa_fd=expt(Migr_FD());
	result_fd=update!(pa_fd)

	f=Misfits.error_squared_euclidean!(nothing, result[2], result_fd[2], nothing, norm_flag=true)

	@test f<1e-15
end

@testset "Testing gradient LS FWI" begin
	expt=x->SeisInvExpt(Fdtd(), x,
			     srcwav=srcwav, ageom=ageom, tgrid=tgrid, modm=model0,
	     		     modm0=model0,
			     modm_obs=model,
			     igrid_interp_scheme=:B2,
			     igrid=broadcast(x->range(x[1],stop=x[end],step=350.),mgrid),
			     parameterization=parameterization,   verbose=false)
	pa=expt(LS())

	update!(pa, bounded_flag=true, solver=:ipopt,
			ipopt_options=[["max_iter", 0],["derivative_test", "first-order"]])
	pa=expt(Migr())

	result=update!(pa)

	pa_fd=expt(Migr_FD());
	result_fd=update!(pa_fd)

	f=Misfits.error_squared_euclidean!(nothing, result[2], result_fd[2], nothing, norm_flag=true)

	@test f<1e-15
end


text:
# Alex Hoffman
module astrofunk

# include the file and then export the functions
include("circularHohmann.jl")
export circHohmann

# include("dirCosMat.jl")
# export DCM
#
# include("jacobiConstant.jl")
# export JC
# export xtest

include("basicEquations.jl")
export JC
export DCM
export CR3BP_EOM
export U_star

end #astrofunk module


text:
@testset "Linear programs ($T)" for T in [Float32, Float64]
    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra

    # Solving LP with AFBA
    #
    #   minimize    c'x             -> f = <c,.>
    #   subject to  Ax = b          -> h = ind{b}, L = A
    #               x >= 0          -> g = ind{>=0}
    #
    # Dual LP
    #
    #   maximize    b'y
    #   subject to  A'y <= c
    #
    # Optimality conditions
    #
    #   x >= 0              [primal feasibility 1]
    #   Ax = b              [primal feasibility 2]
    #   A'y <= c            [dual feasibility]
    #   X'(c - A'y) = 0     [complementarity slackness]
    #

    function assert_lp_solution(c, A, b, x, y, tol)
        # Check and print solution quality measures (for some reason the
        # returned dual iterate is the negative of the dual LP variable y)

        nonneg = -minimum(min.(0.0, x))
        @test nonneg <= tol

        primal_feasibility = norm(A * x - b)
        @test primal_feasibility <= tol

        dual_feasibility = maximum(max.(0.0, -A' * y - c))
        @test dual_feasibility <= tol

        complementarity = abs(dot(c + A' * y, x))
        @test complementarity <= tol
    end

    n = 100 # primal dimension
    m = 80 # dual dimension (i.e. number of linear equalities)
    k = 50 # number of active dual constraints (must be 0 <= k <= n)
    x_star = vcat(rand(T, k), zeros(T, n - k)) # primal optimal point
    s_star = vcat(zeros(T, k), rand(T, n - k)) # dual optimal slack variable
    y_star = randn(T, m) # dual optimal point

    A = randn(T, m, n)
    b = A * x_star
    c = A' * y_star + s_star

    tol = 100 * eps(T)
    maxit = 10_000

    @testset "AFBA" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.AFBA(tol = tol, maxit = maxit)
        x, y, it = solver(x0, y0, f = f, g = g, h = h, L = A)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1000 * tol)

        @test x0 == x0_backup
        @test y0 == y0_backup

    end

    @testset "VuCondat" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.VuCondat(tol = tol, maxit = maxit)
        x, y, it = solver(x0, y0, f = f, g = g, h = h, L = A)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1000 * tol)

        @test x0 == x0_backup
        @test y0 == y0_backup

    end

    # TODO: add Chambolle-Pock (using separable sum of IndPoint and IndNonnegative)

    @testset "DavisYin" begin

        f = IndAffine(A, b)
        g = IndNonnegative()
        h = Linear(c)

        x0 = zeros(T, n)
        x0_backup = copy(x0)

        solver = ProximalAlgorithms.DavisYin(gamma = T(1), tol = tol, maxit = maxit)
        xf, xg, it = solver(x0, f = f, g = g, h = h)

        @test eltype(xf) == T
        @test eltype(xg) == T

        @test it <= maxit

        @assert norm(xf - x_star) <= 1e2 * tol
        @assert norm(xg - x_star) <= 1e2 * tol

        @test x0 == x0_backup

    end

end


text:
#VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Sparrow

using DataFrames, Colors, Cairo
include(Pkg.dir()*"/Processing/src/Processing2D.jl")
const p2d = Processing2D

export DataGroup, PlotFrame, Legend, PageSize
export LinePlot, ScatterPlot, XYZMap, HistPlot
export plot, lineplot, scatterplot, xyzmap, histplot

global dgCount = 0

abstract Graph

include("colors.jl")
include("symbols.jl")
include("plotFrame.jl")
include("datagroup.jl")
include("auxiliary.jl")
include("lineplot.jl")
include("scatterplot.jl")

function plot{T<:Graph}(p::PlotFrame, gs::Array{T})
	findGraphLimits(p, gs)
	drawAxisLabels(p)
	drawTicks(p)
	for g in gs
		plotData(p, g, p.showLegend)
	end

	if p.showLegend
		drawLegend(p, gs)
	end

	p2d.animate(p.wi)
end

function plot(p::PlotFrame, g::Graph)
	findGraphLimits(p, g)
	drawAxisLabels(p)
	drawTicks(p)
	plotData(p, g, p.showLegend)

	if p.showLegend
		drawLegend(p, g)
	end

	p2d.animate(p.wi)
end

function print{T<:Graph}(p::PlotFrame, gs::Array{T}, fn::AbstractString)
	p2d.PDFContext(p.wi, fn)
	p2d.createFont(p.wi, "Arial", 13*(p.wpx/500))
	findGraphLimits(p, gs)
	drawAxisLabels(p)
	drawTicks(p)
	for g in gs
		plotData(p, g, p.showLegend)
	end

	if p.showLegend
		drawLegend(p, gs)
	end

	Cairo.show_page(p2d.contexts[p.wi])
	Cairo.finish(p2d.surfaces[p.wi])
	p2d.popContext(p.wi)
end

function print(p::PlotFrame, g::Graph, fn::AbstractString)
	p2d.PDFContext(p.wi, fn)
	p2d.createFont(p.wi, "Arial", 13*(p.wpx/500))
	findGraphLimits(p, g)
	drawAxisLabels(p)
	drawTicks(p)
	plotData(p, g, p.showLegend)

	if p.showLegend
		drawLegend(p, g)
	end

	Cairo.show_page(p2d.contexts[p.wi])
	Cairo.finish(p2d.surfaces[p.wi])
	p2d.popContext(p.wi)
end

function close(p::PlotFrame)
	p2d.endDrawing(p.wi)
end

end


text:
using Iris
using Test

@testset "BoundaryLayers" begin
    pml1 = PML{1}(.1)
    @test_nowarn pml1.start = 0
    @test_nowarn pml1.stop = pml1.start - pml1.depth
    pml2 = PML{2}(.1)
    @test_nowarn pml2.start = 0
    @test_nowarn pml2.stop = pml2.start + pml2.depth
    @test pml1(-.5) ≈ pml2(.5)
    @test Iris.Common.getside(pml1) == 1
    @test Iris.Common.getside(pml2) == 2

    cpml1 = cPML{1}(.1)
    @test_nowarn cpml1.start = 0
    @test_nowarn cpml1.stop = cpml1.start - cpml1.depth
    cpml2 = cPML{2}(.1)
    @test_nowarn cpml2.start = 0
    @test_nowarn cpml2.stop = cpml2.start + cpml2.depth
    @test cpml1(-.5) ≈ cpml2(.5)
    @test Iris.Common.getside(cpml1) == 1
    @test Iris.Common.getside(cpml2) == 2

    @test pml1(-.5) ≈ -conj(cpml2(.5))

    nbl1 = noBL{1}()
    @test_nowarn nbl1.start = 0
    @test_nowarn nbl1.stop = nbl1.start - nbl1.depth
    nbl2 = noBL{2}(.1)
    @test_nowarn nbl2.start = 0
    @test_nowarn nbl2.stop = nbl2.start + nbl2.depth
    @test nbl1(-.5) ≈ conj(nbl2(-.5))
    @test Iris.Common.getside(nbl1) == 1
    @test Iris.Common.getside(nbl2) == 2
end


text: # NOTE: this file keeps patches for julia itself to make JET.jl test keep passing


text:
### Coefficient

## FIXME: Can't specify both. Update: may be possible now
#@mkapprule Coefficient nargs => 2:3  options => Dict(:Method => :SymPy)

@mkapprule Coefficient options => Dict(:Method => :SymPy)

## TODO Merge this with CoefficientSympy
# @doap function Coefficient(args...; kws...)
#     meth = get(kws,:Method)
#     if meth == :Sympy

@doap function Coefficient(expr, x; kws...)
    res = coefficient0(mx,expr,x)
    res == :nothing ? mx : res
end

@doap function Coefficient(expr, x, pow; kws...)
    pow == 0 && return coefficient_free(mx,expr,x)
    res = coefficient0(mx,expr, mxpr(:Power, x,pow))
    res == :nothing ? mx : res
end

coefficient0(mx,expr,x::Number) = mx
coefficient0(mx,expr,x) = coefficient(expr,x)
coefficient(expr::TimesT, x) = coefficient_times(expr,[x])

function coefficient(expr::TimesT, x::TimesT)
    if isa(x[1],Number)
        fac = x[1]
        res = coefficient_times(expr,margs(x)[2:end])
        res == 0 && return 0
        return mmul(res, mpow(fac,-1))
    else
        coefficient_times(expr,x)
    end
end

function coefficient_times(expr::TimesT, coll)
    expr == coll && return 1
    posns = Symata.oldfindin(expr, coll)
    isempty(posns) && return 0
    (length(posns) != length(coll)) && return 0
    nargs = deleteat!(copy(margs(expr)),posns)
    length(nargs) == 1 && return nargs[1]
    mxpra(mhead(expr), nargs)
end

coefficient(expr,x) = x == expr ? 1 : 0
coefficient(expr::PlusT,x::PlusT) = x == expr ? 1 : 0

function coefficient(expr::PlusT, x)
    nargs = newargs()
    for term in expr
        res = coefficient(term,x)
        if res != 0
            push!(nargs,res)
        end
    end
    isempty(nargs) && return 0
    MPlusA(nargs)
end

coefficient_free(mx,expr,x::Number) = mx
coefficient_free(mx,expr,x) = coefficient_free1(expr,x)

## Very, very strange bug. The use of macro @extomx sometimes causes segfault during compiling, sometimes not,
## even if all other code does not change.
## Eg.
## 1. Check out a particluar commit and compile and run all tests.
## 2. Edit some files
## 3. Check out same commit. Compiling causes segfault.
## Just adding these comments caused the compliation to succeed.
## Adding the previous line "Just..." now causes compilation to fail
function coefficient_free1(expr,x)
    return (freeq(expr,x) && freeq(expr, MPower(x,MBlank())))  ? expr : 0  # this works fine, as well
#    return (freeq(expr,x) && freeq(expr, @extomx( $x^_ ))) ? expr : 0   # This is too tempermental
end

function coefficient_free1(expr::PlusT,x)
    nargs = newargs()
    for term in expr
        res = coefficient_free1(term,x)
        if res != 0
            push!(nargs,res)
        end
    end
    isempty(nargs) && return 0
    MPlusA(nargs)
end


text:
using BinaryBuilder

name = "WCS"
version = v"7.7.0"

# Collection of sources required to build WCS
sources = [
    ArchiveSource("https://cache.julialang.org/ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib-$(version.major).$(version.minor).tar.bz2",
                  "a708e6b4e90e2cd6427434715b591bb9c3d4005c99725ec494b8e0bed2de5355"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/wcslib-*/
if [[ "${target}" == *mingw* ]]; then
    export CFLAGS="${CFLAGS} -DNO_OLDNAMES"
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-fortran --without-cfitsio --without-pgplot --disable-utils
make -j${nproc}
make install

# Remove static library
rm "${prefix}/lib/"libwcs*.a

# Install all license files
install_license COPYING*
"""

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

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

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
abstract type AbstractAcquiredImage{T,N} <: AbstractPooledDenseArray{T,N}
end

"""
    id(image::AcquiredImage)

Return image ID.
"""
id(image::AbstractAcquiredImage) = error("No implementation for $(typeof(image))")

"""
    timestamp(image::AcquiredImage)

Return image timestamp.
"""
timestamp(image::AbstractAcquiredImage) = error("No implementation for $(typeof(image))")

mutable struct AcquiredImage{T,N} <: AbstractAcquiredImage{T,N}
    # Inherits behaviour of AbstractPooledDenseArray, by having the same fields
    array::Array{T,N}
    ref_count::Int
    dispose::Function

    id::Int
    timestamp::Int
    function AcquiredImage(a::Array{T,N}, id, timestamp) where {T,N}
        function dispose(img)
            @debug "Disposing $img"
        end
        new{T,N}(a, 1, dispose, id, timestamp)
    end
end

id(img::AcquiredImage) = img.id
timestamp(img::AcquiredImage) = img.timestamp

function Base.show(io::IO, image::AcquiredImage)
    write(io, "$(nameof(AcquiredImage))(id:$(image.id), timestamp:$(image.timestamp), ref_count:$(image.ref_count), array:$(image.array))")
end


text:
using Documenter

include("make.jl")

deploydocs(repo="github.com/JuliaTime/TimeZones.jl.git")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct BaseMachine <: Machine
        R::Float64
        Xd_p::Float64
        eq_p::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of a Classic Machine: GENCLS in PSSE and PSLF

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: `(0, nothing)`
- `Xd_p::Float64`: Reactance after EMF in machine per unit, validation range: `(0, nothing)`
- `eq_p::Float64`: Fixed EMF behind the impedance, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: BaseMachine has no states
- `n_states::Int`: BaseMachine has no states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct BaseMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in machine per unit"
    Xd_p::Float64
    "Fixed EMF behind the impedance"
    eq_p::Float64
    ext::Dict{String, Any}
    "BaseMachine has no states"
    states::Vector{Symbol}
    "BaseMachine has no states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function BaseMachine(R, Xd_p, eq_p, ext=Dict{String, Any}(), )
    BaseMachine(R, Xd_p, eq_p, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function BaseMachine(; R, Xd_p, eq_p, ext=Dict{String, Any}(), states=Vector{Symbol}(), n_states=0, internal=InfrastructureSystemsInternal(), )
    BaseMachine(R, Xd_p, eq_p, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function BaseMachine(::Nothing)
    BaseMachine(;
        R=0,
        Xd_p=0,
        eq_p=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`BaseMachine`](@ref) `R`."""
get_R(value::BaseMachine) = value.R
"""Get [`BaseMachine`](@ref) `Xd_p`."""
get_Xd_p(value::BaseMachine) = value.Xd_p
"""Get [`BaseMachine`](@ref) `eq_p`."""
get_eq_p(value::BaseMachine) = value.eq_p
"""Get [`BaseMachine`](@ref) `ext`."""
get_ext(value::BaseMachine) = value.ext
"""Get [`BaseMachine`](@ref) `states`."""
get_states(value::BaseMachine) = value.states
"""Get [`BaseMachine`](@ref) `n_states`."""
get_n_states(value::BaseMachine) = value.n_states
"""Get [`BaseMachine`](@ref) `internal`."""
get_internal(value::BaseMachine) = value.internal

"""Set [`BaseMachine`](@ref) `R`."""
set_R!(value::BaseMachine, val) = value.R = val
"""Set [`BaseMachine`](@ref) `Xd_p`."""
set_Xd_p!(value::BaseMachine, val) = value.Xd_p = val
"""Set [`BaseMachine`](@ref) `eq_p`."""
set_eq_p!(value::BaseMachine, val) = value.eq_p = val
"""Set [`BaseMachine`](@ref) `ext`."""
set_ext!(value::BaseMachine, val) = value.ext = val


text:
module SyslogLogging

using Syslogs
using Logging
using Sockets

import Logging: shouldlog, min_enabled_level, catch_exceptions, handle_message
export SyslogLogger

const last_ident = String[""]
function open_syslog(ident::String, facility::Symbol)
    Syslogs.openlog(ident, 0, Syslogs.FACILITIES[facility])
    global last_ident
    last_ident[1] = ident
    nothing
end

function set_syslog_identity(ident::String, facility::Symbol)
    global last_ident
    (last_ident[1] === ident) || open_syslog(ident, facility)
end

"""
Logs messages to a syslog facility.
"""
struct SyslogLogger <: AbstractLogger
    ident::String
    facility::Symbol
    syslog::Syslog
    min_level::LogLevel
    message_limits::Dict{Any,Int}
    lck::Union{Nothing,ReentrantLock}

    function SyslogLogger(ident::String, level=Logging.Info; facility::Symbol=:user, host::Union{Nothing,IPAddr,AbstractString}=nothing, port::Union{Nothing,Int}=nothing, tcp::Bool=false, lck::Union{Nothing,ReentrantLock}=nothing)
        syslog = (host === nothing) ? Syslog(facility) :
                 (port === nothing) ? Syslog(host, facility; tcp=tcp) :
                 Syslog((typeof(host) <: AbstractString) ? getaddrinfo(host) : host, port, facility; tcp=tcp)
        new(ident, facility, syslog, level, Dict{Any,Int}(), lck)
    end
end

idlock(fn, logger::SyslogLogger) = (logger.lck === nothing) ? fn() : lock(fn, logger.lck)

shouldlog(logger::SyslogLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::SyslogLogger) = logger.min_level

catch_exceptions(logger::SyslogLogger) = false

function handle_message(logger::SyslogLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end
    buf = IOBuffer()
    iob = IOContext(buf, logger.syslog)
    levelstr = level == Logging.Warn ? "Warning" : string(level)
    print(iob, levelstr, ": ", chomp(string(message)))
    print(iob, " [", something(_module, "nothing"), " ", something(filepath, "nothing"), ":", something(line, "nothing"), "]")
    if !isempty(kwargs)
        for (key, val) in kwargs
            print(iob, ", [", key, "=", val, "]")
        end
    end
    syslog_level = (level === Logging.Debug) ? :debug :
                   (level === Logging.Info) ? :info :
                   (level === Logging.Warn) ? :warn :
                   :error
    idlock(logger) do
        set_syslog_identity(logger.ident, logger.facility)
        println(logger.syslog, syslog_level, String(take!(buf)))
    end
    nothing
end

end # module


text:
# Tests the version of transit_poly which uses a Type
# for holding all of the variables related to the transit
# computation:

# Defines the transit Type:
#include("transit_structure.jl")
# Defines the function that computes the d_n coefficients:
#include("compute_d_n_struct.jl")

# Defines the function that computes the transit:
include("transit_poly_struct.jl")

# Instantiate a transit with r=0.1, b=0.5, u_n=[0.2,0.3]:

# Transform from u_n to d_n coefficients (this is now added to transit_init function):
#compute_d_n_grad!(trans)

# Load in PyPlot (Julia matplotlib wrapper):
using PyPlot

# Define a function that loops over impact parameter
# and defines arrays to hold the results:

function profile_transit_poly(trans)
nb = 1000001; b=zeros(nb)
flux = zeros(nb); fgrad = zeros(4)
flux_grad = zeros(nb,4)
tic()
for i=0:1000000
  b[i+1] = sqrt(((i-500000.)/500000.*(1.+2.*trans.r))^2)
#  flux[i+1]= transit_poly!(r,b[i+1],u,fgrad)
  trans.b = b[i+1]
  flux[i+1]= transit_poly!(trans)
  flux_grad[i+1,1:2] = trans.dfdrb
  flux_grad[i+1,3:4] = trans.dfdu
end
# Now, plot the flux and derivatives:
toc()
#clf()
plot(b,flux-1,label="flux-1")
plot(b,trans.r*flux_grad[:,1],label="r*dfdr")
plot(b,trans.r*flux_grad[:,2],label="r*dfdb")
plot(b,flux_grad[:,3],label="dfdu1")
plot(b,flux_grad[:,4],label="dfdu2")
legend(loc="lower right")
println("b: ",minimum(b)," ",maximum(b))
println("f: ",minimum(flux)," ",maximum(flux))
return
end

trans = transit_init(0.1,0.5,[0.2,0.3],true)

# Call the function:
profile_transit_poly(trans)

savefig("transit_poly_derivatives.pdf", bbox_inches="tight")


text:
function gc_preserve_call_expr(c, K)
  skip = length(c.args) - 2
  q = Expr(:block, Expr(:meta, :inline))
  for k ∈ 1:K
    push!(c.args, :(@inbounds(args[$k])))
  end
  push!(q.args, gc_preserve_call(c, 1))
  push!(q.args, :A)
  q
end

# TODO: replace this with `LoopVectorization.gc_preserve_vmap!` after making that version switch to `stridedpointers`+ size info
@generated function gc_preserve_map!(
  f::F,
  A::AbstractStrideArray,
  args::Vararg{AbstractArray,K},
) where {F,K}
  gc_preserve_call_expr(:(vmap!(f, A)), K)
end

@inline Base.map!(
  f::F,
  A::AbstractStrideArray,
  arg1::AbstractArray,
  args::Vararg{AbstractArray,K},
) where {F,K} = gc_preserve_map!(f, A, arg1, args...)
# these two definitions are to avoid ambiguities
@inline Base.map!(f::F, A::AbstractStrideArray, arg::AbstractArray) where {F} =
  gc_preserve_map!(f, A, arg)
@inline Base.map!(
  f::F,
  A::AbstractStrideArray,
  arg1::AbstractArray,
  arg2::AbstractArray,
) where {F} = gc_preserve_map!(f, A, arg1, arg2)
@inline Base.map(f::F, A::AbstractStrideArray, args::Vararg{Any,K}) where {F,K} =
  gc_preserve_map!(f, A, args...)
@inline Base.reduce(op::O, A::AbstractStrideArray{<:Number}) where {O} =
  @gc_preserve vreduce(op, A)
@inline Base.reduce(::typeof(vcat), A::AbstractStrideArray{<:Number}) = A
@inline Base.reduce(::typeof(hcat), A::AbstractStrideArray{<:Number}) = A
@inline function gc_preserve_mapreduce(
  f::F,
  op::O,
  A::AbstractStrideArray,
  args::Vararg{AbstractArray,K},
) where {F,O,K}
  gc_preserve_call_expr(:(vmapreduce(f, op, A)), K)
end
@inline Base.mapreduce(
  f::F,
  op::O,
  A::AbstractStrideArray,
  args::Vararg{AbstractArray,K},
) where {F,O,K} = gc_preserve_mapreduce(f, op, A, args...)
@inline function Base.mapreduce(f::F, op::O, A::AbstractStrideArray) where {F,O}
  if (LoopVectorization.check_args(A) && LoopVectorization.all_dense(A))
    @gc_preserve vmapreduce(f, op, A)
  else
    return @gc_preserve Base.invoke(
      Base.mapreduce,
      Tuple{F,O,AbstractArray{eltype(A)}},
      f,
      op,
      A,
    )
  end
end

for (op, r) ∈ ((:+, :sum), (:max, :maximum), (:min, :minimum))
  @eval begin
    @inline Base.reduce(::typeof($op), A::AbstractStrideArray{<:Number}; dims = nothing) =
      @gc_preserve vreduce($op, A, dims = dims)
    @inline Base.$r(A::AbstractStrideArray; dims = nothing) =
      @gc_preserve vreduce($op, A, dims = dims)
  end
end

function Base.copyto!(
  B::AbstractStrideArray{<:Any,<:Any,<:Any,N},
  A::AbstractStrideArray{<:Any,<:Any,<:Any,N},
) where {N}
  @avx for I ∈ eachindex(A, B)
    B[I] = A[I]
  end
  B
end

# why not `vmapreduce`?
@inline function Base.maximum(::typeof(abs), A::AbstractStrideArray{S,D,T}) where {S,D,T}
  s = typemin(T)
  @avx for i ∈ eachindex(A)
    s = max(s, abs(A[i]))
  end
  s
end

@inline function Base.vcat(A::AbstractStrideMatrix, B::AbstractStrideMatrix)
  MA, NA = size(A)
  MB, NB = size(B)
  @assert NA == NB
  TC = promote_type(eltype(A), eltype(B))
  C = StrideArray{TC}(undef, (MA + MB, NA))
  # TODO: Actually handle offsets
  @assert offsets(A) === offsets(B)
  @assert offsets(A) === offsets(C)
  @avx for j ∈ axes(A, 2), i ∈ axes(A, 1)
    C[i, j] = A[i, j]
  end
  @avx for j ∈ axes(B, 2), i ∈ axes(B, 1)
    C[i+MA, j] = B[i, j]
  end
  C
end

@inline function make_stride_dynamic(p::StridedPointer{T,N,C,B,R}) where {T,N,C,B,R}
  si = StrideIndex{N,R,C}(map(Int, strides(p)), offsets(p))
  stridedpointer(pointer(p), si, StaticInt{B}())
end
@inline function make_dynamic(A::PtrArray)
  PtrArray(make_stride_dynamic(stridedpointer(A)), Base.size(A), val_dense_dims(A))
end
@inline function make_dynamic(A::StrideArray)
  StrideArray(make_dynamic(PtrArray(A)), A.data)
end

@inline function Base.exp(A::AbstractStrideArray)
  B = copy(A)
  GC.@preserve C = LinearAlgebra.exp!(B)
  B === C || copyto!(B, C)
  return B
end


text:
# Test that parsing works correctly

# Parsing of terms
@test Const(:atom) == @julog atom
@test Const(1) == @julog 1
@test Const((2,3)) == @julog (2,3)
@test Const("foo") == @julog "foo"
@test Var(:Variable) == @julog Variable
@test Compound(:functor, Term[Const(:a), Var(:B)]) == @julog functor(a, B)

# Parsing of clauses and facts
@test Clause(Const(:head), [Const(:body)]) == @julog head <<= body
@test Clause(Const(:head), [Const(:t1), Const(:t2)]) == @julog head <<= t1 & t2
@test Clause(Const(:head), []) == @julog head <<= true
@test Clause(Const(:head), []) == @julog head'

# Parsing of lists of terms or clauses
@test [Clause(Const(:fact), []), Const(:term)] == @julog [fact', term]
@test [Clause(Const(:a), [Const(:b)])] == @julog Clause[a <<= b]
@test Const[Const(:a), Const(:b)] == @julog [a, b]
@test Term[Const(:a), Const(:b)] == @julog Term[a, b]
@test Const[Const(:a), Const(:b)] == @julog Const[a, b]
@test Var[Var(:A), Var(:B)] == @julog Var[A, B]
@test Compound[Compound(:f, [Const(:arg)])] == @julog Compound[f(arg)]

# Interpolation of constants and Julog expressions
x, y = 2, Const(3)
@test Compound(:even, [Const(x)]) == @julog even($x)
@test Compound(:odd, [y]) == @julog odd(:y)
not_even = Compound(:not, [Compound(:even, [Var(:X)])])
odd_if_not_even = Clause(Compound(:odd, [Var(:X)]), [not_even])
@test odd_if_not_even == @julog odd(X) <<= :not_even


text:
__precompile__(true)
module IPNets

import Base: eltype, lastindex,
length, size, minimum, maximum, extrema, isless,
in, contains, issubset, getindex,
show, string, start, next, done

import Sockets: IPAddr, IPv4, IPv6, @ip_str

export
    # types
    IPNet, IPv4Net, IPv6Net, netmask

IPv4broadcast = typemax(UInt32)
IPv6broadcast = typemax(UInt128)

width(::Type{IPv4}) = UInt8(32)
width(::Type{IPv6}) = UInt8(128)

##################################################
# IPNet
##################################################
abstract type IPNet end

IPNet(ipmask::AbstractString) = ':' in ipmask ? IPv6Net(ipmask) : IPv4Net(ipmask)

##################################################
# Network representations
##################################################

"""Returns the size of an IP network (# of hosts) as a tuple.
"""
function size(net::IPNet)
    numbits = width(typeof(net.netaddr)) - net.netmask
    return (big(2)^numbits, )
end

"""Returns the size of an IP network (# of hosts) as a tuple.
"""
length(net::IPNet) = size(net)[1]

"""String representation of an IP network"""
function string(net::IPNet)
    t = typeof(net)
    s = string("$t(\"")
    s = string(s, net.netaddr, "/", net.netmask, "\")")
    return s
end

show(io::IO, net::IPNet) = print(io, string(net))

# IP Networks are ordered first by starting network address
# and then by network mask. That is, smaller IP nets (with higher
# netmask values) are "less" than larger ones. This corresponds
# to secondary reordering by ending address.
isless(a::T, b::T) where T <: IPNet = a.netaddr == b.netaddr ?
        isless(b.netmask, a.netmask) :
        isless(a.netaddr, b.netaddr)

function issubset(a::T, b::T) where T <: IPNet
    astart, aend = extrema(a)
    bstart, bend = extrema(b)
    return (bstart <= astart <= aend <= bend)
end

"""Membership test for an IP address within an IP network"""
function in(ipaddr::IPAddr, net::IPNet)
    typeof(net.netaddr) == typeof(ipaddr) ||
        error("IPAddr is not the same type as IPNet")
    netstart = net.netaddr.host
    numbits = width(typeof(ipaddr)) - net.netmask
    netend = net.netaddr.host + big(2)^numbits - 1
    return netstart <= ipaddr.host <= netend
end

"""Membership test for an IP address within an IP network"""
contains(net::IPNet, ipaddr::IPAddr) = in(ipaddr, net)

function getindex(net::IPNet, i::Integer)
    t = typeof(net.netaddr)
    ip = t(net.netaddr.host + i - 1)
    ip in net || throw(BoundsError())
    return ip
end

minimum(net::IPNet) = net[1]
maximum(net::IPNet) = net[end]
extrema(net::IPNet) = (minimum(net), maximum(net))
getindex(net::IPNet, r::AbstractRange) = [net[i] for i in r]
start(net::IPNet) = net[1]
next(net::IPNet, s::T) where T <: IPAddr = s, T(s.host + 1)
done(net::IPNet, s::T) where T <: IPAddr = s > net[end]

##################################################
# IPv4
##################################################
"""Type representing an IPv4 network"""
struct IPv4Net <: IPNet
    netaddr::IPv4
    netmask::UInt8
    function IPv4Net(na::IPv4, nmi::Integer)
        (0 <= nmi <= width(IPv4)) || error("Invalid netmask")

        nm = UInt8(nmi)
        mask = _mask2bits(IPv4, nm)
        startip = UInt32(na.host & mask)
        new(IPv4(startip),nm)
    end
end

# "1.2.3.0/24"
function IPv4Net(ipmask::AbstractString)
    if something(findfirst(isequal('/'), ipmask), 0) > 0
        addrstr, netmaskstr = split(ipmask,"/")
        netmask = parse(UInt8, netmaskstr)
    else
        addrstr = ipmask
        netmask = width(IPv4)
    end
    netaddr = IPv4(addrstr)
    return IPv4Net(netaddr,netmask)
end

# "1.2.3.0", "255.255.255.0"
function IPv4Net(netaddr::AbstractString, netmask::AbstractString)
    netaddr = IPv4(netaddr).host
    netmask = _contiguousbitcount(IPv4(netmask).host)
    return IPv4Net(netaddr, netmask)
end

# 123872, 24
IPv4Net(ipaddr::Integer, netmask::Integer) = IPv4Net(IPv4(ipaddr), netmask)

# "(x,y)"
IPv4Net(tuple::Tuple{A,M}) where A where M = IPv4Net(tuple[1],tuple[2])

# "1.2.3.0", 24
IPv4Net(netaddr::AbstractString, netmask::Integer) = IPv4Net(IPv4(netaddr), netmask)

"""Returns the netmask as an IPv4 address"""
netmask(n::IPv4Net) = IPv4(IPv4broadcast-2^(32-n.netmask)+1)

eltype(x::IPv4Net) = IPv4
lastindex(net::IPv4Net) = UInt32(length(net))

##################################################
# IPv6
##################################################

"""Type representing an IPv6 network"""
struct IPv6Net <: IPNet
    # we treat the netmask as a potentially noncontiguous bitmask
    # for speed of calculation and consistency, but RFC2373, section
    # 2 provides for contiguous bitmasks only. We validate this
    # in the internal constructor. This wastes ~15 bytes per addr
    # for the benefit of rapid, consistent computation.
    netaddr::IPv6
    netmask::UInt32

    function IPv6Net(na::IPv6, nmi::Integer)
        (0 <= nmi <= width(IPv6)) || error("Invalid netmask")

        nm = UInt8(nmi)
        mask = _mask2bits(IPv6, nm)
        startip = UInt128(na.host & mask)
        return new(IPv6(startip), nm)
    end
end

# "2001::1/64"
function IPv6Net(ipmask::AbstractString)
    if something(findfirst(isequal('/'), ipmask), 0) > 0
        addrstr, netmaskbits = split(ipmask,"/")
        nmi = parse(Int,netmaskbits)
    else
        addrstr = ipmask
        nmi = width(IPv6)
    end
    netaddr = IPv6(addrstr)
    netmask = nmi
    return IPv6Net(netaddr,netmask)
end

# "2001::1", 64
function IPv6Net(netaddr::AbstractString, netmask::Integer)
    netaddr = IPv6(netaddr)
    return IPv6Net(netaddr, netmask)
end

# 123872, 128
IPv6Net(ipaddr::Integer, netmask::Integer) = IPv6Net(IPv6(ipaddr), netmask)

# (123872, 128)
IPv6Net(t::Tuple{A,M}) where A where M = IPv6Net(t[1],t[2])

eltype(x::IPv6Net) = IPv6
lastindex(net::IPv6Net) = UInt128(length(net))

### Helper functions
function _contiguousbitcount(n::Integer,t=UInt32)
    # takes an integer from 0 to 255 and a type, returns the number
    # of contiguous 1 bits in the number assuming it's of that type,
    # starting from the left.
    # cbc(240,UInt8) == 0x04 ("1111 0000")
    # cbc(252,UInt8) == 0x06 ("1111 1100")
    # cbc(127,UInt8) == error ("0111 1111")

    n = convert(t,n)
    invn = ~n
    bitct = log2(invn + 1)
    isinteger(bitct) || error("noncontiguous bits")

    bitct = floor(Int,bitct)
    return UInt8(sizeof(t)*8 - bitct)
end

function _mask2bits(t::Type, n::Unsigned)
    # takes a number of 1's bits in a
    # netmask and returns an integer representation
    maskbits = Int(width(t)) - Int(n)
    maskbits < 0 && throw(BoundsError())

    return (~(UInt128(2)^maskbits-1))
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 7680525e-633a-11eb-0bb3-6780ff3f200b
using Plots

# ╔═╡ 81dbfe6e-633a-11eb-0681-1b41e42d4faf
import Pkg

# ╔═╡ 85daa704-633a-11eb-0ef7-d1e801b03f0e
Pkg.add("ChemEquations")

# ╔═╡ db0ac6ce-6341-11eb-20a5-49cc6d9715a0
Pkg.add("Plots")

# ╔═╡ 5b0e83f2-6339-11eb-37fc-791c44bf9c27
md" ## _Balancing Chemical Equations_
The purpose of this notebook is to make simplifying chemical equations easier.
"

# ╔═╡ 2092a694-633a-11eb-33ee-17d107ae269b
md" ## Sythesis Reactions
Pattern: *A + B → AB*"

# ╔═╡ d59683e4-633f-11eb-0527-b1eabe8aae85
md" ## Decomposition
Pattern: *AB → A + B *"

# ╔═╡ 44c50086-6342-11eb-19a8-993fa02dd599
md" ## Single-Replacement
Pattern: *AB + C → AC + B *"

# ╔═╡ 46db65e0-6342-11eb-3fde-dffc073cfdfc

# ╔═╡ Cell order:
# ╠═81dbfe6e-633a-11eb-0681-1b41e42d4faf
# ╠═85daa704-633a-11eb-0ef7-d1e801b03f0e
# ╠═db0ac6ce-6341-11eb-20a5-49cc6d9715a0
# ╠═7680525e-633a-11eb-0bb3-6780ff3f200b
# ╠═5b0e83f2-6339-11eb-37fc-791c44bf9c27
# ╠═2092a694-633a-11eb-33ee-17d107ae269b
# ╠═d59683e4-633f-11eb-0527-b1eabe8aae85
# ╠═44c50086-6342-11eb-19a8-993fa02dd599
# ╠═46db65e0-6342-11eb-3fde-dffc073cfdfc


text:
# ## [Пример №3, развернуты, динамическая сетка](@id de_3)
#
# В этом примере:
# - вручную задаём параметры для динамической сетки,
# - решается прямая задача,
# - генерируются экспериментальные данные на динамической сетке.
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: shishkin_mesh;

# ### Заданием параметров для динамической сетки
u_l(t) = -8;                    # ГУ
u_r(t) =  4;                    #
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления, который в обратной
                                ## задаче необходимо определить, но при генерации априорной
                                ## информации мы задаем некоторый коэффициент, который,
                                ## собственно, после имея априорную информацию и будем определять.
ε = 0.01;                       # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.42;                # Область по T
x_tp = 0.12;                    # Положение переходного слоя
M = 500;                        # Кол-во разбиений по X, T
τ = (T-t₀)/M;                   # шаг по T
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
nothing #hide

# ### Создание и решение в случае динамической сетки
#
# Новое поведение будет реализоваться с помощью передачи функции создания сетки внутрь `solve` в качестве аргумента.
# Эта функция должна принимать только один аргумент — положение переходного слоя и формировать соответствующую сетку.
# В пакете есть формирование кусочно-равномерной сетки со сгущениями на границе и на переходном слое.
# Создадим замыкание этой функции, которое будет иметь нужную сигнатуру.
mshfrm(x_tp) = shishkin_mesh(a, b, x_tp, ε, 40, 0.5, 1.0, 1.0, 0.25);
Xₙ = mshfrm(x_tp); ;                            # Сетка по Х
N  = length(Xₙ) - 1                             # Примем за N длину сетки, что получилась в итоге.
qₙ =    qf.(Xₙ);                                # Сеточные значения коэффициента линейного усиления
u₀ = u_init.(Xₙ, ε=ε, x_tp = x_tp);             # Начальные условия

u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, create_mesh = mshfrm);
nothing #hide

# ### Генерация априорной информации
ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Правый вырожденный корень
ϕ       = Φ(ϕl, ϕr, N, M);                                          # Полуразность вырожденных корней
ϕ       = apply_on_dynamic_mesh(ϕ, XX, N, M);                       # Аппроксимация на переменную сетку
ϕl      = apply_on_dynamic_mesh(ϕl, XX, N, M);                      # Аппроксимация на переменную сетку
ϕr      = apply_on_dynamic_mesh(ϕr, XX, N, M);                      # Аппроксимация на переменную сетку
f1_data = f1(ϕ, u, XX, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, XX, N, M);                                 # Значение функции на переходном слое
nothing # hide


text:
function display(d::EMACS_DISPLAY, p::Plots.Plot)
    # TODO: decide which way is best to draw this.
    filename = tempname(ImagesDir()) * ".svg"
    Base.invokelatest(Plots.svg, p, filename)
    println(stdout,"\n<emacs-svg>$filename</emacs-svg>")
end


text:
###############################################################################
#
# Arrhenius corrections
#
###############################################################################
"""
    temperature_correction(
                td_set::AbstractTDParameterSet{FT},
                T::FT
    ) where {FT<:AbstractFloat}

A correction factor based on arrhenius's fitting procedure, given
- `td_set` [`ArrheniusTD`](@ref) or [`ArrheniusPeakTD`](@ref) type struct
- `T` Leaf temperature in `[K]`

The equation used for [`ArrheniusTD`](@ref) is
```math
corr = \\exp \\left( \\dfrac{ΔHa}{R T_0} - \\dfrac{ΔHa}{R T_1} \\right)
```

The equations used for [`ArrheniusPeakTD`](@ref) are
```math
corr = \\exp \\left( \\dfrac{ΔHa}{R T_0} - \\dfrac{ΔHa}{R T_1} \\right)
       \\cdot
       \\dfrac{ 1 + \\exp \\left( \\dfrac{S_v T_0 - H_d}{R T_0} \\right) }
              { 1 + \\exp \\left( \\dfrac{S_v T_1 - H_d}{R T_1} \\right) }
```

The equation used for [`Q10TD`](@ref) is
```math
corr = \\left( \\dfrac{T_1 - T_\\text{REF}}{10} \\right)^{Q_{10}}
```
"""
function temperature_correction(
            td_set::ArrheniusTD{FT},
            T::FT
) where {FT<:AbstractFloat}
    return exp( td_set.ΔHa_to_RT25 - td_set.ΔHa_to_R/T )
end

function temperature_correction(
            td_set::ArrheniusPeakTD{FT},
            T::FT
) where {FT<:AbstractFloat}
    @unpack C, ΔHa_to_RT25, ΔHd_to_R, ΔSv_to_R = td_set;

    # _f_a: activation correction, C/_f_b: de-activation correction
    _f_a::FT = exp( ΔHa_to_RT25 * (1 - T_25(FT)/T) );
    _f_b::FT = 1 + exp(ΔSv_to_R - ΔHd_to_R/T);

    return C /_f_b * _f_a
end

function temperature_correction(
            td_set::Q10TD{FT},
            T::FT
) where {FT<:AbstractFloat}
    return (td_set.Q_10) ^ ( (T - td_set.T_REF) / 10 )
end


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking

ProjDir = @__DIR__

# ### snippet 5.1

println()
df = DataFrame!(CSV.File(rel_path("..", "data", "WaffleDivorce.csv"), delim=';'));
scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
println()

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

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

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         //Priors
  bA ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bA * A;
  D ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5_1s = SampleModel("m5.1", m_5_1);

# Input data for cmdstan

ad_data = Dict(
  "N" => size(df, 1),
  "D" => df[!, :Divorce_s],
  "A" => df[!, :MedianAgeMarriage_s]
);

# Sample using StanSample

rc = stan_sample(m5_1s, data=ad_data);

if success(rc)

  # Describe the draws

  dfa1 = read_samples(m5_1s; output_format=:dataframe)

  # Result rethinking

  rethinking = "
           mean   sd  5.5% 94.5%
    a      0.00 0.10 -0.16  0.16
    bA    -0.57 0.11 -0.74 -0.39
    sigma  0.79 0.08  0.66  0.91
  "

  Particles(dfa1)

end


text:
using YaoBase

export ChainBlock, chain

"""
    ChainBlock{N} <: CompositeBlock{N}

`ChainBlock` is a basic construct tool to create
user defined blocks horizontically. It is a `Vector`
like composite type.
"""
struct ChainBlock{N} <: CompositeBlock{N}
    blocks::Vector{AbstractBlock{N}}
end

ChainBlock(blocks::Vector{<:AbstractBlock{N}}) where N = ChainBlock{N}(blocks)
ChainBlock(blocks::AbstractBlock{N}...) where N = ChainBlock(collect(AbstractBlock{N}, blocks))

"""
    chain(blocks...)

Return a [`ChainBlock`](@ref) which chains a list of blocks with same
[`nqubits`](@ref). If there is lazy evaluated
block in `blocks`, chain can infer the number of qubits and create an
instance itself.
"""
chain(blocks::AbstractBlock{N}...) where N = ChainBlock(blocks...)
chain(blocks::Union{AbstractBlock{N}, Function}...) where N = chain(map(x->parse_block(N, x), blocks)...)

function chain(list::Vector)
    for each in list # check type
        each isa AbstractBlock || error("expect a block, got $(typeof(each))")
    end
    N = nqubits(first(list))
    return ChainBlock(Vector{AbstractBlock{N}}(list))
end

# if not all matrix block, try to put the number of qubits.
chain(n::Int, blocks...) = chain(map(x->parse_block(n, x), blocks)...)
chain(n::Int, itr) = isempty(itr) ? chain(n) : chain(map(x->parse_block(n, x), itr)...)
chain(n::Int, f::Function) = chain(n, parse_block(n, f))
function chain(n::Int, block::AbstractBlock)
    @assert n == nqubits(block) "number of qubits mismatch"
    return ChainBlock(block)
end
chain(blocks::Function...) = @λ(n->chain(n, blocks...))
chain(it) = chain(it...) # forward iterator to vargs, so we could dispatch based on types
chain(blocks...) = @λ(n->chain(n, blocks))

"""
    chain(n)

Return an empty [`ChainBlock`](@ref) which can be used like a list of blocks.
"""
chain(n::Int) = ChainBlock{n}(AbstractBlock{n}[])

"""
    chain()

Return an lambda `n->chain(n)`.
"""
chain() = @λ(n->chain(n))

subblocks(c::ChainBlock) = c.blocks
occupied_locs(c::ChainBlock) =
    unique(Iterators.flatten(occupied_locs(b) for b in subblocks(c)))

chsubblocks(pb::ChainBlock{N}, blocks::Vector{<:AbstractBlock}) where N = length(blocks) == 0 ? ChainBlock{N}([]) : ChainBlock(blocks)
chsubblocks(pb::ChainBlock, it) = chain(it...)

mat(::Type{T}, c::ChainBlock) where T = prod(x->mat(T, x), Iterators.reverse(c.blocks))

function apply!(r::AbstractRegister, c::ChainBlock)
    for each in c.blocks
        apply!(r, each)
    end
    return r
end

cache_key(c::ChainBlock) = Tuple(cache_key(each) for each in c.blocks)

function Base.:(==)(lhs::ChainBlock{N}, rhs::ChainBlock{N}) where {N}
    (length(lhs.blocks) == length(rhs.blocks)) && all(lhs.blocks .== rhs.blocks)
end

Base.copy(c::ChainBlock{N}) where N = ChainBlock{N}(copy(c.blocks))
Base.similar(c::ChainBlock{N}) where {N} = ChainBlock{N}(empty!(similar(c.blocks)))
Base.getindex(c::ChainBlock, index) = getindex(c.blocks, index)
Base.getindex(c::ChainBlock, index::Union{UnitRange, Vector}) = ChainBlock(getindex(c.blocks, index))
Base.setindex!(c::ChainBlock{N}, val::AbstractBlock{N}, index::Integer) where N = (setindex!(c.blocks, val, index); c)
Base.insert!(c::ChainBlock{N}, index::Integer, val::AbstractBlock{N}) where N = (insert!(c.blocks, index, val); c)
Base.adjoint(blk::ChainBlock{N}) where {N} = ChainBlock{N}(map(adjoint, reverse(subblocks(blk))))
Base.lastindex(c::ChainBlock) = lastindex(c.blocks)
## Iterate contained blocks
Base.iterate(c::ChainBlock, st=1) = iterate(c.blocks, st)
Base.length(c::ChainBlock) = length(c.blocks)
Base.eltype(c::ChainBlock) = eltype(c.blocks)
Base.eachindex(c::ChainBlock) = eachindex(c.blocks)
Base.popfirst!(c::ChainBlock) = popfirst!(c.blocks)
Base.pop!(c::ChainBlock) = pop!(c.blocks)
Base.push!(c::ChainBlock{N}, m::AbstractBlock{N}) where N = (push!(c.blocks, m); c)
Base.push!(c::ChainBlock{N}, f::Function) where N = (push!(c.blocks, f(N)); c)
Base.pushfirst!(c::ChainBlock{N}, m::AbstractBlock{N}) where N = (pushfirst!(c.blocks, m); c)
Base.pushfirst!(c::ChainBlock{N}, f::Function) where N = (pushfirst!(c.blocks, f(N)); c)
Base.append!(c::ChainBlock{N}, list::Vector{<:AbstractBlock{N}}) where N = (append!(c.blocks, list); c)
Base.append!(c1::ChainBlock{N}, c2::ChainBlock{N}) where N = (append!(c1.blocks, c2.blocks); c1)
Base.prepend!(c1::ChainBlock{N}, list::Vector{<:AbstractBlock{N}}) where N = (prepend!(c1.blocks, list); c1)
Base.prepend!(c1::ChainBlock{N}, c2::ChainBlock{N}) where N = (prepend!(c1.blocks, c2.blocks); c1)

YaoBase.isunitary(c::ChainBlock) = all(isunitary, c.blocks) || isunitary(mat(c))
YaoBase.isreflexive(c::ChainBlock) = (iscommute(c.blocks...) && all(isreflexive, c.blocks)) || isreflexive(mat(c))
YaoBase.ishermitian(c::ChainBlock) = (all(isreflexive, c.blocks) && iscommute(c.blocks...)) || isreflexive(mat(c))


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Ipopt, Compat
include("convexpenalties.jl")
include("alt_min.jl")

gammaMult=5.0 # Need to vary this to generate the four plots

results_template = DataFrame(n=Int[], r=Int[], p=Real[], y=Real[],
    error_GD=Real[], error_AM=Real[], rank_GD=Int[], rank_AM=Int[], LB_GD=Real[], UB_GD=Real[], UB_AM=Real[])

results_basispursuit=similar(results_template, 0)

pRange=collect(0.1:0.1:1.0)
pSize=size(pRange, 1)
yRange=collect(0.1:0.1:1.0)
ySize=size(yRange, 1)

for ARG in ARGS # 1-100
  for iter in 1:5
  n=100

  array_num = parse(Int, ARG)

  y=yRange[(array_num-1)÷ySize+1]
  p=pRange[(array_num-1)%pSize+1]

  r=Int(round(n*(1.0-sqrt(1.0-y*p))))
  if r==0
      r=1
  end
  k=r #assume we picked it via say cross-validation (can only show two dimensions on the plot).
  M_L=randn(n, r)
  M_R=randn(n, r)
  M=M_L*M_R'
  numEntriesSampled=Int(round(p*n^2))
  gamma=gammaMult/p #setting this to be really high as a good proxy for exact constraints (no noise)

  coords=[repeat(1:n, inner=size(1:n,1)) repeat(1:n, outer=size(1:n,1))];
  entries=sample(1:n^2, Int(numEntriesSampled), replace=false)

  ############## Greedy rounding  #####################################################
  @show k
  X_GD,LB_GD, UB_GD,=getMatrixFrobeniusNorm_greedy(coords, entries, M, n, numEntriesSampled, k, gamma, false)

  u,v,w=svd(X_GD);
  @show rnk_GD=sum(v.>1e-2);
  @show err_GD = norm(X_GD.-M)/norm(M);

  ########### Alternating Minimization ###########################################################
  X_AM, UB_AM=alternatingminimization_matrixcompletion(coords, entries, M, n, numEntriesSampled, k, gamma, 50, X_GD)
  u,v,=svd(X_AM)
  @show rnk_AM=sum(v.>1e-2);
  @show err_AM = norm(X_AM.-M)/norm(M);

  ########### Push results #############################################################

  push!(results_basispursuit, [n, r, p, y, err_GD, err_AM, rnk_GD, rnk_AM, LB_GD, UB_GD, UB_AM])
end
  CSV.write("dualitygapplot_gamma5p0onp.csv", results_basispursuit, append=true)

end

max_stars_repo_path: Matrix Completion/createFigure5.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_path: Matrix Completion/createFigure5.jl

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_path: Matrix Completion/createFigure5.jl

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
"""
Spc fileformat
Described how elementary sectors can be collapsed. Can also contain descriptions
of other bigger airspaces (like ACC...)
See EUROCONTROL NEST Manual for Spc fileformat description

Output:
* name      Airspace name/description   String
* type      Airspace tuype              String  CS, CRSA, AUA, CLUS, NAS, AREA, AUAG, CRAS, REG, UNK
* sectors

Sub-airspace struct with the following elements:
* name      Sub-airspace name
* type      Sub-airspace type

Spc fileformat example
#;AIRSPACE;2;383;20140206;20140305;4082;EAR_P
A;BENELUX;EB/EH;AREA;2
S;EB;NAS
S;EH;NAS
A;BG;GREENLAND (DENMARK);NAS;1
S;BGGLFIR;FIR
A;BI;ICELAND;NAS;3
S;BIFAROER;FIR
S;BIRD;FIR
S;ENJA;FIR
A;BICC;ICELAND AUAG;AUAG;6
S;BIRDCTA;AUA
S;BIRDICTA;AUA
S;BIRDTMA;AUA
S;BIRDTOCA;AUA
S;EKVGCTR;AUA
S;ENJACTR;AUA
"""
module Spc

export read

struct Sector
    name::String
    type::String
end

struct Airspace
    name::String
    type::String
    sectors::Vector{Sector}
end

function read(file)
    airspaces = Dict{String, Airspace}()
    id = ""
    aname = ""
    type = ""
    sectors = Vector{Sector}()
    for line in eachline(file)
        line_elements = split(line, ';')
        if line_elements[1] == "A"
            if id != ""
                airspaces[id] = Airspace(aname, type, sectors)
            end
            id = line_elements[2]
            aname = line_elements[3]
            type = line_elements[4]
            sectors = Vector{Sector}()
        elseif line_elements[1] == "S"
            sname = line_elements[2]
            stype = line_elements[3]
            sector = Sector(sname, stype)
            sectors = vcat(sectors, sector)
        end
    end
    airspaces[id] = Airspace(aname, type, sectors)
    return airspaces
end

end # module


text:
@testset "Manopt.jl Error Measures" begin
    M = Sphere(2)
    N = Power(M,(2,))
    using Random: seed!
    seed!(42)
    w = randomMPoint(M)
    x = randomMPoint(M)
    y = randomMPoint(M)
    z = randomMPoint(M)
    a = PowPoint([w,x])
    b = PowPoint([y,z])
    @test meanSquaredError(M,x,y) == distance(M,x,y)^2
    @test meanSquaredError(N,a,b) == 1/2*distance(N,a,b)^2
    @test meanAverageError(M,x,y) == distance(M,x,y)
    @test meanAverageError(N,a,b) == 1/2*sum( distance.(Ref(M), getValue(a), getValue(b)) )
end


text:
mutable struct Whisker1
    id::Int32
    time::Int32
    len::Int32
    x::Array{Float32,1}
    y::Array{Float32,1}
    thick::Array{Float32,1}
    scores::Array{Float32,1}
end

Whisker1()=Whisker1(0,0,0,Array{Float32,1}(),Array{Float32,1}(),Array{Float32,1}(),Array{Float32,1}())

include("janelia_types.jl")

function Whisker1(w::Whisker2)
    id=w.id
    time=w.time
    len=w.len
    x=unsafe_wrap(Array,w.x,len)
    y=unsafe_wrap(Array,w.y,len)
    thick=unsafe_wrap(Array,w.thick,len)
    scores=unsafe_wrap(Array,w.scores,len)
    Whisker1(id,time,len,x,y,thick,scores)
end

mutable struct Tracker
    data_path::String
    file_name::String
    vid_name::String
    tracking_path::String
    tracking_name::String
    min_length::Int64 #Minimum size of traced element from Janelia tracker, in pixels
    mask::BitArray{2}
    whiskers::Array{Whisker1,1} #Properties of whiskers in active frame
    pad_pos::Tuple{Float32,Float32}
    h::Int64
    w::Int64
    all_whiskers::Array{Array{Whisker1,1},1} #Whiskers on every frame - IS THIS STILL NEEDED?
end

mutable struct discrete_widgets
    win::Gtk.GtkWindowLeaf
    space_button::Gtk.GtkSpinButtonLeaf
    points_button::Gtk.GtkSpinButtonLeaf
    calc_button::Gtk.GtkCheckButtonLeaf
    add_button::Gtk.GtkCheckButtonLeaf
    delete_button::Gtk.GtkButtonLeaf
end

mutable struct janelia_widgets
    win::Gtk.GtkWindowLeaf
    jt_seed_thres_button::Gtk.GtkSpinButtonLeaf
    jt_seed_iterations_button::Gtk.GtkSpinButtonLeaf
end

mutable struct contact_widgets
    win::Gtk.GtkWindowLeaf
    training_num_label::Gtk.GtkLabelLeaf
    fit_button::Gtk.GtkButtonLeaf
    load_predicted_button::Gtk.GtkButtonLeaf
    n_estimators_button::Gtk.GtkSpinButtonLeaf
    forest_depth_button::Gtk.GtkSpinButtonLeaf
    cv_label::Gtk.GtkLabelLeaf
    pred_pole_button::Gtk.GtkCheckButtonLeaf
    pred_pole_position::Gtk.GtkCheckButtonLeaf
    pred_curv::Gtk.GtkCheckButtonLeaf
end

mutable struct classifier
    predictors::Array{Float64,2}
    n_estimators::Int64
    forest_depth::Int64
    clf::PyObject
    cv::Float64
end

classifier()=classifier(zeros(Float64,1,1),100,10,PyObject(1),0.0)

mutable struct image_adjustment_settings
    contrast_min::Int64
    contrast_max::Int64
    sharpen_win ::Int64
    sharpen_reps::Int64
    sharpen_filter::Int64
end

image_adjustment_settings()=image_adjustment_settings(0,255,3,1,1)

mutable struct Normalize_Parameters
    mean_img::Array{Float32,3}
    std_img::Array{Float32,3}
    min_ref::Float32
    max_ref::Float32
    num::Int64 #Number of frames to calculate mean over
end

Normalize_Parameters() = Normalize_Parameters(zeros(Float32,0,0,0),zeros(Float32,0,0,0),0,0,10000)

mutable struct NeuralNetwork
    labels::Array{Float32,4} #Labels for training
    imgs::Array{Float32,4} #Images for training
    norm::Normalize_Parameters #Reference values to scale input images
    hg::StackedHourglass.NN #Deep Learning weights
    epochs::Int64
    losses::Array{Float32,1}
    confidence_thres::Float64
    predict_single::Bool
    normalize_inputs::Bool
    weight_path::String
    features::Int64
    draw_preds::Bool
    use_existing_weights::Bool
    use_existing_labels::Bool
    predicted::Array{Float32,3}
end

NeuralNetwork() = NeuralNetwork(zeros(Float32,0,0,0,0),zeros(Float32,0,0,0,0),Normalize_Parameters(), HG2(64,13,4),10,zeros(Float32,0),0.5,false,
true,pre_train_path,1,false,false,false,zeros(Float32,0,3,0))

mutable struct Save_Paths
    path::String
    temp::String
    images::String
    backup::String
    DLC::String
end

mutable struct DLC_Wrapper
end

function Save_Paths(mypath,make_dirs=true)

    if make_dirs
        mkdir(mypath)
    end

    if is_windows()
        out=Save_Paths(mypath,string(mypath,"\\temp"),string(mypath,"\\images"),string(mypath,"\\backup"),string(mypath,"\\DLC"))
    else
        out=Save_Paths(mypath,string(mypath,"/temp"),string(mypath,"/images"),string(mypath,"/backup"),string(mypath,"/DLC"))
    end

    if make_dirs
        mkdir(out.temp)
        mkdir(out.images)
        mkdir(out.backup)
        mkdir(out.DLC)
    end

    out
end

mutable struct Tracker_Handles
    frame::Int64 #currently active frame number

    b::Gtk.GtkBuilder
    max_frames::Int64
    h::Int64
    w::Int64
    fps::Float64
    frame_loaded::Bool
    requested_frame::Int64
    start_frame::Int64 #First frame for analysis (like taking average)
    end_frame::Int64 #last frame for analysis (like taking average)

    c::Gtk.GtkCanvasLeaf

    plot_frame::Array{UInt32,2}

    current_frame::Array{UInt8,2}
    current_frame2::Array{UInt8,2}

    woi_id::Int64 #Index in array of displayed whiskers which is whisker of interest.
    woi::Dict{Int64,WhiskerTracking.Whisker1} #Dictionary of properties for whisker of interest for every frame
    num_whiskers::Int64
    sw::Int64 #Selected Whisker

    auto_mode::Bool
    erase_mode::Bool

    tracked::Dict{Int64,Bool} #Array of true/false to specify if corresponding frame has been tracked

    combine_mode::Int64
    partial::Whisker1

    draw_mode::Bool

    touch_mode::Bool

    touch_override_mode::Bool
    touch_frames::BitArray{1}
    touch_frames_i::Array{Int64,1}
    wt::Tracker

    im_adj::image_adjustment_settings

    discrete_auto_calc::Bool
    d_spacing::Int64

    frame_list::Array{Int64,1}

    displayed_frame::Int64

    contact_widgets::contact_widgets

    pole_present::Dict{Int64,Bool}
    pole_loc::Dict{Int64,Array{Float32,1}}
    send_frame::Array{UInt8,2}

    selection_mode::Int64 #What the mouse will do when you click

    tracked_whiskers_x::Array{Float64,2}
    tracked_whiskers_y::Array{Float64,2}
    tracked_whiskers_l::BitArray{2}

    show_contact::Bool
    tracked_contact::BitArray{1}
    tracked_pole::Array{Float64,2}

    show_event::Bool
    event_array::BitArray{1}
    save_label_path::String

    class::classifier
    nn::NeuralNetwork

    paths::Save_Paths
    temp_frame::Array{UInt8,2}
end