text:
"""
    CodeType

Struct for holding parms for a custum code type for either the `I` or `Q`
channel

Fields:

- `name::String`: name of code type
- `code_num::Int`: number of codes, including databits, on channel
- `codes::Array{Dict,1}`: array of dictionaries, where each dictionary is a
                          layer of code
- `chipping_rates::Array{Float,1}`: array of chipping rates in Hz for each code
- `include_codes::Array{Bool,1}`: array of `Bool` flags for each code to indicate
                                  whether a given code is being used
- `databits::Bool`: `Bool` flag to specify whether databits used in the set of codes
    * if `true`, databits are the last code in `CodeType.codes` array
- `similar_databits::Bool`: set to true if the databits are the same for all PRNs
"""
struct CodeType{T1,T2,T3,T4,T5}
    name::String           # Name of code type
    code_num::T1           # Number of codes in signal type
                           # Includes primary, secondary, and nav codes
    codes::T2              # Vector of dictionaries
                           # Each dictionary is a primary, secondary, or
                           # or nav code which is specific for each PRN number
                           # The last dictionary in the vector of `codes` is
                           # assumed to be the navbits. This "code" will be
                           # ignored when generating replica signals for
                           # acquisition, fine acquisition, and signal tracking
                           # methods. Primary codes are considered the
                           # first code in `codes`.
    chipping_rates::T3     # Vector of chipping rates of length `code_num`
    code_lengths::T4       # Vector of code lengths of length `code_num`
    include_codes::T5      # Bool vector of length `code_num`
                           # Used to determine if a given code is used for
                           # signal generation.
                           # Will be used in `calc_code_val` method to determine
                           # if a code value will be calculated and XOR`ed to
                           # other code values.
    databits::Bool         # Flag for if there are databits present in code
                           # Databits are assumed to be last entry in the
                           # vector of codes given to definecodetype funciton.
    similar_databits::Bool # True if databits are similar for all PRNs
end

"""
    SignalType

Holds common parms for a signal including its codes for each channel.

Fields:

- `name:String`: name of the signal type
- `I_codes::CodeType`: Custum codes for the I channel of the signal
    * set to `missing` if no `I_codes` given
- `Q_codes::CodeType`: Custum codes for the Q channel of the signal
    * set to `missing` if no `I_codes` given
- `sig_freq::Float64`: carrier frequency of the signal in Hz
- `B_I::Float64`: maximum bandwidth of all I channel codes in Hz
- `B_Q::Float64`: maximum bandwidth of all Q channel codes in Hz
- `include_I::Bool`: set to `true` if `I_codes` is given
- `include_Q::Bool`: set to `true` if `Q_codes` is given
"""
struct SignalType{T1,T2,T3,T4,T5,T6}
    name::T1
    I_codes::T2
    Q_codes::T3
    sig_freq::T4
    B_I::T5
    B_Q::T6
    include_I::Bool
    include_Q::Bool
end

"""
    definecodetype(codes::Vector, chipping_rates::Vector{Float64};
                   databits=missing, name="custom")

Defines a `CodeType` struct that represents the codes of a signal on a
given channel, `I`, `Q`, or `both`.

Required Arguments:

- `codes::Vector`: array of codes, `not including databits`
    * the code is considered the primary code `(NOTE: this must be a dictionary)`
    * dictionary keys are the PRN numbers
    * all other codes after can be either dictionaries or arrays
    * if the code is an array, it is assumed that that layer of code is the
      same for all PRNs
- `chipping_rates::Vector{Float}`: array of chipping rates for codes above,
                                   `not including databits`

Optional Arguments:

- `databits::Vector`: if specified by user, must be a two element
                      array where the first element is the databits
                      as either a dictionary or array and the second
                      element is the chipping rate in Hz of the
                      databits `(default = missing)`
- `name::String`: name of code type `(default is "custom")`

Returns:

- `CodeType` structure
"""
function definecodetype(codes::Vector, chipping_rates::Vector{Float64};
                        databits=missing, name="custom")
    # Check that the first code in `codes`, the primary code, is a Dict
    if ~isa(codes[1], Dict)
        error("Primary code must be a dictionary of codes. The primary code is the first index of `codes`.")
    end
    # Setup for if there is or isn't databits
    if ismissing(databits)
        code_num = length(codes)
        N = code_num
    else
        code_num = length(codes) + 1
        N = code_num - 1
    end
    # Obtain the dictionary keys from the primary code and create Vector of
    # codes of the same dictionary type
    code_keys = collect(keys(codes[1]))
    dict_type = Dict{eltype(code_keys),Vector{eltype(codes[1][code_keys[1]])}}
    signal_codes = Vector{dict_type}(undef, code_num)
    code_lengths = Array{Int}(undef, code_num)
    # Loop through codes and convert non-dict items to dictionaries
    for i in 1:N
        # Check if item in vector is a dictionary.
        # If it is, it is placed into `signal_codes[i]`, otherwise, it is
        # assumed that the item in the `codes` vector is a single code
        # that is applied to all PRNs. It is then copied into a dictionary
        # and associated with all dictionary keys, the PRN numbers.
        if isa(codes[i], Dict)
            signal_codes[i] = codes[i]
            code_lengths[i] = length(codes[i][collect(keys(codes[i]))[1]])
        else
            code = dict_type()
            code_lengths[i] = length(codes[i])
            for code_key in code_keys
                code[code_key] = codes[i]
            end
            signal_codes[i] = code
        end
    end
    # Repeat the above process for databits, if they are given. The same
    # treatment above is used for the databits.
    if ~ismissing(databits)
        push!(chipping_rates, databits[2])
        included_databits = true
        databit_codes = databits[1]
        if isa(databit_codes, Dict)
            similar_databits = false
            signal_codes[code_num] = databit_codes
            code_lengths[code_num] = length(databit_codes[collect(keys(databit_codes))[1]])
        else
            similar_databits = true
            code = dict_type()
            code_lengths[code_num] = length(databit_codes)
            for code_key in code_keys
                code[code_key] = databit_codes
            end
            signal_codes[code_num] = code
        end
    else
        included_databits = false
        similar_databits = false
    end
    include_codes = fill(true, code_num)
    return CodeType(name, code_num, signal_codes, chipping_rates, code_lengths,
                    include_codes, included_databits, similar_databits)
end

"""
    definecodetype(code::Dict, chipping_rate::Float64;
                   databits=missing, name="custom")

Defines a `CodeType` struct that represents the codes of a signal on a
given channel, `I`, `Q`, or `both`.

Required Arguments:

- `code::Dict`: dictionary of primary codes where the keys are the PRN numbers
- `chipping_rate::Float`: chipping rate in Hz for the primary code above

Optional Arguments:

- `databits::Vector`: if specified by user, must be a two element
                      array where the first element is the databits
                      as either a dictionary or array and the second
                      element is the chipping rate in Hz of the
                      databits `(default = missing)`
- `name::String`: name of code type `(default is "custom")`

Returns:

- `CodeType` structure
"""
function definecodetype(code::Dict, chipping_rate::Float64;
                        databits=missing, name="custom")
    # Check that the first code in `codes`, the primary code, is a Dict
    if ~isa(code, Dict)
        error("Code given must be in the form of a dictionary.")
    end
    # Obtain the dictionary keys from the primary code and create Vector of
    # codes of the same dictionary type
    include_code = true
    code_keys = collect(keys(code))
    dict_type = Dict{eltype(code_keys),Vector{eltype(code[code_keys[1]])}}
    code_length = length(code[code_keys[1]])
    # Check if databit codes given is a dictionary or array.
    # If it is, it is placed into `signal_codes[2]`, otherwise, it is
    # assumed that is a single that is applied to all PRNs. It is then copied
    # into a dictionary and associated with all dictionary keys, the PRN numbers.
    if ~ismissing(databits)
        code_num = 2
        signal_codes = Vector{dict_type}(undef, code_num)
        signal_codes[1] = code
        chipping_rate = [chipping_rate, databits[2]]
        included_databits = true
        databit_codes = databits[1]
        if isa(databit_codes, Dict)
            similar_databits = false
            signal_codes[code_num] = databit_codes
            code_length = [code_length, length(databit_codes[code_keys[1]])]
        else
            similar_databits = true
            code = dict_type()
            code_length = [code_length, length(databit_codes)]
            for code_key in code_keys
                code[code_key] = databit_codes
            end
            signal_codes[code_num] = code
        end
    else
        code_num = 1
        included_databits = false
        signal_codes = [code]
        similar_databits = false
    end
    include_codes = fill(true, code_num)
    return CodeType(name, code_num, signal_codes, chipping_rate, code_length,
                    include_codes, included_databits, similar_databits)
end

"""
    definesignaltype(I_codes::CodeType, Q_codes::CodeType, sig_freq;
                     name="custom")

Defines a signal type with `CodeType` structs defined for both the `I` and `Q`
channels.

Required Arguments:

- `I_codes::CodeType`: Custum codes for the I channel of the signal
- `Q_codes::CodeType`: Custum codes for the Q channel of the signal
- `sig_freq::Float64`: carrier frequency of signal in Hz

Optional Arguments:

- `name::String`: name of signal type `(default is "custom")`
- `B`: receiver bandwidth `(default is based off max B in I_codes and Q_codes)`

Returns:

- `SignalType` struct
"""
function definesignaltype(I_codes::CodeType, Q_codes::CodeType, sig_freq;
                          name="custom", B=missing)
    if ismissing(B)
        # Determine maximum bandwidth for I channel codes
        B_I = 2*maximum(I_codes.chipping_rates)
        # Determine maximum bandwidth for Q channel codes
        B_Q = 2*maximum(Q_codes.chipping_rates)
    else
        B_I = B
        B_Q = B
    end
    return SignalType(name, I_codes, Q_codes, sig_freq, B_I, B_Q, true, true)
end

"""
    definesignaltype(codes::CodeType, sig_freq, channel="I"; name="custom")

Defines a signal type with a `CodeType` struct defined for either the `I`, `Q`,
or `both` channels.

Required Arguments:

- `codes::CodeType`: Custum codes for either the `I`, `Q`, or `both` channels
- `sig_freq::Float64`: carrier frequency of signal in Hz

Optional Arguments:

- `channel::String`: set to either `I`, `Q`, or `both` (default is `I`)
- `name::String`: name of signal type (default is `custom`)
- `B`: receiver bandwidth `(default = 2*maximum(codes.chipping_rates))`

Returns:

- `SignalType` struct
"""
function definesignaltype(codes::CodeType, sig_freq, channel="I"; name="custom",
                          B=missing)
    # Determine maximum bandwidth for codes
    if ismissing(B)
        B = 2*maximum(codes.chipping_rates)
    end
    if channel == "both"
        return SignalType(name, codes, codes, sig_freq, B, B, true, true)
    elseif channel == "I"
        return SignalType(name, codes, missing, sig_freq, B, missing, true, false)
    elseif channel == "Q"
        return SignalType(name, missing, codes, sig_freq, missing, B, false, true)
    else
        error("Invalid channel specified.")
    end
end

"""
    GNSSSignal

The most general type of signal type used in `GNSSTools`. Used for both
`GNSSData` and `ReplicaSignal` structs.
"""
abstract type GNSSSignal end

