text:
module ResultFileFunctions
	import HDF5
	using Dates
	import Statistics
	import  ..MasslistFunctions
	import ..InterpolationFunctions

	export MeasurementResult, joinResultsTime, joinResultsMasses, getTraces, getTimetraces, transposeStickCps, getNbrTraces, getTraceSamples, getNbrTraceSamples, findChangingMasses, findVaryingMasses, saturationFromComposition, getIndicesInTimeframe

	#type MeasurementResult
	mutable struct MeasurementResult
	    Times
	    MasslistMasses
	    MasslistElements
	    MasslistElementsMasses
	    MasslistCompositions
	    Traces
	end

	function joinResultsTime(firstResult::MeasurementResult, secondResult::MeasurementResult)
	    if (length(firstResult.Times) > 0) & (length(secondResult.Times) > 0)
		if (firstResult.MasslistMasses == secondResult.MasslistMasses)
		    firstResult.Times = vcat(firstResult.Times, secondResult.Times)
		    firstResult.MasslistCompositions = hcat(firstResult.MasslistCompositions, secondResult.MasslistCompositions)
		    firstResult.Traces = hcat(firstResult.Traces, secondResult.Traces)
		else
		println("Masses did not match, could not merge results!")
		end
	    end
	return firstResult
	end

	function joinResultsMasses(firstResult::MeasurementResult, secondResult::MeasurementResult)
	    if (length(firstResult.MasslistMasses) > 0) & (length(secondResult.MasslistMasses) > 0)
		if (firstResult.Times == secondResult.Times)
		    firstResult.MasslistMasses = vcat(firstResult.MasslistMasses, secondResult.MasslistMasses)
		    firstResult.MasslistCompositions = hcat(firstResult.MasslistCompositions, secondResult.MasslistCompositions)
		    firstResult.Traces = hcat(firstResult.Traces, secondResult.Traces)
		else
		println("Times did not match, could not merge results!")
		end
	    end
	return firstResult
	end

	function loadResults(filename; useAveragesOnly = false, raw = false, startTime::DateTime = DateTime(0), endTime::DateTime = DateTime(3000), massesToLoad=Array{Float64,1}(), massMatchTolerance = 0.00001, masslistOnly = false)
	    println("Loading Times:")
	  if useAveragesOnly
	    timesUnix = HDF5.h5read(filename, "AvgStickCpsTimes")
	  else
	      try
		  timesUnix = HDF5.h5read(filename, "Times")
	      catch
		  println("Could not load high res time, maybe this is an average-only result file?")
		  return []
	      end
	  end

	println("Looking for errors in time axis")

	  ############## Check for Errors in time axis ###################
	#  if ! issorted(timesUnix)
	#    println("Time in result file is not in accending order, selection might not be complete!")
	#  end
	  if (isnan(timesUnix[1])) | (timesUnix[1] < Dates.datetime2unix(DateTime(1975,1,1))) | (timesUnix[1] > Dates.datetime2unix(DateTime(2030,1,1)))
	      timesUnix[1] = timesUnix[1] - (timesUnix[3]-timesUnix[2])
	  end
	  if (isnan(timesUnix[length(timesUnix)]))| (timesUnix[length(timesUnix)] < Dates.datetime2unix(DateTime(1975,1,1))) | (timesUnix[length(timesUnix)] > Dates.datetime2unix(DateTime(2030,1,1)))
	      timesUnix[length(timesUnix)] = timesUnix[length(timesUnix)] + (timesUnix[length(timesUnix)-1]-timesUnix[length(timesUnix)-2])
	  end

	  for i=2:length(timesUnix)-1
	      if (isnan(timesUnix[i])) | (timesUnix[i] < Dates.datetime2unix(DateTime(1975,1,1))) | (timesUnix[i] > Dates.datetime2unix(DateTime(2030,1,1)))
		  timesUnix[i] = 0.5*(timesUnix[i-1]+timesUnix[i+1])
	      end
	  end

	  println("Selecting Times and Masses:")

	  masslistCompositions = HDF5.h5read(filename, "ElementalCompositions")
	  masslistMasses = HDF5.h5read(filename, "MassList")
	  masslistElements = HDF5.h5read(filename, "ElementNames")
	  masslistElementsMasses = HDF5.h5read(filename, "ElementMasses")

	  if masslistOnly
		  return MeasurementResult(Dates.unix2datetime.(timesUnix), masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions, Array{Float64}(undef,0,0))
	  end

	  if length(massesToLoad) > 0
	      selectionMassesIndices = Array{Int,1}()
	      for i=1:length(masslistMasses)
		for j=1:length(massesToLoad)
		  if (isapprox(masslistMasses[i], massesToLoad[j], atol=massMatchTolerance))
		    push!(selectionMassesIndices, i)
		  end
		end
	      end
	  else
	      selectionMassesIndices=1:length(masslistMasses)
	  end

	  selectionTimeIndexStart = searchsortedfirst(timesUnix, Dates.datetime2unix(startTime))
	  println("Start index: $selectionTimeIndexStart")
	  selectionTimeIndexEnd = searchsortedlast(timesUnix, Dates.datetime2unix(endTime))
	  println("End index: $selectionTimeIndexEnd")

	  if (selectionTimeIndexStart > 0) & (selectionTimeIndexEnd > selectionTimeIndexStart)
	    selectedTimesUnix = timesUnix[selectionTimeIndexStart:selectionTimeIndexEnd]
	  else
	    selectedTimesUnix = timesUnix
	    selectionTimeIndexStart = 1
	    selectionTimeIndexEnd = length(timesUnix)
	  end

	  ############# Load Traces here ##################

	  selectedTimes = Dates.unix2datetime.(selectedTimesUnix)

	  selectedMasslistMasses = masslistMasses[selectionMassesIndices]
	  selectedMassesCompositions = masslistCompositions[:,selectionMassesIndices]

	  if length(selectionMassesIndices) == 0
	    return MeasurementResult([], selectedMasslistMasses, masslistElements, masslistElementsMasses, Array{Float64}(undef,0,0), Array{Float64}(undef,0,0))
	  end

	  println("Loading Traces:")

	  traces = getTraces(filename, timeIndexStart=selectionTimeIndexStart, timeIndexEnd=selectionTimeIndexEnd, massIndices=selectionMassesIndices, raw=raw, useAveragesOnly=useAveragesOnly)

	  println("Loaded $(size(traces)) traces")
	  return MeasurementResult(selectedTimes, selectedMasslistMasses, masslistElements, masslistElementsMasses, selectedMassesCompositions, traces)
	end

	function getIndicesInTimeframe(filename, startTime::DateTime, endTime::DateTime)
	  times = HDF5.h5read(filename, "Times")
	  return (1:length(times))[(times.>startTime) & (times.<endTime)]
	end

	function getTraces(filename; timeIndexStart=1, timeIndexEnd=0, massIndices=nothing, raw=false,  useAveragesOnly = false)
	  dsTexists = false

	  fh = HDF5.h5open(filename,"r")
	  if useAveragesOnly
	    if raw
	      ds = fh["AvgStickCps"]
	    else
	      ds = fh["CorrAvgStickCps"]
	    end
	  else
	    if raw
	      ds = fh["StickCps"]
	    else
	      ds = fh["CorrStickCps"]
	      if HDF5.exists(fh, "CorrStickCpsT")
		  dsTexists = true
		  println("A transposed dataset for faster loading is available.")
		  dsT = fh["CorrStickCpsT"]
	      else
		  println("Transposed dataset for faster loading is NOT available. You can create one with ResultFileFunctions.transposeStickCps(filename).")
	      end
	    end
	  end
	  result = 0

	  if timeIndexEnd == 0
	    timeIndexEnd = size(ds)[1]
	    println("No end time given, using all $(size(ds)[1]) elements")
	  end

	  if massIndices == nothing
	    massIndices = ( 1:size(ds)[2])
	    println("No mass selection given, selecting all $(size(ds)[2]) masses.")
	  end

	  if typeof(massIndices) == BitArray{1}
	    massIndices = ( 1:length(massIndices) )[massIndices]
	  end

	  if ((timeIndexStart>0) & (timeIndexEnd<=size(ds)[1])) & (minimum(massIndices)>0) & (maximum(massIndices)<=size(ds)[2])
	    result = similar(ds[1,1],(timeIndexEnd-timeIndexStart+1),length(massIndices))
	    println("Creating result Matrix with size $(size(result))")
	    if length(massIndices) < size(ds)[2]
		println("Subset of masses requested")
		if dsTexists
		    println("Loading from transposed dataset.")
		    for j=1:length(massIndices)
		        result[:,j] = dsT[massIndices[j],timeIndexStart:timeIndexEnd]
		    end
		else
		    for j=1:length(massIndices)
		        result[:,j] = ds[timeIndexStart:timeIndexEnd,massIndices[j]]
		        #println("Loaded $((timeIndexEnd-timeIndexStart)*j) samples")
		    end
		end
	    else
		println("Many masses requested, slicing time first for better performance")
		for j=timeIndexStart:timeIndexEnd
		    tmp = ds[j,:]
		    result[j-timeIndexStart+1,:] = tmp[massIndices]
		    #println("Loaded $(length(massIndices)*j) samples")
		end
	    end
	  end

	  close(fh)
	  return result
	end

	function getTimetraces(filename, indices; raw=false)
	  fh = HDF5.h5open(filename,"r")
	  if raw
	    ds = fh["StickCps"]
	  else
	    ds = fh["CorrStickCps"]
	  end

	  result = 0

	  if typeof(indices) == BitArray{1}
	    indices = ( 1:length(indices) )[indices]
	  end

	  if (typeof(indices) == Array{Int,1}) | (typeof(indices) == Array{Integer})
	    if (minimum(indices)>0) & (maximum(indices)<size(ds)[2])
	      result = similar(ds[1,1],size(ds)[1],length(indices))
	      for i=1:length(indices)
		result[:,i] = ds[:,i]
	      end
	    end
	  end

	  if typeof(indices) == UnitRange
	    if (minimum(indices)>0) & (maximum(indices)<size(ds)[2])
	      result = ds[:,indices]
	    end
	  end

	  close(fh)
	  return result
	end

	function transposeStickCps(filename)
	  fh = HDF5.h5open(filename,"r+")
	  if ! HDF5.exists(fh, "CorrStickCps")
	      println("CorrStickCps not present, skipping transpose.")
	      HDF5.close(fh)
	      return
	  end
	  if HDF5.exists(fh, "CorrStickCpsT")
	      HDF5.o_delete(fh,"CorrStickCpsT")
	  end
	  dsOld = fh["CorrStickCps"]
	  nbrSpectra, nbrMasses = size(dsOld)
	  dsNew = HDF5.d_create(fh, "CorrStickCpsT", HDF5.datatype(Float32), HDF5.dataspace(nbrMasses, nbrSpectra), "chunk", (1,nbrSpectra), "compress", 3)
	  spectraAtOnce = 10000
	  for i=1:Int(floor(nbrSpectra/spectraAtOnce))
	      startIdx = (i-1)*spectraAtOnce + 1
	      endIdx = i*spectraAtOnce
	      print("Transposing Spectra $startIdx to $endIdx for faster loading...")
	      a=transpose(deepcopy(dsOld[startIdx:endIdx,:]))
	      for j=1:nbrMasses
		  dsNew[j,startIdx:endIdx] = a[j,:]
	      end
		  a=nothing
		  GC.gc()
	      println("DONE")
	  end
	  startIdx = (Int(floor(nbrSpectra/spectraAtOnce)))*spectraAtOnce + 1
	  endIdx = nbrSpectra
	  print("Transposing Spectra $startIdx to $endIdx for faster loading...")
	  a=transpose(deepcopy(dsOld[startIdx:endIdx,:]))
	  for j=1:nbrMasses
	      dsNew[j,startIdx:endIdx] = a[j,:]
	  end
	  a=nothing
	  println("DONE")
	  close(fh)
	  GC.gc()
	end

	function getTraceSamples(filename, indices; raw=false)
	  fh = HDF5.h5open(filename,"r")
	  if raw
	    ds = fh["StickCps"]
	  else
	    ds = fh["CorrStickCps"]
	  end
	  result = 0

	  if typeof(indices) == BitArray{1}
	    indices = ( 1:length(indices) )[indices]
	  end

	  if (typeof(indices) == Array{Int,1}) | (typeof(indices) == Array{Integer})
	    if (minimum(indices)>0) & (maximum(indices)<=size(ds)[1])
	      result = similar(ds[1,1],length(indices),size(ds)[2])
	      for i=1:length(indices)
		result[i,:] = ds[i,:]
	      end
	    end
	  end

	  if typeof(indices) == UnitRange{Int}
	    if (minimum(indices)>0) & (maximum(indices)<=size(ds)[1])
	      result = ds[indices,:]
	    end
	  end

	  close(fh)
	  if result == 0
	    println("Failed getting Traces Subset!")
	  end
	  return result
	end

	function getNbrTraces(filename)
	  fh = HDF5.h5open(filename,"r")
	  ds = fh["StickCps"]
	  nbr = size(ds)[2]
	  close(fh)
	  return nbr
	end

	function getNbrTraceSamples(filename)
	  fh = HDF5.h5open(filename,"r")
	  ds = fh["StickCps"]
	  nbr = size(ds)[1]
	  close(fh)
	  return nbr
	end

	function findChangingMasses(masses, compositions, traces, times, bgTimesSelection, signalTimesSelection; minOxygen = 0, sigmaThreshold=3, sorting = "mean", noNitrogen = true, onlySaneMasses = true, filterCrosstalkMasses=true)
	  bgTraces = traces[bgTimesSelection,:]
	  sigTraces = traces[signalTimesSelection,:]
	  meansBG = mean(bgTraces,1)[1,:]
	  meansBG[meansBG.<=0] = 1e-99
	  means = mean(sigTraces,1)[1,:] - meansBG
	  println("Means received: $(size(means))")
	  if size(bgTraces,1) < 3
	    println("Not enough BG points for standard deviation! plotHighTimeRes = true ?")
	  end

	  stderror = std(bgTraces,1)[1,:] ./ sqrt(size(bgTraces,2))

	  selMasses =  (means .> sigmaThreshold*stderror)
	  if filterCrosstalkMasses
	      # Filter out crosstalk masses
	      s = MasslistFunctions.filterMassListByContribution2(masses, means, 5000, 0.05)
	      selMasses = selMasses & s
	  end
	  println("\nRemoving $(length(masses[!s])) masses")

	  if noNitrogen == true
	    selMasses = selMasses & (compositions[5,:] .== 0)
	  end

	  selMasses = selMasses & (compositions[6,:] .>= minOxygen)
	  #selMasses = selMasses| ((masses.>137.1) & (masses.<137.15))

	  ## Filter only sane masses
	  if onlySaneMasses
	    selMasses = selMasses & (compositions[3,:] .> 1*compositions[1,:]) # H:C > 1.3
	    selMasses = selMasses & (compositions[3,:] .< 2.2*compositions[1,:]) # H:C < 2.2
	    selMasses = selMasses & (compositions[6,:] .< 1.5*compositions[1,:]) # O:C < 1.5
	  end
	  println("Selected $(sum(ones(length(masses))[selMasses])) masses!")

	  selIndices = linearindices(masses)[selMasses]
	  if sorting == "mean"
	      println("Sorting masses by mean value!")
	      sorter = sortperm(means[selMasses], rev=true)
	  elseif sorting == "mass"
	      println("Sorting masses!")
	    sorter = sortperm(masses[selMasses], rev=false)
	  else
	      println("Not sorting.")
	    sorter = linearindices(masses)
	  end

	  return selIndices[sorter], means[selIndices[sorter]], stderror[selIndices[sorter]]
	end

	function findVaryingMasses(masses, compositions, traces; minOxygen = 0, sigmaThreshold=3, sorting = "mean", noNitrogen = true, onlySaneMasses = true, filterCrosstalkMasses=true)
	  means = Statistics.mean(traces,dims=1)[1,:]
	  means[means .<= 0] .= 1e-99

	  println("Means received: $(size(means))")
	  if size(traces,1) < 3
	    println("Not enough points for variance! plotHighTimeRes = true ?")
	  end

	  tracesPrime = traces - InterpolationFunctions.smooth(traces, 0.01)
	  stderror = Statistics.std(tracesPrime,dims=1)[1,:] ./ sqrt(size(traces,2))
	  variances = Statistics.var(traces,dims=1)[1,:]
	  lastCount = length(masses)

	  selMasses =  (means .> sigmaThreshold*stderror)
	  println("Removed $(lastCount - sum(selMasses)) masses because mean<sigmaThreshold*stderr")
	  lastCount = sum(selMasses)
	  selMasses = selMasses .& (variances .> means .* sigmaThreshold)

	  println("Removed $(lastCount - sum(selMasses)) masses because variance<sigmaThreshold*mean")
	  lastCount = sum(selMasses)

	  if filterCrosstalkMasses
	      # Filter out crosstalk masses
	      s = MasslistFunctions.filterMassListByContribution2(masses, means, 5000, 0.05)
	      selMasses = selMasses .& s
	  end

	  println("\nRemoved $(lastCount - sum(selMasses)) masses because of crosstalk with neighbors")
	  lastCount = sum(selMasses)

	  println("\nRemoving $(length(masses[.!s])) masses")

	  if noNitrogen == true
	    selMasses = selMasses .& (compositions[5,:] .== 0)
	  end

	  selMasses = selMasses .& (compositions[6,:] .>= minOxygen)
	  #selMasses = selMasses .| ((masses .> 137.1) .& (masses .< 137.15))

	  ## Filter only sane masses
	  if onlySaneMasses
	    selMasses = selMasses .& (compositions[3,:] .> 1*compositions[1,:]) # H:C > 1.3
	    selMasses = selMasses .& (compositions[3,:] .< 2.2*compositions[1,:]) # H:C < 2.2
	    selMasses = selMasses .& (compositions[6,:] .< 1.5*compositions[1,:]) # O:C < 1.5
	  end
	  println("Selected $(sum(ones(length(masses))[selMasses])) masses!")
	  selIndices = (1:length(masses))[selMasses]
	  if sorting == "mean"
	      println("Sorting masses by mean value!")
	      sorter = sortperm(means[selMasses], rev=true)
	  elseif sorting == "mass"
	      println("Sorting masses!")
	    sorter = sortperm(masses[selMasses], rev=false)
	  else
	      println("Not sorting.")
	    sorter = linearindices(masses)
	  end

	  return selIndices[sorter], means[selIndices[sorter]], stderror[selIndices[sorter]]
	end

	function saturationFromComposition(compositions)
	  #return exp(10)*50000*exp.(-compositions[6,:]*5).*exp.(-compositions[1,:])
	  n_C = float(compositions[1,:])
	  n_O = float(compositions[6,:])
	  b_add = n_C
	  b_add[n_C .<= 15] = 0.92  # monomers
	  b_add[n_C .> 15] = 1.15 # dimers

	  return 10 .^((25-n_C)*0.475-n_O.*(2.3-b_add)-2 .*(n_C.*n_O./(n_C+n_O))*(-0.3))
	end
