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

#using Distributions
#include("constants.jl")
#include("newPDST.jl")		#includes a function to calculate if given durations match those observed by Kepler for a given period

# Needed if using full noise model
using LinearAlgebra
using PDMats

#  Starting Section of Observables that are actually used
struct TransitPlanetObs
  # ephem::ephemeris_type     # For now hardcode P and t0, see transit_observation_unused.jl to reinstate
  period::Float64             # days
  t0::Float64                 # days
  depth::Float64              # fractional
  duration::Float64           # days; Full-width, half-max-duration until further notice
  # ingress_duration::Float64   # days;  QUERY:  Will we want to use the ingress/egress duration for anything?
end
TransitPlanetObs() = TransitPlanetObs(0.0,0.0,0.0,0.0)

struct StarObs
  radius::Float64      # in Rsol
  mass::Float64        # in Msol
  id::Int64            # row number in stellar dataframe
end

period(obs::TransitPlanetObs) = obs.period
depth(obs::TransitPlanetObs) = obs.depth
duration(obs::TransitPlanetObs) = obs.duration

period_given_semimajor_axis(a::Float64, M::Float64) = sqrt((4pi^2*a^3)/(grav_const*M))

"""
    semimajor_axis(P, M)
    semimajor_axis(ps, id)

Calculate semimajor axis

# Arguments:
- `P::Float64`: Orbital period (in days)
- `M::Float64`: Planet mass (in M_sol)
- `ps::PlanetarySystemAbstract`: Planetary system object
- `id::Integer`: Planet/orbit index in planetary system object

# Returns:
Semimajor axis of planet (in AU)
"""
semimajor_axis(P::Float64, M::Float64) = (grav_const/(4pi^2)*M*P*P)^(1/3)

function semimajor_axis(ps::PlanetarySystemAbstract, id::Integer)
  M = mass(ps.star) + ps.planet[id].mass   # TODO SCI DETAIL: Replace with Jacobi mass?  Not important unless start including TTVs, even then unlikely to matter
  @assert(M>0.0)
  @assert(ps.orbit[id].P>0.0)
  return semimajor_axis(ps.orbit[id].P,M)
end

"""
    calc_transit_depth(t, s, p)

Calculate fractional transit depth of planet transiting host star (including limb darkening)

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Fractional transit depth at transit midpoint
"""
function calc_transit_depth(t::KeplerTarget, s::Integer, p::Integer)  # WARNING: IMPORTANT: Assumes non-grazing transit
  radius_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  #b = calc_impact_parameter(t.sys[s].planet, p) # If limb darkening should know about which chord the planet takes set b to impact parameter, rather than 0.0.
  depth = depth_at_midpoint(radius_ratio, t.sys[s].star.ld)   # Includes limb darkening
  depth *=  flux(t.sys[s].star)/flux(t)                      # Flux ratio accounts for dilution
end

"""
    calc_transit_duration_central_circ_small_angle_approx(ps, pl)
    calc_transit_duration_central_circ_small_angle_approx(t, s, p)

Calculate transit duration of planet across host star disk if it transits across the center of the stellar disk assuming a circular orbit using the small angle approximation

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
function calc_transit_duration_central_circ_small_angle_approx(ps::PlanetarySystemAbstract, pl::Integer)
  duration = rsol_in_au*ps.star.radius * ps.orbit[pl].P /(pi*semimajor_axis(ps,pl) )
end

calc_transit_duration_central_circ_small_angle_approx(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_circ_small_angle_approx(t.sys[s],p)

"""
    calc_transit_duration_central_circ_with_arcsin(ps, pl)
    calc_transit_duration_central_circ_with_arcsin(t, s, p)

Calculate transit duration of planet across host star disk if it transits across the center of the stellar disk assuming a circular orbit without the small angle approximation

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
function calc_transit_duration_central_circ_with_arcsin(ps::PlanetarySystemAbstract, pl::Integer)
  asin_arg = rsol_in_au*ps.star.radius/semimajor_axis(ps,pl)
  duration = ps.orbit[pl].P/pi * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)