"""
    Data

Structure for holding GNSS signal data.

Fields:

- `file_name::String`: name of data file
- `f_s::Float64`: sampling rate of data in Hz
- `f_if::Float64`: IF frequency of data in Hz
- `t_length::Float64`: length of data loaded in seconds
- `start_data_idx::Int`: starting index in the data file
- `data::Array{Complex{Float64},1}`: loaded data from data file
- `data_type::String`: a `Symbol` type either `:sc4` or `:sc8` for 4 and 8 bit
                       comeplex data, respectively
- `data_start_time::T1`: `Tuple` of length 6
    * format is `(year, month, day, hour, minute, second)`
    * everything other than `second` must be an integer
    * `second` can be an integer or float
    * set to `missing` if no start time given
- `site_loc_lla::T2`: data collection site location
    * format is `(latitude, longitude, height)`
    * `latitude` and `longitude` are in degrees
    * `height` is in meters
    * set to `missing` if no site location given
- `sample_num::Int`: number of data samples in `Data.data`
- `total_data_length::Float64`: total length of data file in seconds
- `nADC::Int`: bit depth of data
    * `sc4` is 4 bit
    * `sc8` is 8 bit
- `start_t::Float64`: the start time of the signal in seconds
"""
struct GNSSData{T1,T2,T3} <: GNSSSignal
    file_name::String
    f_s::Float64
    f_if::Float64
    t_length::Float64
    start_data_idx::Int
    data::Array{T3,1}
    data_type::String
    data_start_time::T1
    site_loc_lla::T2
    sample_num::Int
    total_data_length::Float64
    nADC::Int
    start_t::Float64
end

"""
    GNSSData(file_name, f_s, f_if, t_length, start_data_idx, t,
             data::Array{T3,1}, data_type, data_start_time::T1,
             site_loc_lla::T2, sample_num, total_data_length,
             nADC, start_t) where {T1, T2, T3}

Outer constructor function for `GNSSData` struct initialization.

Returns:

- `GNSSData` struct
"""
function GNSSData(file_name, f_s, f_if, t_length, start_data_idx,
                  data::Array{T3,1}, data_type, data_start_time::T1,
                  site_loc_lla::T2, sample_num, total_data_length,
                  nADC, start_t) where {T1, T2, T3}
    return GNSSData{T1,T2,T3}(file_name, f_s, f_if, t_length, start_data_idx,
                              data, data_type, data_start_time, site_loc_lla,
                              sample_num, total_data_length, nADC, start_t)
end

"""
    ReplicaSignal

A abstract struct for the replica signal structs.

Fields:

- `name::String`: name of replica signal
- `prn::Int`: code PRN to be generated
- `f_s::Float64`: receiver sampling rate of signal in Hz
- `t_length::Float64`: data length in seconds
- `f_if::Float64`: IF frequency in Hz
- `f_d::Float64`: initial carrier Doppler frequency in Hz
- `fd_rate::Float64`: initial carrier Doppler frequency rate in Hz
- `Tsys::Float64`: receiver noise temperature in Kelvin
- `CN0::Float64`: signal carrier to noise ratio (C/N₀)
- `phi::Float64`: initial carrier phase in radians
- `nADC::Int64`: bit depth of receiver
- `code_start_idx::Float64`: index in `ReplicaSignal.data` where all codes in
                             signal start
- `init_code_phases_I::T1`: array of initial code phases for all I channel codes
    * vector length is set to 0 if there are no I channel codes
- `init_code_phases_Q::T2`: array of initial code phases for all Q channel codes
    * vector length is set to 0 if there are no Q channel codes
- `data::Array{Complex{Float64},1}`: where the signal raw I/Q samples are stored
- `include_carrier::Bool`: flag for if signal will be modulated onto a carrier
- `include_adc::Bool`: flag for if signal will undergo ADC quantization
- `include_thermal_noise::Bool`: flag for if thermal noise will be added
- `include_databits_I::Bool`: flag for if databits on the I channel will be put
                              onto signal
- `include_databits_Q::Bool`: flag for if databits on the Q channel will be put
                              onto signal
- `include_phase_noise::Bool`: flag for if phase noise will be added onto signal
- `f_code_d_I::T2`: adjusted chipping rate due to Doppler frequency for I
                    channel codes
    * vector length is set to 0 if there are no I channel codes
- `f_code_dd_I::T3`: adjusted chipping rate rate due to Doppler frequency rate
                     for I channel codes
    * vector length is set to 0 if there are no I channel codes
- `f_code_d_Q::T4`: adjusted chipping rate due to Doppler frequency for Q
                    channel codes
    * vector length is set to 0 if there are no Q channel codes
- `f_code_dd_Q::T5`: adjusted chipping rate rate due to Doppler frequency rate
                     for Q channel codes
    * vector length is set to 0 if there are no Q channel codes
- `sample_num::Int`: number of raw I/Q samples in `ReplicaSignal.data`
- `isreplica::Bool`: flag used to specify that a given `ReplicaSignal` struct
                     will be used as a replica for signal processing
    * if `true`, a `ReplicaSignal` struct will not have noise added onto it
      and will not undergo ADC quantization
    * signal generation will be done using the second method of `generatesignal!`,
      `generatesignal!(signal::ReplicaSignal, isreplica::Bool)`
- `noexp::Bool`: used only if second method of `generatesignal!`, discussed
                 above, is used
    * does not modulate codes onto carrier
- `thermal_noise::Array{Complex{Float64},1}`: contains pre-build thermal noise
                                              vector
    * if `allocate_noise_vectors` is set to `for variable

    endalse` in `definesignal` method,
      the size of this vector is 0, indicating that it will not be used
- `phase_noise::Array{Float64,1}`: contains pre-build phase noise vector
    * if `allocate_noise_vectors` is set to `false` in `definesignal` method,
      the size of this vector is 0, indicating that it will not be used
- `signal_type::T7`: a `SignalType` struct that was used to define the signal
- `receiver_h_parms::Array{Float64,1}`:
- `start_t::Float64`: the start time of the signal in seconds

"""
mutable struct ReplicaSignal{T} <: GNSSSignal
    name::String
    prn::Int
    f_s::Float64
    t_length::Float64
    f_if::Float64
    f_d::Float64
    fd_rate::Float64
    Tsys::Float64
    CN0::Float64
    phi::Float64
    nADC::Int
    code_start_idx::Rational{Int}
    init_code_phases_I::Array{Float64,1}
    init_code_phases_Q::Array{Float64,1}
    data::Array{Complex{Float64},1}
    include_carrier::Bool
    include_adc::Bool
    include_thermal_noise::Bool
    include_databits_I::Bool
    include_databits_Q::Bool
    include_phase_noise::Bool
    f_code_d_I::Array{Float64,1}
    f_code_dd_I::Array{Float64,1}
    f_code_d_Q::Array{Float64,1}
    f_code_dd_Q::Array{Float64,1}
    sample_num::Int
    isreplica::Bool
    noexp::Bool
    thermal_noise::Array{Complex{Float64},1}
    phase_noise::Array{Complex{Float64},1}
    signal_type::T
    receiver_h_parms::Array{Float64,1}
    start_t::Float64
end

"""
    ReplicaSignal(name, prn, f_s, t_length, f_if, f_d, fd_rate, Tsys,
                  CN0, phi, nADC, code_start_idx, init_code_phases_I,
                  init_code_phases_Q, t, data, include_carrier,
                  include_adc, include_thermal_noise, include_databits_I,
                  include_databits_Q, include_phase_noise, f_code_d_I,
                  f_code_dd_I, f_code_d_Q, f_code_dd_Q, sample_num,
                  isreplica, noexp, thermal_noise, phase_noise,
                  signal_type::T, receiver_h_parms, start_t) where {T}

Outer constructor function for `ReplicaSignal` struct initialization.

Returns:

- `ReplicaSignal` struct
"""
function ReplicaSignal(name, prn, f_s, t_length, f_if, f_d, fd_rate, Tsys,
                       CN0, phi, nADC, code_start_idx, init_code_phases_I,
                       init_code_phases_Q, data, include_carrier,
                       include_adc, include_thermal_noise, include_databits_I,
                       include_databits_Q, include_phase_noise, f_code_d_I,
                       f_code_dd_I, f_code_d_Q, f_code_dd_Q, sample_num,
                       isreplica, noexp, thermal_noise, phase_noise,
                       signal_type::T, receiver_h_parms, start_t) where {T}
    return ReplicaSignal{T}(name, prn, f_s, t_length, f_if, f_d, fd_rate, Tsys,
                            CN0, phi, nADC, code_start_idx, init_code_phases_I,
                            init_code_phases_Q, data, include_carrier,
                            include_adc, include_thermal_noise, include_databits_I,
                            include_databits_Q, include_phase_noise, f_code_d_I,
                            f_code_dd_I, f_code_d_Q, f_code_dd_Q, sample_num,
                            isreplica, noexp, thermal_noise, phase_noise,
                            signal_type, receiver_h_parms, start_t)
end

"""
    ReplicaSignals

Struct containing array of `ReplicaSignal` structs to generate a signal that
contains signals from multiple PRNs.

Fields:

- `name::String`:
- `replica_signals::Array{ReplicaSignal,1}`:
- `data::Array{Complex{Float64},1}`:
- `t::Array{Float64,1}`:
- `f_s::Float64`:
- `f_if::Float64`:
- `sample_num::Int`:
- `nADC::Int`:
- `thermal_noise::Array{Complex{Float64},1}`:
- `phase_noise::Array{Float64,1}`:
- `start_t::Float64`: the start time of the signal in seconds
"""
mutable struct ReplicaSignals <: GNSSSignal
    name::String
    replica_signals::Array{ReplicaSignal,1}
    data::Array{Complex{Float64},1}
    t_length::Float64
    f_s::Float64
    f_if::Float64
    B::Float64
    Tsys::Float64
    sample_num::Int
    nADC::Int
    include_carrier::Bool
    include_adc::Bool
    include_thermal_noise::Bool
    include_phase_noise::Bool
    thermal_noise::Array{Complex{Float64},1}
    phase_noise::Array{Complex{Float64},1}
    receiver_h_parms::Array{Float64,1}
    start_t::Float64
end


text:
mutable struct SortingDirectMethod{S<:RatesCacheTrait} <: DecoupledSearchAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
  search_order::Vector{Int}
  search_index::Int
end

RatesCache(::SortingDirectMethod{S}) where S = S()

@inline function initialize!(algorithm::SortingDirectMethod, state, model)
  search_order = algorithm.search_order

  for i in eachindex(search_order)
      @inbounds search_order[i] = i
  end

  algorithm.search_index = 1

  update_jump_rates!(algorithm, state, model)

  return nothing
end

@inline function search_jumps(algorithm::SortingDirectMethod)
  # extract fields
  rates = algorithm.rates
  total_rate = algorithm.total_rate
  search_order = algorithm.search_order

  # define an iterator that searches `rates` along the give order
  itr = (rates[search_order[i]] for i in eachindex(search_order))

  # search the rates using our iterator
  idx = search_jump_rates(RatesCache(algorithm), itr, total_rate)

  # record the index into the search order vector
  algorithm.search_index = idx

  # return the reaction index
  return @inbounds search_order[idx]
end