end


text:
################################################################################
#                                  MODEL QUERIES
################################################################################
# Simple model queries
for op in (:termination_status, :raw_status, :solve_time, :simplex_iterations,
           :barrier_iterations, :node_count, :objective_bound, :relative_gap)
    @eval begin
        @doc """
            JuMP.$($op)(model::InfiniteModel)

        Extend [`JuMP.$($op)`](https://jump.dev/JuMP.jl/v0.21.9/reference/solutions/#JuMP.$($op))
        for `InfiniteModel`s in accordance with that reported by its optimizer
        model. Errors if such a query is not supported or if the optimizer model
        hasn't be solved.
        """
        function JuMP.$op(model::InfiniteModel)
            return JuMP.$op(optimizer_model(model))
        end
    end
end

# Simple result dependent model queries
for op in (:primal_status, :dual_status, :has_values, :has_duals,
           :objective_value, :dual_objective_value)
    @eval begin
        @doc """
            JuMP.$($op)(model::InfiniteModel; [result::Int = 1])

        Extend [`JuMP.$($op)`](https://jump.dev/JuMP.jl/v0.21.9/reference/solutions/#JuMP.$($op))
        for `InfiniteModel`s in accordance with that reported by its optimizer
        model and the result index `result` of the most recent solution obtained.
        Errors if such a query is not supported or if the optimizer model hasn't
        be solved.
        """
        function JuMP.$op(model::InfiniteModel; result::Int = 1)
            return JuMP.$op(optimizer_model(model); result = result)
        end
    end
end

################################################################################
#                                 VALUE QUERIES
################################################################################
"""
    map_value([ref/expr], key::Val{ext_key_name}, result::Int; kwargs...)

Map the value(s) of `ref` to its counterpart in the optimizer model type that is
distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`, `optimizer_model_expression`, and/or
`optimizer_model_constraint`. Such as is the case with reformuations that do not
have a direct mapping between variables and/or constraints in the original
infinite form. Otherwise, `optimizer_model_variable`,
`optimizer_model_expression`, and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are passed on these functions. Here
`result` is the result index used in `value`.
"""
function map_value end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_value(vref::GeneralVariableRef, key, result::Int; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.value(v; result = result), opt_vref)
    else
        return JuMP.value(opt_vref; result = result)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_value(cref::InfOptConstraintRef, key, result::Int; kwargs...)
    func = JuMP.jump_function(JuMP.constraint_object(cref))
    return map_value(func, key, result; kwargs...)
end

# Default method that depends on optimizer_model_expression --> making extensions easier
function map_value(expr::JuMP.AbstractJuMPScalar, key, result::Int; kwargs...)
    opt_expr = optimizer_model_expression(expr, key; kwargs...)
    if opt_expr isa AbstractArray
        return map(e -> JuMP.value(e; result = result), opt_expr)
    else
        return JuMP.value(opt_expr; result = result)
    end
end

## Define dispatch methods to collect value of parameters
# InfiniteParameter
function _get_value(pref, ::Type{<:InfiniteParameterIndex}, result; kwargs...)
    label = get(kwargs, :label, PublicLabel)
    return supports(pref, label = label)
end

# FiniteParameter
function _get_value(pref, ::Type{FiniteParameterIndex}, result; kwargs...)
    return parameter_value(pref)
end

# Others
function _get_value(vref, index_type, result; kwargs...)
    return map_value(vref, Val(optimizer_model_key(JuMP.owner_model(vref))),
                     result; kwargs...)
end