end
calc_transit_duration_central_circ_with_arcsin(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_circ_with_arcsin(t.sys[s],p)

"""
    calc_transit_duration_central_circ(ps, pl)
    calc_transit_duration_central_circ(t, s, p)

Wrapper function to calculate transit duration of planet across host star disk if it transits across the center of the stellar disk assuming a circular orbit

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
#calc_transit_duration_central_circ(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_central_circ_small_angle_approx(ps,pl)
calc_transit_duration_central_circ(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_central_circ_with_arcsin(ps,pl)

calc_transit_duration_central_circ(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_circ(t.sys[s],p)

"""
    calc_transit_duration_central_small_angle_approx(ps, pl)
    calc_transit_duration_central_small_angle_approx(t, s, p)

Calculate transit duration of planet across host star disk if it transits across the center of the stellar disk using the small angle approximation

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
function calc_transit_duration_central_small_angle_approx(ps::PlanetarySystemAbstract, pl::Integer)
  ecc = ps.orbit[pl].ecc
  sqrt_one_minus_ecc_sq = sqrt((1+ecc)*(1-ecc))
  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  vel_fac = sqrt_one_minus_ecc_sq/one_plus_e_sin_w
  duration = calc_transit_duration_central_circ_small_angle_approx(ps,pl) * vel_fac
end
calc_transit_duration_central_small_angle_approx(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_small_angle_approx(t.sys[s],p)

"""
    calc_transit_duration_central_winn2010(ps, pl)
    calc_transit_duration_central_winn2010(t, s, p)

Calculate transit duration of planet across host star disk if it transits across the center of the stellar disk using the formulation from Winn (2010)

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
function calc_transit_duration_central_winn2010(ps::PlanetarySystemAbstract, pl::Integer)
  ecc = ps.orbit[pl].ecc
  sqrt_one_minus_ecc_sq = sqrt((1+ecc)*(1-ecc))
  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  vel_fac = sqrt_one_minus_ecc_sq/one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  asin_arg = rsol_in_au*ps.star.radius/(semimajor_axis(ps,pl))
  # Based on Winn 2010
  duration = ( asin_arg<1.0 ?  asin(asin_arg) : 1.0 ) * ps.orbit[pl].P*radial_separation_over_a/(pi*sqrt_one_minus_ecc_sq)
end
calc_transit_duration_central_winn2010(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_winn2010(t.sys[s],p)

"""
    calc_transit_duration_central_kipping2010(ps, pl)
    calc_transit_duration_central_kipping2010(t, s, p)

Calculate transit duration of planet across host star disk if it transits across the center of the stellar disk using the formulation from Kipping (2010)

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
function calc_transit_duration_central_kipping2010(ps::PlanetarySystemAbstract, pl::Integer)
  ecc = ps.orbit[pl].ecc
  sqrt_one_minus_ecc_sq = sqrt((1+ecc)*(1-ecc))
  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  vel_fac = sqrt_one_minus_ecc_sq/one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  # Based on pasting cos i = 0 into Eqn 15 from Kipping 2010
  asin_arg = rsol_in_au*ps.star.radius/(semimajor_axis(ps,pl)* radial_separation_over_a)
  duration = ps.orbit[pl].P*radial_separation_over_a^2/(pi*sqrt_one_minus_ecc_sq) * ( asin_arg<1.0 ?  asin(asin_arg) : 1.0 )
end
calc_transit_duration_central_kipping2010(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_kipping2010(t.sys[s],p)

"""
    calc_transit_duration_central(ps, pl)
    calc_transit_duration_central(t, s, p)

Wrapper function to calculate transit duration of planet across host star disk if it transits across the center of the stellar disk

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days) if planet transits across center of stellar disk
"""
#calc_transit_duration_central(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_central_small_angle_approx(ps,pl)
#calc_transit_duration_central(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_central_winn2010(ps,pl)
calc_transit_duration_central(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_central_kipping2010(ps,pl)

calc_transit_duration_central(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central(t.sys[s],p)
calc_transit_duration_eff_central(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central(t.sys[s],p)

"""
    calc_transit_duration_factor_for_impact_parameter_b(b, p)

Calculate scaling factor to convert central transit duration to transit duration given an orbit with impact parameter b
NOTE: Used for single observations of transit duration

# Arguments:
- `b::Real`: Impact parameter
- `p::Real`: Planet-to-star radius ratio

# Returns:
Transit duration scaling factor
"""
function calc_transit_duration_factor_for_impact_parameter_b(b::T, p::T)  where T <:Real
    @assert(zero(b)<=b)         # b = Impact Parameter
    @assert(zero(p)<=p<one(p))  # p = R_p/R_star
    if b < 1-p
          duration_ratio = sqrt((1-b)*(1+b))  # Approximation to (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2, which is itself an approximation
    else
          return zero(b)
    end
end

"""
    calc_transit_duration_eff_factor_for_impact_parameter_b(b, p)

Calculate scaling factor to convert central transit duration to transit duration given an orbit with impact parameter b
NOTE: Used for single observations of effective transit duration for use in SNR calculations.

# Arguments:
- `b::Real`: Impact parameter
- `p::Real`: Planet-to-star radius ratio

# Returns:
Transit duration scaling factor
"""
function calc_transit_duration_eff_factor_for_impact_parameter_b(b::T, p::T)  where T <:Real
  @assert(zero(b)<=b)         # b = Impact Parameter
  @assert(zero(p)<=p<one(p))  # p = R_p/R_star
  if b < 1-3p                 # Far enough from grazing for approximation
        #duration_ratio = sqrt(1-b^2)  # Approximation to (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2, which is itself an approximation
        duration_ratio = sqrt((1-b)*(1+b))  # Approximation to (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2, which is itself an approximation
    elseif b < 1-p            # Planet is fully inscribed at mid-transit
        #duration_ratio = (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2 # Average of full and flat transit durations approximates to duration between center of planet being over limb of star
        duration_ratio = (sqrt(((1+p)+b)*((1+p)-b))+sqrt(((1-p)+b)*((1-p)-b)))/2 # Average of full and flat transit durations approximates to duration between center of planet being over limb of star
    elseif b < 1+p            # Planet never fully inscribed by star
        #duration_ratio = sqrt((1+p)^2-b^2)/2  # /2 since now triangular
        duration_ratio = sqrt(((1+p)+b)*((1+p)-b))/2  # /2 since now triangular
    else                      # There's no transit
        duration_ratio = zero(b)
    end
    return duration_ratio
end

"""
    calc_effective_transit_duration_factor_for_impact_parameter_b(b, p)

Calculate scaling factor to convert central transit duration to transit duration given an orbit with impact parameter b
NOTE: Used for sky-averaging observations where multiple transit durations are sampled for different values of b.

# Arguments:
- `b::Real`: Impact parameter
- `p::Real`: Planet-to-star radius ratio

# Returns:
Transit duration scaling factor
"""
function calc_effective_transit_duration_factor_for_impact_parameter_b(b::T, p::T)  where T <:Real
  @assert(zero(b)<=b)         # b = Impact Parameter
  @assert(zero(p)<=p<one(p))  # p = R_p/R_star
  if b < 1-3p                 # Far enough from grazing for approximation
        duration_ratio = sqrt((1+b)*(1-b))  # Approximation to (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2, which is itself an approximation
        area_ratio = one(p)
    elseif b < 1-p            # Planet is fully inscribed at mid-transit
        #duration_ratio = (sqrt((1+p)^2-b^2)+sqrt((1-p)^2-b^2))/2 # Average of full and flat transit durations approximates to duration between center of planet being over limb of star
        duration_ratio = (sqrt(((1+p)+b)*((1+p)-b))+sqrt(((1-p)+b)*((1-p)-b)))/2 # Average of full and flat transit durations approximates to duration between center of planet being over limb of star
        area_ratio = one(p)
    elseif b < 1+p            # Planet never fully inscribed by star
        #duration_ratio = sqrt((1+p)^2-b^2)/2  # /2 since now triangular
        duration_ratio = sqrt(((1+p)+b)*((1+p)-b))/2  # /2 since now triangular
        #area_ratio = (p^2*acos((b^2+p^2-1)/(2*b*p))+acos((b^2+1-p^2)/(2b))-0.5*sqrt((1+p-b)*(p+b-1)*(1-p+b)*(1+p+b))) / (pi*p^2)
        #area_ratio = (p^2*acos((b^2-(1-p)*(1+p))/(2*b*p))+acos((b^2+(1+p)*(1-p))/(2b))-0.5*sqrt((1+p-b)*(p+b-1)*(1-p+b)*(1+p+b))) / (pi*p^2)
        acos_arg1 = max(-1.0,min(1.0,(b^2-(1-p)*(1+p))/(2*b*p)))
        acos_arg2 = max(-1.0,min(1.0,(b^2+(1+p)*(1-p))/(2b)))
        sqrt_arg = max(0.0,(1+p-b)*(p+b-1)*(1-p+b)*(1+p+b))
        area_ratio = (p^2*acos(acos_arg1)+acos(acos_arg2)-0.5*sqrt(sqrt_arg)) / (pi*p^2)
    else                      # There's no transit
        duration_ratio = zero(b)
        area_ratio = zero(p)
    end
    return duration_ratio*area_ratio
end

"""
    calc_depth_correction_for_grazing_transit(b, p)

Calculate scaling factor to correct transit depth to account for grazing transits when only part of the planet disk blocks starlight at mid-transit. Source: Eqn 14 of http://mathworld.wolfram.com/Circle-CircleIntersection.html
NOTE: Assumes uniform surface brightness star

# Arguments:
- `b::Real`: Impact parameter
- `p::Real`: Planet-to-star radius ratio

# Returns:
Transit depth scaling factor
"""
function calc_depth_correction_for_grazing_transit(b::T, p::T)  where T <:Real
  @assert(zero(b)<=b)         # b = Impact Parameter
  @assert(zero(p)<=p<one(p))  # p = R_p/R_star
  if b < 1-p                  # Planet fully inscribed by star
        area_ratio = one(p)
    elseif b < 1+p            # Planet never fully inscribed by star
        #area_ratio = (p^2*acos((b^2+p^2-1)/(2*b*p))+acos((b^2+1-p^2)/(2b))-0.5*sqrt((1+p-b)*(p+b-1)*(1-p+b)*(1+p+b))) / (pi*p^2)
        acos_arg1 = max(-1.0,min(1.0,(b^2-(1-p)*(1+p))/(2*b*p)))
        acos_arg2 = max(-1.0,min(1.0,(b^2+(1+p)*(1-p))/(2b)))
        sqrt_arg = max(0.0,(1+p-b)*(p+b-1)*(1-p+b)*(1+p+b))
        area_ratio = (p^2*acos(acos_arg1)+acos(acos_arg2)-0.5*sqrt(sqrt_arg)) / (pi*p^2)
    else                      # There's no transit
        area_ratio = zero(p)
    end
    return area_ratio
end

# Transit durations to be used for observations of transit duration

"""
    calc_transit_duration_small_angle_approx(ps, pl)
    calc_transit_duration_small_angle_approx(t, s, p)

Calculate transit duration of planet across host star disk using the small angle approximation

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days)
"""
function calc_transit_duration_small_angle_approx(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)-size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  duration_ratio_for_impact_parameter = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w

  duration = duration_central_circ * duration_ratio_for_impact_parameter * vel_fac
end
calc_transit_duration_small_angle_approx(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_small_angle_approx(t.sys[s],p)

"""
    calc_transit_duration_central_winn2010(ps, pl)
    calc_transit_duration_central_winn2010(t, s, p)

Calculate transit duration of planet across host star disk using the formulation from Winn (2010)

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days)
"""
function calc_transit_duration_winn2010(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)-size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  arcsin_circ_central = pi/ps.orbit[pl].P*duration_central_circ

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  duration_ratio_for_impact_parameter = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)

  # WARNING: This is technically an approximation.  It avoids small angle for non-grazing transits, but does use a variant of the small angle approximation for nearly grazing transits.
  asin_arg = (arcsin_circ_central * duration_ratio_for_impact_parameter)
  duration = ps.orbit[pl].P/pi * radial_separation_over_a/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)
  duration = duration_central_cric * radial_separation_over_a/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)/arcsin_circ_central
end
calc_transit_duration_winn2010(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_winn2010(t.sys[s],p)

"""
    calc_transit_duration_central_kipping2010(ps, pl)
    calc_transit_duration_central_kipping2010(t, s, p)

Calculate transit duration of planet across host star disk using the formulation from Kipping (2010)

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days)
"""
function calc_transit_duration_kipping2010(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)-size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  arcsin_circ_central = pi/ps.orbit[pl].P*duration_central_circ

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  duration_ratio_for_impact_parameter = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)

  # WARNING: This is technically an approximation (see Kipping 2010 Eqn 15).  It avoids small angle for non-grazing transits, but does use a variant of the small angle approximation for nearly grazing transits.
  asin_arg = (arcsin_circ_central * duration_ratio_for_impact_parameter/radial_separation_over_a)
  duration = ps.orbit[pl].P/pi * radial_separation_over_a^2/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)
  #duration = duration_central_circ * radial_separation_over_a^2/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)/arcsin_circ_central
end
calc_transit_duration_kipping2010(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_kipping2010(t.sys[s],p)

"""
    calc_transit_duration(ps, pl)
    calc_transit_duration(t, s, p)

Wrapper function to calculate transit duration of planet across host star disk

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Transit duration (in days)
"""
#calc_transit_duration(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_small_angle_approx(ps,pl)
#calc_transit_duration(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_winn2010(ps,pl)
calc_transit_duration(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_kipping2010(ps,pl)
calc_transit_duration(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration(t.sys[s],p)

# Effective transit durations to be used for SNR calculations

"""
    calc_transit_duration_eff_small_angle_approx(ps, pl)
    calc_transit_duration_eff_small_angle_approx(t, s, p)

Calculate effective transit duration of planet across host star disk using the small angle approximation
NOTE: To be used for SNR calculations

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Effective transit duration (in days)
"""
function calc_transit_duration_eff_small_angle_approx(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)+size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  duration_ratio_for_impact_parameter = calc_transit_duration_eff_factor_for_impact_parameter_b(b,size_ratio)

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w

  duration = duration_central_circ * duration_ratio_for_impact_parameter * vel_fac
end
calc_transit_duration_eff_small_angle_approx(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_eff_small_angle_approx(t.sys[s],p)

"""
    calc_transit_duration_eff_winn2010(ps, pl)
    calc_transit_duration_eff_winn2010(t, s, p)

Calculate effective transit duration of planet across host star disk using the formulation from Winn (2010)
NOTE: To be used for SNR calculations

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Effective transit duration (in days)
"""
function calc_transit_duration_eff_winn2010(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)+size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  arcsin_circ_central = pi/ps.orbit[pl].P*duration_central_circ

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  duration_ratio_for_impact_parameter = calc_transit_duration_eff_factor_for_impact_parameter_b(b,size_ratio)

  # WARNING: This is technically an approximation.  It avoids small angle for non-grazing transits, but does use a variant of the small angle approximation for nearly and grazing transits.
  asin_arg = (arcsin_circ_central * duration_ratio_for_impact_parameter)
  duration = ps.orbit[pl].P/pi * radial_separation_over_a/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)
  duration = duration_central_cric * radial_separation_over_a/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)/arcsin_circ_central
end
calc_transit_duration_eff_winn2010(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_eff_winn2010(t.sys[s],p)

"""
    calc_transit_duration_eff_kipping2010(ps, pl)
    calc_transit_duration_eff_kipping2010(t, s, p)

Calculate effective transit duration of planet across host star disk using the formulation from Kipping (2010)
NOTE: To be used for SNR calculations

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Effective transit duration (in days)
"""
function calc_transit_duration_eff_kipping2010(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = calc_impact_parameter(ps, pl)
  size_ratio = ps.planet[pl].radius/ps.star.radius
  @assert !isnan(b)
  @assert zero(b)<=b
  if b>one(b)+size_ratio
     return zero(b)
  end
  duration_central_circ = calc_transit_duration_central_circ(ps,pl)
  arcsin_circ_central = pi/ps.orbit[pl].P*duration_central_circ

  one_plus_e_sin_w = 1+ecc*sin(ps.orbit[pl].omega)
  sqrt_one_minus_e_sq = sqrt((1+ecc)*(1-ecc))
  vel_fac = sqrt_one_minus_e_sq / one_plus_e_sin_w
  radial_separation_over_a = (1+ecc)*(1-ecc)/one_plus_e_sin_w
  duration_ratio_for_impact_parameter = calc_transit_duration_eff_factor_for_impact_parameter_b(b,size_ratio)

  # WARNING: This is technically an approximation (see Kipping 2010 Eqn 15).  It avoids small angle for non-grazing transits, but does use a variant of the small angle approximation for nearly and grazing transits.
  asin_arg = (arcsin_circ_central * duration_ratio_for_impact_parameter/radial_separation_over_a)
  duration = ps.orbit[pl].P/pi * radial_separation_over_a^2/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)
  #duration = duration_central_circ * radial_separation_over_a^2/sqrt_one_minus_e_sq * (asin_arg < 1.0 ? asin(asin_arg) : 1.0)/arcsin_circ_central
end
calc_transit_duration_eff_kipping2010(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_eff_kipping2010(t.sys[s],p)

"""
    calc_transit_duration_eff(ps, pl)
    calc_transit_duration_eff(t, s, p)

Wrapper function to calculate effective transit duration of planet across host star disk to be used for SNR calculations

# Arguments:
- `ps::PlanetarySystemAbstract`: Planetary system object
- `pl::Integer`: Planet/orbit index in planetary system object
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target

# Returns:
Effective transit duration (in days)
"""
#calc_transit_duration_eff(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_eff_small_angle_approx(ps,pl)
#calc_transit_duration_eff(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_eff_winn2010(ps,pl)
calc_transit_duration_eff(ps::PlanetarySystemAbstract, pl::Integer) = calc_transit_duration_eff_kipping2010(ps,pl)
calc_transit_duration_eff(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration_eff(t.sys[s],p)

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

Calculates (expected) number of transits observed for planet

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target
- `sim_param::SimParam`: Simulation parameter object

# Returns:
Effective (expected) number of transits
"""
function calc_expected_num_transits(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
 period = t.sys[s].orbit[p].P
 exp_num_transits = t.duty_cycle * t.data_span/period
 return exp_num_transits
end

"""
    get_durations_searched_Kepler(period, duration)

Returns appropriate transit duration for use with 1-sigma depth function (or CDPP) given the range of transit durations searched by the Kepler pipeline.

# Arguments:
- `period::Float64`: Orbital period (in days)
- `duration::Float64`: Transit duration (in days)

# Returns:
Valid transit duration searched by Kepler pipeline that is closest to the actual transit duration input to the function.
"""
#function get_legal_durations(period::Float64,duration::Float64)
function get_durations_searched_Kepler(period::Float64,duration::Float64)
  num_dur = length(cdpp_durations) # 14
  min_duration = 0.0
  max_duration = 0.0
  min_periods = [0.5, 0.52, 0.6517, 0.7824, 0.912, 1.178, 1.3056, 1.567, 1.952, 2.343, 2.75, 3.14, 3.257, 3.91]
  max_periods = [50.045, 118.626, 231.69, 400.359, 635.76, 725, 725, 725, 725, 725, 725, 725, 725, 725]
  i = 1
  #determine what maximum and minimum durations were searched for this period
  while min_duration == 0.0 || max_duration == 0.0
    if i > 14
      @warn "No duration searched match this period and duration." period duration
      return (min_duration > 0.0) ? min_duration/24 : 0.0
    end
    if period <= max_periods[i] && min_duration == 0.0
      min_duration = cdpp_durations[i]
    end
    if period >= min_periods[num_dur+1-i] && max_duration == 0.0
      max_duration = cdpp_durations[num_dur+1-i]
    end
    i+=1
  end
  if duration<=max_duration/24 && duration>=min_duration/24
    return duration
  elseif duration>=max_duration/24
    return max_duration/24
  elseif duration <=min_duration/24
    return min_duration/24
  end
end

include("transit_detection_model.jl")
include("transit_prob_geometric.jl")

const has_sc_bit_array_size = 7*8        # WARNING: Must be big enough given value of num_quarters (assumed to be <=17)

mutable struct KeplerTargetObs                        # QUERY:  Do we want to make this type depend on whether the catalog is based on simulated or real data?
  obs::Vector{TransitPlanetObs}
  sigma::Vector{TransitPlanetObs}           # Simplistic approach to uncertainties for now.  QUERY: Should estimated uncertainties be part of Observations type?
  # phys_id::Vector{Tuple{Int32,Int32}}     # So we can lookup the system's properties # Commented out since Not used

  prob_detect::SystemDetectionProbsAbstract  # QUERY: Specialize type of prob_detect depending on whether for simulated or real data?

  has_sc::BitArray{1}                        # Note: Changed from Array{Bool}.  Alternatively, we could try StaticArray{Bool} so fixed size?  Do we even need to keep this?

  star::StarObs
end
#KeplerTargetObs(n::Integer) = KeplerTargetObs( fill(TransitPlanetObs(),n), fill(TransitPlanetObs(),n), fill(tuple(0,0),n),  ObservedSystemDetectionProbsEmpty(),  fill(false,num_quarters), StarObs(0.0,0.0) )
KeplerTargetObs(n::Integer) = KeplerTargetObs( fill(TransitPlanetObs(),n), fill(TransitPlanetObs(),n), ObservedSystemDetectionProbsEmpty(),  falses(has_sc_bit_array_size), StarObs(0.0,0.0,0) )
num_planets(t::KeplerTargetObs) = length(t.obs)

"""
    calc_target_obs_sky_ave(t, sim_param)

Simulate observation of Kepler target (with associated planets) and return observed quantities and detection probabilities of all planets averaging over all viewing angles.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `sim_param::SimParam`: Simulation parameter object, must contain the following values:
    - max_tranets_in_sys: Maximum number of transiting planets per Kepler target
    - transit_noise_model: Function to use for simulating noise in the observed properties of the transiting planets

# Returns:
Kepler observable object containing observed properties and detection probabilities for every planet
"""
function calc_target_obs_sky_ave(t::KeplerTarget, sim_param::SimParam)
   max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
   transit_noise_model = get_function(sim_param,"transit_noise_model")
   min_detect_prob_to_be_included = 0.0  # get_real(sim_param,"min_detect_prob_to_be_included")
   num_observer_samples = 1 # get_int(sim_param,"num_viewing_geometry_samples")
  vetting_efficiency = get_function(sim_param,"vetting_efficiency")

  np = num_planets(t)
  obs = Array{TransitPlanetObs}(undef,np)
  sigma = Array{TransitPlanetObs}(undef,np)
  #id = Array{Tuple{Int32,Int32}}(np)
  #id = Array{Tuple{Int32,Int32}}(np)
  ns = length(t.sys)
  sdp_sys = Array{SystemDetectionProbsAbstract}(undef,ns)
  i = 1
  for (s,sys) in enumerate(t.sys)
    pdet = zeros(num_planets(sys))
    for (p,planet) in enumerate(sys.planet)
      if get(sim_param,"verbose",false)
         println("# s=",s, " p=",p," num_sys= ",length(t.sys), " num_pl= ",num_planets(sys) )
      end
        period = sys.orbit[p].P
        duration_central = calc_transit_duration_central(t,s,p)
        size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
        depth = calc_transit_depth(t,s,p)
        ntr = calc_expected_num_transits(t, s, p, sim_param)

        # cdpp_central = interpolate_cdpp_to_duration(t, duration_central)
        # snr_central = calc_snr_if_transit_cdpp(t, depth, duration_central, cdpp_central, sim_param, num_transit=ntr)
        # pdet_ave = calc_ave_prob_detect_if_transit_from_snr_cdpp(t, snr_central, period, duration_central, size_ratio, cdpp_central, sim_param, num_transit=ntr)

        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
        snr_central = calc_snr_if_transit(t, depth, duration_central, osd_central, sim_param, num_transit=ntr)
        pdet_ave = calc_ave_prob_detect_if_transit_from_snr(t, snr_central, period, duration_central, size_ratio, osd_central, sim_param, num_transit=ntr)

	add_to_catalog = pdet_ave > min_detect_prob_to_be_included  # Include all planets with sufficient detection probability

	if add_to_catalog
           pdet_central = calc_prob_detect_if_transit(t, snr_central, period, duration_central, sim_param, num_transit=ntr)
           threshold_pdet_ratio = rand()
	   hard_max_num_b_tries = 100
	   max_num_b_tries = min_detect_prob_to_be_included == 0. ? hard_max_num_b_tries : min(hard_max_num_b_tries,convert(Int64,1/min_detect_prob_to_be_included))
           # We compute measurement noise based on a single value of b.  We draw from a uniform distribution for b and then using rejection sampling to reduce probability of higher impact parameters
           pdet_this_b = 0.0
           for j in 1:max_num_b_tries
              b = rand()
              transit_duration_factor = calc_effective_transit_duration_factor_for_impact_parameter_b(b,size_ratio)

	          duration = duration_central * transit_duration_factor   # WARNING:  Technically, this duration may be slightly reduced for grazing cases to account for reduction in SNR due to planet not being completely inscribed by star at mid-transit.  But this will be a smaller effect than limb-darkening for grazing transits.  Also, makes a variant of the small angle approximation

              # cdpp = interpolate_cdpp_to_duration(t, duration)
              # snr = snr_central * (cdpp_central/cdpp) * sqrt(transit_duration_factor)

              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
              snr = snr_central * (osd_central/osd)

              pdet_this_b = calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=ntr)
              pvet = vetting_efficiency(t.sys[s].planet[p].radius, period)

              if pdet_this_b >= threshold_pdet_ratio * pdet_central
                  #println("# Adding pdet_this_b = ", pdet_this_b, " pdet_c = ", pdet_central, " snr= ",snr, " cdpp= ",cdpp, " duration= ",duration, " b=",b, " u01= ", threshold_pdet_ratio)
	                 pdet[p] = pdet_ave*pvet
#####
                     duration_central = calc_transit_duration_central(t,s,p)
                     transit_duration_factor = calc_effective_transit_duration_factor_for_impact_parameter_b(b,size_ratio)
                     duration = duration_central * transit_duration_factor   # WARNING:  Makes a variant of the small angle approximation
####
                 obs[i], sigma[i] = transit_noise_model(t, s, p, depth, duration, snr, ntr, b=b)
      	         i += 1
                 break
              end # if add to obs and sigma lists
           end # for j
	else # add_to_catalog
	   # Do anything for planets that are extremely unlikely to be detected even if they were to transit?
	end
    end
    resize!(obs,i-1)
    resize!(sigma,i-1)
    sdp_sys[s] = calc_simulated_system_detection_probs(sys, pdet, max_tranets_in_sys=max_tranets_in_sys, min_detect_prob_to_be_included=min_detect_prob_to_be_included, num_samples=num_observer_samples)
  end
  # TODO SCI DETAIL: Combine sdp_sys to allow for target to have multiple planetary systems
  s1 = findfirst(x->num_planets(x)>0,sdp_sys)  # WARNING IMPORTANT: For now just take first system with planets
  if s1 == nothing
     s1 = 1
  end
  sdp_target = sdp_sys[s1]

  has_no_sc = falses(has_sc_bit_array_size)
  star_obs = StarObs( t.sys[1].star.radius, t.sys[1].star.mass, t.sys[1].star.id )  # NOTE:  This sets the observed star properties to be those in the stellar catalog.  If want to incorporate uncertainty in stellar properties, that would be done elsewhere when translating depths into planet radii.
  #return KeplerTargetObs(obs, sigma, id, sdp_target, has_no_sc, star_obs )
  return KeplerTargetObs(obs, sigma, sdp_target, has_no_sc, star_obs )
end

"""
    calc_target_obs_single_obs(t, sim_param)

Simulate observation of Kepler target (with associated planets) and return observed quantities and detection probabilities of all planets assuming a single viewpoint.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `sim_param::SimParam`: Simulation parameter object, must contain the following values:
    - max_tranets_in_sys: Maximum number of transiting planets per Kepler target
    - transit_noise_model: Function to use for simulating noise in the observed properties of the transiting planets

# Returns:
Kepler observable object containing observed properties and detection probabilities for every planet
"""
function calc_target_obs_single_obs(t::KeplerTarget, sim_param::SimParam)
  # max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
   transit_noise_model = get_function(sim_param,"transit_noise_model")
   min_detect_prob_to_be_included = 0.0  # get_real(sim_param,"min_detect_prob_to_be_included")
  transit_noise_model = get_function(sim_param,"transit_noise_model")

  np = num_planets(t)
  obs = Array{TransitPlanetObs}(undef,np)
  sigma = Array{TransitPlanetObs}(undef,np)
  ns = length(t.sys)
  sdp_sys = Array{ObservedSystemDetectionProbs}(undef,ns)
  i = 1
  cuantos = 1000			#indicator for testing OSD interpolator.
  for (s,sys) in enumerate(t.sys)
    pdet = zeros(num_planets(sys))
    for (p,planet) in enumerate(sys.planet)
      if get(sim_param,"verbose",false)
         println("# s=",s, " p=",p," num_sys= ",length(t.sys), " num_pl= ",num_planets(sys) )
      end
        duration = calc_transit_duration_eff(t,s,p)
	if duration <= 0.
	   continue
	end
        period = sys.orbit[p].P
        ntr = calc_expected_num_transits(t, s, p, sim_param)
        depth = calc_transit_depth(t,s,p)
        # Apply correction to snr if grazing transit
        size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
        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)
        # snr = calc_snr_if_transit_cdpp(t, depth, duration, cdpp, sim_param, num_transit=ntr)

        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
        snr = calc_snr_if_transit(t, depth, duration, osd, sim_param, num_transit=ntr)

        pdet[p] = calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=ntr)

	if pdet[p] > min_detect_prob_to_be_included
           pvet = vetting_efficiency(t.sys[s].planet[p].radius, period)
           pdet[p] *= pvet
           duration = calc_transit_duration(t,s,p)
            obs[i], sigma[i] = transit_noise_model(t, s, p, depth, duration, snr, ntr)
      	   i += 1
	end
    end
    resize!(obs,i-1)
    resize!(sigma,i-1)
    sdp_sys[s] = ObservedSystemDetectionProbs(pdet)
  end
  # TODO SCI DETAIL: Combine sdp_sys to allow for target to have multiple planetary systems
  s1 = findfirst(x->num_planets(x)>0,sdp_sys)  # WARNING: For now just take first system with planets, assumes not two stars wht planets in one target
  if s1 == nothing
     s1 = 1
  end
  sdp_target = sdp_sys[s1]

  has_no_sc = falses(3*num_quarters)
  star_obs = StarObs( t.sys[1].star.radius, t.sys[1].star.mass, t.sys[1].star.id )  # NOTE: This just copies star properties directly
  return KeplerTargetObs(obs, sigma, sdp_target, has_no_sc, star_obs )
end

function test_transit_observations(sim_param::SimParam; verbose::Bool=false)  # TODO TEST: Add more tests
  #transit_param = TransitParameter( EphemerisLinear(10.0, 0.0), TransitShape(0.01, 3.0/24.0, 0.5) )
  generate_kepler_target = get_function(sim_param,"generate_kepler_target")
   max_it = 100000
  local obs
  for i in 1:max_it
    target = generate_kepler_target(sim_param)::KeplerTarget
    while num_planets(target) == 0
      target = generate_kepler_target(sim_param)::KeplerTarget
    end
    #calc_transit_prob_single_planet_one_obs(target,1,1)
    calc_transit_prob_single_planet_obs_ave(target,1,1)
    obs = calc_target_obs_single_obs(target,sim_param)
    obs = calc_target_obs_sky_ave(target,sim_param)
    if verbose && (num_planets(obs) > 0)
      println("# i= ",string(i)," np= ",num_planets(obs), " obs= ", obs )
      break
    end
  end
  return obs
end

randtn() = rand(truncated(Normal(0.0,1.0),-0.999,0.999))

"""
    transit_noise_model_no_noise(t, s, p, depth, duration, snr, num_tr, b = 0.0)

Transit noise model that assumes no additional noise introduced to observations.

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target
- `depth::Float64`: (Fractional) transit depth
- `duration::Float64`: Transit duration (in days)
- `snr::Float64`: Expected multiple event statistic (signal-to-noise ratio) of transit
- `num_tr::Float64`: Expected number of transits observed
- `b::Float64 = 0.0`: Impact parameter of planet's orbit

# Returns
- TransitPlanetObs object containing simulated period, t0, transit depth, and transit duration
- TransitPlanetObs object containing simulated uncertainties for period, t0, transit depth, and transit duration
"""
function transit_noise_model_no_noise(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Float64; b::Float64 = 0.0)
  period = t.sys[s].orbit[p].P
  t0 = period*rand()    # WARNING: Not being calculated from orbit
  sigma_period = 0.0
  sigma_t0 = 0.0
  sigma_depth =  0.0
  sigma_duration =  0.0
  sigma = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration )
  obs = TransitPlanetObs( period, t0, depth,duration)
  return obs, sigma
end

"""
    transit_noise_model_no_noise(t, s, p, depth, duration, snr, num_tr, b = 0.0)

Transit noise model that assumes fixed amounts of uncertainty on observed transit properties. Uncertainties are as follows:
- Period = 10^-6
- t0 = 10^-4
- Depth = 10^-1
- Duration = 10^-2

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target
- `depth::Float64`: (Fractional) transit depth
- `duration::Float64`: Transit duration (in days)
- `snr::Float64`: Expected multiple event statistic (signal-to-noise ratio) of transit
- `num_tr::Float64`: Expected number of transits observed
- `b::Float64 = 0.0`: Impact parameter of planet's orbit

# Returns
- TransitPlanetObs object containing simulated period, t0, transit depth, and transit duration
- TransitPlanetObs object containing simulated uncertainties for period, t0, transit depth, and transit duration
"""
function transit_noise_model_fixed_noise(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Float64; b::Float64 = 0.0)
  period = t.sys[s].orbit[p].P
  t0 = period*rand()    # WARNING: Not being calculated from orbit

  sigma_period = 1e-6
  sigma_t0 = 1e-4
  sigma_depth =  0.1
  sigma_duration =  0.01

  sigma = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration)
  #obs = TransitPlanetObs( period, t0, depth, duration)
  obs = TransitPlanetObs( period*(1.0+sigma.period*randtn()), t0*(1.0+sigma.period*randtn()), depth*(1.0+sigma.depth*randtn()),duration*(1.0+sigma.duration*randtn()))
  return obs, sigma
end

#make_matrix_pos_def_count = 0
function make_matrix_pos_def(A::Union{AbstractArray{T1,2},Symmetric{AbstractArray{T1,2}}}; verbose::Bool = false) where {T1<:Real}
    @assert size(A,1) == size(A,2)
    #global make_matrix_pos_def_count
    A = (typeof(A) <: Symmetric) ? A : Symmetric(A)
    smallest_eigval = eigvals(A,1:1)[1]
    if smallest_eigval > 0.0
        return PDMat(A)
    else
        #make_matrix_pos_def_count += 1
        ridge = 1.01 * abs(smallest_eigval)
        if verbose
            println("# Warning: Adding ridge (",ridge,") to matrix w/ eigenvalue ", smallest_eigval," (#", make_matrix_pos_def_count,").")
        end
        return PDMat(A + Diagonal(ridge*ones(size(A,1))))
    end
end

"""
    transit_noise_model_diagonal(t, s, p, depth, duration, snr, num_tr, b = calc_impact_parameter(t.sys[s], p))

Transit noise model that uses the Fisher information matrix formulation assuming finite integration from Price & Rogers (2014).

# Arguments:
- `t::KeplerTarget`: Kepler target object
- `s::Integer`: Index of star in Kepler target
- `p::Integer`: Index of planet in Kepler target
- `depth::Float64`: (Fractional) transit depth
- `duration::Float64`: Transit duration (in days)
- `snr::Float64`: Expected multiple event statistic (signal-to-noise ratio) of transit
- `num_tr::Float64`: Expected number of transits observed
- `b::Float64 = calc_impact_parameter(t.sys[s], p)`: Impact parameter of planet's orbit

# Returns
- TransitPlanetObs object containing simulated period, t0, transit depth, and transit duration
- TransitPlanetObs object containing simulated uncertainties for period, t0, transit depth, and transit duration
"""
function transit_noise_model_diagonal(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Float64; b::Float64 = calc_impact_parameter(t.sys[s],p) )
    transit_noise_model_price_rogers(t, s, p, depth, duration, snr, num_tr; b=b, diagonal=true )
end

function transit_noise_model_price_rogers(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Float64; b::Float64 = calc_impact_parameter(t.sys[s],p), diagonal::Bool = false )
  period = t.sys[s].orbit[p].P
  t0 = period*rand()    # WARNING: Not being calculated from orbit

	# Use variable names from Price & Rogers
	one_minus_e2 = (1-t.sys[s].orbit[p].ecc)*(1+t.sys[s].orbit[p].ecc)
	a_semimajor_axis = semimajor_axis(t.sys[s],p)

	tau0 = rsol_in_au*t.sys[s].star.radius*period/(a_semimajor_axis*2pi)
	tau0 *= sqrt(one_minus_e2)/(1+t.sys[s].orbit[p].ecc*sin(t.sys[s].orbit[p].omega))
	r = t.sys[s].planet[p].radius/t.sys[s].star.radius
	sqrt_one_minus_b2 = (0.0<=b<1.0) ? sqrt((1-b)*(1+b)) : 0.0
	@assert(sqrt_one_minus_b2>=0.0)
    if(b<1)   # trapezoidal transit shape
 	   T = 2*tau0*sqrt_one_minus_b2
	   tau = 2*tau0*r/sqrt_one_minus_b2
	   delta = depth
   else      # triangular transit shape, TODO: SCI DETAIL: Could improve treatment, but so rare this should be good enough for most purposes not involving EBs
       @assert b<=1+r
       T = 2*tau0*sqrt((1+r+b)*(1+r-b))
       tau = T/2
	   delta = depth/2
    end

	Ttot = period
	I = LC_integration_time      # WARNING: Assumes LC only
	Lambda_eff = LC_rate * num_tr # calc_expected_num_transits(t, s, p, sim_param)
	sigma = interpolate_cdpp_to_duration(t, duration)

	# Price & Rogers Eqn A8 & Table 1 # Thanks to Danley for finding typeos.
	tau3 = tau^3
	I3 = I^3
	a1 = (10*tau3+2*I^3-5*tau*I^2)/tau3
	a2 = (5*tau3+I3-5*tau*tau*I)/tau3
	a3 = (9*I^5*Ttot-40*tau3*I*I*Ttot+120*tau^4*I*(3*Ttot-2*tau))/tau^6
	a4 = (a3*tau^5+I^4*(54*tau-35*Ttot)-12*tau*I3*(4*tau+Ttot)+360*tau^4*(tau-Ttot))/tau^5
	a5 = (a2*(24T*T*(I-3*tau)-24*T*Ttot*(I-3*tau))+tau3*a4)/tau3
	a6 = (3*tau*tau+T*(I-3*tau))/(tau*tau)
	a7 = (-60*tau^4+12*a2*tau3*T-9*I^4+8*tau*I3+40*tau3*I)/(tau^4)
	a8 = (2T-Ttot)/tau
	a9 = (-3*tau*tau*I*(-10*T*T+10*T*Ttot+I*(2*I+5*Ttot))-I^4*Ttot+8*tau*I3*Ttot)/(tau^5)
	a10 = ((a9+60)*tau*tau+10*(-9*T*T+9*T*Ttot+I*(3*I+Ttot))-75*tau*Ttot)/(tau*tau)
	a11 = (I*Ttot-3*tau*(Ttot-2*tau))/(tau*tau)
	a12 = (-360*tau^5-24*a2*tau3*T*(I-3*tau)+9*I^5-35*tau*I^4-12*tau*tau*I3-40*tau3*I*I+360*tau^4*I)/(tau^5)
	a13 = (-3*I3*(8*T*T-8*T*Ttot+3*I*Ttot)+120*tau*tau*T*I*(T-Ttot)+8*tau*I3*Ttot)/tau^5
	a14 = (a13*tau*tau+40*(-3*T*T+3*T*Ttot+I*Ttot)-60*tau*Ttot)/(tau*tau)
	a15 = (2*I-6*tau)/tau
	b1  = (6*I*I-3*I*Ttot+tau*Ttot)/(I*I)
	b2  = (tau*T+3*I*(I-T))/(I*I)
	b3 = (tau3-12*T*I*I+8*I3+20*tau*I*I-8*tau*tau*I)/I3
	b4 = (6*T*T-6*T*Ttot+I*(5*Ttot-4*I))/(I*I)
	b5 = (10*I-3*tau)/I
	b6 = (12*b4*I3+4*tau*(-6*T*T+6T*Ttot+I*(13*Ttot-30*I)))/I3
	b7 = (b6*I^5+4*tau*tau*I*I*(12*I-11*Ttot)+tau3*I*(11*Ttot-6*I)-tau^4*Ttot)/I^5
	b8 = (3T*T-3*T*Ttot+I*Ttot)/(I*I)
	b9 = (8*b8*I^4+20*tau*I*I*Ttot-8*tau*tau*I*Ttot+tau3*Ttot)/I^4
	b10 =  (-tau^4+24*T*I*I*(tau-3I)+60*I^4+52*tau*I3-44*tau*tau*I*I+11*tau3*I)/I^4
	b11 =  (-15*b4*I3+10*b8*tau*I*I+15*tau*tau*(2*I-Ttot))/I3
	b12 =  (b11*I^5+2*tau3*I*(4*Ttot-3*I)-tau^4*Ttot)/I^5
	b13 =  (Ttot-2*T)/I
	b14 =  (6*I-2*tau)/I

        Q = snr/sqrt(num_tr)
        sigma_t0 = tau>=I ?  sqrt(0.5*tau*T/(1-I/(3*tau)))/Q : sqrt(0.5*I*T/(1-tau/(3*I)))/Q
        sigma_period = sigma_t0/sqrt(num_tr)
        sigma_duration = tau>=I ? sigma*sqrt(abs(6*tau*a14/(delta*delta*a5)) /Lambda_eff )  : sigma*sqrt(abs(6*I*b9/(delta*delta*b7)) / Lambda_eff)
        sigma_depth = tau>=I ? sigma*sqrt(abs(-24*a11*a2/(tau*a5)) / Lambda_eff)  : sigma*sqrt(abs(24*b1/(I*b7)) / Lambda_eff)

        sigma_obs = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration )

        local obs
        if diagonal     # Assume uncertainties uncorrelated (Diagonal)
  	         obs = TransitPlanetObs( period*(1.0+sigma_obs.period*randtn()), t0*(1.0+sigma_obs.period*randtn()), depth*(1.0+sigma_obs.depth*randtn()),duration*(1.0+sigma_obs.duration*randtn()))
        else  # TODO WARNING TEST: Should test before using full covariance matrix
            cov = zeros(4,4)
            if tau>=I
	           # cov[0,0] = -3*tau/(delta*delta*a15)
               cov[1,1] = 24*tau*a10/(delta*delta*a5)
               cov[1,2] = cov[2,1] = 36*a8*tau*a1/(delta*delta*a5)
	           cov[1,3] = cov[3,1] = -12*a11*a1/(delta*a5)
	           cov[1,4] = cov[4,1] = -12*a6*a1/(delta*a5)
	           cov[2,2] = 6*tau*a14/(delta*delta*a5)
	           cov[2,3] = cov[3,2] = 72*a8*a2/(delta*a5)
	           cov[2,4] = cov[4,2] = 6*a7/(delta*a5)
	           cov[3,3] = -24*a11*a2/(tau*a5)
	           cov[3,4] = cov[4,3] = -24*a6*a2/(tau*a5)
	           cov[4,4] = a12/(tau*a5)
           else
	           # cov[0,0] = 3*I/(delta*delta*b14)
	           cov[1,1] = -24*I*I*b12/(delta*delta*tau*b7)
	           cov[1,2] = cov[2,1] = 36*I*b13*b5/(delta*delta*b7)
	           cov[1,3] = cov[3,1] = 12*b5*b1/(delta*b7)
	           cov[1,4] = cov[4,1] = 12*b5*b2/(delta*b7)
	           cov[2,2] = 6*I*b9/(delta*delta*b7)
	           cov[2,3] = cov[3,2] = 72*b13/(delta*b7)
	           cov[2,4] = cov[4,2] = 6*b3/(delta*b7)
	           cov[3,3] = 24*b1/(I*b7)
	           cov[3,4] = cov[4,3] = 24*b2/(I*b7)
	           cov[4,4] = b10/(I*b7)
	       end
	       cov .*= sigma*sigma/Lambda_eff
	       cov = make_matrix_pos_def(cov)
           obs_dist = MvNormal(zeros(4),cov)

	       local obs_period, obs_duration, obs_depth, sigma_period, sigma_duration, sigma_depth
           isvalid = false
	       while !isvalid
	          obs_vec = rand(obs_dist)
              obs_duration = duration + obs_vec[2]
              obs_depth = depth + obs_vec[3]
              if (obs_duration>0.0) && (obs_depth>0.0)
                  isvalid = true
	          end
           end # while
     	   obs = TransitPlanetObs( period*(1.0+sigma_obs.period*randn()), t0*(1.0+sigma_obs.t0*randn()), obs_depth, obs_duration)
        end
  	return obs, sigma_obs
end


text:
using IncrementalInference

"""
    ScalarPointToMultiPoint

Factor modeling an uncertain association between a point and a set of other
points. Measurement model operates on the difference between the measured
field value at the point and the interpolated field value (from the point
set).
"""
struct ScalarPointToMultiPoint{T<:SamplableBelief}
    # nothing
end

function IIF.getSample(s::CalcFactor{<:ScalarPointToMultiPoint}, N::Int=1)
    # generate 0-mean 1D  samples
    return (reshape(rand(s.factor.Z, N), 1, N), )
end

function (s::CalcFactor{<:ScalarPointToMultiPoint})(point, points)
    # calc difference between field value at point and
    # interpolated field value from points

    h = Interpolations.LinearInterpolation((points[1]..., points[2]...),points[3]...)

    return point[3] - h(point[1],point[2])
end

ext: jl


text:
using Random

"Minibatched data"
mutable struct Data; x; y; batchsize; length; partial; indices; shuffle; xsize; ysize; xtype; ytype; end

"""
    minibatch(x, [y], batchsize; shuffle, partial, xtype, ytype, xsize, ysize)

Return an iterable of minibatches [(xi,yi)...] given data tensors x, y
and batchsize.  The last dimension of x and y should match and give
the number of instances. `y` is optional.  Keyword arguments:

- `shuffle=false`: Shuffle the instances before minibatching.
- `partial=false`: If true include the last partial minibatch < batchsize.
- `xtype=typeof(x)`: Convert xi in minibatches to this type.
- `ytype=typeof(y)`: Convert yi in minibatches to this type.
- `xsize=size(x)`: Convert xi in minibatches to this shape.
- `ysize=size(y)`: Convert yi in minibatches to this shape.
"""
function minibatch(x,y,batchsize; shuffle=false,partial=false,xtype=typeof(x),ytype=typeof(y),xsize=size(x), ysize=size(y))
    nx = size(x)[end]
    if nx != size(y)[end]; throw(DimensionMismatch()); end
    x2 = reshape(x, :, nx)
    y2 = reshape(y, :, nx)
    # indices = shuffle ? randperm(nx) : 1:nx --> do this at the beginning of for loop
    Data(x2,y2,batchsize,nx,partial,1:nx,shuffle,xsize,ysize,xtype,ytype)
end

function minibatch(x,batchsize; shuffle=false,partial=false,xtype=typeof(x),xsize=size(x))
    nx = size(x)[end]
    x2 = reshape(x, :, nx)
    # indices = shuffle ? randperm(nx) : 1:nx --> do this at the beginning of for loop
    Data(x2,nothing,batchsize,nx,partial,1:nx,shuffle,xsize,nothing,xtype,nothing)
end

function Base.iterate(d::Data,i=0)     # returns data in d.indices[i+1:i+batchsize]
    if i == 0 && d.shuffle; d.indices = randperm(d.length); end
    if i >= d.length || (!d.partial && i + d.batchsize > d.length); return nothing; end
    j=min(i+d.batchsize, d.length)
    ids = d.indices[i+1:j]
    xbatch = convert(d.xtype, reshape(d.x[:,ids],d.xsize[1:end-1]...,length(ids)))
    if d.y == nothing
        return (xbatch,j)
    else
        ybatch = convert(d.ytype, reshape(d.y[:,ids],d.ysize[1:end-1]...,length(ids)))
        return ((xbatch,ybatch),j)
    end
end

function Base.length(d::Data)
    n = d.length / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

function Random.rand(d::Data)
    i = rand(0:(d.length-d.batchsize))
    return next(d, i)[1]
end


text:
import Base: AbstractCmd, CmdRedirect
using OrderedCollections: OrderedDict

const BATCH_ENVS = (
    "AWS_BATCH_JOB_ID" => "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9",
    "AWS_BATCH_JQ_NAME" => "HighPriority"
)

const SUBMIT_JOB_RESP = Dict(
    "jobName" => "example",
    "jobId" => "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9",
)

const REGISTER_JOB_DEF_RESP = Dict(
    "jobDefinitionName" => "sleep60",
    "jobDefinitionArn" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1",
    "revision"=>1,
)

const DESCRIBE_JOBS_DEF_RESP = Dict(
    "jobDefinitions" => [
        Dict(
            "type" => "container",
            "containerProperties" => Dict(
                "command" => [
                    "sleep",
                    "60"
                ],
                "environment" => [

                ],
                "image" => "busybox",
                "memory" => 128,
                "mountPoints" => [

                ],
                "ulimits" => [

                ],
                "vcpus" => 1,
                "volumes" => [

                ],
                "jobRoleArn" => "arn:aws:iam::012345678910:role/sleep60",
            ),
            "jobDefinitionArn" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1",
            "jobDefinitionName" => "sleep60",
            "revision" => 1,
            "status" => "ACTIVE"
        )
    ]
)

const DESCRIBE_JOBS_RESP = Dict(
    "jobs" => [
        Dict(
            "container" => Dict(
                "command" => [
                    "sleep",
                    "60"
                ],
                "containerInstanceArn" => "arn:aws:ecs:us-east-1:012345678910:container-instance/5406d7cd-58bd-4b8f-9936-48d7c6b1526c",
                "environment" => [

                ],
                "exitCode" => 0,
                "image" => "busybox",
                "memory" => 128,
                "mountPoints" => [

                ],
                "ulimits" => [

                ],
                "vcpus" => 1,
                "volumes" => [

                ],
                "jobRoleArn" => "arn:aws:iam::012345678910:role/sleep60",
            ),
            "createdAt" => 1480460782010,
            "dependsOn" => [

            ],
            "jobDefinition" => "sleep60",
            "jobId" => "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9",
            "jobName" => "example",
            "jobQueue" => "arn:aws:batch:us-east-1:012345678910:job-queue/HighPriority",
            "parameters" => Dict(

            ),
            "startedAt" => 1480460816500,
            "status" => "SUCCEEDED",
            "stoppedAt" => 1480460880699
        )
    ]
)

function describe_compute_environments_patch(output::Vector=[])
    @patch function AWSBatch.describe_compute_environments(d::Dict)
        compute_envs = d["computeEnvironments"]
        @assert length(compute_envs) == 1
        ce = first(compute_envs)

        key = startswith(ce, "arn:") ? "computeEnvironmentArn" : "computeEnvironmentName"
        results = filter(d -> d[key] == ce, output)
        OrderedDict("computeEnvironments" => results)
    end
end

function describe_compute_environments_patch(output::OrderedDict)
    describe_compute_environments_patch([output])
end

function describe_job_queues_patch(output::Vector=[])
    @patch function AWSBatch.describe_job_queues(d::Dict)
        queues = d["jobQueues"]
        @assert length(queues) == 1
        queue = first(queues)

        key = startswith(queue, "arn:") ? "jobQueueArn" : "jobQueueName"
        results = filter(d -> d[key] == queue, output)
        OrderedDict("jobQueues" => output)
    end
end

function describe_job_queues_patch(output::OrderedDict)
    describe_job_queues_patch([output])
end

function log_events_patches(; log_stream_name="mock_stream", events=[], exception=nothing)
    job_descriptions = if log_stream_name === nothing
        Dict("jobs" => [Dict("container" => Dict())])
    else
        Dict("jobs" => [Dict("container" => Dict("logStreamName" => log_stream_name))])
    end

    get_log_events_patch = if exception !== nothing
        @patch get_log_events(; kwargs...) = throw(exception)
    else
        @patch function get_log_events(; kwargs...)
            if get(kwargs, :nextToken, nothing) === nothing
                Dict("events" => events, "nextForwardToken" => "0")
            else
                Dict("events" => [], "nextForwardToken" => "0")
            end
        end
    end

    return [
        @patch describe_jobs(args...; kwargs...) = job_descriptions
        get_log_events_patch
    ]
end


text:
module CUDNN

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, libraryPropertyType
using ..CUDA: libcudnn, @retry_reclaim, isdebug

using CEnum

using Memoize

using DataStructures

import NNlib

# core library
include("libcudnn_common.jl")
include("error.jl")
include("libcudnn.jl")
include("libcudnn_deprecated.jl")

# low-level wrappers
include("util.jl")
include("base.jl")
include("tensor.jl")
include("conv.jl")
include("pooling.jl")
include("activation.jl")
include("filter.jl")
include("softmax.jl")
include("batchnorm.jl")
include("dropout.jl")
include("rnn.jl")

# high-level integrations
include("nnlib.jl")

include("compat.jl")

function math_mode(mode=CUDA.math_mode())
    if mode == CUDA.PEDANTIC_MATH
        # don't use tensor cores.
        # on A100, only use them for TF32
        CUDNN_DEFAULT_MATH
    elseif mode == CUDA.DEFAULT_MATH
        # allow tensor core usage
        CUDNN_TENSOR_OP_MATH
    elseif mode == CUDA.FAST_MATH
        # also downcast inputs
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION
    end
end

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,cudnnHandle_t}}()

# cache for created, but unused handles
const old_handles = DefaultDict{CuContext,Vector{cudnnHandle_t}}(()->cudnnHandle_t[])

function handle()
    CUDA.detect_state_changes()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUDNN, ctx)) do
            handle = if isempty(old_handles[ctx])
                cudnnCreate()
            else
                pop!(old_handles[ctx])
            end

            finalizer(current_task()) do task
                push!(old_handles[ctx], handle)
            end
            # TODO: cudnnDestroy to preserve memory, or at exit?

            cudnnSetStream(handle, stream())

            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

@inline function set_stream(stream::CuStream)
    ctx = context()
    tls = task_local_storage()
    handle = get(tls, (:CUDNN, ctx), nothing)
    if handle !== nothing
        cudnnSetStream(handle, stream)
    end
    return
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end
end

function log_message(sev, udata, dbg_ptr, cstr)
    # "Each line of this message is terminated by \0, and the end of the message is
    # terminated by \0\0"
    len = 0
    while true
        if unsafe_load(cstr, len+1) == '\0' && unsafe_load(cstr, len+2)
            break
        end
        len += 1
    end
    str = unsafe_string(cstr, len)
    lines = split(str, '\0')
    msg = join(str, '\n')

    # TODO: inspect `sev` to generate an appropriate message (@debug, @info, etc)
    dbg = unsafe_load(dbg_ptr)
    println(msg)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUDNN
    # FIXME: this doesn't work, and the mask remains 0 (as observed with cudnnGetCallback)
    if isdebug(:init, CUDNN)
        callback = @cfunction(log_message, Nothing,
                              (cudnnSeverity_t, Ptr{Cvoid}, Ptr{cudnnDebug_t}, Cstring))
        cudnnSetCallback(typemax(UInt32), C_NULL, callback)
    end
end

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Heatex_NTU
	Heatex_NTU()=begin
		new(
			Heatex_Basic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Type of Heat Exchanger",
				:Valid=>["Counter Flow","Cocurrent Flow", "Shell and Tube"],
				:Default=>"Cocurrent Flow"
			)),
			NTU_Basic(Dict{Symbol,Any}(
				:Brief=>"NTU Method of Calculation",
				:Symbol=>" "
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Fraction in Hot Side",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar Fraction in Hot Side",
				:Hidden=>true
			)),(NComp)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Fraction in Hot Side",
				:Hidden=>true
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Fraction in Cold Side",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar Fraction in Cold Side",
				:Hidden=>true
			)),(NComp)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Fraction in Cold Side",
				:Hidden=>true
			)),
			[
				:([vh, xh, yh] = _base_1.PP.Flash(_base_1.InletHot.T, _base_1.InletHot.P, _base_1.InletHot.z)),
				:([vc, xc, yc] = _base_1.PP.Flash(_base_1.InletCold.T, _base_1.InletCold.P, _base_1.InletCold.z)),
				:(Method.NTU*Method.Cmin = U*A),
				:(Method.Cmin = min([Method.Ch,Method.Cc])),
				:(Method.Cmax = max([Method.Ch,Method.Cc])),
				:(Method.Cr = Method.Cmin/Method.Cmax),
				:(Q = Method.Eft*Method.Cmin*(_base_1.InletHot.T-_base_1.InletCold.T)),
				:(Method.Ch = _base_1.InletHot.F*((1-_base_1.InletHot.v)*_base_1.PP.LiquidCp(0.5*_base_1.InletHot.T+0.5*_base_1.OutletHot.T,0.5*_base_1.InletHot.P+0.5*_base_1.OutletHot.P,xh)+ _base_1.InletHot.v*_base_1.PP.VapourCp(0.5*_base_1.InletHot.T+0.5*_base_1.OutletHot.T,0.5*_base_1.InletHot.P+0.5*_base_1.OutletHot.P,yh))),
				:(Method.Cc = _base_1.InletCold.F*((1-_base_1.InletCold.v)*_base_1.PP.LiquidCp(0.5*_base_1.InletCold.T+0.5*_base_1.OutletCold.T,0.5*_base_1.InletCold.P+0.5*_base_1.OutletCold.P,xc)+ _base_1.InletCold.v*_base_1.PP.VapourCp(0.5*_base_1.InletCold.T+0.5*_base_1.OutletCold.T,0.5*_base_1.InletCold.P+0.5*_base_1.OutletCold.P,yc))),
				:(Method.Eft1 = 1),
				:(Method.Eft = 1-exp(-Method.NTU)),
				:(Method.Eft = (1-exp(-Method.NTU*(1+Method.Cr)))/(1+Method.Cr)),
				:(Method.Eft = Method.NTU/(1+Method.NTU)),
				:(Method.Eft = (1-exp(-Method.NTU*(1-Method.Cr)))/(1-Method.Cr*exp(-Method.NTU*(1-Method.Cr)))),
				:(Method.Eft = 2*(1+Method.Cr+sqrt(1+Method.Cr^2)*((1+exp(-Method.NTU*sqrt(1+Method.Cr^2)))/(1-exp(-Method.NTU*sqrt(1+Method.Cr^2)))) )^(-1)),
			],
			[
				"Flash Calculation in Hot Side","Flash Calculation in Cold Side","Number of Units Transference","Minimum Heat Capacity","Maximum Heat Capacity","Thermal Capacity Ratio","Duty","Hot Stream Average Heat Capacity","Cold Stream Average Heat Capacity","Effectiveness Correction","Effectiveness","Effectiveness in Cocurrent Flow","Effectiveness in Counter Flow","Effectiveness in Counter Flow","TEMA E Shell Effectiveness",
			],
			[:ExchangerType,],
			[:Method,:xh,:yh,:vh,:xc,:yc,:vc,]
		)
	end
	_base_1::Heatex_Basic
	ExchangerType::DanaSwitcher
	Method::NTU_Basic
	xh::Array{fraction}
	yh::Array{fraction}
	vh::fraction
	xc::Array{fraction}
	yc::Array{fraction}
	vc::fraction
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Heatex_NTU
function setEquationFlow(in::Heatex_NTU)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	if Method.Cr == 0
		addEquation(11)
	else
		let switch=ExchangerType
			if switch=="Cocurrent Flow"
				addEquation(12)
			elseif switch=="Counter Flow"
				if Method.Cr == 1
					addEquation(13)
				else
					addEquation(14)
				end
			elseif switch=="Shell and Tube"
				addEquation(15)
			end
		end
	end
end
function atributes(in::Heatex_NTU,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Heatex"
	fields[:Brief]="Simplified model for Heat Exchangers"
	fields[:Info]="This model perform material and heat balance using the NTU-Effectiveness Approach.
This shortcut calculation does not require exchanger configuration or geometry data.

== Assumptions ==
* Steady-State operation;
* No heat loss to the surroundings.

== Specify ==
* The Inlet streams: Hot and Cold.

== References ==
[1] E.A.D. Saunders, Heat Exchangers: Selection, Design and
 Construction, Longman, Harlow, 1988.

"
	drive!(fields,_)
	return fields
end
Heatex_NTU(_::Dict{Symbol,Any})=begin
	newModel=Heatex_NTU()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
mutable struct OperationsProblemTemplate
    transmission::Type{<:PM.AbstractPowerModel}
    devices::Dict{Symbol, DeviceModel}
    branches::Dict{Symbol, DeviceModel}
    services::Dict{Symbol, ServiceModel}
end

"""
    OperationsProblemTemplate(::Type{T}) where {T<:PM.AbstractPowerFormulation}
Creates a model reference of the Power Formulation, devices, branches, and services.
# Arguments
- `model::Type{T<:PM.AbstractPowerFormulation}`:
- `devices::Dict{Symbol, DeviceModel}`: device dictionary
- `branches::Dict{Symbol, BranchModel}`: branch dictionary
- `services::Dict{Symbol, ServiceModel}`: service dictionary
# Example
```julia
template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
```
"""
function OperationsProblemTemplate(::Type{T}) where {T <: PM.AbstractPowerModel}
    return OperationsProblemTemplate(
        T,
        Dict{Symbol, DeviceModel}(),
        Dict{Symbol, DeviceModel}(),
        Dict{Symbol, ServiceModel}(),
    )
end

OperationsProblemTemplate() = OperationsProblemTemplate(PM.AbstractPowerModel)


text:
export ELFSections, ELFSection, ELFSection32, ELFSection64, ELFSectionRef

# Export certain extensions to the base API
export section_type

"""
    ELFSections

ELF Section header table type, containing information about the number of
sections within the ELF object, the location of the section headers, etc...
"""
struct ELFSections{H <: ELFHandle} <: Sections{ELFHandle}
    handle::H
end

Sections(oh::ELFHandle) = ELFSections(oh)

# Implement Sections API
handle(sections::ELFSections) = sections.handle
lastindex(sections::ELFSections) = header(handle(sections)).e_shnum

"""
    ELFSection

ELF Section type, containing information about a `Section` within the ELF
object, such as the `Section`'s name, its size, etc...
"""
abstract type ELFSection{H <: ELFHandle} <: Section{H} end
@io struct ELFSection32{H <: ELFHandle} <: ELFSection{H}
    sh_name::UInt32
    sh_type::UInt32
    sh_flags::UInt32
    sh_addr::UInt32
    sh_offset::UInt32
    sh_size::UInt32
    sh_link::UInt32
    sh_info::UInt32
    sh_addralign::UInt32
    sh_entsize::UInt32
end
@io struct ELFSection64{H <: ELFHandle} <: ELFSection{H}
    sh_name::UInt32
    sh_type::UInt32
    sh_flags::UInt64
    sh_addr::UInt64
    sh_offset::UInt64
    sh_size::UInt64
    sh_link::UInt32
    sh_info::UInt32
    sh_addralign::UInt64
    sh_entsize::UInt64
end

# Implement Section API
function read(oh::ELFHandle, ::Type{ST}) where {ST <: ELFSection}
    type_func = section_header_type
    size_func = section_header_size
    return read_struct(oh, type_func, size_func, "Section Header")
end

section_name(s::ELFSection) = string("strtab@", s.sh_name)
section_type(s::ELFSection) = s.sh_type
section_size(s::ELFSection) = s.sh_size
section_offset(s::ELFSection) = s.sh_offset
section_address(s::ELFSection) = s.sh_addr

"""
    ELFSectionRef

ELF `SectionRef` type, containing an `ELFSection` and important metadata.
"""
struct ELFSectionRef{H<:ELFHandle} <: SectionRef{H}
    sections::ELFSections{H}
    section::ELFSection{H}
    idx::UInt32
end

function SectionRef(sections::ELFSections, s::ELFSection, idx)
    return ELFSectionRef(sections, s, UInt32(idx))
end
deref(s::ELFSectionRef) = s.section
sections(s::ELFSectionRef) = s.sections
handle(s::ELFSectionRef) = handle(sections(s))
section_number(s::ELFSectionRef) = s.idx
@derefmethod section_type(s::ELFSectionRef)

function section_name(section::ELFSectionRef)
    # Actually do the strtab lookup to get the name of this guy
    strtab = StrTab(handle(section))
    return strtab_lookup(strtab, deref(section).sh_name)
end

"""
    section_type_string(s::ELFSection)

Return the given `ELFSection`'s section type as a string.
"""
function section_type_string(s::ELFSection)
    global SHT_TYPES

    sh_type = section_type(s)
    if haskey(SHT_TYPES, sh_type)
        return SHT_TYPES[sh_type]
    end
    return string("Unknown Section Type (0x", string(sh_type, base=16), ")")
end
@derefmethod section_type_string(s::ELFSectionRef)


text:
#=
An example of Type unions
=#

mutable struct Point2D
    x::Float64
    y::Float64
end

mutable struct Vector2D
    x::Float64
    y::Float64
end

VecOrPoint2D = Union{Vector2D,Point2D}

import Base.+

function +(u::VecOrPoint2D, v::VecOrPoint2D)
    @show typeof(u),typeof(v)
    if typeof(u) == Vector2D && typeof(v) == Vector2D
        Vector2D(u.x+v.x,u.y+v.y)
    else
        Point2D(u.x+v.x,u.y+v.y)
    end
end

p = Point2D(2,4)
q = Point2D(4,3)
v = Vector2D(3,2)
w = Vector2D(-1,2)

@show typeof(p+q)
@show typeof(p+v)
@show typeof(w+p)
@show typeof(v+w)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using SailRoute, Test

@testset "Min angle" begin
    @test SailRoute.min_angle(10.0, 20.0) ≈ 10.0
    @test SailRoute.min_angle(-10.0, 30.0) ≈ 40.0
    @test SailRoute.min_angle(-30.0, 10.0) ≈ 40.0
end

@testset "Resolve vectors" begin
    ws = SailRoute.w_c(0.0, 10.0, 0.0, 0.0)
    @test ws[1] ≈ 10.0
    @test ws[2] ≈ 0.0
end


text:
using Statistics:quantile

"""
    classify_percentile(threshold)

Create a percentile-based classifiction function that converts `scores_train::Scores` and `scores_test::Scores` to an
array of classes with `"normal"` indicating normal data and `"outlier"` indicating outliers. The conversion is based on
percentiles of the training data, i.e. all datapoints above the `threshold` percentile are considered outliers.

Parameters
----------
    threshold::Real
The score threshold (number between 0 and 1) used to classify the samples into inliers and outliers.

    scores_train::AbstractVector{<:Real}
A vector of training scores, typically the result of [`fit`](@ref) with a detector.

    scores_test::AbstractVector{<:Real}
A vector of test scores, typically the result of [`score`](@ref) using a previously fitted detector.

Returns
----------
    classes::AbstractVector{String}
The vector of classes consisting of `"outlier"` and `"normal"` elements.

Examples
----------
    classify = classify_percentile(0.9)
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    classify(scores_train, scores_train) # [1, 1, -1]
    classify(scores_train, scores_test) # [-1, -1, 1, 1, 1]
"""
function classify_percentile(threshold::Real)
    function percentile(scores_train::Scores, scores_test::Scores)::Tuple{Labels, Labels}
        @assert 0 < threshold < 1
        t = quantile(scores_train, threshold)
        f = scores -> ifelse.(scores .> t, CLASS_OUTLIER, CLASS_NORMAL)
        f(scores_train), f(scores_test)
    end
    return percentile
end


text:
module PropertyGraphs
export AbstractPropertyGraph, PropertyGraph, SymmetricPropertyGraph,
  gprops, vprops, eprops, get_gprop, get_vprop, get_eprop,
  set_gprop!, set_vprop!, set_eprop!, set_gprops!, set_vprops!, set_eprops!

using ...Present, ...CSetDataStructures, ..BasicGraphs
using ..BasicGraphs: TheoryGraph, TheorySymmetricGraph
import ..BasicGraphs: nv, ne, src, tgt, inv, edges, vertices,
  has_edge, has_vertex, add_edge!, add_edges!, add_vertex!, add_vertices!

# Data types
############

""" Abstract type for graph with properties.

Concrete types are [`PropertyGraph`](@ref) and [`SymmetricPropertyGraph`](@ref).
"""
abstract type AbstractPropertyGraph{T} end

@present TheoryPropertyGraph <: TheoryGraph begin
  Props::AttrType
  vprops::Attr(V,Props)
  eprops::Attr(E,Props)
end

@abstract_acset_type __AbstractPropertyGraph <: HasGraph

const _AbstractPropertyGraph{T} = __AbstractPropertyGraph{S, Tuple{Dict{Symbol,T}}} where {S}

@acset_type __PropertyGraph(TheoryPropertyGraph, index=[:src,:tgt]) <: __AbstractPropertyGraph

const _PropertyGraph{T} = __PropertyGraph{Dict{Symbol,T}}

""" Graph with properties.

"Property graphs" are graphs with arbitrary named properties on the graph,
vertices, and edges. They are intended for applications with a large number of
ad-hoc properties. If you have a small number of known properties, it is better
and more efficient to create a specialized C-set type using `@acset_type`.

See also: [`SymmetricPropertyGraph`](@ref).
"""
struct PropertyGraph{T,G<:_AbstractPropertyGraph{T}} <: AbstractPropertyGraph{T}
  graph::G
  gprops::Dict{Symbol,T}
end

PropertyGraph{T,G}(; kw...) where {T,G<:_AbstractPropertyGraph{T}} =
  PropertyGraph(G(), Dict{Symbol,T}(kw...))
PropertyGraph{T}(; kw...) where T = PropertyGraph{T,_PropertyGraph{T}}(; kw...)

@present TheorySymmetricPropertyGraph <: TheorySymmetricGraph begin
  Props::AttrType
  vprops::Attr(V,Props)
  eprops::Attr(E,Props)

  compose(inv,eprops) == eprops # Edge involution preserves edge properties.
end

@abstract_acset_type __AbstractSymmetricPropertyGraph <: HasGraph

const _AbstractSymmetricPropertyGraph{T} = __AbstractSymmetricPropertyGraph{S, Tuple{Dict{Symbol,T}}} where {S}

@acset_type __SymmetricPropertyGraph(TheorySymmetricPropertyGraph, index=[:src,:tgt]) <:
  __AbstractSymmetricPropertyGraph

const _SymmetricPropertyGraph{T} = __SymmetricPropertyGraph{Dict{Symbol,T}}

""" Symmetric graphs with properties.

The edge properties are preserved under the edge involution, so these can be
interpreted as "undirected" property (multi)graphs.

See also: [`PropertyGraph`](@ref).
"""
struct SymmetricPropertyGraph{T,G<:_AbstractSymmetricPropertyGraph{T}} <:
    AbstractPropertyGraph{T}
  graph::G
  gprops::Dict{Symbol,T}
end

SymmetricPropertyGraph{T,G}(; kw...) where {T,G<:_AbstractSymmetricPropertyGraph{T}} =
  SymmetricPropertyGraph(G(), Dict{Symbol,T}(kw...))
SymmetricPropertyGraph{T}(; kw...) where T =
  SymmetricPropertyGraph{T,_SymmetricPropertyGraph{T}}(; kw...)

# Accessors and mutators
########################

""" Graph-level properties of a property graph.
"""
gprops(g::AbstractPropertyGraph) = g.gprops

""" Properties of vertex in a property graph.
"""
vprops(g::AbstractPropertyGraph, v) = subpart(g.graph, v, :vprops)

""" Properties of edge in a property graph.
"""
eprops(g::AbstractPropertyGraph, e) = subpart(g.graph, e, :eprops)

""" Get graph-level property of a property graph.
"""
get_gprop(g::AbstractPropertyGraph, key::Symbol) = gprops(g)[key]

""" Get property of vertex or vertices in a property graph.
"""
get_vprop(g::AbstractPropertyGraph, v, key::Symbol) =
  broadcast(v) do v; vprops(g,v)[key] end

""" Get property of edge or edges in a property graph.
"""
get_eprop(g::AbstractPropertyGraph, e, key::Symbol) =
  broadcast(e) do e; eprops(g,e)[key] end

""" Set graph-level property in a property graph.
"""
set_gprop!(g::AbstractPropertyGraph, key::Symbol, value) =
  (gprops(g)[key] = value)

""" Set property of vertex or vertices in a property graph.
"""
set_vprop!(g::AbstractPropertyGraph, v, key::Symbol, value) =
  broadcast(v, value) do v, value; vprops(g,v)[key] = value end

""" Set property of edge or edges in a property graph.
"""
set_eprop!(g::AbstractPropertyGraph, e, key::Symbol, value) =
  broadcast(e, value) do e, value; eprops(g,e)[key] = value end

""" Set multiple graph-level properties in a property graph.
"""
set_gprops!(g::AbstractPropertyGraph; kw...) = merge!(gprops(g), kw)
set_gprops!(g::AbstractPropertyGraph, d::AbstractDict) = merge!(gprops(g), d)

""" Set multiple properties of a vertex in a property graph.
"""
set_vprops!(g::AbstractPropertyGraph, v::Int; kw...) = merge!(vprops(g,v), kw)
set_vprops!(g::AbstractPropertyGraph, v::Int, d::AbstractDict) =
  merge!(vprops(g,v), d)

""" Set multiple properties of an edge in a property graph.
"""
set_eprops!(g::AbstractPropertyGraph, e::Int; kw...) = merge!(eprops(g,e), kw)
set_eprops!(g::AbstractPropertyGraph, e::Int, d::AbstractDict) =
  merge!(eprops(g,e), d)

@inline nv(g::AbstractPropertyGraph) = nv(g.graph)
@inline ne(g::AbstractPropertyGraph) = ne(g.graph)
@inline src(g::AbstractPropertyGraph, args...) = src(g.graph, args...)
@inline tgt(g::AbstractPropertyGraph, args...) = tgt(g.graph, args...)
@inline inv(g::SymmetricPropertyGraph, args...) = inv(g.graph, args...)
@inline vertices(g::AbstractPropertyGraph) = vertices(g.graph)
@inline edges(g::AbstractPropertyGraph) = edges(g.graph)
@inline has_vertex(g::AbstractPropertyGraph, v::Int) = has_vertex(g.graph, v)
@inline has_edge(g::AbstractPropertyGraph, e::Int) = has_edge(g.graph, e)

add_vertex!(g::AbstractPropertyGraph{T}; kw...) where T =
  add_vertex!(g, Dict{Symbol,T}(kw...))
add_vertex!(g::AbstractPropertyGraph{T}, d::Dict{Symbol,T}) where T =
  add_part!(g.graph, :V, vprops=d)

add_vertices!(g::AbstractPropertyGraph{T}, n::Int) where T =
  add_parts!(g.graph, :V, n, vprops=[Dict{Symbol,T}() for i=1:n])

add_edge!(g::AbstractPropertyGraph{T}, src::Int, tgt::Int; kw...) where T =
  add_edge!(g, src, tgt, Dict{Symbol,T}(kw...))

# Non-symmetric case.

add_edge!(g::PropertyGraph{T}, src::Int, tgt::Int, d::Dict{Symbol,T}) where T =
  add_part!(g.graph, :E, src=src, tgt=tgt, eprops=d)

function add_edges!(g::PropertyGraph{T}, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}, eprops=nothing) where T
  @assert (n = length(srcs)) == length(tgts)
  if isnothing(eprops)
    eprops = [Dict{Symbol,T}() for i=1:n]
  end
  add_parts!(g.graph, :E, n, src=srcs, tgt=tgts, eprops=eprops)
end

# Symmetric case.

add_edge!(g::SymmetricPropertyGraph{T}, src::Int, tgt::Int,
          d::Dict{Symbol,T}) where T =
 add_edges!(g, src:src, tgt:tgt, [d])

function add_edges!(g::SymmetricPropertyGraph{T}, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}, eprops=nothing) where T
  @assert (n = length(srcs)) == length(tgts)
  if isnothing(eprops)
    eprops = [ Dict{Symbol,T}() for i=1:n ]
  end
  invs = nparts(g.graph, :E) .+ [(n+1):2n; 1:n]
  eprops = [eprops; eprops] # Share dictionaries to ensure equal properties.
  add_parts!(g.graph, :E, 2n, src=[srcs; tgts], tgt=[tgts; srcs],
             inv=invs, eprops=eprops)
end

# Constructors from graphs
##########################

function PropertyGraph{T}(g::HasGraph, make_vprops, make_eprops;
                          gprops...) where T
  pg = PropertyGraph{T}(; gprops...)
  add_vertices!(pg, nv(g))
  add_edges!(pg, src(g), tgt(g))
  for v in vertices(g)
    set_vprops!(pg, v, make_vprops(v))
  end
  for e in edges(g)
    set_eprops!(pg, e, make_eprops(e))
  end
  pg
end

PropertyGraph{T}(g::HasGraph; gprops...) where T =
  PropertyGraph{T}(g, v->Dict(), e->Dict(); gprops...)

function SymmetricPropertyGraph{T}(g::HasGraph, make_vprops, make_eprops;
                                   gprops...) where T
  pg = SymmetricPropertyGraph{T}(; gprops...)
  add_vertices!(pg, nv(g))
  for v in vertices(g)
    set_vprops!(pg, v, make_vprops(v))
  end
  for e in edges(g)
    if !has_subpart(g, :inv) || e <= inv(g,e)
      e1, e2 = add_edge!(pg, src(g,e), tgt(g,e))
      set_eprops!(pg, e1, make_eprops(e))
    end
  end
  pg
end

SymmetricPropertyGraph{T}(g::HasGraph; gprops...) where T =
  SymmetricPropertyGraph{T}(g, v->Dict(), e->Dict(); gprops...)

end


text:
using PyPlot, Bridge, Colors
import PyPlot: PyObject
using Bridge: runmean

#simname = "sumgamma"
oneplot = false
withbeta = transdim
bstr = ["", "b"][1 + withbeta]
params = readdlm(joinpath("output", simname, "params.txt"), Float64; skipstart=1)[:,2:end];

withrho1 = N > 0 && rhosigma[1] > 0

function PyObject(t::Color)
    trgb = convert(RGB, t)
    ctup = map(float, (red(trgb), green(trgb), blue(trgb)))
    o = PyPlot.PyObject(ctup)
    return o
end
cd = 40
cd = 30
trqs(x) = 1-sqrt(1-x)

i = 1
dcol1 = PyObject(RGB(trqs.(Bridge._viridis[cd*i])...))
col1 = PyObject(RGB((Bridge._viridis[cd*i])...))
lcol1 = PyObject(RGB(sqrt.(Bridge._viridis[cd*i])...))
i = 2
dcol2 = PyObject(RGB(trqs.(Bridge._viridis[cd*i])...))
col2 = PyObject(RGB((Bridge._viridis[cd*i])...))
lcol2 = PyObject(RGB(sqrt.(Bridge._viridis[cd*i])...))
i = 3
dcol3 = PyObject(RGB(trqs.(Bridge._viridis[cd*i])...))
col3 = PyObject(RGB((Bridge._viridis[cd*i])...))
lcol3 = PyObject(RGB(sqrt.(Bridge._viridis[cd*i])...))

dcol, col, lcol = dcol2, col2, lcol3

if simid == 2
    alpha1 = 2.0
    alpha2 = alpha1/10
    beta1 = 0.4
    beta2 = beta1/10
end
#=
beta = beta1 + beta2

alpha = (beta1*alpha1 + beta2*alpha2)/(beta1 + beta2)
@show alpha
=#
######## alpha and beta

if oneplot

    figure(figsize=(8,5))
    subplot(231)
else
    if withbeta
        figure(figsize=(8,5))
        subplot(121)
    else
        figure(figsize=(4,5))
    end
end

@show n = size(params, 1)
N = div(size(params, 2)-2,2)
skip = max(iterations÷1_000, 1)

plot(params[:,1], color=lcol , lw = 0.5, label=latexstring("\\alpha"))
abar = runmean(params[:,1])
plot(abar, color=dcol , lw = 0.5, label=latexstring("\\bar\\alpha"))

#=
annotate(latexstring("\\alpha"),
    xy=[div(2n, 3); params[1+div(2n, 3), 1]],
    xytext=[0, -10],
    textcoords="offset points",
    fontsize=10.0,
    ha="right",
    va="bottom")
annotate(latexstring("\\bar\\alpha"),
    xy=[n; abar[end]],
    xytext=[10,0],
    textcoords="offset points",
    fontsize=10.0,
    ha="right",
    va="bottom")=#

plot(1:n, fill(alpha0, n), ":", color=:darkorange , lw = 1.0, label=L"\varnothing\alpha")
legend()
grid(linestyle=":", axis="y")

if withbeta
    if oneplot
        subplot(234)
    else
        subplot(122)
    end

    plot(params[:,2], color=lcol , lw = 0.5, label=latexstring("\\beta"))
    bbar = runmean(params[:,2])
    plot(bbar, color=dcol , lw = 0.5, label=latexstring("\\bar\\beta"))

    #=
    annotate(latexstring("\\beta"),
        xy=[div(2n, 3); params[1+div(2n, 3), 2]],
        xytext=[0, -10],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
    annotate(latexstring("\\bar\\beta"),
        xy=[n; bbar[end]],
        xytext=[10,0],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
    =#
    if simid == 2
        plot(1:n, fill(beta0, n), ":", color=:darkorange , lw = 1.0, label=L"\beta_1+\beta_2")
    elseif simid ==3
        plot(1:n, fill(beta0, n), ":", color=:darkorange , lw = 1.0, label=L"\hat\beta")
    end
    legend()
    grid(linestyle=":", axis="y")

end

if oneplot
    subplot(132)
else
    savefig(joinpath("output", simname, "traceplot1$(bstr).pdf"))
    figure()
    subplot(121)
end

######## theta and rho (added or not)

for i in 1:N
    dcol = PyObject(RGB(trqs.(Bridge._viridis[cd*i])...))
    col = PyObject(RGB((Bridge._viridis[cd*i])...))
    lcol = PyObject(RGB(sqrt.(Bridge._viridis[cd*i])...))

    plot(params[:, 2 + i], color=lcol, lw = 0.2)
    thbar = runmean(params[:, 2 + i])
    plot(thbar, color=col , lw = 0.5)
    annotate(latexstring("\\theta_$i"),
        xy=[div(2n, 3); params[div(2n, 3), 2 + i]],
        xytext=[0, -10],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
    annotate(latexstring("\\bar\\theta_$i"),
        xy=[n; thbar[end]],
        xytext=[10,0],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
end

#legend()
grid(linestyle=":", axis="y")

if oneplot
    subplot(133)
else
    subplot(122)
end

for i in (2-withrho1):N
    dcol = PyObject(RGB(trqs.(Bridge._viridis[cd*i])...))
    col = PyObject(RGB((Bridge._viridis[cd*i])...))
    lcol = PyObject(RGB(sqrt.(Bridge._viridis[cd*i])...))

    plot(params[:, 2+N + i], color = lcol, lw = 0.2)
    rbar = runmean(params[:, 2+N + i])
    plot(rbar, color = col, lw = 0.5)
    annotate(latexstring("\\rho_$i"),
        xy=[div(2n, 3); params[div(2n, 3), 2+N+i]],
        xytext=[0, -10],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
    annotate(latexstring("\\bar\\rho_$i"),
        xy=[n; rbar[end]],
        xytext=[0,0],
        textcoords="offset points",
        fontsize=10.0,
        ha="right",
        va="bottom")
end
grid(linestyle=":", axis="y")

savefig(joinpath("output", simname, "traceplot2$(bstr).pdf"))

truth = readdlm(joinpath("output", simname, "truth.txt"), header=true)

b = truth[1][6:6+N-1] # bin boundaries

function θ(x, params, b)
    N = length(b)
    N == 0 && return 0.0
    x <= b[1] && return 0.0
    x > b[N] && return params[2 + N] * x + params[end]
    k = first(searchsorted(b, x)) - 1
    params[2 + k] * x + params[2 + N + k]
end

theta0(x, alpha1, beta1, alpha2, beta2) = -log((beta1*exp(-alpha1*x) + beta2*exp(-alpha2*x))/(beta1+beta2)) - (beta1*alpha1 + beta2*alpha2)/(beta1 + beta2)*x
phi0(x, alpha1, beta1, alpha2, beta2) = -log((beta1*exp(-alpha1*x) + beta2*exp(-alpha2*x))/(beta1+beta2))

if N > 0
xx = (0.002:0.01:2)*b[end]
else
    xx = (0.002:0.01:8)
end

thetahat(xx, p) = [θ(x, p, b) for x in xx]
logxvhat(xx, p) = [θ(x, p, b) + p[1]*x - log(p[2]) for x in xx]
logxv0(x, alpha1, beta1, alpha2, beta2) = -log((beta1*exp(-alpha1*x) + beta2*exp(-alpha2*x)))
logxv0(x, alpha1, beta1) = -log(beta1) + alpha1*x

#theta2(x) = θ(x, phat + theta[1] - alpha, b) + alpha*x

#println([quantile(params[:,1] + params[:,4], q) for q in [0.05, 0.5, 0.95]])
n1 = div(n, 4)
#n1 = 20000
n2 = n
phat = median(params[n1:n2,:], 1)

A = zeros(length(xx), length(n1:skip:n2))
i_ = 1
for i in n1:skip:n2
    A[:, i_] = logxvhat(xx, params[i, :])
    i_ += 1
end

upper = mapslices(v-> quantile(v,0.975), A, 2)
med = median(A, 2)
lower = mapslices(v-> quantile(v,0.025), A, 2)

figure(figsize=(8,5))
#subplot(121)
fill_between(xx, upper[:], lower[:], edgecolor=lcol2, facecolor=(lcol2..., 0.2), hatch="X", label="95% cred. band")
#plot(xx, med, label = "marg. post. median", color=:blue)
if simid == 2
    plot(xx, logxv0.(xx, alpha1, beta1, alpha2, beta2), label=L"-\log(x v_0(x))", color=:darkorange)
elseif simid == 3
    sa, sb = [0.0376041, 0.1005920]*1.96

    plot(xx, logxv0.(xx,  0.5889420, 1.8127392 *365/7 ), label=L"-\log(x \hat v(x))", color=:darkorange)
  #  fill_between(xx,  logxv0.(xx, 0.5889420+sa,(1.8127392 -sb)*365/7), logxv0.(xx, 0.5889420-sa ,(1.8127392 + sb)*365/7), edgecolor=:orange, facecolor=:None, label="freq. band")
    #plot(xx, logxv0.(xx, (phat[1] + phat[end-N]), (phat[2] * exp(-phat[end]))))
end
legend()
savefig(joinpath("output", simname, "bands$(bstr)1.pdf"))
#savefig(joinpath("output", simname, "bands$(bstr)1.svg"))

z = sort(diff(X.yy))
figure();plot(z, range(0, stop=1, length=length(z)))
#plot(z, cdf.(Gamma( var(z)/mean(z), mean(z)^2/var(z)), z))
plot(z, cdf.(Gamma(beta0*dt, 1/alpha), z), label="Gamma ML")
plot(z, cdf.(Gamma(phat[2]*dt, 1/phat[1]), z), label="Gamma(beta,alpha)")
#plot(z, cdf.(Gamma(params[end,2]*dt, 1/params[end,1]), z), label="Gamma(beta,alpha)")
N > 1 && try
    plot(z, cdf.(Gamma((phat[2] * exp(-phat[end-1]))*dt, 1/(phat[1] + phat[end-N-1])), z), label="Gamma mid")
catch
end
plot(z, cdf.(Gamma((phat[2] * exp(-phat[end]))*dt, 1/(phat[1] + phat[end-N])), z), label="Gamma end")
simid==2 && axis([0.0, 6.0, 0.4, 1.0])
legend()
#=
figure();plot(z, linspace(0,1,length(z))-cdf.(Gamma(  var(z)/mean(z), mean(z)^2/var(z)), z))
=#

#=
betahat = var(z)/mean(z)/dt
alphahat = 1/(mean(z)^2/var(z))

smooth(z, hw, D = Normal()) = conv(z, normalize(pdf.(D, linspace(-std(D)*3,std(D)*3, 2*hw)), 1))[hw:end-hw]
sz = smooth(z, 20)[1:end-20]
figure();plot(z[2:end-20], 1/length(z)./diff(sz), linewidth=0.2)
G = Gamma(  var(z)/mean(z), mean(z)^2/var(z))
plot(z, pdf.(G, z))
=#
#=
Plot Levy density of end piece
plot(xx, -llevyx.(GammaProcess((phat[2] * exp(-phat[end])), (phat[1] + phat[end-N])), xx))
=#

z = sort(diff(X.yy))
figure();plot(z, range(0, stop=1, length=length(z)))
plot(z, cdf.(Gamma(beta0*dt, 1/alpha), z), label="Gamma ML")
plot(z, cdf.(Gamma(phat[2]*dt, 1/phat[1]), z), label="Gamma(beta,alpha)")
if N > 0
    i = searchsortedfirst(z, b[1]/2)
else
    i = length(z)÷2
end
axis([0.0, z[i], 0, i/length(z)])
legend()

ps = mean(diff(params,1).!=0,1)
ps *= 5
ps[2] /= 4
ps


text:
include("..\\src\\MathProgComplex.jl")
using Test

using Ipopt

rosenbrock = Problem()

x = Variable("x", Real)
y = Variable("y", Real)
set_objective!(rosenbrock, Base.power_by_squaring((1-x),2) + 100*Base.power_by_squaring((y-Base.power_by_squaring(x,2)),2))

mysolver = IpoptSolver(print_level = 0)
m, variables_jump, ctr_jump, ctr_exp = get_JuMP_cartesian_model(rosenbrock, mysolver)
print(m)
solve(m)

sol = get_JuMP_solution(m, variables_jump, rosenbrock)

@test sol[x] ≈ 1.0 atol=1e-5
@test sol[y] ≈ 1.0 atol=1e-5
@test get_objective(rosenbrock, sol) ≈ 0.0 atol=1e-5