@inline function update!(algorithm::SortingDirectMethod, state, model, j)
  # update dependent propensities
  update_jump_rates!(algorithm, state, model, j)

  search_order = algorithm.search_order
  search_index = algorithm.search_index

  # perform a sorting step
  if search_index > 1
      @inbounds tmp = search_order[search_index]

      @inbounds search_order[search_index]     = search_order[search_index - 1]
      @inbounds search_order[search_index - 1] = tmp
  end

  return nothing
end

ext: jl


text:
module OshiroBoss

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Oshiro Boss" => Ahorn.EntityPlacement(
        Maple.FriendlyGhost
    )
)

sprite = "characters/oshiro/boss13.png"

function Ahorn.selection(entity::Maple.FriendlyGhost)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.FriendlyGhost, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end


text:
using BenchmarkTools, BioMotifInference, BioBackgroundModels, BioSequences
println(Threads.nthreads())
source_pwm = [.7 .1 .1 .1
.1 .1 .1 .7
.1 .1 .7 .1]
source_pwm_2 = [.6 .1 .1 .2
.2 .1 .1 .6
.1 .2 .6 .1]
sources = [(log.(source_pwm), 1),(log.(source_pwm_2), 1)]

bg_scores = log.(fill(.1, (150,10000)))

obs=zeros(UInt8,151,10000)
obs[1:150,1:end].=1
obsl=[findfirst(iszero,obs[:,o])-1 for o in 1:size(obs)[2]]
mix=trues(10000,2)

println(median(@benchmark (BioMotifInference.IPM_likelihood($sources, $obs, $obsl, $bg_scores, $mix)) samples=5 evals=20))

println(median(@benchmark (BioMotifInference.dev_likelihood($sources, $obs, $obsl, $bg_scores, $BitMatrix(transpose(mix)))) samples=5 evals=20))


text:
@testset "Operations" begin
    for p in permutations(primitives_array, 2)
        for op in operations_array
            op(p...)
        end
    end
end


text:
module CommunalHelperDreamFallingBlock

using ..Ahorn, Maple
using Ahorn.Cairo
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/DreamFallingBlock" DreamFallingBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    featherMode::Bool=false,
    oneUse::Bool=false,
    refillCount::Integer=-1,
    noCollide::Bool=false,
    below::Bool=false,
    fallDistance::Integer=64,
    centeredChain::Bool=false,
    chainOutline::Bool=true,
    indicator::Bool=false,
    indicatorAtStart::Bool=false,
    chained::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Dream Falling Block (Communal Helper)" => Ahorn.EntityPlacement(
        DreamFallingBlock,
        "rectangle",
    ),
    "Dream Falling Block (Chained) (Communal Helper)" => Ahorn.EntityPlacement(
        DreamFallingBlock,
        "rectangle",
        Dict{String, Any}(
            "chained" => true,
        ),
    ),
)

Ahorn.minimumSize(entity::DreamFallingBlock) = 8, 8
Ahorn.resizable(entity::DreamFallingBlock) = true, true

Ahorn.selection(entity::DreamFallingBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.editingIgnored(entity::DreamFallingBlock, multiple::Bool=false) =
    multiple ?
    String["x", "y", "width", "height", "chained"] : (Bool(get(entity.data, "chained", false)) ?
    String["chained"] : String["chained", "fallDistance", "centeredChain", "chainOutline", "indicator", "indicatorAtStart"])

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamFallingBlock)
    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if Bool(get(entity.data, "chained", false))
        save(ctx)

        set_dash(ctx, [0.6, 0.2])
        set_antialias(ctx, 1)
        set_line_width(ctx, 1)

        # width & height extend up and left for some reason.
        Ahorn.drawRectangle(ctx, 0, 0, width, height + Integer(get(entity.data, "fallDistance", 64)), (0.0, 0.0, 0.0, 0.0), (1.0, 1.0, 1.0, 0.5))

        restore(ctx)
    end

    renderDreamBlock(ctx, 0, 0, width, height, entity.data)
end

end

max_stars_repo_name: vitellaryjr/CommunalHelper

max_issues_repo_name: vitellaryjr/CommunalHelper

max_forks_repo_name: vitellaryjr/CommunalHelper


text:
# Assertion (B.19)

export @cuassert

"""
    @assert cond [text]

Signal assertion failure to the CUDA driver if `cond` is `false`. Preferred syntax for
writing assertions, mimicking `Base.@assert`. Message `text` is optionally displayed upon
assertion failure.

!!! warning
    A failed assertion will crash the GPU, so use sparingly as a debugging tool.
    Furthermore, the assertion might be disabled at various optimization levels, and thus
    should not cause any side-effects.
"""
macro cuassert(ex, msgs...)
    # message handling copied from Base.@assert
    msg = isempty(msgs) ? ex : msgs[1]
    if isa(msg, AbstractString)
        msg = msg # pass-through
    elseif !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        # message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    elseif isdefined(Main, :Base) && isdefined(Main.Base, :string) && applicable(Main.Base.string, msg)
        msg = Main.Base.string(msg)
    else
        # string() might not be defined during bootstrap
        msg = :(Main.Base.string($(Expr(:quote,msg))))
    end

    return :($(esc(ex)) ? $(nothing)
                        : cuassert_fail($(Val(Symbol(msg))),
                                        $(Val(__source__.file)),
                                        $(Val(__source__.line))))
end

assert_counter = 0

@generated function cuassert_fail(::Val{msg}, ::Val{file}, ::Val{line}) where
                                 {msg, file, line}
    JuliaContext() do ctx
        T_void = LLVM.VoidType(ctx)
        T_int32 = LLVM.Int32Type(ctx)
        T_pint8 = LLVM.PointerType(LLVM.Int8Type(ctx))

        # create function
        llvm_f, _ = create_function(T_void)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            global assert_counter
            assert_counter += 1

            message = globalstring_ptr!(builder, String(msg), "assert_message_$(assert_counter)")
            file = globalstring_ptr!(builder, String(file), "assert_file_$(assert_counter)")
            line = ConstantInt(T_int32, line)
            func = globalstring_ptr!(builder, "unknown", "assert_function_$(assert_counter)")
            charSize = ConstantInt(Csize_t(1), ctx)

            # invoke __assertfail and return
            # NOTE: we don't mark noreturn since that control flow might confuse ptxas
            assertfail_typ =
                LLVM.FunctionType(T_void,
                                [T_pint8, T_pint8, T_int32, T_pint8, llvmtype(charSize)])
            assertfail = LLVM.Function(mod, "__assertfail", assertfail_typ)
            call!(builder, assertfail, [message, file, line, func, charSize])

            ret!(builder)
        end

        call_function(llvm_f, Nothing, Tuple{})
    end
end


text:
module TestSetindex
using Accessors
using Test

"""
    ==ₜ(x, y)

Check that _type_ and value of `x` and `y` are equal.
"""
==ₜ(_, _) = false
==ₜ(x::T, y::T) where T = x == y

@testset "==ₜ" begin
    @test 1 ==ₜ 1
    @test !(1.0 ==ₜ 1)
end

function ref_alloc_test()
    ref = Ref((UInt(10), 100))
    ref2 = @set ref[][1] *= -1

    ref2[]
end

@testset "setindex" begin
    arr = [1,2,3]
    @test_throws MethodError Base.setindex(arr, 10, 1)
    @test Accessors.setindex(arr, 10, 1) == [10, 2, 3]
    @test arr == [1,2,3]
    @test @set(arr[1] = 10) == [10, 2, 3]
    @test arr == [1,2,3]
    @test Accessors.setindex(arr, 10.0, 1) ==ₜ Float64[10.0, 2.0, 3.0]

    d = Dict(:a => 1, :b => 2)
    @test_throws MethodError Base.setindex(d, 10, :a)
    @test Accessors.setindex(d, 10, :a) == Dict(:a=>10, :b=>2)
    @test d == Dict(:a => 1, :b => 2)
    @test @set(d[:a] = 10) == Dict(:a=>10, :b=>2)
    @test d == Dict(:a => 1, :b => 2)
    @test Accessors.setindex(d, 30, "c") ==ₜ Dict(:a=>1, :b=>2, "c"=>30)
    @test Accessors.setindex(d, 10.0, :a) ==ₜ Dict(:a=>10.0, :b=>2.0)

    ref = Ref((; a = 1, b = 2, c = (; aa = 3)))
    @test @set(ref[].a = 90)[] == (; a = 90, b = 2, c = (; aa = 3))
    @test @set(ref[].b = "2")[] ==ₜ (; a = 1, b = "2", c = (; aa = 3))
    @test @set(ref[].c.aa += 2)[] == (; a = 1, b = 2, c = (; aa = 5))

    ref = Ref(1::Int)
    @set ref[] = "no mutation"
    @test ref[] === 1
    @test typeof(ref) == Base.RefValue{Int}

    if VERSION >= v"1.5.0"
        _ = ref_alloc_test()
        @test @allocated(ref_alloc_test()) == 0
    end
end

end


text:
__precompile__()

module GenomicFeatures

export
    Strand,
    STRAND_NA,
    STRAND_POS,
    STRAND_NEG,
    STRAND_BOTH,

    Interval,
    seqname,
    leftposition,
    rightposition,
    strand,
    metadata,
    isoverlapping,

    IntervalCollection,
    eachoverlap,
    coverage,

    BED,
    GFF3,
    BigWig,
    BigBed,
    isfilled,
    hasseqname,
    hasleftposition,
    hasrightposition

import BGZFStreams
import BioCore: BioCore, seqname, leftposition, rightposition, isoverlapping, isfilled, hasseqname, hasleftposition, hasrightposition, metadata
import DataStructures
import IntervalTrees
import Base.@propagate_inbounds

include("strand.jl")
include("interval.jl")
include("intervalcollection.jl")
include("queue.jl")
include("overlap.jl")
include("coverage.jl")
include("indexes/indexes.jl")
include("bed/bed.jl")
include("gff3/gff3.jl")
include("bbi/bbi.jl")
include("bigwig/bigwig.jl")
include("bigbed/bigbed.jl")

end # module


text:
Dict(
     "T5" => "T5 - en_Latn_US",
    )


text:
abstract type AbstractFilteringAlgorithm{T1<:TimeType, T2<:TimeType} end
abstract type AbstractFilterState end

"""
    update!(filter_state, filter_algo, obs, dt) --> filter_state

Updates the filter state by performing one forward step of the model and then assimilating the observation.
"""
function update!(
    filter_state::AbstractFilterState,
    filter_algo::AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime},
    obs,
    dt
) end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# ========================================
# Projection operations
# ========================================

function _project_oa(R::AbstractLazyReachSet, vars, t; remove_zero_generators=true)
    return Project(R, vars)
end

function _project_oa(R::AbstractTaylorModelReachSet, vars, t; remove_zero_generators=true)
    Z = overapproximate(R, Zonotope, t, remove_zero_generators=remove_zero_generators)
    return project(set(Z), vars; remove_zero_generators=remove_zero_generators)
end

# ========================================
# Decomposition operations
# ========================================

# decompose a set into the Cartesian product of intervals
using LazySets.Arrays: SingleEntryVector