"""
    JuMP.value(vref::GeneralVariableRef; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Extend `JuMP.value` to return the value(s) of `vref` in accordance with its
reformulation variable(s) stored in the optimizer model and the result index
`result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check
if a result exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

To provide context for the
results it may be helpful to also query the variable's `parameter_refs` and
`supports` which will have a one-to-one correspondence with the value(s).
It may also be helpful to query via [`optimizer_model_variable`](@ref) to
retrieve the variables(s) that these values are based on. These functions should
all be called with the same keyword arugments for consistency.

For extensions, this only works if
[`optimizer_model_variable`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> value(z)
42.0
```
"""
function JuMP.value(vref::GeneralVariableRef; result::Int = 1, kwargs...)
    return _get_value(vref, _index_type(vref), result; kwargs...)
end

"""
    JuMP.value(cref::InfOptConstraintRef; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Extend `JuMP.value` to return the value(s) of `cref` in accordance with its
reformulation constraint(s) stored in the optimizer model and the result index
`result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check if a result
exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

To provide context for
the results it may be helpful to also query the constraint's `parameter_refs`
and `supports` which will have a one-to-one correspondence with the value(s).
It may also be helpful to query via [`optimizer_model_constraint`](@ref) to
retrieve the constraint(s) that these values are based on. By default, only the
values corresponding to public supports are returned. These functions should
all be called with the same keyword arugments for consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> value(c1)
4-element Array{Float64,1}:
 -0.0
 20.9
 20.9
 20.9
```
"""
function JuMP.value(cref::InfOptConstraintRef; result::Int = 1,
                    kwargs...)
    return map_value(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                     result; kwargs...)
end

"""
    JuMP.value(expr::JuMP.AbstractJuMPScalar; [result::Int = 1,
               label::Type{<:AbstractSupportLabel} = PublicLabel,
               ndarray::Bool = false, kwargs...])

Return the value(s) of `expr` in accordance with the optimized variable values
the result index `result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check if a result
exists before asking for values.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the values associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the values of infinite
expressions are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the expression has multiple
infinite parameter dependencies.

To provide context for the results it may be
helpful to also query the expression's `parameter_refs` and `supports` which
will have a one-to-one correspondence with the value(s). It may also be helpful
to query via [`optimizer_model_expression`](@ref) to retrieve the expression(s)
that these values are based on. These should use the same keyword arguments for
consistency.

For extensions, this only works if [`optimizer_model_expression`](@ref) has been
extended correctly and/or [`map_value`](@ref) has been extended for expressions.

**Example**
```julia-repl
julia> value(my_finite_expr)
23.34

julia> value(my_infinite_expr)
4-element Array{Float64,1}:
 -0.0
 20.9
 20.9
 20.9
```
"""
function JuMP.value(
    expr::AbstractInfOptExpr;
    result::Int = 1,
    kwargs...
    )
    # get the model
    model = _model_from_expr(expr)
    # if no model then the expression only contains a constant
    if model === nothing
        expr isa NLPExpr && error("Cannot evaluate the value of `$expr`,",
                                  "because it doesn't have variables.")
        return JuMP.constant(expr)
    # otherwise let's call map_value
    else
        key = optimizer_model_key(model)
        return map_value(expr, Val(key), result; kwargs...)
    end
end

################################################################################
#                                 REDUCED COST
################################################################################
"""
    map_reduced_cost(vref::GeneralVariableRef, key::Val{ext_key_name},
                      result::Int; kwargs...)

Map the reduced cost(s) of `vref` to its counterpart in the optimizer model type that is
distininguished by its extension key `key` as type `Val{ext_key_name}`.
This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`. Such as is the case with reformulations
that do not have a direct mapping between variables in the original
infinite form. Otherwise, `optimizer_model_variable`, is used to make
these mappings by default where `kwargs` are passed on these functions. Here
`result` is the result index used in `value`.
"""
function map_reduced_cost end

# Default definition for when optimizer_model_variable is defined
function map_reduced_cost(vref::GeneralVariableRef, key; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.reduced_cost(v), opt_vref)
    else
        return JuMP.reduced_cost(opt_vref)
    end
end

"""
    JuMP.reduced_cost(vref::GeneralVariableRef)

Extend `JuMP.reduced_cost`. This returns the reduced cost(s) of a variable. This
will be a vector of scalar values for an infinite variable or will be a scalar
value for finite variables.

**Example**
```julia-repl
julia> reduced_cost(x)
12.81
```
"""
function JuMP.reduced_cost(vref::GeneralVariableRef; kwargs...)
    return map_reduced_cost(vref, Val(optimizer_model_key(JuMP.owner_model(vref)));
                            kwargs...)
end

################################################################################
#                             OPTIMIZER INDEX QUERIES
################################################################################
"""
    map_optimizer_index(ref, key::Val{ext_key_name}; kwargs...)

Map the `MathOptInterface` index(es) of `ref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are passed on as well.
"""
function map_optimizer_index end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_optimizer_index(vref::GeneralVariableRef, key; kwargs...)
    opt_vref = optimizer_model_variable(vref, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> JuMP.optimizer_index(v), opt_vref)
    else
        return JuMP.optimizer_index(opt_vref)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_optimizer_index(cref::InfOptConstraintRef, key; kwargs...)
    if JuMP.jump_function(JuMP.constraint_object(cref)) isa NLPExpr
        error("`optimizer_index` not defined for general nonlinear constraints.")
    end
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> JuMP.optimizer_index(c), opt_cref)
    else
        return JuMP.optimizer_index(opt_cref)
    end
end

"""
    JuMP.optimizer_index(vref::GeneralVariableRef;
                         [label::Type{<:AbstractSupportLabel} = PublicLabel,
                         ndarray::Bool = false, kwargs...])

Extend `JuMP.optimizer_index` to return the `MathOptInterface` index(es) of
`vref` in accordance with its reformulation variable(s) stored in the optimizer
model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the optimizer indices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the indices of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_variable`](@ref) to
retrieve the variables(s) that these indices are based on. These should use the
same keyword arguments for consistency.

For extensions, this
only works if [`optimizer_model_variable`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> optimizer_index(x)
4-element Array{MathOptInterface.VariableIndex,1}:
 MathOptInterface.VariableIndex(2)
 MathOptInterface.VariableIndex(3)
 MathOptInterface.VariableIndex(4)
 MathOptInterface.VariableIndex(5)
```
"""
function JuMP.optimizer_index(vref::GeneralVariableRef; kwargs...)
    return map_optimizer_index(vref, Val(optimizer_model_key(JuMP.owner_model(vref)));
                               kwargs...)
end

"""
    JuMP.optimizer_index(cref::InfOptConstraintRef;
                         [label::Type{<:AbstractSupportLabel} = PublicLabel,
                         ndarray::Bool = false, kwargs...])

Extend `JuMP.optimizer_index` to return the `MathOptInterface` index(es) of
`cref` in accordance with its reformulation constraints(s) stored in the
optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the optimizer indices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the indices of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_constraint`](@ref) to
retrieve the constraints(s) that these indices are based on. The same keyword
arguments should be used for consistency.

For extensions, this
only works if [`optimizer_model_constraint`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> optimizer_index(c1)
4-element Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},1}:
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(1)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(2)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(3)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(4)
```
"""
function JuMP.optimizer_index(cref::InfOptConstraintRef; kwargs...)
    return map_optimizer_index(cref, Val(optimizer_model_key(JuMP.owner_model(cref)));
                               kwargs...)
end

################################################################################
#                                  DUAL QUERIES
################################################################################
"""
    map_dual(cref::InfOptConstraintRef, key::Val{ext_key_name}, result::Int;
             kwargs...)

Map the dual(s) of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default where `kwargs` are also pass on to. Here `result` is
the result index that is used in `dual`.
"""
function map_dual end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_dual(cref::InfOptConstraintRef, key, result::Int; kwargs...)
    opt_cref = optimizer_model_constraint(cref, key; kwargs...)
    if opt_cref isa AbstractArray && first(opt_cref) isa JuMP.NonlinearConstraintRef
        return map(c -> JuMP.dual(c), opt_cref)
    elseif opt_cref isa AbstractArray
        return map(c -> JuMP.dual(c; result = result), opt_cref)
    elseif opt_cref isa JuMP.NonlinearConstraintRef
        return JuMP.dual(opt_cref)
    else
        return JuMP.dual(opt_cref; result = result)
    end
end

"""
    JuMP.dual(cref::InfOptConstraintRef; [result::Int = 1,
              label::Type{<:AbstractSupportLabel} = PublicLabel,
              ndarray::Bool = false, kwargs...])

Extend `JuMP.dual` to return the dual(s) of `cref` in accordance with its
reformulation constraint(s) stored in the optimizer model and the result index
`result` of the most recent solution obtained. Use
[`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel)) to check if a result
exists before asking for duals.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the duals associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the duals of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to
query via [`optimizer_model_constraint`](@ref) to retrieve the constraint(s)
that these duals are based on. Calling `parameter_refs` and `supports` may also
be insightful. Be sure to use the same keyword arguments for consistency.

For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_dual`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> dual(c1)
4-element Array{Float64,1}:
 -42.0
 -42.0
 32.3
 0.0
```
"""
function JuMP.dual(cref::InfOptConstraintRef; result::Int = 1, kwargs...)
    return map_dual(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                    result; kwargs...)
end

# Error redriect for variable call
function JuMP.dual(vref::GeneralVariableRef; kwargs...)
    return JuMP.dual(JuMP.VariableRef(JuMP.Model(), MOI.VariableIndex(1)))
end

################################################################################
#                              SHADOW PRICE QUERIES
################################################################################
# Dispatch functions for computing the shadow price
function _process_shadow_price(::MOI.LessThan, sense, duals)
    if sense == MOI.MAX_SENSE
        return -duals
    elseif sense == MOI.MIN_SENSE
        return duals
    else
        error(
            "The shadow price is not available because the objective sense " *
            "$sense is not minimization or maximization.",
        )
    end
end
function _process_shadow_price(::MOI.GreaterThan, sense, duals)
    if sense == MOI.MAX_SENSE
        return duals
    elseif sense == MOI.MIN_SENSE
        return -duals
    else
        error(
            "The shadow price is not available because the objective sense " *
            "$sense is not minimization or maximization.",
        )
    end
end
function _process_shadow_price(::MOI.EqualTo, sense, duals)
    if sense == MOI.MAX_SENSE
        return map(d -> d > 0 ? d : -d, duals)
    elseif sense == MOI.MIN_SENSE
        return map(d -> d > 0 ? -d : d, duals)
    else
        error(
            "The shadow price is not available because the objective sense " *
            "$sense is not minimization or maximization.",
        )
    end
end

"""
    JuMP.shadow_price(cref::InfOptConstraintRef;
                      [label::Type{<:AbstractSupportLabel} = PublicLabel,
                      ndarray::Bool = false, kwargs...])

Extend `JuMP.shadow_price` to return the shadow price(s) of `cref` in accordance
with its reformulation constraint(s) stored in the optimizer model. Use
[`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel)) to check if a result
exists before asking for the shadow price (it uses the duals).

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ.

By default only the shadow prices associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, the prices of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

It may also be helpful to query via [`optimizer_model_constraint`](@ref) to
retrieve the constraint(s) that these shadow prices are based on. Calling
`parameter_refs` and `supports` may also be insightful. Be sure to use the same
keyword arguments for consistency.

For extensions, this only works if [`optimizer_model_constraint`](@ref) has been
extended correctly and/or [`map_dual`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> shadow_price(c1)
4-element Array{Float64,1}:
 42.0
 42.0
 -32.3
 -0.0
```
"""
function JuMP.shadow_price(cref::InfOptConstraintRef; kwargs...)
    model = JuMP.owner_model(cref)
    set = JuMP.moi_set(JuMP.constraint_object(cref))
    sense = JuMP.objective_sense(model)
    if !JuMP.has_duals(model)
        error("The shadow price is not available because no dual result is " *
              "available.")
    end
    duals = map_dual(cref, Val(optimizer_model_key(model)), 1; kwargs...)
    return _process_shadow_price(set, sense, duals)