function _decompose_1D(X0::LazySet{N}) where {N}
    n = dim(X0)
    out = Vector{Interval{Float64, IA.Interval{Float64}}}(undef, n)

    @inbounds for i in 1:n
        eᵢ = SingleEntryVector(i, n, one(N))
        out[i] = Interval(-ρ(-eᵢ, X0), ρ(eᵢ, X0))
    end
    return CartesianProductArray(out)
end

# ==================
# Overapproximation
# ==================

function LazySets.overapproximate(X::UnionSetArray, ::Type{Interval})
    return overapproximate(ConvexHullArray(array(X)), Interval)
end

# ==============================================================
# Construct the initial states for the continuous post-operator
# ==============================================================

using ReachabilityAnalysis: TaylorModel1, TaylorModelN, fp_rpa, zeroBox, symBox

abstract type AbstractReconstructionMethod end

struct CartesianProductReconstructor <: AbstractReconstructionMethod end

function _reconstruct(method::CartesianProductReconstructor, P₀::LazySet, U₀::LazySet, R, ti)
    Q₀ = P₀ × U₀
    return Q₀
end

@with_kw struct TaylorModelReconstructor <: AbstractReconstructionMethod
    box::Bool = true
end

# if no Taylor model is available => use the given set P₀
function _reconstruct(method::TaylorModelReconstructor, P₀::LazySet, U₀, R::Nothing, ti) where {N}
    return _reconstruct(CartesianProductReconstructor(), P₀, U₀, R, ti)
end

function _reconstruct(method::TaylorModelReconstructor, P₀::LazySet, U₀::Vector{<:LazySet}, R::TaylorModelReachSet{N}, ti) where {N}
    @assert length(U₀) == 1 "expected the length of U₀ to be 1, got $(lenght(U₀))"
    return _reconstruct(method, P₀, first(U₀), R, ti)
end

function _reconstruct(method::TaylorModelReconstructor, P₀::LazySet, U₀::LazySet, R::TaylorModelReachSet{N}, ti) where {N}
    # evaluate R at the final time of the reach-set
    S = set(R)
    tn = sup(domain(R)) # assume that the reach set spans the whole period (checked outside this method)
    X_Δt = evaluate(S, tn)

    n = dim(P₀) # number of state variables = dim(P₀) - dim(U₀)
    m = dim(U₀) # control variables

    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n + m)

    # construct state variables
    orderT = get_order(first(S))
    orderQ = get_order(X_Δt[1])

    zeroI = interval(zero(N), zero(N))
    Δtn = zeroI
    @inbounds for i in 1:n
        W = TaylorModelN(X_Δt[i], zeroI, zeroBox(n + m), symBox(n + m))
        Ŵ = fp_rpa(W)
        p = Taylor1(TaylorN(polynomial(Ŵ)), orderT)
        rem = remainder(Ŵ)
        vTM[i] = TaylorModel1(p, rem, zeroI, Δtn)
    end

    # fill the components for the inputs
    if method.box || isa(U₀, AbstractHyperrectangle)
        B₀ = convert(IntervalBox, box_approximation(U₀))
        @inbounds for i in 1:m
            I = B₀[i]
            p = mid(I) + zero(TaylorN(n+m, order=orderQ))
            d = diam(I) / 2
            rem = interval(-d, d)
            vTM[n+i] = TaylorModel1(Taylor1(p, orderT), rem, zeroI, Δtn)
        end
    else
        Z₀ = ReachabilityAnalysis._convert_or_overapproximate(U₀, Zonotope)
        Z₀ = ReachabilityAnalysis._reduce_order(Z₀, 2, force_reduction=true)
        # NOTE if we used _overapproximate_structured directly :
        # Utm₀ = set(ReachabilityAnalysis._overapproximate_structured(Z₀, TaylorModelReachSet, orderT=orderT, orderQ=orderQ))
        # @inbounds for i in 1:m
        #     vTM[n+i] = Utm₀[i]
        # end

        x = set_variables("x", numvars=n+m, order=orderQ)
        xc = view(x, n+1:n+m)
        G = Z₀.generators
        c = Z₀.center
        @assert size(G) == (m, 2m)
        M = view(G, :, 1:m)
        D = view(G, :, m+1:2m)
        @assert isdiag(D)
        @inbounds for i in 1:m
            p = c[i] + sum(view(M, i, :) .* xc)
            d = abs(D[i, i])
            rem = interval(-d, d)
            vTM[n+i] = TaylorModel1(Taylor1(p, orderT), rem, zeroI, Δtn)
        end
    end
    return TaylorModelReachSet(vTM, Δtn)
end


text:
using PyPlot
using ADSeismic
using ADCME
using JLD2

mpi_init()

M = N = 5
n = 100
@info M
param = MPIAcousticPropagatorParams(NX = n*M, NY = n*N, n = n,
     Rcoef=0.2, DELTAX=10, DELTAY=10, DELTAT=0.05, NSTEP = 2000)
compute_PML_Params!(param; check_mpi_size = false)

Vals = Array{Float64, 3}[]
for i = 1:M
    for j = 1:N
        r = (i-1)*N + j - 1
        @load "data/dat$(r)-$(M*N).jld2" u
        push!(Vals, u)
    end
end

p = visualize_wavefield(Vals, param)
saveanim(p, "wavefield.gif")


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

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON <: SwaggerModel

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

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

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON)
    true
end

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


text:
module TestUtils
using AugmentedGPLikelihoods
const AGPL = AugmentedGPLikelihoods
using AugmentedGPLikelihoods.SpecialDistributions
using Distributions
using GPLikelihoods: AbstractLikelihood
using LinearAlgebra
using MeasureBase
using MeasureTheory: For
using Random
using Test
using TupleVectors
# Test API for augmented likelihood
remove_ntdist_wrapper(d::NTDist) = d.d
remove_ntdist_wrapper(d) = d

function test_auglik(
    lik::AbstractLikelihood;
    n=10,
    f=randn(n),
    qf=Normal.(randn(n), 1.0),
    rng::AbstractRNG=Random.GLOBAL_RNG,
)
    y = rand.(rng, lik.(f))
    nf = nlatent(lik)
    # Testing sampling
    @testset "Sampling" begin
        Ω = init_aux_variables(lik, n)
        @test Ω isa TupleVector
        @test first(Ω) isa NamedTuple
        @test length(Ω) == n
        Ω = aux_sample!(rng, Ω, lik, y, f)
        @test Ω isa TupleVector
        new_Ω = aux_sample(rng, lik, y, f)
        @test new_Ω isa TupleVector
        @test length(Ω) == n

        βs = auglik_potential(lik, Ω, y)
        γs = auglik_precision(lik, Ω, y)
        β2, γ2 = auglik_potential_and_precision(lik, Ω, y)
        @test length(γs) == length(βs) == nf # Check that there are n latent vectors
        @test first(βs) isa AbstractVector
        @test first(γs) isa AbstractVector
        @test all(map(≈, βs, β2))
        @test all(map(≈, γs, γ2))
        @test all(x -> all(>=(0), x), γs) # Check that the variance is positive

        sumlogtilt = logtilt(lik, Ω, y, f)
        @test sumlogtilt isa Real
        @test logtilt(lik, AGPL.aux_field(lik, first(Ω)), first(y), first(f)) isa Real
        sumlogtilt_alt = mapreduce(+, AGPL.aux_field(lik, Ω), y, f) do ωᵢ, yᵢ, fᵢ
            logtilt(lik, ωᵢ, yᵢ, fᵢ)
        end
        @test sumlogtilt ≈ sumlogtilt_alt

        @test aug_loglik(lik, Ω, y, f) isa Real

        pΩ = aux_prior(lik, y)
        @test logdensity_def(pΩ, Ω) isa Real
        pω = aux_prior(lik, first(y)) # Scalar version
        @test pω == remove_ntdist_wrapper(first(marginals(pΩ)))

        # Test that the full conditional is correct
        @testset "Full conditional Ω" begin
            pcondΩ = aux_full_conditional(lik, y, f) # Compute the full conditional of Ω
            Ω₁ = tvrand(rng, pcondΩ) # Sample a set of aux. variables
            Ω₂ = tvrand(rng, pcondΩ) # Sample another set of aux. variables
            # We compute p(f, y) by doing C = p(f,y) = p(y|Ω,f)p(Ω)/p(Ω|y,f)
            # This should be the same no matter what Ω is
            logC₁ = logtilt(lik, Ω₁, y, f) + logdensity_def(pΩ, Ω₁) - logdensity_def(pcondΩ, Ω₁)
            logC₂ = logtilt(lik, Ω₂, y, f) + logdensity_def(pΩ, Ω₂) - logdensity_def(pcondΩ, Ω₂)
            @test logC₁ ≈ logC₂ atol = 1e-5
        end

        @testset "Full conditional f" begin
            pcondΩ = aux_full_conditional(lik, y, f) # Compute the full conditional of Ω
            Ω = tvrand(rng, pcondΩ) # Sample a set of aux. variables
            K = (x -> x * x')(rand(n, n)) # Prior Covariance matrix
            S = inv(Symmetric(inv(K) + Diagonal(only(auglik_precision(lik, Ω, y)))))
            m = S * (only(auglik_potential(lik, Ω, y)))
            qF = MvNormal(m, S)
            pF = MvNormal(K)
            f₁ = rand(rng, qF)
            f₂ = rand(rng, qF)
            logC₁ = logtilt(lik, Ω, y, f₁) + logpdf(pF, f₁) - logpdf(qF, f₁)
            logC₂ = logtilt(lik, Ω, y, f₂) + logpdf(pF, f₂) - logpdf(qF, f₂)
            @test logC₁ ≈ logC₂ atol = 1e-5
        end
    end

    #Testing variational inference
    @testset "Variational Inference" begin
        qΩ = init_aux_posterior(lik, n)
        @test qΩ isa For
        @test length(qΩ) == n
        qΩ = aux_posterior!(qΩ, lik, y, qf)
        @test qΩ isa For
        new_qΩ = aux_posterior(lik, y, qf)
        @test new_qΩ isa For
        @test length(new_qΩ) == n

        βs = expected_auglik_potential(lik, qΩ, y)
        γs = expected_auglik_precision(lik, qΩ, y)
        β2, γ2 = expected_auglik_potential_and_precision(lik, qΩ, y)
        @test length(γs) == length(βs) == nf # Check that there are n latent vectors
        @test first(βs) isa AbstractVector
        @test first(γs) isa AbstractVector
        @test all(map(≈, βs, β2))
        @test all(map(≈, γs, γ2))

        @test all(x -> all(>=(0), x), γs) # Check that the variance is positive

        # TODO test that aux_posterior parameters return the minimizing
        φ = TupleVectors.unwrap(only(aux_posterior(lik, y, qf).inds)) # TupleVector
        φ_opt = vcat(values(φ)...)
        s = keys(φ)
        n_var = length(s)
        function loss(φ)
            q = For(
                qΩ.f,
                TupleVector(
                    NamedTuple{s}(
                        collect(φ[((j - 1) * n_var + 1):(j * n_var)] for j in 1:n_var)
                    ),
                ),
            )
            return -expected_logtilt(lik, q, y, qf) + aux_kldivergence(lik, q, y)
        end
        ϵ = 1e-2
        # Test that by perturbing the value in random directions, the loss does not decrease
        for i in n_var * n
            (lik isa PoissonLikelihood && i <= n) && continue # We do not want to vary y
            Δ = zeros(n_var * n)
            Δ[i] = ϵ # We try one element at a time
            @test loss(φ_opt) <= loss(φ_opt + Δ)
            @test loss(φ_opt) <= loss(φ_opt - Δ)
        end
        # Optim.optimize(loss, φ_opt)
        # values of the ELBO
        pΩ = aux_prior(lik, y)
        @test pΩ isa For
        @test kldivergence(first(marginals(qΩ)), first(marginals(pΩ))) isa Real
        @test expected_logtilt(lik, qΩ, y, qf) isa Real
        @test aux_kldivergence(lik, qΩ, pΩ) isa Real
    end
end
end


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

function plot_draw_prior_distribution(df,N,beta,sigma,y_l,x_lab,n,m)
               Pkn_NGG_FK_ = smooth_pk(df.Pkn_FK[1:N])
               Pkn_NGG_SB_ = smooth_pk(df.Pkn_SB[1:N])
               R"p = ggplot(data.frame(k = 1:$N,
                                Pkn_5_numeric = $(df.Pkn_numeric[1:N]),
                                Pkn_4_order2 = $(df.Pkn_order2[1:N]),
                                Pkn_3_NGG_sM = $(df.Pkn_NGGM[1:N]),
                                Pkn_2_NGG_FK = $Pkn_NGG_FK_,
                                Pkn_1_NGG_SB = $Pkn_NGG_SB_
                            ) %>%
                    gather(Process_type, density, Pkn_5_numeric:Pkn_1_NGG_SB),
               aes(x=k, y = density, colour = Process_type)) +
                geom_line(aes(linetype =Process_type) ) + xlab($x_lab) + scale_linetype_manual(values=c('solid','solid','solid','solid','dashed')) +
                ylab('') + ggtitle(TeX(sprintf('$\\tau =%2.f$, $\\sigma = %.2f$,$\\n = %3.f$',$beta,$sigma,$n)))+
                scale_color_viridis(discrete=TRUE, direction = -1)+ theme_minimal()+ylim(0,$y_l)+xlim(1,$N)+scale_x_continuous(limits = c(1, $N), expand = c(0, 0),breaks= c(1,seq(0,$N,length=5)[2:5]))+
               theme(plot.title = element_text(hjust = 0.5,size = 10), axis.text.x = element_text(size=10), plot.margin = unit(c(1,$m, 0, 0),'pt'))"
    return R"p"
end

DF = load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_3.jld")
DF_100_1000 = DF["DF_100_1000"]

n=1000
β= 1.0
ntr=250
sigma_vec= [0.25,0.75]

N_plot = [40, 100]
y_l = [0.2,0.1]
P_all_approx_1000 =Array{RObject{VecSxp}}(undef,length(sigma_vec))
x_lab=[" ","k"]
for i in (1:length(sigma_vec))
             P_all_approx_1000[i]= plot_draw_prior_distribution(DF_100_1000[1][i],N_plot[i],β,sigma_vec[i],y_l[i],x_lab[i],n,0)

end

β= 10.0
N_plot = [40, 100]
y_l = [0.2,0.1]
P_all_approx_10_100 =Array{RObject{VecSxp}}(undef,length(sigma_vec))
x_lab=[" ","k"]
for i in (1:length(sigma_vec))
             P_all_approx_10_1000[i]= plot_draw_prior_distribution(DF_100_1000[2][i],N_plot[i],β,sigma_vec[i],y_l[i],x_lab[i],n,0)

end

R"
m1=as.list($P_all_approx_1000)
m2 = as.list($P_all_approx_10_1000)
prow <- plot_grid(
  m1[[1]] + theme(legend.position='none'),
  m2[[1]] + theme(legend.position='none'),
  m1[[2]] + theme(legend.position='none'),
  m2[[2]]+ theme(legend.position='none'),
  nrow = 2)
legend_b <- get_legend(m1[[1]]+theme(legend.position ='top'))
p <- plot_grid(prow,ncol = 1,rel_heights = c(10, 1))
#ggsave(file = '~/Documents/GitHub/GibbsTypePriors/test/comparison/Plots_sigma_all_approximation_100_n.pdf', width= 6, height = 4,p)
#save(m1,file ='~/Documents/GitHub/GibbsTypePriors/P_100_df_b1.Rdata')
#save(m2,file ='~/Documents/GitHub/GibbsTypePriors/P_100_df_b10.Rdata')
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_3.pdf')
plot(p)
dev.off()
"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
using Test, SegmentIntersections

@testset "Test results" begin
    segments = [
        Segment(4.9, 3.56, 2.86, 1.62),
        Segment(3.63, 3.88, 5.45, 2.32),
        Segment(2.99, 3.29, 4.47, 1.8),
        Segment(2.66, 2.85, 3.81, 1.44),
    ];
    expected_intersections = [Point(4.48, 3.16), Point(3.78, 2.5), Point(3.31, 2.05)];
    intersections = find_intersections_brute(segments)
    for expected in expected_intersections
        match = false
        for intersection in intersections
            if isapprox(intersection, expected, rtol = 1e-2)
                match = true
            end
        end
        @test match == true
    end
end


text:
using Base.Test
import WaterData
import WaterData.testhelpers: between, noerror

# quis custodiet ipsos custodes?
@testset "Test helper functions" begin
    @testset "`between`" begin
        @test between(3,5)(4)
        @test !between(7,9)(7)
        @test !between(2,3)(3)
        @test !between(-3,5)(-4)
        @test !between(-7,0)(2)
        @test_throws AssertionError between(2,2)
        @test_throws AssertionError between(3,-2)
    end

    @testset "`noerror`" begin
        @test noerror(1 + 1)
    end
end


text:
module QuantumTimeSteppers

using ApproxFun, FastTransforms, Plots, LinearAlgebra, FastGaussQuadrature
import ApproxFun: points, plan_transform, plan_itransform, TransformPlan, ITransformPlan, domain, canonicaldomain, spacescompatible, evaluate
import ApproxFun: Multiplication, ConcreteMultiplication, Derivative, ConcreteDerivative, space, rangespace, bandwidths
import Base: *, first, last, getindex
export HermiteFSE
export Lie_Trotter_evolve, Strang_evolve, Zassenhaus2_evolve

include("splitting_methods.jl")
include("HermiteFSE.jl")

end


text:
function bfgs_trace!(tr::OptimizationTrace,
                     x_old::Vector,
                     f_x::Real,
                     iteration::Integer,
                     gradient_new::Vector,
                     B::Matrix,
                     store_trace::Bool,
                     show_trace::Bool)
    dt = Dict()
    dt["g(x_new)"] = copy(gradient_new)
    dt["~inv(H)"] = copy(B)
    os = OptimizationState(copy(x_old), f_x, iteration, dt)
    if store_trace
        push!(tr, os)
    end
    if show_trace
        println(os)
    end
end

function bfgs(d::DifferentiableFunction,
              initial_x::Vector;
              initial_B::Matrix = eye(length(initial_x)),
              tolerance::Real = 1e-8,
              iterations::Integer = 1_000,
              store_trace::Bool = false,
              show_trace::Bool = false)

    # Keep track of the number of iterations
    iteration = 0

    # Count function and gradient calls
    f_calls = 0
    g_calls = 0

    # Keep a record of our current position
    x_old = copy(initial_x)
    x_new = copy(initial_x)

    # Determine number of parameters
    n = length(initial_x)
    gradient_old = Array(Float64, n)
    gradient_new = Array(Float64, n)

    # Keep a record of the current gradient
    f_calls += 1
    g_calls += 1
    f_x = d.fg!(x_old, gradient_old)
    copy!(gradient_new, gradient_old)

    # Initialize our approximate inverse Hessian
    B = copy(initial_B)

    # Iterate until convergence
    converged = false

    # Show state of the system
    tr = OptimizationTrace()
    if store_trace || show_trace
        bfgs_trace!(tr, x_old, f_x, iteration, gradient_new, B, store_trace, show_trace)
    end

    # Maintain arrays for position and gradient changes
    s = Array(Float64, n)
    y = Array(Float64, n)

    # Reuse arrays during every line search
    ls_x = Array(Float64, n)
    ls_gradient = Array(Float64, n)

    while !converged && iteration < iterations
        # Increment the iteration counter
        iteration += 1

        # Set the search direction
        p = -B * gradient_new

        # Calculate a step-size
        # TODO: Clean up backtracking_line_search
        alpha, f_update, g_update =
          backtracking_line_search!(d, x_new, p, ls_x, ls_gradient)
        f_calls += f_update
        g_calls += g_update

        # Update our position
        copy!(x_old, x_new)
        for i in 1:n
            s[i] = alpha * p[i]
            x_new[i] = x_old[i] + s[i]
        end

        # Update the gradient
        copy!(gradient_old, gradient_new)
        f_x = d.fg!(x_new, gradient_new)
        f_calls += 1
        g_calls += 1

        # Update the change in the gradient
        for i in 1:n
            y[i] = gradient_new[i] - gradient_old[i]
        end

        # Update the inverse Hessian approximation
        rho = 1.0 / dot(y, s)
        if isinf(rho)
           @printf "Cannot decrease the objective function along the current search direction\n"
           break
        end
        # TODO: Reuse storage here
        v = eye(size(B, 1)) - rho * y * s'
        if iteration == 1
            B = v'v * dot(y, s) / dot(y, y)
        else
            B = v' * B * v + rho * s * s'
        end

        # Show state of the system
        if store_trace || show_trace
            bfgs_trace!(tr, x_old, f_x, iteration, gradient_new, B, store_trace, show_trace)
        end

        # Assess convergence
        if norm(gradient_new, Inf) <= tolerance
           converged = true
        end
    end

    OptimizationResults("BFGS",
                        initial_x,
                        x_new,
                        f_x,
                        iteration,
                        converged,
                        tr,
                        f_calls,
                        g_calls)
end


text:
struct FileLogger <: AbstractLogger
    logger::SimpleLogger
    always_flush::Bool
end

"""
    FileLogger(path::AbstractString; append=false, always_flush=true)

Create a logger sink that write messages to a file specified with `path`.
To append to the file (rather than truncating the file first), use `append=true`.
If `always_flush=true` the stream is flushed after every handled log message.
"""
function FileLogger(path; append=false, kwargs...)
    filehandle = open(path, append ? "a" : "w")
    FileLogger(filehandle; kwargs...)
end

"""
    FileLogger(io::IOStream; always_flush=true)

Create a logger sink that write messages to the `io::IOStream`. The stream
is expected to be open and writeable.
If `always_flush=true` the stream is flushed after every handled log message.

# Examples
```julia
io = open("path/to/file.log", "a") # append to the file
logger = FileLogger(io)
```
"""
function FileLogger(filehandle::IOStream; always_flush=true)
    FileLogger(SimpleLogger(filehandle, BelowMinLevel), always_flush)
end

function handle_message(filelogger::FileLogger, args...; kwargs...)
    handle_message(filelogger.logger, args...; kwargs...)
    filelogger.always_flush && flush(filelogger.logger.stream)
end
shouldlog(filelogger::FileLogger, arg...) = true
min_enabled_level(filelogger::FileLogger) = BelowMinLevel
catch_exceptions(filelogger::FileLogger) = catch_exceptions(filelogger.logger)