end

################################################################################
#                           LP SENSITIVITY ANALYSIS
################################################################################
"""
    InfOptSensitivityReport

A wrapper `DataType` for `JuMP.SensitivityReport`s in `InfiniteOpt`.
These are generated based on the optimizer model and should be made via the use of
[`lp_sensitivity_report`](@ref JuMP.lp_sensitivity_report(::InfiniteModel)). Once
made these can be indexed to get the sensitivies with respect to variables and/or
constraints. The indexing syntax for these is:
```julia
report[ref::[GeneralVariableRef/InfOptConstraintRef];
       [label::Type{<:AbstractSupportLabel} = PublicLabel,
       ndarray::Bool = false, kwargs...]]
```

This is enabled in user-defined optimizer model extensions by appropriately
extending [`optimizer_model_variable`](@ref) and [`optimizer_model_constraint`](@ref).

**Fields**
- `opt_report::JuMP.SensitivityReport`: The LP sensitivity captured from the optimizer model.
"""
struct InfOptSensitivityReport
    opt_report::JuMP.SensitivityReport
end

# Extend Base.getindex for variables on InfOptSensitivityReport
function Base.getindex(s::InfOptSensitivityReport, v::GeneralVariableRef; kwargs...)
    key = Val(optimizer_model_key(JuMP.owner_model(v)))
    opt_vref = optimizer_model_variable(v, key; kwargs...)
    if opt_vref isa AbstractArray
        return map(v -> s.opt_report[v], opt_vref)
    else
        return s.opt_report[opt_vref]
    end
end

# Extend Base.getindex for constraints on InfOptSensitivityReport
function Base.getindex(s::InfOptSensitivityReport, c::InfOptConstraintRef; kwargs...)
    key = Val(optimizer_model_key(JuMP.owner_model(c)))
    opt_cref = optimizer_model_constraint(c, key; kwargs...)
    if opt_cref isa AbstractArray
        return map(c -> s.opt_report[c], opt_cref)
    else
        return s.opt_report[opt_cref]
    end
end

"""
    JuMP.lp_sensitivity_report(model::InfiniteModel;
                               [atol::Float64 = 1e-8])::InfOptSensitivityReport

Extends `JuMP.lp_sensitivity_report` to generate and return an LP sensitivity
report in accordance with the optimizer model. See
[`InfOptSensitivityReport`](@ref) for syntax details on how to query it. `atol`
denotes the optimality tolerance and should match that used by the solver to
compute the basis. Please refer to `JuMP`'s documentation for more technical
information on interpretting the output of the report.

**Example**
```julia-repl
julia> report = lp_sensitivity_report(model);

julia> report[x]
(0.0, 0.5)
```
"""
function JuMP.lp_sensitivity_report(
    model::InfiniteModel;
    atol::Float64 = 1e-8
    )::InfOptSensitivityReport
    opt_model = optimizer_model(model)
    opt_report = JuMP.lp_sensitivity_report(opt_model, atol = atol)
    return InfOptSensitivityReport(opt_report)
end

ext: jl


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("car", "SLID"), x="Wages", color="Language", Stat.histogram, Geom.line)


text:
"""
Google Cloud Stackdriver Logging API
"""
module _logging

export logging

using ..api
using ...root

"""
Google Cloud Stackdriver Logging API root.
"""
logging = APIRoot(
    "https://logging.googleapis.com/v2",
    Dict(
        "cloud-platform" => "View and manage your data across Google Cloud Platform services",
        "cloud-platform.read-only" => "View your data across Google Cloud Platform services",
        "logging.admin" => "Administrate log data for your projects",
        "logging.read" => "View log data for your projects",
        "logging.write" => "Submit log data for your projects",
    );
    Entry=APIResource("entries";
        list=APIMethod(:POST, ":list", "Lists log entries."),
        write=APIMethod(:POST, ":write", "Writes log entries to Stackdriver Logging."),
    ),
    MonitoredResourceDescriptor=APIResource("monitoredResourceDescriptors";
        list=APIMethod(:GET, "", "Lists the monitored resource descriptors used by Stackdriver Logging."),
    ),
    OrganizationLog=APIResource("organizations/{organization}/logs";
        delete=APIMethod(:DELETE, "{log}", "Deletes a log and all its log entries."),
    ),
    OrganizationSink=APIResource("organizations/{organization}/sinks";
        create=APIMethod(:POST, "", "Creates a sink."),
        delete=APIMethod(:DELETE, "{sink}", "Deletes a sink."),
        get=APIMethod(:GET, "{sink}", "Gets a sink."),
        list=APIMethod(:GET, "", "Lists sinks."),
        update=APIMethod(:PUT, "{sink}", "Updates or creates a sink."),
    ),
    ProjectLog=APIResource("projects/{project}/logs";
        delete=APIMethod(:DELETE, "{log}", "Deletes a log and all its log entries."),
    ),
    ProjectMetric=APIResource("projects/{project}/metrics";
        create=APIMethod(:POST, "", "Creates a logs-based metric."),
        delete=APIMethod(:DELETE, "{metric}", "Deletes a logs-based metric."),
        get=APIMethod(:GET, "{metric}", "Gets a logs-based metric."),
        list=APIMethod(:GET, "", "Lists logs-based metrics."),
        update=APIMethod(:PUT, "{metric}", "Creates or updates a logs-based metric."),
    ),
    ProjectSink=APIResource("projects/{project}/sinks";
        create=APIMethod(:POST, "", "Creates a sink."),
        delete=APIMethod(:DELETE, "{sink}", "Deletes a sink."),
        get=APIMethod(:GET, "{sink}", "Gets a sink."),
        list=APIMethod(:GET, "", "Lists sinks."),
        update=APIMethod(:PUT, "{sink}", "Updates or creates a sink."),
    ),
)

end


text:
using Documenter, InverseLaplace

makedocs(
    format = :html,
    sitename = "InverseLaplace.jl",
    modules = [InverseLaplace],
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo = "github.com/jlapeyre/InverseLaplace.jl.git",
    target = "build",
    julia  = "0.7",
    osname = "linux",
    deps = nothing,
    make = nothing
)

# using Documenter, InverseLaplace

# makedocs(
#     modules = InverseLaplace,
# #    clean   = false,
# )

# # deploydocs(
# #     deps = Deps.pip("pygments", "mkdocs", "mkdocs-material"),
# #     repo = "github.com/MichaelHatherly/PrivateModules.jl.git",
# # )


text:
module QuakePAK

const _fileid = 0x5041434b      # The chars "PACK" used as a signature of PAK archives

struct ReadableFile <: IO
  _io::IO
  filename:: String
  offset:: Int
  size:: Int
end

Base.position(rf::ReadableFile) = position(rf._io) - rf.offset
Base.eof(rf.ReadableFile) = position(rf._io) < rf.offset + rf.size

Base.show(io::IO,  p::PakFileEntry) =
  print(io, "$(p.pak_filename)/$(p.filename) [offset: $(p.offset), size: $(p.size) byte]")

function read_file_entry(io::IO)
  pos = position(io)
  name = readuntil(io, '\0')
  seek(io, pos+56)
  offset = read(io, Int32)
  size = read(io, Int32)
  return name, offset, size
end

function open_pak(filename::String)
  fileentries = []
  open(filename) do f
    # Check if the file is a PAK-File
    @assert read(f, UInt32) == _fileid
    # Read meta data
    offset = read(f, Int32)      # Where the file entries start
    size = read(f, Int32)        # Size of file entry table
    nbentries = size÷64          # Each file entry is 64 bytes long
    # Iterate over all file entries
    seek(f, offset)
    for _ in 1:nbentries
      name, offset, size = read_file_entry(f)
      push!(fileentries, PakFileEntry(filename, name, offset, size))
    end
  end
  return fileentries
end

export open_pak

end # module


text:
using HDF5;

function read_dataset(data_path::String, train::Bool)

    # Read train/test dataset
    if train
        split = "train"
    else
        split = "test"
    end
    folders = open("$data_path/2018LA_Seg_Training Set/$split.list", "r") do file
        return readlines(file)
    end

    # Read all h5 images & masks
    images = []
    labels = []
    for folder in folders
        img, label = h5open("$data_path/2018LA_Seg_Training Set/$folder/mri_norm2.h5", "r") do file
            read(file, "image"), read(file, "label")
        end
        push!(images, img)
        push!(labels, label)
    end

    return images, labels
end

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
immutable IS <: InferenceAlgorithm
  n_samples :: Int
end

type ImportanceSampler{IS} <: Sampler{IS}
  alg :: IS
  model :: Function
  samples :: Array{Dict{Symbol,Any}}
  logweights :: Array{Float64}
  logevidence :: Float64
  predicts :: Dict{Symbol,Any}
  function ImportanceSampler(alg :: IS, model :: Function)
    samples = Array{Dict{Symbol,Any}}(alg.n_samples)
    for i = 1:alg.n_samples
      samples[i] = Dict{Symbol,Any}()
    end
    logweights = zeros(Float64, alg.n_samples)
    logevidence = 0
    predicts = Dict{Symbol,Any}()
    new(alg, model, samples, logweights, logevidence, predicts)
  end
end

function Base.run(spl :: Sampler{IS})
  n = spl.alg.n_samples
  for i = 1:n
    consume(Task(spl.model))
    spl.samples[i] = spl.predicts
    spl.logweights[i] = spl.logevidence
    spl.logevidence = 0
    spl.predicts = Dict{Symbol,Any}()
  end
  spl.logevidence = logsum(spl.logweights) - log(n)
  results = Dict{Symbol,Any}()
  results[:logevidence] = spl.logevidence
  results[:logweights] = spl.logweights
  results[:samples] = spl.samples
  return results
end

function assume(spl :: ImportanceSampler{IS}, distribution :: Distribution)
  return rand(distribution)
end

function param(spl :: ImportanceSampler{IS}, distribution :: Distribution)
  return assume(spl, distribution)
end

function observe(spl :: ImportanceSampler{IS}, score :: Float64)
  spl.logevidence += score
end

function predict(spl :: ImportanceSampler{IS}, name :: Symbol, value)
  spl.predicts[name] = value
end

sample(model :: Function, alg :: IS) =
  (global sampler = ImportanceSampler{IS}(alg, model); run(sampler))


text:
# function clust!(d::Dict, p::Any, l::Any)::Dict

dc = Dict("chr1::150::+" => [Int16(100),Int16(200)], "chr2::160::-" => [Int16(300),Int16(400)])
k1 = "chr1::150::+"
k2 = "chr2::160::-"
k3 = "chr3::170::."
vl1 = Int16(300)
vl2 = Int16(400)
vl3 = Int16(600)

@test Dict("chr1::150::+" => [100,200,300], "chr2::160::-" => [300,400]) == clust!(dc, k1, vl1)
@test Dict("chr1::150::+" => [100,200,300], "chr2::160::-" => [300,400,400]) == clust!(dc, k2, vl2)
@test Dict("chr1::150::+" => [100,200,300], "chr2::160::-" => [300,400,400], "chr3::170::." => [600] ) == clust!(dc, k3, vl3)

# function cluster!(d::Dict, chr::String, pos::Int, l::Array{Int16,1}, str::String, k::Int)::Dict