text:
## ExoplanetsSysSim/src/transit_detection_model.jl
## (c) 2015 Eric B. Ford

# Several functions below based on https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
# That follows the procedure outlined in Burke et al.(2015).
# However we don't currently interpolate the mesthreshold to the relevant duration
import SpecialFunctions.lgamma

function real_log_choose(m::Float64, n::Float64)::Float64
  lgamma(m+1)-lgamma(n+1)-lgamma(m-n+1.0)
end

function real_binom(k::Float64, BigM::Float64, f::Float64)::Float64
    F1 = real_log_choose(BigM,k)
    F2 = k*log(f)
    F3 = (BigM-k)*log(1.0-f)
    x = exp(F1+F2+F3)
    return x
end

"""
    kepler_window_function_binomial_model(exp_num_transits_no_gaps, duty_cycle; min_transits = 3.0)

Binomial window function model for use with Kepler taken from Burke et al. (2015).

# Arguments:
- `exp_num_transits_no_gaps::Float64`: Expected number of transits assuming no gaps in the observation
- `duty_cycle::Float64`: Kepler duty cycle (i.e. fraction of observation time with usable data)
- `min_transits::Float64 = 3.0`: Minimum number of observed transits required for consideration as planet candidate

# Returns:
Probability of detecting at least min_transits given the provided properties of the Kepler target star and planet.
"""
function kepler_window_function_binomial_model(exp_num_transits_no_gaps::Float64, duty_cycle::Float64; min_transits::Float64 = 3.0)::Float64
  if exp_num_transits_no_gaps < min_transits
     return 0.0
  else
     return max(1.0 - real_binom(min_transits,exp_num_transits_no_gaps,duty_cycle), 0.0)
  end
end

"""
    kepler_window_function_binomial_model(t, exp_num_transits_no_gaps, period, duration; min_transits = 3.0)

Binomial window function model for use with Kepler taken from Burke et al. (2015).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `exp_num_transits_no_gaps::Float64`: Expected number of transits assuming no gaps in the observation
- `period::Float64`: Orbital period of planet (in days)
- `duration::Float64`: Transit duration of planet (in days)
- `min_transits::Float64 = 3.0`: Minimum number of observed transits required for consideration as planet candidate

# Returns:
Probability of detecting at least min_transits given the provided properties of the Kepler target star and planet.
"""
function kepler_window_function_binomial_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64; min_transits::Float64 = 3.0)::Float64
   kepler_window_function_binomial_model(exp_num_transits_no_gaps, t.duty_cycle, min_transits=min_transits)
end

"""
    kepler_window_function_dr25_model(t, exp_num_transits_no_gaps, period, duration)

Window function model for use with Kepler that were fit per target for DR25 (see KSCI-19101-002).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `exp_num_transits_no_gaps::Float64`: Expected number of transits assuming no gaps in the observation
- `period::Float64`: Orbital period of planet (in days)
- `duration::Float64`: Transit duration of planet (in days)

# Returns:
Probability of detecting at least 3 transits given the provided properties of the Kepler target star and planet.
"""
function kepler_window_function_dr25_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64)::Float64
   ExoplanetsSysSim.WindowFunction.eval_window_function(t.window_function_id, Duration=duration, Period=period)
end

# WARNING: Hard-coded function variable for computational efficiency,
#          replace with appropriate window function model.
#kepler_window_function = kepler_window_function_binomial_model
kepler_window_function = kepler_window_function_dr25_model

"""
    frac_depth_to_tps_depth(frac_depth)

Converts fractional transit depth to depth used by Kepler Transiting Planet Search (TPS) module. (See Christiansen et al. (2015))
NOTE: This assumes a linear limbdarkening coefficient of 0.6

# Arguments:
- `frac_depth::Float64`: Fractional transit depth of planet at transit center.

# Returns:
Kepler TPS depth.
"""
function frac_depth_to_tps_depth(frac_depth::Float64)
    alp = 1.0874
    bet = 1.0187
    REALDEPTH2TPSSQUARE = 1.0  # WARNING: Waiting for this to be confirmed
    k = sqrt(frac_depth)
    tps_depth = min( (alp-bet*k) * frac_depth* REALDEPTH2TPSSQUARE, 1.0)   # NOTE: I added the max based on common sense
    return tps_depth::Float64
end

"""
    detection_efficiency_theory(mes, expected_num_transits; min_pdet_nonzero = 0.0)

Detection efficiency model assuming perfect theoretical 7.1-sigma error function response.

# Arguments:
- `mes::Float64`: Estimated multiple event statistic (or signal-to-noise ratio) for planet
- `expected_num_transits::Float64`: Expected number of planet transits (assuming no observation gaps)
- `min_pdet_nonzero::Float64 = 0.0`: Minimum probability of detection (if transiting) to be treated as identifiable

# Results:
Probability of detection (if transiting) for planet
"""
function detection_efficiency_theory(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)
   muoffset =  0.0
   sig =  1.0
   mesthresh = 7.1
   mes *= 1.003
   if mes > (9.0 - mesthresh - muoffset)
      return 1.0
   else
      pdet = 0.5 + 0.5*erf((mes - mesthresh - muoffset) / sqrt(2.0*sig*sig))
      pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
      return pdet
   end
end

"""
    detection_efficiency_fressin2013(mes, expected_num_transits)

Detection efficiency model using the linear ramp from Fressin (2013).

# Arguments:
- `mes::Float64`: Estimated multiple event statistic (or signal-to-noise ratio) for planet
- `expected_num_transits::Float64`: Expected number of planet transits (assuming no observation gaps)

# Results:
Probability of detection (if transiting) for planet
"""
function detection_efficiency_fressin2013(mes::Float64, expected_num_transits::Float64)
    mesmin =  6.0
    mesmax =  16.0
   if mes <= mesmin
      return 0.0
   elseif mes >= mesmax
      return 1.0
   else
     return (mes - mesmin) / (mesmax - mesmin)
  end
end

"""
    detection_efficiency_christiansen2015(mes, expected_num_transits; mes_threshold = 7.1, min_pdet_nonzero = 0.0)

Detection efficiency model using a Gamma function taken from Christiansen (2015).

# Arguments:
- `mes::Float64`: Estimated multiple event statistic (or signal-to-noise ratio) for planet
- `expected_num_transits::Float64`: Expected number of planet transits (assuming no observation gaps)
- `mes_threshold::Float64 = 7.1`: Minimum multiple event statistic value for signal to have been considered a potential planet candidate signal
- `min_pdet_nonzero::Float64 = 0.0`: Minimum probability of detection (if transiting) to be treated as identifiable

# Results:
Probability of detection (if transiting) for planet
"""
function detection_efficiency_christiansen2015(mes::Float64, expected_num_transits::Float64; mes_threshold::Float64 = 7.1, min_pdet_nonzero::Float64 = 0.0)
    a =  4.65  # from code for detection_efficiency(...) at https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
   # b =  1.05
   # a =  4.35 # from arxiv abstract. Informal testing showed it didn't matter
    b =  0.98
   mes *= 1.003
   usemes = max(0.0,mes - 7.1 - (mes_threshold - 7.1))
   pdet = cdf(Gamma(a,b), usemes)
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

"""
    detection_efficiency_dr25_simple(mes, expected_num_transits; min_pdet_nonzero = 0.0)

Detection efficiency model using a Gamma function fit over a FGK sample of DR25 targets, taken from KSCI-19110-001.

# Arguments:
- `mes::Float64`: Estimated multiple event statistic (or signal-to-noise ratio) for planet
- `expected_num_transits::Float64`: Expected number of planet transits (assuming no observation gaps)
- `min_pdet_nonzero::Float64 = 0.0`: Minimum probability of detection (if transiting) to be treated as identifiable

# Results:
Probability of detection (if transiting) for planet
"""
function detection_efficiency_dr25_simple(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)::Float64
   a = 30.87  # from pg 16 of https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf
   b = 0.271
   c = 0.940
   mes *= 1.003
   dist = Gamma(a,b)
   pdet::Float64 = c*cdf(dist, mes)::Float64
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

"""
    get_param_for_detection_and_vetting_efficiency_depending_on_num_transits(num_tr)

Determine appropriate parameters to use in Gamma function model for detection efficiency and vetting of Kepler planet candidates.  Fit using the simulated transit injection tests on Kepler DR25 targets.

# Arguments:
- `num_tr::Integer`: (Expected) number of observed transits for planet

# Returns:
alpha, beta, and C, the three parameters to be used for the Gamma detection efficiency function
"""
function get_param_for_detection_and_vetting_efficiency_depending_on_num_transits(num_tr::Integer)
    if num_tr <= 3
        return (33.3884, 0.264472, 0.699093)
    elseif num_tr <= 4
        return  (32.886, 0.269577, 0.768366)
    elseif num_tr <= 5
        return (31.5196, 0.282741, 0.833673)
    elseif num_tr <= 6
        return (30.9919, 0.286979, 0.859865)
    elseif num_tr <= 9
        return (30.1906, 0.294688, 0.875042)
    elseif num_tr <= 18
        return (31.6342, 0.279425, 0.886144)
    elseif num_tr <= 36
        return (32.6448, 0.268898, 0.889724)
    else
        return (27.8185, 0.32432, 0.945075)
    end
end

"""
    detection_and_vetting_efficiency_model_v1(mes, expected_num_transits; min_pdet_nonzero = 0.0)

Detection efficiency model using a Gamma function whose parameters were determined using a fit over the simulated transit injection tests of a sample of FGK DR25 targets (see Hsu et al. (2019)).
# WARNING: Combined detection and vetting efficiency model - do NOT include additional vetting efficiency

# Arguments:
- `mes::Float64`: Estimated multiple event statistic (or signal-to-noise ratio) for planet
- `expected_num_transits::Float64`: Expected number of planet transits (assuming no observation gaps)
- `min_pdet_nonzero::Float64 = 0.0`: Minimum probability of detection (if transiting) to be treated as identifiable

# Returns:
Probability of detection (if transiting) for planet
"""
function detection_and_vetting_efficiency_model_v1(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)::Float64
    mes *= 1.003
    num_transit_int = convert(Int64,floor(expected_num_transits))
    num_transit_int += rand() < expected_num_transits-num_transit_int ? 1 : 0
    a, b, c = get_param_for_detection_and_vetting_efficiency_depending_on_num_transits(num_transit_int)
    dist = Gamma(a,b)
    pdet::Float64 = c*cdf(dist, mes)::Float64
    pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
    return pdet
end

# WARNING: Hardcoded choice of transit detection efficiency here for speed and so as to not have it hardcoded in multiple places
#detection_efficiency_model = detection_efficiency_christiansen2015
# detection_efficiency_model = detection_efficiency_dr25_simple
detection_efficiency_model = detection_and_vetting_efficiency_model_v1

"""
    vetting_efficiency_none(R_p, P)

Vetting efficiency model assuming all identified planet candidate signals are true positives

# Arguments:
- `R_p::Real`: Planet radius
- `P::Real`: Planet orbital period

# Returns:
1
"""
function vetting_efficiency_none(R_p::Real, P::Real)
    return 1.0
end

"""
    vetting_efficiency_dr25_mulders(R_p, P)

Vetting efficiency model from Mulders et al. (2018) using parameters fit for Gaia DR2 stellar properties applied to Kepler DR25 targets independent of reliability score (through e-mail communication).

# Arguments:
- `R_p::Real`: Planet radius (in R_sol)
- `P::Real`: Planet orbital period (in days)

# Returns:
Probability of planet with given properties to be vetted as a planet candidate.
"""
function vetting_efficiency_dr25_mulders(R_p::Real, P::Real)
    c = 0.93
    a_R = -0.03
    P_break = 205.
    a_P = 0.00
    b_P = -0.24
    pvet = c*(R_p/earth_radius)^a_R
    if P < P_break
        pvet *= (P/P_break)^a_P
    else
        pvet *= (P/P_break)^b_P
    end
    return pvet
end

"""
    vetting_efficiency_dr25_mulders_score_cut(R_p, P)

Vetting efficiency model from Mulders et al. (2018) using parameters fit for Gaia DR2 stellar properties applied to Kepler DR25 targets selected via reliability score (arXiv 1805.08211).

# Arguments:
- `R_p::Real`: Planet radius (in R_sol)
- `P::Real`: Planet orbital period (in days)

# Returns:
Probability of planet with given properties to be vetted as a planet candidate.
"""
function vetting_efficiency_dr25_mulders_score_cut(R_p::Real, P::Real)
    c = 0.63
    a_R = 0.19
    P_break = 53.
    a_P = -0.07
    b_P = -0.39
    pvet = c*(R_p/earth_radius)^a_R
    if P < P_break
        pvet *= (P/P_break)^a_P
    else
        pvet *= (P/P_break)^b_P
    end
    return pvet
end

# WARNING: Hardcoded choice of planet vetting efficiency here for speed and so as to not have it hardcoded in multiple places
vetting_efficiency = vetting_efficiency_none

"""
    interpolate_cdpp_to_duration_use_target_cdpp(t, duration)

Linearly interpolated CDPP using list of provided CDPP at different transit durations stored in Kepler target object (DEPRECATED).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `duration::Float64`: Transit duration (in days)

# Results:
Estimated CDPP at given transit duration for Kepler target.
"""
function interpolate_cdpp_to_duration_use_target_cdpp(t::KeplerTarget, duration::Float64)::Float64
   duration_in_hours = duration *24.0
   dur_idx = searchsortedlast(cdpp_durations,duration_in_hours)   # cdpp_durations is defined in constants.jl
   if dur_idx <= 0
      cdpp = t.cdpp[1,1]
   elseif dur_idx==length(cdpp_durations) && (duration_in_hours >= cdpp_durations[end]) # Should be 15 cdpp_durations.
      cdpp = t.cdpp[length(cdpp_durations),1]
   else
      w = ((duration_in_hours)-cdpp_durations[dur_idx]) / (cdpp_durations[dur_idx+1]-cdpp_durations[dur_idx])
      cdpp = w*t.cdpp[dur_idx+1,1] + (1-w)*t.cdpp[dur_idx,1]
   end
   return cdpp
end

"""
    interpolate_cdpp_to_duration_lookup_cdpp(t, duration)

Linearly interpolated CDPP using list of provided CDPP at different transit durations stored in stellar catalog dataframe.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `duration::Float64`: Transit duration (in days)

# Results:
Estimated CDPP at given transit duration for Kepler target.
"""
function interpolate_cdpp_to_duration_lookup_cdpp(t::KeplerTarget, duration::Float64)::Float64
   duration_in_hours = duration *24.0
   dur_idx = searchsortedlast(cdpp_durations,duration_in_hours)   # cdpp_durations is defined in constants.jl
   get_cdpp(i::Integer) = 1.0e-6*sqrt(cdpp_durations[i]/24.0/LC_duration)*star_table(t,duration_symbols[i])

   if dur_idx <= 0
      cdpp = get_cdpp(1)
   elseif dur_idx==length(cdpp_durations) && (duration_in_hours >= cdpp_durations[end]) # Should be 15 cdpp_durations.
      cdpp = get_cdpp(length(cdpp_durations))
   else
      w = ((duration_in_hours)-cdpp_durations[dur_idx]) / (cdpp_durations[dur_idx+1]-cdpp_durations[dur_idx])
      cdpp = w*get_cdpp(dur_idx+1) + (1-w)*get_cdpp(dur_idx)
   end
   return cdpp
end

#interpolate_cdpp_to_duration = interpolate_cdpp_to_duration_use_target_cdpp
interpolate_cdpp_to_duration = interpolate_cdpp_to_duration_lookup_cdpp

"""
    calc_snr_if_transit_cdpp(t, depth, duration, cdpp, sim_param; num_transit = 1)

Calculate the expected multiple event statistic (signal-to-noise ratio) for planet around Kepler target star using CDPP.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `depth::Float64`: (Fractional) transit depth
- `duration::Float64`: Transit duration (in days)
- `cdpp::Float64`: CDPP for target star given transit duration
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
(Expected) multiple event statistic (signal-to-noise ratio)
"""
function calc_snr_if_transit_cdpp(t::KeplerTarget, depth::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  # depth_tps = frac_depth_to_tps_depth(depth)                  # TODO SCI:  WARNING: Hardcoded this conversion.  Remove once depth calculated using limb darkening model
  # snr = depth_tps*sqrt(num_transit*duration*LC_rate)/cdpp     # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
  snr = depth*sqrt(num_transit*duration*LC_rate)/cdpp     # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
end

"""
    calc_snr_if_transit(t, depth, duration, osd, sim_param; num_transit = 1)

Calculate the expected multiple event statistic (signal-to-noise ratio) for planet around Kepler target star using 1-sigma depth function (OSD).
# NOTE: Assumes OSD functions have already been read in for all relevant Kepler targets.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `depth::Real`: (Fractional) transit depth
- `duration::Real`: Transit duration (in days)
- `osd::Real`: OSD for target star given transit duration and period
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Real = 1`: Expected number of transits observed

# Returns:
(Expected) multiple event statistic (signal-to-noise ratio)
"""
function calc_snr_if_transit(t::KeplerTarget, depth::Real, duration::Real, osd::Real, sim_param::SimParam; num_transit::Real = 1)
  # depth_tps = frac_depth_to_tps_depth(depth)                  # WARNING: Hardcoded this conversion
  # snr = depth_tps/osd*1.0e6 # osd is in ppm
  snr = depth/osd*1.0e6 # osd is in ppm
end

"""
    calc_snr_if_transit_central(t, s, p, sim_param)

Calculate the expected multiple event statistic (signal-to-noise ratio) for planet around Kepler target star using 1-sigma depth function (OSD) at transit center.
# NOTE: Assumes OSD functions have already been read in for all relevant Kepler targets.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
(Expected) multiple event statistic (signal-to-noise ratio) at transit center
"""
function calc_snr_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_eff_central(t,s,p)
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration_central = get_durations_searched_Kepler(period,duration_central)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd_central = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration_central)
  if osd_duration_central > duration_central				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd_central = osd_central*osd_duration_central/duration_central
  end
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit(t,depth,duration_central,osd_central, sim_param,num_transit=num_transit)
end

"""
    calc_snr_if_transit_central_cdpp(t, s, p, sim_param)

Calculate the expected multiple event statistic (signal-to-noise ratio) for planet around Kepler target star using CDPP at transit center.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
(Expected) multiple event statistic (signal-to-noise ratio) at transit center
"""
function calc_snr_if_transit_central_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_eff_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit_cdpp(t,depth,duration_central,cdpp, sim_param,num_transit=num_transit)
end

"""
    calc_prob_detect_if_transit(t, snr, period, duration, sim_param; num_transit = 1)
    calc_prob_detect_if_transit(t, depth, period, duration, osd, sim_param; num_transit = 1)

Calculate probability of detecting planet signal (if planet transits) using 1-sigma depth (OSD) function.
# NOTE: Assumes OSD functions already read in.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `snr::Float64`: Expected multiple event statistic (signal-to-noise ratio)
- `depth::Float64`: (Fractional) transit depth
- `period::Float64`: Orbital period (in days)
- `duration::Float64`: Transit duration (in days)
- `osd::Float64`: OSD for target star given transit duration and period
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
Probability of detecting planet (if it transits)
"""
function calc_prob_detect_if_transit(t::KeplerTarget, snr::Float64, period::Float64, duration::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  min_pdet_nonzero = 1.0e-4                                                # TODO OPT: Consider raising threshold to prevent a plethora of planets that are very unlikely to be detected due to using 0.0 or other small value here
  wf = kepler_window_function(t, num_transit, period, duration)
  return wf*detection_efficiency_model(snr, num_transit, min_pdet_nonzero=min_pdet_nonzero)
end

function calc_prob_detect_if_transit(t::KeplerTarget, depth::Float64, period::Float64, duration::Float64, osd::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  snr = calc_snr_if_transit(t,depth,duration,osd, sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=num_transit)
end

"""
    calc_prob_detect_if_transit_cdpp(t, depth, period, duration, cdpp, sim_param; num_transit = 1)

Calculate probability of detecting planet signal (if planet transits) using CDPP.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `depth::Float64`: (Fractional) transit depth
- `period::Float64`: Orbital period (in days)
- `duration::Float64`: Transit duration (in days)
- `cdpp::Float64`: CDPP for target star given transit duration
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
Probability of detecting planet (if it transits)
"""
function calc_prob_detect_if_transit_cdpp(t::KeplerTarget, depth::Float64, period::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  snr = calc_snr_if_transit_cdpp(t,depth,duration,cdpp, sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=num_transit)
end

"""
    calc_prob_detect_if_transit_central(t, s, p, sim_param)

Calculate probability of detecting planet signal (if planet transits) at transit center using 1-sigma depth (OSD) functions.
# NOTE: Assumes OSD functions already read in.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
Probability of detecting planet (if it transits) at transit center
"""
function calc_prob_detect_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_eff_central(t,s,p)
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration_central = get_durations_searched_Kepler(period,duration_central)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd_central = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration_central)
  if osd_duration_central > duration_central				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd_central = osd_central*osd_duration_central/duration_central
  end
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration_central,osd_central, sim_param, num_transit=ntr)
end

"""
    calc_prob_detect_if_transit_central_cdpp(t, s, p, sim_param)

Calculate probability of detecting planet signal (if planet transits) at transit center using CDPP

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
Probability of detecting planet (if it transits) at transit center
"""
function calc_prob_detect_if_transit_central_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_eff_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit_cdpp(t,depth,period,duration_central,cdpp, sim_param, num_transit=ntr)
end

"""
    calc_prob_detect_if_transit_with_actual_b(t, s, p, sim_param)

Calculate probability of detecting planet signal (if planet transits) at transit center using 1-sigma depth (OSD) functions and the impact parameter of the provided orbit.
# NOTE: Assumes OSD functions already read in.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
Probability of detecting planet (if it transits) given the impact parameter of the planet's orbit
"""
function calc_prob_detect_if_transit_with_actual_b(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration_eff(t,s,p)
  b = calc_impact_parameter(t.sys[s],p)
  snr_correction = calc_depth_correction_for_grazing_transit(b,size_ratio)
  depth *= snr_correction
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration = get_durations_searched_Kepler(period,duration)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration)
  if osd_duration > duration				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd = osd*osd_duration/duration
  end
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration,osd, sim_param, num_transit=ntr)
end