dict = Dict("chr1::10::+" => Dict(10 => [2,400,120]))
chr1 = "chr1"
pos1 = 8
l1 = [Int64(80)]
str1 = "+"
k1 = 5

@test Dict("chr1::10::+" => Dict(10=>[2,400,120], 8=>[1,80])) == cluster!(dict, chr1, pos1, l1, str1, k1)

cluster!(dict, chr1, 20, l1, str1, k1)
test_dict1 = Dict("chr1::10::+" => Dict(10=>[2,400,120], 8=>[1,80]), "chr1::20::+" => Dict(20=>[1,80]))

for key in sort(collect(keys(dict)))
    @test dict[key] == test_dict1[key]
end

cluster!(dict, chr1, 20, l1, str1, k1)
test_dict2 = Dict("chr1::10::+" => Dict(10=>[2,400,120], 8=>[1,80]), "chr1::20::+" => Dict(20=>[2,80,80]))

for key in sort(collect(keys(dict)))
    @test dict[key] == test_dict2[key]
end

# when center changes with a median.
for i in 1:10
    cluster!(dict, chr1, 20+i, l1, str1, k1)
    cluster!(dict, chr1, 20-i, l1, str1, k1)
end

test_dict3 = Dict(14=>[1,80], 10=>[3,400,120,80], 13=>[1,80], 11=>[1,80],
                  8=>[1,80], 12=>[1,80])
for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict3[i]
end

test_dict4 = Dict(16=>[1,80], 21=>[1,80], 25=>[1,80], 19=>[1,80], 17=>[1,80],
                  22=>[1,80], 20=>[2,80,80], 15=>[1,80], 18=>[1,80])
for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict4[i]
end

test_dict5 = Dict(23=>[1,80], 24=>[1,80], 25=>[1,80], 26=>[1,80], 29=>[1,80], 27=>[1,80],
                  28=>[1,80], 30=>[1,80])
for i in keys(dict["chr1::26::+"])
    @test dict["chr1::26::+"][i] == test_dict5[i]
end

# when center changes with frequency.
for i in 1:10
    cluster!(dict, chr1, 27, l1, str1, k1)
end

test_dict6 = Dict(14=>[1,80], 10=>[3, 400, 120,80], 13=>[1,80], 11=>[1,80],
                  8=>[1,80], 12=>[1,80])
for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict6[i]
end

test_dict7 = Dict(16=>[1,80], 21=>[1,80], 25=>[1,80], 19=>[1,80], 17=>[1,80],
                  22=>[1,80], 24=>[1,80], 20=>[2,80,80], 23=>[1,80], 15=>[1,80],
                  18=>[1,80])
for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict7[i]
end

test_dict8 = Dict(24=>[1,80], 25=>[1,80], 27=>[11,80,80,80,80,80,80,80,80,80,80,80],
                  29=>[1,80], 26=>[1,80], 28=>[1,80], 30=>[1,80])
for i in keys(dict["chr1::27::+"])
    @test dict["chr1::27::+"][i] == test_dict8[i]
end

# when center changes and recursion follows for discarded TS on the left.
for i in 1:5
    cluster!(dict, chr1, 40+i, l1, str1, k1)
end

cluster!(dict, chr1, 46, l1, str1, k1)

for i in 1:5
    cluster!(dict, chr1, 47, l1, str1, k1)
end

for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict6[i]
end

for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict7[i]
end

for i in keys(dict["chr1::27::+"])
    @test dict["chr1::27::+"][i] == test_dict8[i]
end

test_dict9 = Dict(43=>[1,80], 41=>[1,80], 42=>[1,80], 44=>[1,80])
for i in keys(dict["chr1::42::+"])
    @test dict["chr1::42::+"][i] == test_dict9[i]
end

test_dict10 = Dict(46=>[1,80], 47=>[5,80,80,80,80,80], 45=>[1,80])
for i in keys(dict["chr1::47::+"])
    @test dict["chr1::47::+"][i] == test_dict10[i]
end

# when center changes and recursion follows for discarded TS on the rigth.

for i in 1:5
    cluster!(dict, chr1, 71-i, l1, str1, k1)
end

for i in 1:2
    cluster!(dict, chr1, 65, l1, str1, k1)
end

cluster!(dict, chr1, 72, l1, str1, k1)

for i in 1:5
    cluster!(dict, chr1, 61, l1, str1, k1)
end

for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict6[i]
end

for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict7[i]
end

for i in keys(dict["chr1::27::+"])
    @test dict["chr1::27::+"][i] == test_dict8[i]
end

for i in keys(dict["chr1::42::+"])
    @test dict["chr1::42::+"][i] == test_dict9[i]
end

for i in keys(dict["chr1::47::+"])
    @test dict["chr1::47::+"][i] == test_dict10[i]
end

test_dict11 = Dict(61=>[5,80,80,80,80,80], 65=>[2,80,80], 66=>[1,80], 67=>[1,80])

for i in keys(dict["chr1::61::+"])
    @test dict["chr1::61::+"][i] == test_dict11[i]
end

test_dict12 = Dict(68=>[1,80], 69=>[1,80], 70=>[1,80], 72=>[1,80])
for i in keys(dict["chr1::70::+"])
    @test dict["chr1::70::+"][i] == test_dict12[i]
end

# when center changes and adjacent left cluster falls in.
for i in 1:5
    cluster!(dict, chr1, 99+i, l1, str1, k1)
end

for i in 1:2
    cluster!(dict, chr1, 105, l1, str1, k1)
end

cluster!(dict, chr1, 98, l1, str1, k1)

for i in 1:5
    cluster!(dict, chr1, 101, l1, str1, k1)
end

for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict6[i]
end

for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict7[i]
end

for i in keys(dict["chr1::27::+"])
    @test dict["chr1::27::+"][i] == test_dict8[i]
end

for i in keys(dict["chr1::42::+"])
    @test dict["chr1::42::+"][i] == test_dict9[i]
end

for i in keys(dict["chr1::47::+"])
    @test dict["chr1::47::+"][i] == test_dict10[i]
end

for i in keys(dict["chr1::61::+"])
    @test dict["chr1::61::+"][i] == test_dict11[i]
end

for i in keys(dict["chr1::70::+"])
    @test dict["chr1::70::+"][i] == test_dict12[i]
end

test_dict13 = Dict(100=>[1,80],102=>[1,80],98=>[1,80],101=>[6,80,80,80,80,80,80],
                   103=>[1,80],104=>[1,80],105=>[2,80,80])
for i in keys(dict["chr1::101::+"])
    @test dict["chr1::101::+"][i] == test_dict13[i]
end

# when center changes and adjacent right cluster falls in.
for i in 1:5
    cluster!(dict, chr1, 199+i, l1, str1, k1)
end

for i in 1:2
    cluster!(dict, chr1, 199, l1, str1, k1)
end

cluster!(dict, chr1, 205, l1, str1, k1)

for i in 1:5
    cluster!(dict, chr1, 205, l1, str1, k1)
end

for i in keys(dict["chr1::10::+"])
    @test dict["chr1::10::+"][i] == test_dict6[i]
end

for i in keys(dict["chr1::20::+"])
    @test dict["chr1::20::+"][i] == test_dict7[i]
end

for i in keys(dict["chr1::27::+"])
    @test dict["chr1::27::+"][i] == test_dict8[i]
end

for i in keys(dict["chr1::42::+"])
    @test dict["chr1::42::+"][i] == test_dict9[i]
end

for i in keys(dict["chr1::47::+"])
    @test dict["chr1::47::+"][i] == test_dict10[i]
end

for i in keys(dict["chr1::61::+"])
    @test dict["chr1::61::+"][i] == test_dict11[i]
end

for i in keys(dict["chr1::70::+"])
    @test dict["chr1::70::+"][i] == test_dict12[i]
end

for i in keys(dict["chr1::101::+"])
    @test dict["chr1::101::+"][i] == test_dict13[i]
end

test_dict14 = Dict(102=>[1,80], 103=>[1,80], 104=>[1,80], 105=>[2,80,80])
for i in keys(dict["chr1::105::+"])
    @test dict["chr1::105::+"][i] == test_dict14[i]
end

test_dict15 = Dict(199=>[2,80,80], 200=>[1,80], 201=>[1,80], 202=>[1,80])
for i in keys(dict["chr1::199::+"])
    @test dict["chr1::199::+"][i] == test_dict15[i]
end

test_dict16 = Dict(203=>[1,80], 204=>[1,80], 205=>[6,80,80,80,80,80,80])
for i in keys(dict["chr1::205::+"])
    @test dict["chr1::205::+"][i] == test_dict16[i]
end

# when center changes and adjecent cluster falls in with greater frequency of the center.
for i in 1:5
    cluster!(dict, chr1, 300+i, l1, str1, k1)
end

for i in 1:3
    cluster!(dict, chr1, 302, l1, str1, k1)
end

cluster!(dict, chr1, 307, l1, str1, k1)

for i in 1:5
    cluster!(dict, chr1, 306, l1, str1, k1)
end

test_dict17 = Dict(301=>[1,80], 302=>[4,80,80,80,80], 303=>[1,80], 304=>[1,80],
                   305=>[1,80], 306=>[5,80,80,80,80,80], 307=>[1,80])
for i in keys(dict["chr1::306::+"])
    @test dict["chr1::306::+"][i] == test_dict17[i]
end

# test all positions.

collection_of_positions = Int64[]
for i in keys(dict)
    for position in keys(dict[i])
        push!(collection_of_positions, position)
    end
end

test_collection_of_positions = [10, 8, 20, 21, 19, 22, 18, 23, 17, 24, 16, 25,
                                15, 26, 14, 27, 13, 28, 12, 29, 11, 30, 41, 42,
                                43, 44, 45, 46, 47, 70, 69, 68, 67, 66, 65, 61,
                                100, 101, 102, 103, 104, 105, 98, 199, 200, 201,
                                202, 203, 204, 205, 72, 301, 302, 303, 304, 305,
                                306, 307]

@test sort(test_collection_of_positions) == sort(collection_of_positions)

# function merge_adj_clusters!(d::Dict, m=Int64; verbose::Bool=false)::Dict

merge_adj_clusters!(dict, 1)

test_dict18 = Dict("chr1::205::+"=>Dict(200=>[1,80], 203=>[1,80], 199=>[2,80,80],
                                        204=>[1,80], 201=>[1,80], 205=>[6,80,80,80,80,80,80],
                                        202=>[1,80]),
                   "chr1::27::+"=>Dict(30=>[1,80], 16=>[1,80], 11=>[1,80], 21=>[1,80],
                                       26=>[1,80], 25=>[1,80], 10=>[3,400,120,80],
                                       29=>[1,80], 19=>[1,80], 17=>[1,80], 8=>[1,80],
                                       22=>[1,80], 24=>[1,80], 28=>[1,80], 20=>[2,80,80],
                                       23=>[1,80], 14=>[1,80], 13=>[1,80],
                                       27=>[11,80,80,80,80,80,80,80,80,80,80,80],
                                       15=>[1,80], 12=>[1,80], 18=>[1,80]),
                   "chr1::65::+"=>Dict(68=>[1,80], 66=>[1,80], 69=>[1,80], 67=>[1,80],
                                       72=>[1,80], 65=>[2,80,80], 70=>[1,80]),
                   "chr1::61::+"=>Dict(61=>[5,80,80,80,80,80]),
                   "chr1::101::+"=>Dict(100=>[1,80], 102=>[1,80], 98=>[1,80],
                                        101=>[6,80,80,80,80,80,80], 103=>[1,80],
                                        104=>[1,80], 105=>[2,80,80]),
                   "chr1::47::+"=>Dict(43=>[1,80], 47=>[5,80,80,80,80,80], 45=>[1,80],
                                       42=>[1,80], 41=>[1,80], 44=>[1,80], 46=>[1,80]),
                   "chr1::306::+"=>Dict(305=>[1,80], 306=>[5,80,80,80,80,80],
                                        304=>[1,80], 307=>[1,80], 301=>[1,80],
                                        303=>[1,80], 302=>[4,80,80,80,80]))