"""
    calc_prob_detect_if_transit_with_actual_b_cdpp(t, s, p, sim_param)

Calculate probability of detecting planet signal (if planet transits) at transit center using CDPP and the impact parameter of the provided orbit.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters

# Returns:
Probability of detecting planet (if it transits) given the impact parameter of the planet's orbit
"""
function calc_prob_detect_if_transit_with_actual_b_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration_eff(t,s,p)
  b = calc_impact_parameter(t.sys[s],p)
  snr_correction = calc_depth_correction_for_grazing_transit(b,size_ratio)
  depth *= snr_correction
  cdpp = interpolate_cdpp_to_duration(t, duration)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit_cdpp(t,depth,period,duration,cdpp, sim_param, num_transit=ntr)
end

"""
    calc_ave_prob_detect_if_transit_from_snr(t, snr_central, period, duration_central, size_ratio, osd_central, sim_param; num_transit = 1)

Calculate probability of detecting planet signal (if planet transits) using 1-sigma depth (OSD) function and averaged over impact parameters b~U[0,1).
# NOTE: Assumes OSD functions already read in.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `snr_central::Float64`: Expected multiple event statistic (signal-to-noise ratio) at transit center
- `period::Float64`: Orbital period (in days)
- `duration_central::Float64`: Transit duration (in days) at transit center
- `size_ratio::Float64`: Ratio of planet-to-star radii
- `osd_central::Float64`: OSD for target star given transit duration at transit center and period
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
Probability of detecting planet (if it transits) averaged over impact parameter
"""
function calc_ave_prob_detect_if_transit_from_snr(t::KeplerTarget, snr_central::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, osd_central::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  min_pdet_nonzero = 1.0e-4
  wf = kepler_window_function(t, num_transit, period, duration_central)

  detection_efficiency_central = detection_efficiency_model(snr_central, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  if wf*detection_efficiency_central <= min_pdet_nonzero
     return 0.
  end

  # Breaking integral into two sections [0,1-b_boundary) and [1-b_boundary,1], so need at least 5 points to evaluate integral via trapezoid rule
  num_impact_param_low_b =  7                            # Number of points to evaluate integral over [0,1-b_boundary) via trapezoid rule
  num_impact_param_high_b = 5 # (size_ratio<=0.05) ? 5 : 11  # Number of points to evaluate integral over [1-b_boudnary,1) via trapezoid rule.  If using 2*size_ratio for bondary for small planets, then keep this odd, so one point lands on 1-size_ratio.
  @assert(num_impact_param_low_b >= 5)
  @assert(num_impact_param_high_b >= 3)
  num_impact_param = num_impact_param_low_b+num_impact_param_high_b-1 # One point is shared
  b_boundary = (size_ratio <= 0.15) ? 2*size_ratio : min(max(0.3,size_ratio),0.5)
  b = Array{Float64}(undef,num_impact_param)
  weight = Array{Float64}(undef,num_impact_param)
  b[1:num_impact_param_low_b] = range(0.0,stop=1-b_boundary,length=num_impact_param_low_b)
  b[num_impact_param_low_b:num_impact_param] .= range(1-b_boundary,stop=1.0,length=num_impact_param_high_b)
  weight[1:num_impact_param_low_b] .= (1-b_boundary)/(num_impact_param_low_b-1)  # Points for first integral
  weight[1] *= 0.5                        # Lower endpoint of first integral
  weight[num_impact_param_low_b] *= 0.5   # Upper endpoint of first integral
  weight[num_impact_param_low_b] += 0.5*(b_boundary)/(num_impact_param_high_b-1) # Also lower endpoint of second integral
  weight[(num_impact_param_low_b+1):num_impact_param] .= b_boundary/(num_impact_param_high_b-1)
  weight[num_impact_param] *= 0.5         # Upper endpoint of second integral
  @assert isapprox(sum(weight),1.0)

  function integrand(b::Float64)::Float64
     depth_factor = calc_depth_correction_for_grazing_transit(b,size_ratio)
     duration_factor = calc_transit_duration_eff_factor_for_impact_parameter_b(b,size_ratio)
     kepid = StellarTable.star_table(t.sys[1].star.id, :kepid)
     osd_duration = get_durations_searched_Kepler(period,duration_central*duration_factor)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
     osd = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration)
     if osd_duration > duration_central*duration_factor				#use a correcting factor if this duration is lower than the minimum searched for this period.
	osd = osd*osd_duration/(duration_central*duration_factor)
     end
     snr_factor = depth_factor*(osd_central/osd)
     detection_efficiency_model(snr_central*snr_factor, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  end

  ave_detection_efficiency = sum(weight .* map(integrand,b)::Vector{Float64} )

  return wf*ave_detection_efficiency
end

"""
    calc_ave_prob_detect_if_transit_from_snr_cdpp(t, snr_central, period, duration_central, size_ratio, osd_central, sim_param; num_transit = 1)

Calculate probability of detecting planet signal (if planet transits) using CDPP and averaged over impact parameters b~U[0,1).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `snr_central::Float64`: Expected multiple event statistic (signal-to-noise ratio) at transit center
- `period::Float64`: Orbital period (in days)
- `duration_central::Float64`: Transit duration (in days) at transit center
- `size_ratio::Float64`: Ratio of planet-to-star radii
- `cdpp_central::Float64`: CDPP for target star given transit duration at transit center
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
Probability of detecting planet (if it transits) averaged over impact parameter
"""
function calc_ave_prob_detect_if_transit_from_snr_cdpp(t::KeplerTarget, snr_central::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, cdpp_central::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  min_pdet_nonzero = 1.0e-4
  wf = kepler_window_function(t, num_transit, period, duration_central)

  detection_efficiency_central = detection_efficiency_model(snr_central, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  if wf*detection_efficiency_central <= min_pdet_nonzero
     return 0.
  end

  # Breaking integral into two sections [0,1-b_boundary) and [1-b_boundary,1], so need at least 5 points to evaluate integral via trapezoid rule
  num_impact_param_low_b =  7                            # Number of points to evaluate integral over [0,1-b_boundary) via trapezoid rule
  num_impact_param_high_b = 5 # (size_ratio<=0.05) ? 5 : 11  # Number of points to evaluate integral over [1-b_boudnary,1) via trapezoid rule.  If using 2*size_ratio for bondary for small planets, then keep this odd, so one point lands on 1-size_ratio.
  @assert(num_impact_param_low_b >= 5)
  @assert(num_impact_param_high_b >= 3)
  num_impact_param = num_impact_param_low_b+num_impact_param_high_b-1 # One point is shared
  b_boundary = (size_ratio <= 0.15) ? 2*size_ratio : min(max(0.3,size_ratio),0.5)
  b = Array{Float64}(undef,num_impact_param)
  weight = Array{Float64}(undef,num_impact_param)
  b[1:num_impact_param_low_b] = range(0.0,stop=1-b_boundary,length=num_impact_param_low_b)
  b[num_impact_param_low_b:num_impact_param] .= range(1-b_boundary,stop=1.0,length=num_impact_param_high_b)
  weight[1:num_impact_param_low_b] .= (1-b_boundary)/(num_impact_param_low_b-1)  # Points for first integral
  weight[1] *= 0.5                        # Lower endpoint of first integral
  weight[num_impact_param_low_b] *= 0.5   # Upper endpoint of first integral
  weight[num_impact_param_low_b] += 0.5*(b_boundary)/(num_impact_param_high_b-1) # Also lower endpoint of second integral
  weight[(num_impact_param_low_b+1):num_impact_param] .= b_boundary/(num_impact_param_high_b-1)
  weight[num_impact_param] *= 0.5         # Upper endpoint of second integral
  @assert isapprox(sum(weight),1.0)

  function integrand(b::Float64)::Float64
     depth_factor = calc_depth_correction_for_grazing_transit(b,size_ratio)
     duration_factor = calc_transit_duration_eff_factor_for_impact_parameter_b(b,size_ratio)
     cdpp = interpolate_cdpp_to_duration(t,duration_central*duration_factor)
     snr_factor = depth_factor*sqrt(duration_factor)*(cdpp_central/cdpp)
     detection_efficiency_model(snr_central*snr_factor, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  end

  ave_detection_efficiency = sum(weight .* map(integrand,b)::Vector{Float64} )

  return wf*ave_detection_efficiency
end

"""
    calc_ave_prob_detect_if_transit_cdpp(t, depth, period, duration_central, size_ratio, sim_param; num_transit = 1)
    calc_ave_prob_detect_if_transit_cdpp(t, s, p, sim_param)

Calculate probability of detecting planet signal (if planet transits) using CDPP and averaged over impact parameters b~U[0,1).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `depth::Float64`: (Fractional) transit depth
- `period::Float64`: Orbital period (in days)
- `duration_central::Float64`: Transit duration (in days) at transit center
- `size_ratio::Float64`: Ratio of planet-to-star radii
- `s::Integer`: Star index within Kepler target
- `p::Integer`: Planet index within Kepler target
- `sim_param::SimParam`: Simulation parameters
- `num_transit::Float64 = 1`: Expected number of transits observed

# Returns:
Probability of detecting planet (if it transits) averaged over impact parameter
"""
function calc_ave_prob_detect_if_transit_cdpp(t::KeplerTarget, depth::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  cdpp_central = interpolate_cdpp_to_duration(t, duration_central)
  snr_central = calc_snr_if_transit(t,depth,duration_central,cdpp_central, sim_param, num_transit=num_transit)
  return calc_ave_prob_detect_if_transit_from_snr_cdpp(t, snr_central, period, duration_central, size_ratio, cdpp_central, sim_param, num_transit=num_transit)
end

function calc_ave_prob_detect_if_transit_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  period = t.sys[s].orbit[p].P
  duration_central = calc_transit_duration_eff_central(t,s,p)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_ave_prob_detect_if_transit_cdpp(t,depth,period,duration_central, size_ratio, sim_param, num_transit=ntr)
end


text:
__precompile__()

module Drs

include("DrsMathProgSolverInterface.jl")
using .DrsMathProgSolverInterface

export DrsMathProgSolver

end

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
@testset "Plotting" begin
  if visualtests
    @testset "h-scatter" begin
      @plottest begin
        sdata = readgeotable(joinpath(datadir,"samples2D.tsv"), delim='\t', coordnames=[:x,:y])
        hscatter(sdata, :value, lags=[0.,1.,2.,3.], layout=(2,2), size=(600,600))
      end joinpath(datadir,"HScatter.png") !istravis
    end

    # UNCOMMENT WHEN GR ADDS SUPPORT TO POLAR PLOTS
    # @testset "varplane" begin
      # imgdata = readdlm(joinpath(datadir,"anisotropic.tsv"))
      # geodata = RegularGridData{Float64}(Dict(:z => imgdata))

      # function plot_varplane(fname)
        # varplane(geodata, :z, maxlag=50., size=(500,500))
        # png(fname)
      # end
      # refimg = joinpath(datadir,"VarPlane.png")
      # @test test_images(VisualTest(plot_varplane, refimg), popup=!istravis) |> success
    # end
  end
end