for k in keys(dict)
    @test test_dict18[k] == sort(dict[k])
end

# function stats_poly_a(d::Dict{String,Array{Int16,1}})::DataFrame

@test DataFrame(Chrmosome=String["chr1","chr2","chr3"],
             Position=Int32[150,160,170],
             Strand=String["+","-","."],
             Median=Float16[200,400,600],
             Minimum=Int16[100,300,600],
             Maximum=Int16[300,400,600],
             Counts=Int32[3,3,1]) == sort!(stats_poly_a(dc))

# function parserecord(r::GFF3.Record)::String

record = GFF3.Record(b"chr1\tHAVANA\texon\t43930132\t43931077\t.\t+\t.\tID=exon:ENST00000643283.1:12;Parent=ENST00000643283.1;gene_id=ENSG00000126091.20;transcript_id=ENST00000643283.1;gene_type=protein_coding;gene_name=ST3GAL3;transcript_type=nonsense_mediated_decay;transcript_name=RP11-7O11.1-079;exon_number=12;exon_id=ENSE00003818117.1;level=2;protein_id=ENSP00000494746.1;tag=RNA_Seq_supported_only;havana_gene=OTTHUMG00000007561.18;havana_transcript=OTTHUMT00000494570.1")
record2 = GFF3.Record(b"chr1\tHAVANA\texon\t43930132\t43931077\t.\t+\t.\tParent=ENST00000643283.1;gene_id=ENSG00000126091.20;transcript_id=ENST00000643283.1;gene_type=protein_coding;gene_name=ST3GAL3;transcript_type=nonsense_mediated_decay;transcript_name=RP11-7O11.1-079;exon_number=12;exon_id=ENSE00003818117.1;level=2;protein_id=ENSP00000494746.1;tag=RNA_Seq_supported_only;havana_gene=OTTHUMG00000007561.18;havana_transcript=OTTHUMT00000494570.1")
record3 = GFF3.Record(b"chr1\tHAVANA\texon\t43930132\t43931077\t.\t+\t.\tgene_id=ENSG00000126091.20;transcript_id=ENST00000643283.1;gene_type=protein_coding;gene_name=ST3GAL3;transcript_type=nonsense_mediated_decay;transcript_name=RP11-7O11.1-079;exon_number=12;exon_id=ENSE00003818117.1;level=2;protein_id=ENSP00000494746.1;tag=RNA_Seq_supported_only;havana_gene=OTTHUMG00000007561.18;havana_transcript=OTTHUMT00000494570.1")
record4 = GFF3.Record(b"chr1\tHAVANA\texon\t43930132\t43931077\t.\t+\t.\t")
record5 = GFF3.Record(b"chr1\tHAVANA\texon\t43930132\t43931077\t.\t+\t.\tgene_id=ENSG00000126091.20")

@test String("ID=exon:ENST00000643283.1:12;"*
             "Parent=ENST00000643283.1;"*
             "gene_id=ENSG00000126091.20;"*
             "transcript_id=ENST00000643283.1;"*
             "gene_type=protein_coding;"*
             "gene_name=ST3GAL3;"*
             "ftype=exon") == parserecord(record)
@test String("ID=;"*
             "Parent=ENST00000643283.1;"*
             "gene_id=ENSG00000126091.20;"*
             "transcript_id=ENST00000643283.1;"*
             "gene_type=protein_coding;"*
             "gene_name=ST3GAL3;"*
             "ftype=exon") == parserecord(record2)
@test String("ID=;"*
             "Parent=;"*
             "gene_id=ENSG00000126091.20;"*
             "transcript_id=ENST00000643283.1;"*
             "gene_type=protein_coding;"*
             "gene_name=ST3GAL3;"*
             "ftype=exon") == parserecord(record3)
@test String("ID=;"*
             "Parent=;"*
             "gene_id=;"*
             "transcript_id=;"*
             "gene_type=;"*
             "gene_name=;"*
             "ftype=exon") == parserecord(record4)
@test String("ID=;"*
             "Parent=;"*
             "gene_id=ENSG00000126091.20;"*
             "transcript_id=;"*
             "gene_type=;"*
             "gene_name=;"*
             "ftype=exon") == parserecord(record5)

# function getintervals(dframe::DataFrame)::IntervalCollection{String}

tdf = DataFrame(Chrmosome=String["chr1","chr2","chr3"],
                Position=Int32[150,160,170],
                Strand=String["+","-","."],
                Median=Float16[200,400,600],
                Minimum=Int16[100,300,600],
                Maximum=Int16[300,400,600],
                Counts=Int32[3,3,1],
                Smaple=String["A","A","A"],
                TotalReads=Int64[10000,10000,10000],
                PassedReads=Int64[1000,1000,1000],
                PolyAReads=Int64[500,500,500]
                )

@test IntervalCollection([Interval("chr1",150,150,Strand('+'),"Median=200.0;Min=100;Max=300;Counts=3")
                          Interval("chr2",160,160,Strand('-'),"Median=400.0;Min=300;Max=400;Counts=3")
                          Interval("chr3",170,170,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                          ]) == getintervals(tdf)

# function annotate_polya_sites(a::IntervalCollection,  b::Dict{String, IntervalCollection{String}})::DataFrame

intcol1 = IntervalCollection([Interval("chr1",150,150,Strand('+'),"Median=200.0;Min=100;Max=300;Counts=3")
                          Interval("chr2",160,160,Strand('-'),"Median=400.0;Min=300;Max=400;Counts=3")
                          Interval("chr3",170,170,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                          Interval("chr3",3000,3000,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                          ])

intcol2 = Dict("chr1:1:10000" => IntervalCollection([Interval("chr1",150,150,Strand('+'),
                                                          String("ID=exon:ENST00000643283.1:12;"*
                                                                 "Parent=ENST00000643283.1;"*
                                                                 "gene_id=ENSG00000126091.20;"*
                                                                 "transcript_id=ENST00000643283.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL3;"*
                                                                 "ftype=exon"))
                                                         Interval("chr1",160,160,Strand('-'),
                                                          String("ID=exon:ENST00000643282.1:12;"*
                                                                 "Parent=ENST00000643282.1;"*
                                                                 "gene_id=ENSG00000126092.20;"*
                                                                 "transcript_id=ENST00000643282.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL32;"*
                                                                 "ftype=intron"))
                                                         Interval("chr1",170,170,Strand('.'),
                                                          String("ID=exon:ENST00000643283.4:12;"*
                                                                 "Parent=ENST00000643283.4;"*
                                                                 "gene_id=ENSG00000126094.20;"*
                                                                 "transcript_id=ENST00000643284.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL43;"*
                                                                 "ftype=exon"))
                                                         ]),
               "chr2:1:10000" => IntervalCollection([Interval("chr2",150,150,Strand('+'),
                                                          String("ID=exon:ENST00000643283.1:12;"*
                                                                 "Parent=ENST00000643283.1;"*
                                                                 "gene_id=ENSG00000126091.20;"*
                                                                 "transcript_id=ENST00000643283.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL3;"*
                                                                 "ftype=exon"))
                                                         Interval("chr2",160,160,Strand('-'),
                                                          String("ID=exon:ENST00000643282.1:12;"*
                                                                 "Parent=ENST00000643282.1;"*
                                                                 "gene_id=ENSG00000126092.20;"*
                                                                 "transcript_id=ENST00000643282.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL32;"*
                                                                 "ftype=intron"))
                                                         Interval("chr2",170,170,Strand('.'),
                                                          String("ID=exon:ENST00000643283.4:12;"*
                                                                 "Parent=ENST00000643283.4;"*
                                                                 "gene_id=ENSG00000126094.20;"*
                                                                 "transcript_id=ENST00000643284.1;"*
                                                                 "gene_type=protein_coding;"*
                                                                 "gene_name=ST3GAL43;"*
                                                                 "ftype=exon"))
                                                         ]),
               "chr3:1:10000" => IntervalCollection([Interval("chr3",150,150,Strand('+'),
                                                      String("ID=exon:ENST00000643283.1:12;"*
                                                             "Parent=ENST00000643283.1;"*
                                                             "gene_id=ENSG00000126091.20;"*
                                                             "transcript_id=ENST00000643283.1;"*
                                                             "gene_type=protein_coding;"*
                                                             "gene_name=ST3GAL3;"*
                                                             "ftype=exon"))
                                                     Interval("chr3",160,160,Strand('-'),
                                                      String("ID=exon:ENST00000643282.1:12;"*
                                                             "Parent=ENST00000643282.1;"*
                                                             "gene_id=ENSG00000126092.20;"*
                                                             "transcript_id=ENST00000643282.1;"*
                                                             "gene_type=protein_coding;"*
                                                             "gene_name=ST3GAL32;"*
                                                             "ftype=intron"))
                                                     Interval("chr3",170,170,Strand('.'),
                                                      String("ID=exon:ENST00000643283.4:12;"*
                                                             "Parent=ENST00000643283.4;"*
                                                             "gene_id=ENSG00000126094.20;"*
                                                             "transcript_id=ENST00000643284.1;"*
                                                             "gene_type=protein_coding;"*
                                                             "gene_name=ST3GAL43;"*
                                                             "ftype=exon"))
                                                     ])
               )

@test DataFrame(Chr=String["chr1","chr2","chr3","chr3"],
                Start=Int64[149,159,169,2999],
                End=Int64[150,160,170,3000],
                Name=String["ST3GAL3.1","ST3GAL32.1","ST3GAL43.1","NA.1"],
                Counts=Int32[3,3,1,1],
                Strand=String["+","-",".","."],
                Feature=String["exon","intron","exon","NA"],
                Median=Float32[200.0,400.0,600.0,600.0],
                Min=Int16[100,300,600,600],
                Max=Int16[300,400,600,600],
                Biotype=String["protein_coding","protein_coding","protein_coding","NA"]
                ) == sort!(annotate_polya_sites(intcol1, intcol2), [:Chr, :Start])

intcol1 = IntervalCollection([Interval("chr1",1500,1500,Strand('+'),"Median=200.0;Min=100;Max=300;Counts=3")
                              Interval("chr1",3000,3000,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                              Interval("chr1",10000,10000,Strand('-'),"Median=400.0;Min=300;Max=400;Counts=3")
                              Interval("chr1",10001,10001,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                              Interval("chr1",17000,17000,Strand('.'),"Median=600.0;Min=600;Max=600;Counts=1")
                              ])

intcol2 = Dict("chr1:1:10000" => IntervalCollection([Interval("chr1",1,1000,Strand('+'),
                                                              String("ID=exon:ENST00000643283.1:12;"*
                                                                     "Parent=ENST00000643283.1;"*
                                                                     "gene_id=ENSG00000126091.20;"*
                                                                     "transcript_id=ENST00000643283.1;"*
                                                                     "gene_type=protein_coding;"*
                                                                     "gene_name=1;"*
                                                                     "ftype=exon"))
                                                     Interval("chr1",1,20000,Strand('.'),
                                                              String("ID=exon:ENST00000643283.4:12;"*
                                                                     "Parent=ENST00000643283.4;"*
                                                                     "gene_id=ENSG00000126094.20;"*
                                                                     "transcript_id=ENST00000643284.1;"*
                                                                     "gene_type=protein_coding;"*
                                                                     "gene_name=2;"*
                                                                     "ftype=exon"))
                                                     Interval("chr1",1000,10000,Strand('-'),
                                                              String("ID=exon:ENST00000643282.1:12;"*
                                                                     "Parent=ENST00000643282.1;"*
                                                                     "gene_id=ENSG00000126092.20;"*
                                                                     "transcript_id=ENST00000643282.1;"*
                                                                     "gene_type=protein_coding;"*
                                                                     "gene_name=3;"*
                                                                     "ftype=intron"))
                                                         ]),
                "chr1:10001:20000" => IntervalCollection([Interval("chr1",1,20000,Strand('.'),
                                                                   String("ID=exon:ENST00000643283.4:12;"*
                                                                          "Parent=ENST00000643283.4;"*
                                                                          "gene_id=ENSG00000126094.20;"*
                                                                          "transcript_id=ENST00000643284.1;"*
                                                                          "gene_type=three_prime_UTR;"*
                                                                          "gene_name=4;"*
                                                                          "ftype=exon"))
                                                          Interval("chr1",150,150,Strand('+'),
                                                                   String("ID=exon:ENST00000643283.1:12;"*
                                                                          "Parent=ENST00000643283.1;"*
                                                                          "gene_id=ENSG00000126091.20;"*
                                                                          "transcript_id=ENST00000643283.1;"*
                                                                          "gene_type=protein_coding;"*
                                                                          "gene_name=5;"*
                                                                          "ftype=exon"))
                                                          Interval("chr1",160,160,Strand('-'),
                                                                   String("ID=exon:ENST00000643282.1:12;"*
                                                                          "Parent=ENST00000643282.1;"*
                                                                          "gene_id=ENSG00000126092.20;"*
                                                                          "transcript_id=ENST00000643282.1;"*
                                                                          "gene_type=protein_coding;"*
                                                                          "gene_name=6;"*
                                                                          "ftype=intron"))
                                                      ])
                  )

@test DataFrame(Chr=String["chr1","chr1","chr1","chr1","chr1"],
              Start=Int64[1499,2999,9999,10000,16999],
              End=Int64[1500,3000,10000,10001,17000],
              Name=String["NA.1","2.1","3.1","4.1","4.2"],
              Counts=Int32[3,1,3,1,1],
              Strand=String["+",".","-",".","."],
              Feature=String["NA","exon","intron","exon","exon"],
              Median=Float32[200.0,600.0,400.0,600.0,600.0],
              Min=Int16[100,600,300,600,600],
              Max=Int16[300,600,400,600,600],
              Biotype=String["NA","protein_coding","protein_coding","three_prime_UTR","three_prime_UTR"]
              ) == sort!(annotate_polya_sites(intcol1, intcol2), [:Chr, :Start])

# function rmdups(dframe::DataFrame)::DataFrame

tdf = DataFrame(Chr=String["chr1","chr1","chr2","chr2","chr2","chr2"],
                Start=Int64[150,160,120,120,120,180],
                End=Int64[150,160,120,120,120,180],
                Name=String["A.1","A.2","B.1","B.1","B.1","C.1"],
                Counts=Int32[1,2,4,4,4,5],
                Strand=String["+","-","+","+","+","+"],
                Feature=String["exon","transcript","gene","transcript","CDS","gene"],
                Median=Float32[200.0,400.0,600.0,600.0,600.0,100.0],
                Min=Int16[100,300,100,100,100,200],
                Max=Int16[300,400,800,800,800,400],
                Biotype=String["AC","AC","BC","BC","BC","GC"]
                )

tdf2 = DataFrame(Chr=String["chr1","chr1","chr2","chr2"],
                 Start=Int64[150,160,120,180],
                 End=Int64[150,160,120,180],
                 Name=String["A.1","A.2","B.1","C.1"],
                 Counts=Int32[1,2,4,5],
                 Strand=String["+","-","+","+"],
                 Feature=String["exon","transcript","CDS","gene"],
                 Median=Float32[200.0,400.0,600.0,100.0],
                 Min=Int16[100,300,100,200],
                 Max=Int16[300,400,800,400],
                 Biotype=String["AC","AC","BC","GC"]
                 )

tdf4 = DataFrame(Chr=String["chr1","chr1","chr2","chr2","chr2","chr2"],
                 Start=Int64[150,160,120,120,120,180],
                 End=Int64[150,160,120,120,120,180],
                 Name=String["A.1","A.2","B.1","B.2.1","B.3.1","C.1"],
                 Counts=Int32[1,2,4,4,4,5],
                 Strand=String["+","-","+","+","+","+"],
                 Feature=String["exon","transcript","gene","transcript","CDS","gene"],
                 Median=Float32[200.0,400.0,600.0,600.0,600.0,100.0],
                 Min=Int16[100,300,100,100,100,200],
                 Max=Int16[300,400,800,800,800,400],
                 Biotype=String["AC","AC","BC","BC","BC","GC"]
                 )

tdf5 = DataFrame(Chr=String["chr1","chr1","chr2","chr2","chr2","chr2"],
                 Start=Int64[150,160,120,120,120,180],
                 End=Int64[150,160,120,120,120,180],
                 Name=String["A.1","A.2","B.1","B.2.1","B.3.1","C.1"],
                 Counts=Int32[1,2,4,4,4,5],
                 Strand=String["+","-","+","+","+","+"],
                 Feature=String["exon","transcript","gene","transcript","CDS","gene"],
                 Median=Float32[200.0,400.0,600.0,600.0,600.0,100.0],
                 Min=Int16[100,300,100,100,100,200],
                 Max=Int16[300,400,800,800,800,400],
                 Biotype=String["AC","AC","BC","BC","BC","GC"]
                 )

@test tdf2 == rmdups(tdf)
@test tdf5 == rmdups(tdf4)

# function enumeratenames!(df::DataFrame)::DataFrame

tdf3 = DataFrame(Chr=String["chr1","chr1","chr2","chr2","chr2","chr2"],
                 Start=Int64[150,160,120,120,120,180],
                 End=Int64[150,160,120,120,120,180],
                 Name=String["A.1.1","A.2.1","B.1.1","B.1.2","B.1.3","C.1.1"],
                 Counts=Int32[1,2,4,4,4,5],
                 Strand=String["+","-","+","+","+","+"],
                 Feature=String["exon","transcript","gene","transcript","CDS","gene"],
                 Median=Float32[200.0,400.0,600.0,600.0,600.0,100.0],
                 Min=Int16[100,300,100,100,100,200],
                 Max=Int16[300,400,800,800,800,400],
                 Biotype=String["AC","AC","BC","BC","BC","GC"]
                 )

@test tdf3 == enumeratenames!(tdf)

# function get_split_key(chr::String, x::Int; step::Int64=10000)::String

@test ["chr1:1:10"] == get_split_key("chr1",2,4; step=10)
@test ["chr2:1291:1300",
       "chr2:1301:1310",
       "chr2:1311:1320"] == get_split_key("chr2",1300,1315; step=10)
@test ["chr2:1291:1300",
      "chr2:1301:1310",
      "chr2:1311:1320",
      "chr2:1321:1330"] == get_split_key("chr2",1300,1325; step=10)
@test ["chr4:1:10000"] == get_split_key("chr4",1300,1300; step=10000)
@test ["chr4:1:10000"] == get_split_key("chr4",10000,10000; step=10000)
@test ["chr4:10001:20000"] == get_split_key("chr4",10290,10290; step=10000)
@test ["chr4:10001:20000"] == get_split_key("chr4",10590,10590; step=10000)
@test ["chr4:500001:510000"] == get_split_key("chr4",500590,500590; step=10000)
@test ["chr21:45070001:45080000"] == get_split_key("chr21",45073853,45073853)
@test ["chr21:46680001:46690000"] == get_split_key("chr21", 46689753,46689753)
@test ["chr1:11:20",
       "chr1:21:30",
       "chr1:31:40"] == get_split_key("chr1",12,35; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",12,25; step=10)
@test ["chr1:11:20"] == get_split_key("chr1",11,20; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",11,21; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",20,30; step=10)
@test ["chr1:11:20",
       "chr1:21:30",
       "chr1:31:40"] == get_split_key("chr1",20,31; step=10)
@test ["chr1:11:20"] == get_split_key("chr1",12,20; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",12,21; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",11,22; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",11,30; step=10)
@test ["chr1:11:20",
       "chr1:21:30",
       "chr1:31:40"] == get_split_key("chr1",11,31; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",12,30; step=10)
@test ["chr1:11:20",
       "chr1:21:30",
       "chr1:31:40"] == get_split_key("chr1",11,32; step=10)
@test ["chr1:11:20",
       "chr1:21:30",
       "chr1:31:40"] == get_split_key("chr1",12,32; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",11,29; step=10)
@test ["chr1:11:20",
       "chr1:21:30"] == get_split_key("chr1",12,29; step=10)
@test ["chr1:11:20"] == get_split_key("chr1",11,19; step=10)
@test ["chr1:11:20"] == get_split_key("chr1",12,19; step=10)

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
module CrossDecomposition
using StatsBase
import StatsBase.coef

include("canoncor.jl")

export canoncor, coef, scores
end # module


text:
## this script shows how to do IV estimation, in the GMM form, using
#  the well-known Card data on returns to schooling

##
using Econometrics, CSV, DataFrames, DataFramesMeta, Statistics
cd(@__DIR__)
card = CSV.read("../Julia/cooked.csv", DataFrame)   # this is the prepared data
                                                    # from the script BasicDataAnalysis.jl
display(first(card,6))

## first, OLS
y = Matrix(@select(card, :lnwage))
x = Matrix(@select(card, :educ, :exper, :expsq, :black, :smsa, :south))
x = [ones(size(y,1)) x]
names = ["const", "educ", "exp", "expsq", "black", "smsa", "south"]
βols, Vols, junk = ols(y,x,names=names);

## define instruments
w = Matrix(@select(card, :nearc4, :age, :agesq, :black, :smsa, :south))
w, junk = stnorm(w) # this type of scaling the istruments is good for numeric accuracy
w = [ones(size(w,1)) w]

## define moments
moments = θ -> w.*(y-x*θ)

## CUE
βcue, junk, Vcue, junk =  gmmresults(moments, βols, "", "CUE GMM", names, true);

## two-step
βgmm1, junk =  gmm(moments, βols, eye(7))
ms = moments(βgmm1)
W = inv(cov(ms))
βgmm, junk =  gmmresults(moments, βgmm1, W, "two step GMM", names, true);

## Note: this is all just identified, so the following should
# all be the same. Without the scaling instruments, step1 will
# be different, which means it didn't actually find the minimizer
prettyprint([βcue βgmm1 βgmm], ["cue", "step1", "step2"])

## Are the instruments ok? Look at first stage regressions
R = [zeros(3) eye(3) zeros(3,3)]
r = zeros(3)
y = Matrix{Float64}(@select(card, :educ))
F1 = (TestStatistics(y, w, R, r)[1]) # first return of TestStatistics is F
y = Matrix(@select(card, :exper))
F2 = (TestStatistics(y, w, R, r)[1])
y = Matrix(@select(card, :expsq))
F3 = (TestStatistics(y, w, R, r)[1])
println("F test of instrument strength for the 3 endog variables")
prettyprint([F1 F2 F3],["educ", "age", "agesq"])
println("according to the simple rule that F should be ≥ 10")
println("there is concern that education is not well-instrumented")

## Hausman test: compare CUE and OLS
# this is assuming that OLS is efficient, which
# is a bit of a stretch (possible HET, at least)
using LinearAlgebra, Distributions
e = βcue-βols
H = dot(e, inv(Vcue-Vols), e)
df = size(e,1)
pval = 1. - cdf(Chisq(df),H)
println("Hausman test statistic: $H, degrees of freedom: $df, p-value: $pval")

##
# to summarize, given the weakness of the nearc4 instrument,
# and the Hausman test, which gives some support to OLS,
# there's not too much reason to favor the IV results over
# the OLS results

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
module Changepoints

using Distributions
using Distributions: Sampleable
using Base.Meta
using RecipesBase
using Statistics
using Statistics: median
import Base.rand

export PELT, @PELT, BS, @BS, CROPS, @segment_cost, NormalVarSegment, NormalMeanSegment, NormalMeanVarSegment, PoissonSegment, BetaSegment, ExponentialSegment, ChangepointSampler, GammaShapeSegment, GammaRateSegment, NonparametricSegment, OLSSegment, @changepoint_sampler, WBS, @WBS, get_WBS_changepoints, MOSUM, @MOSUM, get_cps_mosum, @MOSUM_multi_scale, sSIC, MOSUM_multi_scale

include("segment_costs.jl")
include("PELT.jl")
include("CROPS.jl")
include("BS.jl")
include("sim_type.jl")
include("macros.jl")
include("WBS.jl")
include("MOSUM.jl")
include("plots.jl")

end # module


text:
function var_param_map_and_boundschecks(;
    nexpressions::Int,
    variables::Vector{Symbol},
    parameters::Vector{Symbol},
    t::Union{Symbol,Nothing} = nothing,
    # taylor = false,
    U::Bool = false,
    u::Bool = false,
)
    var_map = Dict{Symbol,Union{Symbol,Expr}}()
    for (i, v) in enumerate(variables)
        var_map[v] = :(x[$i])
    end
    for (i, v) in enumerate(parameters)
        var_map[v] = :(p[$i])
    end
    if t !== nothing
        var_map[t] = :(t)
    end
    m = length(variables)
    l = length(parameters)
    checks = Expr[]
    if U
        push!(checks, :(@boundscheck checkbounds(U, 1:$nexpressions, 1:$m)))
    end
    if u
        push!(checks, :(@boundscheck u === nothing || checkbounds(u, 1:$nexpressions)))
    end
    push!(checks, :(@boundscheck checkbounds(x, 1:$m)))
    push!(checks, :(@boundscheck p === nothing || checkbounds(p, 1:$l)))

    var_map, Expr(:block, checks...)
end

function unroll_pow(var, n)
    n == 0 && return :(one($var))
    n == 1 && return var
    n == 2 && return :(sqr($var))
    n == 3 && return :(@fastmath $var * sqr($var))
    n == 4 && return :(sqr(sqr($var)))
    n == 5 && return :(@fastmath $var * sqr(sqr($var)))

    # base to expansion shows up which power it is needed to compute
    d = digits(n, base = 2)
    x = :x
    exprs = [:(local $(Symbol(x, 1)) = $var)]
    for i = 2:length(d)
        push!(exprs, :(local $(Symbol(x, 1 << (i - 1))) = sqr($(Symbol(x, 1 << (i - 2))))))
    end
    prods = Symbol[]
    for (i, di) in enumerate(d)
        if !is_zero(di)
            push!(prods, Symbol(x, 1 << (i - 1)))
        end
    end
    if length(prods) > 1
        push!(exprs, :(@fastmath $(Expr(:call, :*, prods...))))
    end
    quote
        let
            $(exprs...)
        end
    end
end

function generate_evaluate_impl(
    list::InstructionList,
    out::Vector;
    variables::Vector{Symbol},
    t = nothing,
    parameters = Symbol[],
)
    var_param_map, boundschecks = var_param_map_and_boundschecks(
        nexpressions = length(out),
        variables = variables,
        parameters = parameters,
        t = t,
        u = true,
        U = false,
    )
    slp = to_expr(list; variables = var_param_map)
    assign = map(1:length(out), out) do i, a
        :(u[$i] = $(to_expr_arg(get(var_param_map, a, a))))
    end
    quote
        $boundschecks
        @inbounds begin
            $slp
            $(assign...)
        end
        u
    end
end

function generate_evaluate_and_jacobian_impl(
    list::InstructionList,
    out::Vector;
    variables::Vector{Symbol},
    t = nothing,
    parameters = Symbol[],
)
    grad_list, out, jac_out = ModelKit.gradient(list, out, variables)

    var_param_map, boundschecks = var_param_map_and_boundschecks(
        nexpressions = length(out),
        variables = variables,
        parameters = parameters,
        t = t,
        u = true,
        U = true,
    )
    slp = to_expr(grad_list; variables = var_param_map)
    jac_assign = Expr[]
    all_assigned = true
    for l = 1:length(variables), k = 1:length(out)
        a = jac_out[k, l]
        if a != 0
            push!(jac_assign, :(U[$k, $l] = $(to_expr_arg(get(var_param_map, a, a)))))
        else
            all_assigned = false
        end
    end

    quote
        $boundschecks
        $((!all_assigned ? (:(fill!(U, zero(eltype(U)))),) : ())...)
        @inbounds begin
            $slp
            if !(u isa Nothing)
                $(
                    (
                        :(u[$i] = $(to_expr_arg(get(var_param_map, a, a)))) for
                        (i, a) in enumerate(out)
                    )...
                )
            end
            $(jac_assign...)
        end
        nothing
    end
end

function generate_taylor_impl(
    list::InstructionList,
    out::Vector;
    variables::Vector{Symbol},
    t = nothing,
    parameters = Symbol[],
)
    var_param_map, boundschecks = var_param_map_and_boundschecks(
        nexpressions = length(out),
        variables = variables,
        parameters = parameters,
        t = t,
        u = true,
        U = false,
    )
    slp = to_expr(list; variables = var_param_map) do op, args
        Expr(:call, taylor_call_op(op), :(Val(K)), args...)
    end
    assign_taylor = map(1:length(out), out) do i, a
        :(u[$i] = $(to_expr_arg(get(var_param_map, a, a))))
    end
    assign_highest = map(1:length(out), out) do i, a
        :(u[$i] = last($(to_expr_arg(get(var_param_map, a, a)))))
    end
    if !isnothing(t)
        quote
            $boundschecks
            @inbounds let t = (t, one(t))
                $slp
                if u isa TaylorVector
                    $(assign_taylor...)
                else
                    $(assign_highest...)
                end
            end
            u
        end
    else
        quote
            $boundschecks
            @inbounds begin
                $slp
                if u isa TaylorVector
                    $(assign_taylor...)
                else
                    $(assign_highest...)
                end
            end
            u
        end
    end
end


text:
#=
fld-read.jl
Jeff Fessler
=#

export fld_header, fld_read

using AVSfldIO # temporary
import FileIO # temporary

@deprecate fld_header(args... ; kwargs...) AVSfldIO.fld_header(args... ; kwargs...)
@deprecate fld_read(args... ; kwargs...) FileIO.load(args... ; kwargs...)


text:
using TeachSigmoid
using Test

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


text:
s = LightGraphs.sample!([1:10;], 3)
@test length(s) == 3
for  e in s
    @test 1 <= e <= 10
end

s = LightGraphs.sample!([1:10;], 6, exclude=[1,2])
@test length(s) == 6
for  e in s
    @test 3 <= e <= 10
end


text:
################################################################################
# Animations
################################################################################

function animate_drone!(vis::Visualizer, anim::MeshCat.Animation, model::AbstractGameModel,
		s::AbstractVector; name::Symbol=:Robot)
	default_background!(vis)
	for t in 1:length(s)
		MeshCat.atframe(anim, t) do
			set_drone!(vis, model, s[t], name=name)
		end
	end
	setanimation!(vis, anim)
	return nothing
end

function visualize_drone!(vis::Visualizer, model::AbstractGameModel, s::AbstractVector;
		h=0.01, α=1.0,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=:Robot)

	build_drone!(vis, model, name=name, α=α)
	animate_drone!(vis, anim, model, s, name=name)
	return anim
end

function visualize_drone!(vis::Visualizer, model::AbstractGameModel, traj::Algames.Traj;
		sample=max(1, Int(floor(traj[1].dt*length(traj) / 100))), h=traj[1].dt*sample, α=1.0,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=:Robot)

	visualize_drone!(vis, model, Algames.state.(traj); anim=anim, name=name, h=h, α=α)
	return anim
end


text:
using Distributions
import Distributions: logpdf, pdf, sample

using ConjugatePriors
import ConjugatePriors: NormalInverseGamma, NormalWishart, NormalInverseWishart, posterior

include("transform.jl")

import Base.LinAlg: Cholesky
import Base.Random: rand
import Base: getindex, mean, scale

export NormalInverseGamma, NormalWishart, NormalInverseWishart, posterior

flatten(x :: Float64, d :: Distribution) = ( [link(d,x)], a -> invlink(d,a[1]) )
flatten(xs :: Array{Float64}, d :: Distribution) = ( link(d,xs), a -> invlink(d,a) )

function flatten(xs :: Matrix{Float64}, d :: Distribution)
  ys = link(d,xs)
  shape = size(ys)
  rebuild(zs) = invlink(d, reshape(zs, shape))
  return (ys[:], rebuild)
end

function flatten(x :: Tuple{Real,Real}, d :: NormalInverseGamma)
  y = link(d,x)
  z = [y[1],y[2]]
  return (z, a -> invlink(d, (a[1],a[2])))
end


text:
using Random, ProgressMeter

import Base: copy

copy(strat::Strategy) = Strategy(strat.universe, strat.indicator, strat.rules)

#TODO: parallel processing
function optimize(strat::Strategy; samples::Int=0, seed::Int=0, verbose::Bool=true, summary_fun::Function=cum_pnl, args...)::Matrix
    original = copy(strat)
    if samples > 0
        seed >= 0 ? Random.seed!(seed) : nothing
        sample_index = rand(collect(1:samples), samples)
    else
        samples = count_runs(strat.indicator.paramset)
        sample_index = collect(1:samples)
    end
    combos = generate_combinations(strat.indicator.paramset)[sample_index,:]
    optimization = zeros(samples, 1)
    verbose ? progress = Progress(length(sample_index), 1, "Optimizing Backtest") : nothing
    for (i, combo) in enumerate(sample_index)
        verbose ? next!(progress) : nothing
        strat.indicator.paramset.arg_defaults = combos[i,:]
        generate_trades!(strat, verbose=false)
        backtest!(strat, verbose=false; args...)
        optimization[i] = summary_fun(strat.backtest)
    end
    # prevent out-of-scope alteration of strat object
    strat = original
    return [combos optimization]
end

function optimize!(strat::Strategy; samples::Int=0, seed::Int=0, verbose::Bool=true, summary_fun::Function=cum_pnl, args...)::Nothing
    optimization = optimize(strat, samples=samples, seed=seed, verbose=verbose, summary_fun=summary_fun; args...)
    strat.backtest.optimization = optimization
    return nothing
end