text:
verificationstats=[:cor, :betamb, :spear, #=:kend, =# :N]
verificationlabels=Dict{Symbol, Any}(
  :cor=>"Pearson Cor",
    :cort=>"t-stat",
    :corhac=>"HAC z",
    :cormbbz=>"MBB z",
    :cormbbp=>"MBB p",
    :corar1h0z=>"AR(1) \$H_0\$ z",
    :corar1h0p=>"AR(1) \$H_0\$ p",
  :spear=>"\\midrule Spearman Cor",
    :speart=>"t-stat",
    :spearmbbz=>"MBB z",
    :spearmbbp=>"MBB p",
    :spearar1h0z=>"AR(1) \$H_0\$ z",
    :spearar1h0p=>"AR(1) \$H_0\$ p",
  :kend=>"\\midrule Kendall's Tau",
    :kendz=>"z-stat",
    :kendmbbz=>"MBB z",
    :kendmbbp=>"MBB p",
    :kendar1h0z=>"AR(1) \$H_0\$ z",
    :kendar1h0p=>"AR(1) \$H_0\$ p",
  :betabm=>"\\midrule Beta", #bm=benchmark v measure
    :betabmmw=>"HC z",
    :betabmhac=>"HAC z",
    :betabmmbbz=>"MBB z",
    :betabmmbbp=>"MBB p",
    :betabmar1h0z=>"AR(1) \$H_0\$ z",
    :betabmar1h0p=>"AR(1) \$H_0\$ p",
  :betamb=>"\\midrule Beta", #bm=benchmark v measure
    :betambmw=>"HC z",
    :betambhac=>"HAC z",
    :betambmbbz=>"MBB z",
    :betambmbbp=>"MBB p",
    :betambar1h0z=>"AR(1) \$H_0\$ z",
    :betambar1h0p=>"AR(1) \$H_0\$ p",
  :N=>"\\midrule Months",
  :comom=>Dict(
    :lG_WXplmomentum_LP366d=>"Mom", #industry adjusted
    :lGVP_WXplmomentum_LP366d=>"Mom Participation",
    :lGP_WXplmomentum_LP366d=>"Mom Pt-Weighted",
    :lG_WXplmomentum121s_LP366d=>"Mom", #non-industry adjusted
    :lGVP_WXplmomentum121s_LP366d=>"Mom Participation",
    :lGP_WXplmomentum121s_LP366d=>"Mom Pt-Weighted",
    :M_WXplmomentum121s_LP366d=>"\$Mom/|ME|\$",
    :GlG_WXplmomentum121s_LP366d=>"\$\\Delta\$ Mom",
    :GlGP_WXplmomentum121s_LP366d=>"\$\\Delta\$ Mom Pt-Weighted",
    :GM_WXplmomentum121s_LP366d=>"\$\\Delta Mom/|ME|\$",
  ),
  :measureadj=>Dict(
#    :lG_A_aumequityXw_WXmomentum121_LP366d=>"\$w=\\beta\$",
    :lG_A_aumequityXabsw_WXmomentum121s_LP366d=>"\$w=|\\beta_i|\$",
    :lG_A_aumequityXabscorw_WXmomentum121s_LP366d=>"\$w=|\\beta_i|_{>0.1}\$",
  ),
  :umdadj=>Dict(
  #  :lG_A_aumequityXw_F_ff3m_umd_LP366d=>"\$w=\\beta\$",
    :lG_A_aumequityXabsw_F_ff3m_umd_LP366d=>"\$w=|\\beta_i|\$",
    :lG_A_aumequityXabscorw_F_ff3m_umd_LP366d=>"\$w=|\\beta_i|_{>0.1}\$",
  ),
)

verificationsubstats=Dict{Symbol, Vector{Symbol}}(
  :cor=>[#=:cort, =#:corhac, :cormbbz, :cormbbp, :corar1h0z, :corar1h0p],
  :spear=>[#=:speart,=# :spearmbbz, :spearmbbp, :spearar1h0z, :spearar1h0p],
  :kend=>[#=:kendz,=# :kendmbbz, :kendmbbp, :kendar1h0z, :kendar1h0p],
  :betabm=>[#=:betamw,=# :betabmhac, :betabmmbbz, :betabmmbbp, :betabmar1h0z, :betabmar1h0p],
  :betamb=>[#=:betamw,=# :betambhac, :betambmbbz, :betambmbbp, :betambar1h0z, :betambar1h0p],
  :N=>Symbol[])

verificationformats = Dict(
  :cor=>(;scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:"),
  :cort=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :corhac=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :cormbbz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :cormbbp=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),
  :corar1h0z=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :corar1h0p=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),

  :spear=>(;scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:"),
  :speart=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :spearmbbz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :spearmbbp=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}") ,
  :spearar1h0z=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :spearar1h0p=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),

  :kend=>(;scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:"),
  :kendz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :kendmbbz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :kendmbbp=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),
  :kendar1h0z=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :kendar1h0p=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),

  :betabm=>(;scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:"),
  :betabmmw=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betabmhac=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betabmmbbz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betabmmbbp=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),
  :betabmar1h0z=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betabmar1h0p=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),

  :betamb=>(;scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:"),
  :betambmw=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betambhac=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betambmbbz=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betambmbbp=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),
  :betambar1h0z=>(;scalefactor=1.0, decimals=2, prefix="\\mathit{", suffix="}\\:\\:"),
  :betambar1h0p=>(;scalefactor=1.0, decimals=3, prefix="\\mathit{", suffix="}"),
  :N=>(;scalefactor=1.0, decimals=0, prefix="", suffix=""),
)

verificationrows = [(s->[s;verificationsubstats[s];]).(verificationstats)...;]

#=mfbenchmarkcols = [
  :lG_A_aumequityXw_WXmomentum121_LP366d
  :lG_A_aumequityXw_F_ff3m_umd_LP366d
  :lG_A_aumequityXabsw_WXmomentum121_LP366d,
  :lG_A_aumequityXabscorw_WXmomentum121_LP366d,
]=#

#=comombenchmarkcols = [ #NOTE- these are the industry adjusted versions
  :lG_WXplmomentum_LP366d,
  :lGVP_WXplmomentum_LP366d,
  :lGP_WXplmomentum_LP366d,]=#
stdcomombenchmarkcols = [
  :lG_WXplmomentum121_LP366d,
  :lGP_WXplmomentum121_LP366d,
  :M_WXplmomentum121_LP366d,
  ]
Gcomombenchmarkcols = [
  :GlG_WXplmomentum121_LP366d,
  :GlGP_WXplmomentum121_LP366d,
  :GM_WXplmomentum121_LP366d,
  ]

fundbenchmarkcols = Dict(
  :umdadj=>[
    #:lG_A_aumequityXw_F_ff3m_umd_LP366d,
    :lG_A_aumequityXabsw_F_ff3m_umd_LP366d,
    :lG_A_aumequityXabscorw_F_ff3m_umd_LP366d,],
  :measureadj => [
    #:lG_A_aumequityXw_WXmomentum121_LP366d,
    :lG_A_aumequityXabsw_WXmomentum121s_LP366d,
    :lG_A_aumequityXabscorw_WXmomentum121s_LP366d
    ])

fundmeasurecol = :lG_WXmomentum121s_LP366d

#create a table of the comomomentum results
function comomentumtables()
  comomentumtable(
    comomtablenamesuffix="",
    comombenchmarkcols=stdcomombenchmarkcols,
    benchmark=:comom)
  comomentumtable(
    comomtablenamesuffix="G",
    comombenchmarkcols=Gcomombenchmarkcols,
    benchmark=:Gcomom)
end
function comomentumtable(;
    tablepath = PARAM[:tablepath],
    measurefilename = PARAM[:tabcomommeasurefilename],
    rid = PARAM[:tabcomomrid],
    resultsname = PARAM[:tabcomomresultsname],
    analysispath = PARAM[:analysispath],
    comomtablename = PARAM[:tabcomomtablename],
    comomtablenamesuffix,
    comombenchmarkcols,
    benchmark)

  comom = CSV.File("$analysispath\\$measurefilename\\$rid\\$(rid)_$resultsname.csv") |> DataFrame
  comom.measure = comom.measure .|> Symbol
  comom.benchmark = comom.benchmark .|> Symbol
  gcomom = groupby(comom, [:measure, :benchmark]) #this effectivelly keys the df

  #get the headers
  colnames::Vector{Vector{String}} = [comombenchmarkcols .|> c->verificationlabels[:comom][c]]
  descrownames::Vector{String} = verificationrows .|> c->verificationlabels[c]

  Ncols = maximum(length.(colnames))
  Nrows = length(descrownames)
  desccontent = [Vector{String}(undef, Ncols) for i ∈ 1:Nrows]

  #extracts the values for a column of the output table from a row of the input
  function verrow2str(r::DataFrameRow)
    outcol = Vector{String}()
    for s ∈ verificationrows
      f = verificationformats[s]
      valstring = num2str(r[s], f.decimals, scalefactor=f.scalefactor)
      push!(outcol, "$(f.prefix)$(valstring)$(f.suffix)")
    end
    return outcol
  end
  #  (s->num2str(r[s], verificationdecimals[s])).(verificationrows)

  function verrow2str(r::AbstractDataFrame)
    @assert nrow(r) == 1
    verrow2str(r[1,:])
  end

  #fill in the content, double checking each row and column entry
  for (c,col) ∈ enumerate(comombenchmarkcols)
    @assert colnames[1][c] === verificationlabels[:comom][col]
    colcontent = verrow2str(gcomom[(;measure=col, benchmark=benchmark)])
    for (r, valstring) ∈ enumerate(colcontent)
      @assert descrownames[r] === verificationlabels[verificationrows[r]]
      desccontent[r][c] = valstring
    end
  end

  #now compile the content
  t = textable(;colnames, descrownames, desccontent)

  writetextable(t, path=tablepath, outname="$(comomtablenamesuffix)$(comomtablename).tex")

  return nothing
end

function fundtables()

  fundtable(fundadjtype=:umdadj)
  fundtable(fundadjtype=:measureadj)
end

#this creates a table of verification results
#I seperate by the weighting adjustment, either French's UMD or
  #returns of the actual momentum measure
function fundtable(;
    tablepath = PARAM[:tablepath],
    measurefilename = PARAM[:tabfundmeasurefilename],
    rid = PARAM[:tabfundrid],
    analysispath = PARAM[:analysispath],
    fundmeasurecol::Symbol=fundmeasurecol,
    resultsname::String=PARAM[:tabfundresultsname],
    fundadjtype::Symbol,
    fundtablename::String = PARAM[:tabfundtablename],
    fundbenchmarkcols::Vector{Symbol} = fundbenchmarkcols[fundadjtype],
    )

  #load the results
  fund = CSV.File("$analysispath\\$measurefilename\\$rid\\$(rid)_$resultsname.csv") |> DataFrame
  fund.label = fund.label .|> Symbol
  fund.measure = fund.measure .|> Symbol
  fund.benchmark = fund.benchmark .|> Symbol
  gfund = groupby(fund, [:label, :measure, :benchmark]) #this effectivelly keys the df

  #get the headers
  colnames::Vector{Vector{String}} = [["\$\\beta\$-Weighted MF AUM Growth",
    "\$\\beta\$-Weighted HF AUM Growth"],
    [fundbenchmarkcols .|> c->verificationlabels[fundadjtype][c];
      fundbenchmarkcols .|> c->verificationlabels[fundadjtype][c]]]
  widthcolnames = [[length(fundbenchmarkcols), length(fundbenchmarkcols)],
    ones(Int, length(fundbenchmarkcols) * 2)]
  alignmentcolnames = [["c" for c ∈ 1:length(colnames[1])],
    ["r" for c ∈ 1:length(colnames[2])]]
  descrownames::Vector{String} = verificationrows .|> c->verificationlabels[c]
  alignmentstring= "l | rr | rr"

  Ncols = maximum(length.(colnames))
  Nfundcols = length(fundbenchmarkcols)
  Nrows = length(descrownames)
  desccontent = [Vector{String}(undef, Ncols) for i ∈ 1:Nrows]

  #extracts the values for a column of the output table from a row of the input
  function verrow2str(r::DataFrameRow)
    outcol = Vector{String}()
    for s ∈ verificationrows
      f = verificationformats[s]
      valstring = num2str(r[s], f.decimals, scalefactor=f.scalefactor)
      push!(outcol, "$(f.prefix)$(valstring)$(f.suffix)")
    end
    return outcol
  end
  #  (s->num2str(r[s], verificationdecimals[s])).(verificationrows)

  function verrow2str(r::AbstractDataFrame)
    @assert nrow(r) == 1
    verrow2str(r[1,:])
  end

  #fill in the content, double checking each row and column entry
  #for this table, we are populating both mutual funds and hedge funds
  #but each fund type uses the same benchmark column adjustment
  for (cmf,col) ∈ enumerate(fundbenchmarkcols)
    chf = cmf + Nfundcols
    @assert colnames[2][cmf] === verificationlabels[fundadjtype][col]
    @assert colnames[2][chf] === verificationlabels[fundadjtype][col]
    colmfcontent = verrow2str(gfund[(;label=:mf, measure=fundmeasurecol, benchmark=col)])
    colhfcontent = verrow2str(gfund[(;label=:hf, measure=fundmeasurecol, benchmark=col)])

    #now iterate over the columns of content
    for (r, (mfvalstring, hfvalstring)) ∈ enumerate(zip(colmfcontent, colhfcontent))
      @assert descrownames[r] === verificationlabels[verificationrows[r]]
      desccontent[r][cmf] = mfvalstring
      desccontent[r][chf] = hfvalstring
    end
  end

  #now compile the content
  t = textable(;colnames, descrownames, desccontent,
    widthcolnames, alignmentstring, alignmentcolnames)

  writetextable(t, path=tablepath, outname="$fundtablename-$fundadjtype.tex")
end

#function univariatelltable(;measurecols, returncols)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Tilings
using Luxor

@draw begin
	tile = PolyTile([10*n*Point(real(exp(im*2*pi*x/n)), imag(exp(im*2*pi*x/n))) for x in 0:(n-1)], fill=false, stroke=false)
	draw(tile)
end

ext: jl


text:
# blpapi_defs.h

const BLPAPI_DATETIME_YEAR_PART         = 0x01
const BLPAPI_DATETIME_MONTH_PART        = 0x02
const BLPAPI_DATETIME_DAY_PART          = 0x04
const BLPAPI_DATETIME_OFFSET_PART       = 0x08
const BLPAPI_DATETIME_HOURS_PART        = 0x10
const BLPAPI_DATETIME_MINUTES_PART      = 0x20
const BLPAPI_DATETIME_SECONDS_PART      = 0x40
const BLPAPI_DATETIME_MILLISECONDS_PART = 0x80
const BLPAPI_DATETIME_FRACSECONDS_PART  = 0x80
const BLPAPI_DATETIME_DATE_PART         = BLPAPI_DATETIME_YEAR_PART | BLPAPI_DATETIME_MONTH_PART| BLPAPI_DATETIME_DAY_PART

const BLPAPI_DATETIME_TIME_PART         = BLPAPI_DATETIME_HOURS_PART | BLPAPI_DATETIME_MINUTES_PART | BLPAPI_DATETIME_SECONDS_PART
const BLPAPI_DATETIME_TIMEMILLI_PART    = BLPAPI_DATETIME_TIME_PART  | BLPAPI_DATETIME_MILLISECONDS_PART
const BLPAPI_DATETIME_TIMEFRACSECONDS_PART = BLPAPI_DATETIME_TIME_PART | BLPAPI_DATETIME_FRACSECONDS_PART

const DATETIME_WITH_TIMEZONE = BLPAPI_DATETIME_DATE_PART | BLPAPI_DATETIME_TIMEMILLI_PART | BLPAPI_DATETIME_OFFSET_PART

has_part(blp_datetime::BLPDateTime, part::UInt8) = blp_datetime.parts & part == part

function Dates.Date(blp_datetime::BLPDateTime)
    # checks that parts mask implies a pure Date value
    @assert blp_datetime.parts == BLPAPI_DATETIME_DATE_PART "BLPDateTime is not a Date. (parts mask = $P)."
    return Dates.Date(blp_datetime.year, blp_datetime.month, blp_datetime.day)
end

function Dates.Time(blp_datetime::BLPDateTime)
    @assert blp_datetime.parts & (BLPAPI_DATETIME_DATE_PART + BLPAPI_DATETIME_OFFSET_PART) == 0
    return Dates.Time(blp_datetime.hours, blp_datetime.minutes, blp_datetime.seconds, blp_datetime.milliSeconds)
end

function Base.show(io::IO, dt::BLPDateTime)

    @printf(io, "%04d-%02d-%02dT%02d:%02d:%02d.%03d", dt.year, dt.month, dt.day, dt.hours, dt.minutes, dt.seconds, dt.milliSeconds)

    if dt.offset >= 0
        print(io, '+')
    else
        print(io, '-')
    end

    abs_offset = abs(dt.offset)
    hours = div(abs_offset, 60)
    minutes = abs_offset - hours*60
    @printf(io, "%02d:%02d", hours, minutes)
end

function BLPDateTime(dt::Dates.DateTime, offset::Integer)
    return BLPDateTime(
            DATETIME_WITH_TIMEZONE,
            Dates.hour(dt),
            Dates.minute(dt),
            Dates.second(dt),
            Dates.millisecond(dt),
            Dates.month(dt),
            Dates.day(dt),
            Dates.year(dt),
            Int16(offset)
        )
end

function BLPDateTime(dt::Dates.DateTime)
    return BLPDateTime(
            BLPAPI_DATETIME_DATE_PART | BLPAPI_DATETIME_TIMEMILLI_PART,
            Dates.hour(dt),
            Dates.minute(dt),
            Dates.second(dt),
            Dates.millisecond(dt),
            Dates.month(dt),
            Dates.day(dt),
            Dates.year(dt),
            0
        )
end

function BLPDateTime(dt::Dates.Date)
    return BLPDateTime(
            BLPAPI_DATETIME_DATE_PART,
            0,
            0,
            0,
            0,
            Dates.month(dt),
            Dates.day(dt),
            Dates.year(dt),
            0
        )
end

BLPDateTime(x::BLPDateTime) = x

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
# forward models of he diffusion
# define constants

const Rjoules = 8.314 # gas constant
const sec_in_yrs = 3.1558e7
const lambda_f = 8.46e-17/sec_in_yrs
const lambda_38 = 1.55125 * 1e-10/sec_in_yrs
const lambda_35 = 9.8485*1e-10/sec_in_yrs # Steiger Jaeger values
const lambda_32 = 4.948*1e-11/sec_in_yrs # see discussion in Schoene, 10
const τ38 = 1.0./lambda_38
const τ35 = 1.0./lambda_35
const τ32 = 1.0./lambda_32

const Na = 6.02214e23  # Avogadro's number
const atomic_mass_U235 = 235.04393 # g/mol
const atomic_mass_U238 = 238.050791  # g/mol
const atomic_mass_Th232 = 232.0380553  # g/mol, CIAAW value

"""

Calculate radioactive ingrowth and diffusion of He4 using modified eqns
 from Meesters & Dunai, 2002 for single crystal.

# Arguments:
- `alpha::Float64`: α -  1/5 fit factors from
  [Flowers et al. 2009](https://doi.org/10.1016/j.gca.2009.01.015) and
  [Ketcham et al. (2007)](https://doi.org/10.2138/am.2007.2281)
- `c0`: 2/5 fit factor: Fanning Curvilinear fit
- `c1`: 3/5 fit factor: Fanning Curvilinear fit
- `c2`: 4/5 fit factor: Fanning Curvilinear fit
- `c3`: 5/5 fit factor: Fanning Curvilinear fit
- `rmr0`: reduced mean length of B2 apatite, see Ketcham et al., 2007
- `eta_q`: product η*q; η - fission-track etching efficiency, q - scaling factor

"""
function rdaam_forward_diffusion(alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,psi,omega,Etrap,
  R,Ea,logD0_a2,n_iter,U238,U238_V,U235,U235_V,Th232,Th232_V,L,
  tT...)
  tTcopy = collect(tT)
  times = tTcopy[1:ceil(Int,length(tT)/2)] # type unstable but fastest so far
  T = tTcopy[ceil(Int,length(tT)/2)+1:end]

  kappa = 1.04 - rmr0
  rdaam2nd_root = 0.5274 # cutoff bellow which rho_r values can be negative
  rdaam2nd_root_cutoff = (rdaam2nd_root^(1/kappa))*(1-rmr0)+rmr0 # also prevents
  # square root term from being negative

  # diffusivity mod

  # preallocate
  (N_t_segs,e_rho_s,rcb2,rho_r,D0,F,zeta,dzeta,dfdchi) = preallocate_diff_data(T) # 1st arg is only for rad dam

  n_times = length(times)
  dt_equiv = 0.0

U238_0 = U238*(exp(lambda_38*times[1])) # U238 is measured at present
U235_0 = U235*(exp(lambda_35*times[1])) # U238 is measured at present
Th232_0 = Th232*(exp(lambda_32*times[1])) # U238 is measured at present

"""
Damage calculation

"""
for ind in eachindex(times)
  if ind>1
      t_dam = times[n_times-ind+1]-times[n_times-ind+2]+dt_equiv
      rcb2_base = (c0 + c1*(log(t_dam)-c2)/(log((1.0/T[n_times-ind+1]))-c3))
      rcb2_fill_val = ((rcb2_base^(1.0/alpha))+1)^(-1)
      rcb2[ind,2] = rcb2_fill_val
      for jj in 1:(ind-1)
          if rcb2[ind,2]>=rdaam2nd_root_cutoff
            rho_r[n_times-ind+jj,n_times-jj] = ((rcb2_fill_val-rmr0)/(1.0-rmr0))^kappa
            rho_r_copy = rho_r[n_times-ind+jj,n_times-jj] # otherwise if AND elseif

            if rho_r_copy>=0.765
                rho_r[n_times-ind+jj,n_times-jj] = 1.6*rho_r[n_times-ind+jj,n_times-jj] - 0.6
            elseif rho_r_copy<0.765
                rho_r[n_times-ind+jj,n_times-jj] = 9.205*rho_r[n_times-ind+jj,n_times-jj]^2 - 9.157*rho_r[n_times-ind+jj,n_times-jj] + 2.269
            end
          end
      end
      # use concept of equivalent time
    if ind<n_times
      dt_equiv = exp((((((1.0/rcb2[ind,2])-1.0)^alpha-c0)*(log(1.0/T[n_times-ind])-c3))/c1)+c2)
    end
  end
end

w_rho_r = -diff(times)
prepend!(w_rho_r,w_rho_r[1])
append!(w_rho_r,1.0)
  ## prep diffusivity input
  for ind in eachindex(times)

        """
        Damage calculation

        """
      if ind>1
         update_damage_rdaam!(times,eta_q,L_dist,U238_V,U235_V,Th232_V,e_rho_s,ind,rho_r,w_rho_r)
      end

      # from Meesters & Dunai 02 formulation
      F[ind] = 8*(U238_0)*(-exp(-times[ind]/τ38))
      F[ind] += 7*(U235_0)*(-exp(-times[ind]/τ35))
      F[ind] += 6*(Th232_0)*(-exp(-times[ind]/τ32))
      if ind>1
        damage_mod_term = ((psi*e_rho_s[ind] + omega*e_rho_s[ind]^3)*exp(Etrap/(R*T[ind-1])))+1
        D0[ind] = arrhenius(L,T[ind-1],R,Ea,logD0_a2)/damage_mod_term
        dzeta_val = times[ind-1]-times[ind]
        dzeta[ind] = (D0[ind-1]+D0[ind])*(times[ind-1]-times[ind])/2
        zeta[ind] = zeta[ind-1]+dzeta[ind]

      else
        damage_mod_term = ((psi*e_rho_s[ind] + omega*e_rho_s[ind]^3)*exp(Etrap/(R*T[1])))+1
        D0[ind] = arrhenius(L,T[1],R,Ea,logD0_a2)/damage_mod_term
      end

      if (ind<N_t_segs+1) & (ind>1)
         dfdchi[ind-1]=(F[ind-1]-F[ind])/(dzeta[ind])
      end

    end

   zeta_end = zeta[end]

    uterm = fill_u_term(L,n_iter,N_t_segs,F,dfdchi,dzeta,zeta,zeta_end)

    uF = uterm*(8)/pi
    uF = (uF/(pi*4/3))

    return uF

end
# end rdaam

## calc new damage
function update_damage_rdaam!(times,eta_q,L_dist,U238_V,U235_V,Th232_V,e_rho_s,ind,rho_r_final,w_rho_r)

  dt = times[ind-1]-times[ind]

  rho_v = (8/8)*U238_V*(exp(lambda_38*times[ind-1])-exp(lambda_38*times[ind]))
  rho_v += (7/8)*U235_V*(exp(lambda_35*times[ind-1])-exp(lambda_35*times[ind]))
  rho_v += (6/8)*Th232_V*(exp(lambda_32*times[ind-1])-exp(lambda_32*times[ind]))

  e_rho_s[ind]= rho_v*sum(rho_r_final[end-(ind-1):(end-1),ind].*(w_rho_r[2:ind]./w_rho_r[ind]))*eta_q*L_dist*lambda_f/lambda_38

#./w_rho_r[end-(ind-2):(end)]
  return nothing
end

function zrdaam_forward_diffusion(n_iter,c0,c1,c2,c3,alpha,Ba,interconnection,SV,lint_lattice,D0,Ea,R,D0N17,EaN17,U238,U235,Th232,
  L,tT...)
  equivalent_total_anneal_len=0.36/1.25+0.2

  tTcopy = collect(tT)
  times = tTcopy[1:ceil(Int,length(tT)/2)] # type unstable but fastest so far
  T = tTcopy[ceil(Int,length(tT)/2)+1:end]

  # preallocate
  (N_t_segs,e_rho_s,rcb2,rho_r,D,F,zeta,dzeta,dfdchi) = preallocate_diff_data(T) # 1st arg is only for rad dam

  n_times = length(times)
  dt_equiv = 0.0

  U238_0 = U238*(exp(lambda_38*times[1])) # U238 is measured at present
  U235_0 = U235*(exp(lambda_35*times[1])) # U238 is measured at present
  Th232_0 = Th232*(exp(lambda_32*times[1])) # U238 is measured at present

"""
Damage calculation

"""
for ind in eachindex(times)
  if ind>1
    t_dam = times[n_times-ind+1]-times[n_times-ind+2]+dt_equiv
    rcb2_base = (c0 + c1*(log(t_dam)-c2)/(log((1.0/T[n_times-ind+1]))-c3))
    rcb2_fill_val = ((rcb2_base^(1.0/alpha))+1)^(-1)
    rcb2[ind,2] = rcb2_fill_val
    for jj in 1:(ind-1)
      if !(isapprox(rcb2_fill_val, 0, atol=1e-9))
        rcb2_new = 1.25*(rcb2_fill_val-0.2)
      else
        rcb2_new = 0.0
      end
      if rcb2_new<equivalent_total_anneal_len
        rho_r[n_times-ind+jj,n_times-jj] = 0.0
      else
        rho_r[n_times-ind+jj,n_times-jj] = rcb2_new
      end
    end
      # use concept of equivalent time
    if ind<n_times
      dt_equiv = exp((((((1.0/rcb2[ind,2])-1.0)^alpha-c0)*(log(1.0/T[n_times-ind])-c3))/c1)+c2)
     end
  end
end

# weighting based on segment length
w_rho_r = -diff(times)
prepend!(w_rho_r,w_rho_r[1])
append!(w_rho_r,1.0)

for ind in eachindex(times)

    """
    Damage calculation

    """
  if ind>1
    update_damage_zrdaam!(U238*Na,U235*Na,Th232*Na,times,e_rho_s,rho_r,w_rho_r,ind)
  end

  # from Meesters & Dunai 02 formulation
  F[ind] = 8*(U238_0)*(-exp(-times[ind]/τ38))
  F[ind] += 7*(U235_0)*(-exp(-times[ind]/τ35))
  F[ind] += 6*(Th232_0)*(-exp(-times[ind]/τ32))
  scaling_factor = 10000

  # if ind==1
  #   e_rho_s[ind] = 1e9
  #   fa=1.0-exp(-Ba*e_rho_s[ind]) # using naming conventions from DAAM Guenthner, 2021
  #   DI=1.0-exp(-Ba*e_rho_s[ind]*interconnection)
  #   lint=(4.2/(fa*SV))-2.5
  #   tortuosity=(lint_lattice/lint)^2
  #   Dtort=(1/tortuosity)*D0*exp(-Ea/(R*T[1])) #cm2/s
  #   Dtort_a2=Dtort/(L*(1.0/scaling_factor)*(1.0-DI))^2 #1/s
  #   DN17=D0N17*exp(-EaN17/(R*T[1])) #cm2/s
  #   DN17a2=DN17/(L*(1.0/scaling_factor)*DI)^2 #1/s
  #   D[ind] = (DI/DN17a2+(1.0-DI)/Dtort_a2)^(-1) # cm2/s
  #   D[ind] = D[ind]*L^2*1e-8 # microns -> centimeters

  # end

  if ind>1
    fa=1.0-exp(-Ba*e_rho_s[ind]) # using naming conventions from DAAM Guenthner, 2021
    DI=1.0-exp(-Ba*e_rho_s[ind]*interconnection)
    lint=(4.2/(fa*SV))-2.5
    tortuosity=(lint_lattice/lint)^2
    Dtort=(1.0/tortuosity)*D0*exp(-Ea/(R*T[ind-1])) #cm2/s

    Dtort_a2=Dtort/(L*(1.0/scaling_factor)*(1.0-DI))^2 #1/s
    DN17=D0N17*exp(-EaN17/(R*T[ind-1])) #cm2/s
    DN17a2=DN17/(L*(1/scaling_factor)*DI)^2 #1/s
    D[ind]= (DI/DN17a2+(1.0-DI)/Dtort_a2)^(-1) # cm2/s
    D[ind] = D[ind]*L^2*1e-12 # microns -> centimeters
    if ind>2
      dzeta[ind] = (D[ind-1]+D[ind])*(times[ind-1]-times[ind])/2
      zeta[ind] = zeta[ind-1]+dzeta[ind]
    end

  end

  if (ind<N_t_segs+1) & (ind>1)
    dfdchi[ind-1]=(F[ind-1]-F[ind])/(dzeta[ind])
  end

end

zeta = zeta[2:end-1]
dzeta = dzeta[2:end-1]
dfdchi = dfdchi[2:end-1]

zeta_end = zeta[end]

uterm = fill_u_term(L,n_iter,N_t_segs-2,F,dfdchi,dzeta,zeta,zeta_end)

uF = uterm*(8)/pi
uF = (uF/(pi*4/3))

  return uF

end

function update_damage_zrdaam!(U238_V,U235_V,Th232_V,times,e_rho_s,rho_r,w_rho_r,ind)

  rho_v = (8)*U238_V*(exp(lambda_38*times[ind-1])-exp(lambda_38*times[ind]))
  rho_v += (7)*U235_V*(exp(lambda_35*times[ind-1])-exp(lambda_35*times[ind]))
  rho_v += (6)*Th232_V*(exp(lambda_32*times[ind-1])-exp(lambda_32*times[ind]))
  e_rho_s[ind]= rho_v*sum(rho_r[end-(ind-1):(end-1),ind].*(w_rho_r[2:ind]./w_rho_r[ind]))
end

function arrhenius(L,T,R,Ea,log10D0_a2)

  D = 10^(log10D0_a2)*L^2.0.*exp(-Ea./(R*T))

  return D
end

function fill_u_term(L,n_iter,N_t_segs,F,dfdchi,dzeta,zeta,zeta_end)
  uterm = 0.0
  mu_n = 0.0
  for n in 1:n_iter
          utermbase = 0.0
          mu_n = (pi/(L))^2*n^2
          for j in 1:(N_t_segs-1)

             utermbase+= dfdchi[j]*(exp(-mu_n*(zeta_end-zeta[j+1]))-exp(-mu_n*(zeta_end-zeta[j])))

          end
          uterm += (1/mu_n)*(1/(n*n))*utermbase

  end
  return uterm
end

function preallocate_diff_data(T)

  Tvv = eltype(T)
  N_t_segs = length(T)+1
  e_rho_s= zeros(Tvv,N_t_segs)

  D0 = zeros(Tvv,N_t_segs)
  F = zeros(Tvv,N_t_segs)
  zeta = zeros(Tvv,N_t_segs)
  dzeta = zeros(Tvv,N_t_segs)
  dfdchi = zeros(Tvv,(N_t_segs-1))
  rcb2 = zeros(Tvv,(N_t_segs,N_t_segs))
  rho_r = zeros(Tvv,(N_t_segs,N_t_segs))

  return N_t_segs,e_rho_s,rcb2,rho_r,D0,F,zeta,dzeta,dfdchi
end

function mod_constraints(T::Tvv...) where {Tvv<:Real}

n_smooth = length(T)-2
smooth_tot = 0.0
for i in 1:n_smooth
  smooth_tot += ((T[i]-2*T[i+1]+T[i+2])^2)*1e0
end

return smooth_tot

end


text:
mutable struct BarplotGraphics{R<:Number} <: GraphicsArea
    bars::Vector{R}
    color::Symbol
    width::Int
    max_freq::R
    max_len::Int
    symb::String

    function BarplotGraphics{R}(
            bars::AbstractVector{R},
            width::Int,
            color::Symbol,
            symb::String) where R
        width = max(width, 5)
        max_freq = maximum(bars)
        max_len = length(string(max_freq))
        new{R}(bars, color, width, max_freq, max_len, symb)
    end
end

nrows(c::BarplotGraphics) = length(c.bars)
ncols(c::BarplotGraphics) = c.width

function BarplotGraphics(
        bars::AbstractVector{R},
        width::Int;
        color::Symbol = :blue,
        symb::String = "▪") where {R <: Number}
    BarplotGraphics{R}(bars, width, color, symb)
end

function addrow!(c::BarplotGraphics{R}, bars::AbstractVector{R}) where {R <: Number}
    append!(c.bars, bars)
    c.max_freq = maximum(c.bars)
    c.max_len = length(string(c.max_freq))
    c
end

function addrow!(c::BarplotGraphics{R}, bar::Number) where {R <: Number}
    push!(c.bars, R(bar))
    c.max_freq = max(c.max_freq, R(bar))
    c.max_len = length(string(c.max_freq))
    c
end

function printrow(io::IO, c::BarplotGraphics, row::Int)
    0 < row <= nrows(c) || throw(ArgumentError("Argument row out of bounds: $row"))
    bar = c.bars[row]
    max_bar_width = max(c.width - 2 - c.max_len, 1)
    bar_len = c.max_freq > 0 ? round(Int, bar/c.max_freq * max_bar_width, RoundNearestTiesUp) : 0
    bar_str = c.max_freq > 0 ? repeat(c.symb, bar_len) : ""
    bar_lbl = string(bar)
    printstyled(io, bar_str; color = c.color)
    printstyled(io, " ", bar_lbl; color = :white)
    pan_len = max(max_bar_width + 1 + c.max_len - bar_len - length(bar_lbl), 0)
    pad = repeat(" ", round(Int, pan_len))
    print(io, pad)
end


text:
export sum

import Base: sum
import Base: show

"""
    sum(; from = nothing)

Creates a sum operator, which applies a sum accumulator function
over the source Observable, and returns the accumulated result when the source completes,
given an optional initial value.

The `sum` operator is similar to `reduce(T, T, +)` (see [`reduce`](@ref)).

# Arguments
- `from`: optional initial accumulation value, if nothing first value will be used instead

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ i for i in 1:42 ])
subscribe!(source |> sum(), logger())
;

# output

[LogActor] Data: 903
[LogActor] Completed

```

```jldoctest
using Rocket

source = from([ i for i in 1:42 ])
subscribe!(source |> sum(from = 97), logger())
;

# output

[LogActor] Data: 1000
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`reduce`](@ref), [`logger`](@ref)
"""
sum(; from::D = nothing) where D = SumOperator{D}(from)

struct SumOperator{D} <: InferableOperator
    from :: D
end

function on_call!(::Type{L}, ::Type{Union{L, Nothing}}, operator::SumOperator, source) where L
    return proxy(Union{L, Nothing}, source, SumProxy(operator.from !== nothing ? convert(L, operator.from) : nothing))
end

operator_right(operator::SumOperator, ::Type{L}) where L = Union{L, Nothing}

struct SumProxy{D} <: ActorProxy
    from :: D
end

actor_proxy!(::Type{Union{L, Nothing}}, proxy::SumProxy, actor::A) where { L, A } = SumActor{L, A}(actor, proxy.from)

mutable struct SumActor{L, A} <: Actor{L}
    actor   :: A
    current :: Union{L, Nothing}
end

getcurrent(actor::SumActor)         = actor.current
setcurrent!(actor::SumActor, value) = actor.current = value

function on_next!(actor::SumActor{L}, data::L) where L
    current = getcurrent(actor)
    if current === nothing
        setcurrent!(actor, data)
    else
        setcurrent!(actor, current + data)
    end
end

function on_error!(actor::SumActor, err)
    error!(actor.actor, err)
end

function on_complete!(actor::SumActor)
    next!(actor.actor, getcurrent(actor))
    complete!(actor.actor)
end

Base.show(io::IO, ::SumOperator)         = print(io, "SumOperator()")
Base.show(io::IO, ::SumProxy)            = print(io, "SumProxy()")
Base.show(io::IO, ::SumActor{L}) where L = print(io, "SumActor($L)")


text:
struct Morris <: GSAMethod
    p_steps::Array{Int,1}
    relative_scale::Bool
    num_trajectory::Int
    total_num_trajectory::Int
    len_design_mat::Int
end

Morris(;p_steps::Array{Int,1}=Int[], relative_scale::Bool=false, num_trajectory::Int=10,
total_num_trajectory::Int=5*num_trajectory, len_design_mat::Int=10) = Morris(p_steps, relative_scale, num_trajectory, total_num_trajectory, len_design_mat)

struct MatSpread{T1,T2}
    mat::T1
    spread::T2
end

struct MorrisResult{T1,T2}
    means::T1
    means_star::T1
    variances::T1
    elementary_effects::T2
end

function generate_design_matrix(p_range, p_steps, rng;len_design_mat = 10)
    ps = [range(p_range[i][1], stop=p_range[i][2], length=p_steps[i]) for i in 1:length(p_range)]
    indices = [rand(rng, 1:i) for i in p_steps]
    all_idxs = Vector{typeof(indices)}(undef,len_design_mat)

    for i in 1:len_design_mat
        j = rand(rng, 1:length(p_range))
        indices[j] += (rand(rng) < 0.5 ? -1 : 1)
        if indices[j] > p_steps[j]
            indices[j] -= 2
        elseif indices[j] < 1.0
            indices[j] += 2
        end
        all_idxs[i] = copy(indices)
    end

    B = Array{Array{Float64}}(undef,len_design_mat)
    for j in 1:len_design_mat
        cur_p = [ps[u][(all_idxs[j][u])] for u in 1:length(p_range)]
        B[j] = cur_p
    end
    reduce(hcat, B)
end

function calculate_spread(matrix)
    spread = 0.0
    for i in 2:size(matrix,2)
        spread += sqrt(sum(abs2.(matrix[:,i] - matrix[:,i-1])))
    end
    spread
end

function sample_matrices(p_range,p_steps, rng;num_trajectory=10,total_num_trajectory=5*num_trajectory,len_design_mat=10)
    matrix_array = []
    if total_num_trajectory<num_trajectory
        error("total_num_trajectory should be greater than num_trajectory preferably atleast 3-4 times higher")
    end
    for i in 1:total_num_trajectory
        mat = generate_design_matrix(p_range, p_steps, rng;len_design_mat = len_design_mat)
        spread = calculate_spread(mat)
        push!(matrix_array,MatSpread(mat,spread))
    end
    sort!(matrix_array,by = x -> x.spread,rev=true)
    matrices = [i.mat for i in matrix_array[1:num_trajectory]]
    reduce(hcat,matrices)
end

function gsa(f, method::Morris, p_range::AbstractVector; batch=false, rng::AbstractRNG = Random.default_rng(), kwargs...)
    @unpack p_steps, relative_scale, num_trajectory, total_num_trajectory, len_design_mat  = method
    if !(length(p_steps) == length(p_range))
        for i in 1:length(p_range)-length(p_steps)
            push!(p_steps,100)
        end
    end

    design_matrices = sample_matrices(p_range, p_steps, rng;num_trajectory=num_trajectory,
                                        total_num_trajectory=total_num_trajectory,len_design_mat=len_design_mat)

    multioutput = false
    desol = false
    local y_size
    if batch
        all_y = f(design_matrices)
        multioutput = all_y isa AbstractMatrix
    else
        _y = [f(design_matrices[:,i]) for i in 1:size(design_matrices,2)]
        multioutput = !(eltype(_y) <: Number)
        if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
            y_size = size(_y[1])
            _y = vec.(_y)
            desol = true
        end
        all_y = multioutput ? reduce(hcat,_y) : _y
    end

    effects = []
    for i in 1:num_trajectory
        y1 = multioutput ? all_y[:,(i-1)*len_design_mat+1] : all_y[(i-1)*len_design_mat+1]
        for j in (i-1)*len_design_mat+1:(i*len_design_mat)-1
            y2 = y1
            del = design_matrices[:,j+1] - design_matrices[:,j]
            change_index = 0
            for k in 1:length(del)
                if abs(del[k]) > 0
                    change_index = k
                    break
                end
            end
            del = sum(del)
            y1 = multioutput ? all_y[:,j+1] : all_y[j+1]
            if relative_scale == false
                effect = @. (y1-y2)/(del)
                elem_effect = typeof(y1) <: Number ? effect : mean(effect, dims = 2)
            else
                if del > 0
                    effect = @. (y1-y2)/(y2*del)
                    elem_effect = typeof(y1) <: Number ? effect : mean(effect, dims = 2)
                else
                    effect = @. (y1-y2)/(y1*del)
                    elem_effect = typeof(y1) <: Number ? effect : mean(effect, dims = 2)
                end
            end
            if length(effects) >= change_index && change_index > 0
                push!(effects[change_index],elem_effect)
            elseif change_index > 0
                while(length(effects) < change_index-1)
                    push!(effects,typeof(elem_effect)[])
                end
                push!(effects,[elem_effect])
            end
        end
    end
    means = eltype(effects[1])[]
    means_star = eltype(effects[1])[]
    variances = eltype(effects[1])[]
    for k in effects
        if !isempty(k)
            push!(means, mean(k))
            push!(means_star, mean(x -> abs.(x), k))
            push!(variances, var(k))
        else
            push!(means, zero(effects[1][1]))
            push!(means_star, zero(effects[1][1]))
            push!(variances, zero(effects[1][1]))
        end
    end
    if desol
        f_shape = x -> [reshape(x[:,i],y_size) for i in 1:size(x,2)]
        means = map(f_shape,means)
        means_star = map(f_shape,means_star)
        variances = map(f_shape,variances)
    end
    MorrisResult(reduce(hcat, means),reduce(hcat, means_star),reduce(hcat, variances),effects)
end


text:
using StatisticalRethinking
#gr(size=(600,600));

N = 20
p_grid = range( 0 , stop=1 , length=N )

prior = ones( 20 );

likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
likelihood[1:5]

unstd_posterior = likelihood .* prior;

posterior = unstd_posterior  ./ sum(unstd_posterior);

p1 = plot( p_grid , posterior ,
    xlab="probability of water" , ylab="posterior probability",
    lab = "interpolated", title="20 points" )
p2 = scatter!( p1, p_grid , posterior, lab="computed" )

prior1 = [p < 0.5 ? 0 : 1 for p in p_grid]
prior2 = [exp( -5*abs( p - 0.5 ) ) for p in p_grid]

p3 = plot(p_grid, prior1,
  xlab="probability of water" , ylab="posterior probability",
  lab = "semi_uniform", title="Other priors" )
scatter!(p3, p_grid, prior1, lab = "semi_uniform grid point")
plot!(p3, p_grid, prior2,  lab = "double_exponential" )
scatter!(p3, p_grid, prior2,  lab = "double_exponential grid point" )

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module TestDicts

using BangBang: merge!!
using MicroCollections
using Test

@testset "empty" begin
    @test EmptyDict() === emptyshim(Dict)
    @test eltype(EmptyDict()) === Pair{Union{},Union{}}
    @test EmptyDict(Dict{Symbol,Int}()) === EmptyDict{Symbol,Int}()
    @test EmptyDict(Union{}[]) === EmptyDict{Union{},Union{}}()
    @test EmptyDict{Char,Nothing}(Dict{Symbol,Int}()) === EmptyDict{Char,Nothing}()
    @test EmptyDict{Char,Nothing}(Union{}[]) === EmptyDict{Char,Nothing}()
    @test_throws ArgumentError EmptyDict(Dict(:a => 1))
    @test_throws ArgumentError EmptyDict{Char,Nothing}(Dict(:a => 1))
end

@testset "singleton" begin
    @test SingletonDict(:a => 0) === singletonshim(Dict, :a => 0)
    @test SingletonDict(:a => 0)::AbstractDict{Symbol,Int} == Dict(:a => 0)
    @test SingletonDict(Dict(:a => 0)) === singletonshim(Dict, :a => 0)
    @test SingletonDict{Symbol,Float64}(Dict(:a => 0)) === singletonshim(Dict, :a => 0.0)
    @test_throws ArgumentError SingletonDict(Dict())
    @test_throws ArgumentError SingletonDict(Dict(:a => 0, :b => 1))
end

@testset "$_merge" for _merge in [merge, merge!!]
    @test _merge(SingletonDict(:a => 0), Dict(:b => 0.5))::Dict{Symbol,Float64} ==
          Dict(:a => 0.0, :b => 0.5)
    @test _merge(SingletonDict(:a => 0), SingletonDict(:b => 0.5))::Dict{Symbol,Float64} ==
          Dict(:a => 0.0, :b => 0.5)
    @test _merge(SingletonDict(:a => 0), EmptyDict())::Dict{Symbol,Int} == Dict(:a => 0)
    @test _merge(EmptyDict(), SingletonDict(:a => 0))::Dict{Symbol,Int} == Dict(:a => 0)
    @test _merge(SingletonDict(:a => 0), EmptyDict())::Dict{Symbol,Int} == Dict(:a => 0)

    @test _merge(+, SingletonDict(:a => 1), Dict(:a => 0.5))::Dict{Symbol,Float64} ==
          Dict(:a => 1.5)
    @test _merge(
        +,
        SingletonDict(:a => 1),
        SingletonDict(:a => 0.5),
    )::Dict{Symbol,Float64} == Dict(:a => 1.5)
    @test _merge(+, EmptyDict(), SingletonDict(:a => 1))::Dict{Symbol,Int} == Dict(:a => 1)
    @test _merge(+, SingletonDict(:a => 1), EmptyDict())::Dict{Symbol,Int} == Dict(:a => 1)
end

end  # module


text:
module DistributedOperations

using Distributed

struct TypeFutures{T} f::Dict{Int,Future} end
Base.getindex(a::TypeFutures, i::Int) = getindex(a.f, i)
Base.setindex!(a::TypeFutures, f::Future, i::Int) = setindex!(a.f, f, i)
Base.keys(a::TypeFutures) = keys(a.f)

"""
    x = TypeFutures(T, f, pids, fargs...)

Construt a `x::TypeFutures` of type `T` on workers defined by the process id's `pids`.  On each
worker `pid`, `f` is evaluated, and a future for what is returned by `f` is stored.

# Example
```
using Distributed
addprocs(2)
@everywhere using DistributedOperations
@everywhere struct MyStruct
    x::Vector{Float64}
    y::Vector{Float64}
end
@everywhere foo() = MyStruct(rand(10), rand(10))
x = TypeFutures(MyStruct, foo, procs())
@show remotecall_fetch(localpart, workers()[1], x)
rmprocs(workers())
```
"""
function TypeFutures(_T::Type{T}, f::Function, pids::AbstractArray, fargs::Vararg) where {T}
    pids[1] == myid() || error("expected myid()==pids[1], got pids[1]=$(pids[1]) where-as myid()=$(myid())")

    futures = Dict()
    @sync for pid in pids
        futures[pid] = remotecall(f, pid, fargs...)
        @async remotecall_fetch(wait, pid, futures[pid])
    end
    TypeFutures{T}(futures)
end
TypeFutures(_T::Type{T}, f::Function, fargs::Vararg) where {T} = TypeFutures(_T, f, procs(), fargs...)

"""
    x = TypeFutures(y::T, f[, pids=procs()], fargs...)

Construt a `x::TypeFutures` from `y::T` on workers defined by the process id's `pids`.  On
each worker `pid`, `f` is evaluated, and a future for what is returned by `f` is stored.

# Example
```
using Distributed
addprocs(2)
@everywhere using DistributedOperations
@everywhere struct MyStruct
    x::Vector{Float64}
    y::Vector{Float64}
end
@everywhere foo() = MyStruct(rand(10), rand(10))
x = foo()
x = TypeFutures(x, foo, procs())
@show remotecall_fetch(localpart, workers()[1], x)
rmprocs(workers())
```
"""
function TypeFutures(x::T, f::Function, pids::AbstractArray, fargs::Vararg) where {T}
    futures = Dict()
    @sync for pid in pids
        if pid == myid()
            futures[pid] = Future()
            put!(futures[pid], x)
        else
            futures[pid] = remotecall(f, pid, fargs...)
            @async remotecall_fetch(wait, pid, futures[pid])
        end
    end
    TypeFutures{T}(futures)
end
TypeFutures(x::T, f::Function, fargs::Vararg) where {T} = TypeFutures(x, f, procs(), fargs...)

ArrayFutures{T,N} = TypeFutures{Array{T,N}}

"""
    x = ArrayFutures(T, n::NTuple{N,Int}[, pids=procs()])

Create `x::TypeFutures`, and where each proccess id (pid) in `pids` is
assigned `zeros(T,n)`.

# Example
```
using Distributed
addprocs(2)
@everywhere using DistributedOperations
x = ArrayFutures(Float32, (10,20), procs())
localpart(x)
rmprocs(workers())
```
"""
ArrayFutures(_T::Type{T}, n::NTuple{N,I}, pids=procs()) where {T,N,I<:Integer} = TypeFutures(Array{T,N}, zeros, pids, T, n)

"""
    x = ArrayFutures(x::Array[, pids=procs()])

Create `x::TypeFutures`, and where `myid()` is assigned `x`, and all other processes
are assigned `zeros(eltype(x), size(x))`.
"""
ArrayFutures(x::Array{T,N}, pids=procs()) where {T,N} = TypeFutures(x, zeros, pids, T, size(x))

"""
    bcast!(x::TypeFutures, pids)

Broadcast an existing `x::TypeFutures` to `pids`.  This is useful for
elastic computing where the cluster may grow after the construction
and broadcast of `x::TypeFuture`.
"""
function bcast!(x::TypeFutures{T}, pids) where {T}
    _pids = 1 ∈ pids ? pids : [1;pids]
    _x = bcast(localpart(x)::T, _pids)
    merge!(x.f, _x.f)
    x
end

"""
    bcast(x[, pids=procs()])

Broadcast `x` to `pids`.

# Example
```
using Distributed
addprocs(2)
@everywhere using DistributedOperations
x = rand(10)
_x = bcast(x)
y = remotecall_fetch(localpart, workers()[1], _x)
y ≈ x  # true
rmprocs(workers())
```
"""
function bcast(x::T, pids=procs()) where {T}
    pids[1] == myid() || error("expected myid()==pids[1], got pids[1]=$(pids[1]) where-as myid()=$(myid())")

    M = length(pids)
    L = round(Int,log2(prevpow(2,M)))

    _f(x) = x
    futures = Dict(pids[1]=>remotecall(_f, myid(), x))

    for l = 1:L
        m = 2^(l-1)
        @sync for i = 1:m
            futures[pids[i+m]] = remotecall(fetch, pids[i+m], futures[pids[i]])
            @async remotecall_fetch(wait, pids[i+m], futures[pids[i+m]])
        end
    end

    m = 2^L
    R = M - m

    if R != 0
        @sync for i = 1:R
            futures[pids[i+m]] = remotecall(fetch, pids[i+m], futures[pids[i]])
            @async remotecall_fetch(wait, pids[i+m], futures[pids[i+m]])
        end
    end

    TypeFutures{T}(futures)
end

@inline paralleloperations_reduce!(y, x) = begin y .+= x; nothing end

"""
    y = reduce!(x::TypeFutures[, reducemethod!=DistributedOperations.paralleloperations_reduce!])

Parallel reduction of `x::TypeFutures` using `reducemethod!`.  By default, the reduction
is a mutating in-place element-wise addition, such that `y=localpart(x)`.

# Example
```
using Distributed
addprocs(2)
@everywhere using DistributedOperations
x = ArrayFutures(Float64, (3,))
fill!(x, 1, workers())
y = reduce!(x)
y ≈ [2.0,2.0,2.0] # true
localpart(x) ≈ [2.0,2.0,2.0] # true
rmprocs(workers())
```
"""
function reduce!(futures::TypeFutures{T}, reducemethod!::Function=paralleloperations_reduce!) where {T}
    function _reduce!(future_mine, future_theirs, reducemethod!, _T::Type{T}) where {T}
        x = remotecall_fetch(fetch, future_theirs.where, future_theirs)::T
        y = fetch(future_mine)::T
        reducemethod!(y, x)
        nothing
    end

    pids = sort(collect(keys(futures)))
    M = length(pids)
    L = round(Int,log2(prevpow(2,M)))
    m = 2^L
    R = M - m

    if R != 0
        @sync for i = 1:R
            @async remotecall_fetch(_reduce!, pids[i], futures[pids[i]], futures[pids[m+i]], reducemethod!, T)
        end
    end

    for l = L:-1:1
        m = 2^(l-1)
        @sync for i = 1:m
            @async remotecall_fetch(_reduce!, pids[i], futures[pids[i]], futures[pids[m+i]], reducemethod!, T)
        end
    end
    fetch(futures[myid()])::T
end

@inline paralleloperations_copy!(y, x) = begin y .= x; nothing end

"""
    copy!(to::TypeFutures, from::TypeFutures[, copymethod!=DistributedOperations.paralleloperations_copy!, pids])

Copy `from` into `to` using `copymethod!`.
"""
function Base.copy!(to::TypeFutures{T}, from::TypeFutures{T}, copymethod!::Function, pids::AbstractArray=Int[]) where {T}
    pids = isempty(pids) ? keys(to) : pids
    function _copy!(future_to, future_from, copymethod!, _T::Type{T}) where {T}
        y = fetch(future_to)::T
        x = fetch(future_from)::T
        copymethod!(y, x)
        nothing
    end
    @sync for pid in pids
        @async pid ∈ keys(to) && pid ∈ keys(from) && remotecall_fetch(_copy!, pid, to[pid], from[pid], copymethod!, T)
    end
end
Base.copy!(to::TypeFutures, from::TypeFutures, pids::AbstractArray=procs()) = copy!(to, from, paralleloperations_copy!, pids)

@inline paralleloperations_fill!(x, a) = begin x .= a; nothing end

"""
    fill!(x::TypeFutures, a[, fillmethod!=DistributedOperations.fillmethod!, pids])

Fill `x` with `a::Number` using the `fillmethod!::Function`.
"""
function Base.fill!(futures::TypeFutures{T}, a::Number, fillmethod!::Function, pids::AbstractArray=Int[]) where {T}
    pids = isempty(pids) ? keys(futures) : pids
    function _fill!(future, a, fillmethod!, _T::Type{T}) where {T}
        x = fetch(future)::T
        fillmethod!(x,a)
        nothing
    end
    @sync for pid in pids
        @async remotecall_fetch(_fill!, pid, futures[pid], a, fillmethod!, T)
    end
end
Base.fill!(futures::TypeFutures, a::Number, pids::AbstractArray=Int[]) = fill!(futures, a, paralleloperations_fill!, pids)

using DistributedArrays
import DistributedArrays.localpart

"""
    localpart(x::TypeFutures)

Get the piece of `x::TypeFutures` that is local to `myid()`.
"""
localpart(futures::TypeFutures{T}) where {T} = fetch(futures[myid()])::T

Base.show(io::IO, futures::TypeFutures) = write(io, "TypeFutures with pids=$(keys(futures)) and type $(typeof(localpart(futures)))")
Base.show(io::IO, futures::ArrayFutures) = write(io, "ArrayFutures with pids=$(keys(futures)) and type $(size(localpart(futures)))")

export ArrayFutures, TypeFutures, bcast, bcast!, localpart, reduce!

end


text:
################################################################################
# same with LLL
################################################################################

function single_env{T}(c::ClassGrpCtx{T}, I::Hecke.SmallLLLRelationsCtx, nb::Int, rat::Float64, max_good::Int = 2)
  bad_norm = 0
  rk = rank(c.M)
  good = 0
  while true
    e = class_group_small_lll_elements_relation_next(I)
    n = norm_div(e, norm(I.A), nb)
    if I.cnt > length(I.b)^2
      break
    end
    if nbits(numerator(n)) > nb - 25
      bad_norm += 1
      if I.cnt > 100  && bad_norm / I.cnt > 0.1
        @vprint :ClassGroup 2 "norm too large, $(I.cnt) has size $(nbits(numerator(n))) should be <= $(nb - 20) $bad_norm $(I.cnt)\n"
        break
      end
      continue
    end
    fl = class_group_add_relation(c, e, n, norm(I.A), integral = true)
    if !fl  && I.cnt/(good+1) > 2*c.expect
      @vprint :ClassGroup 2 "not enough progress $(I.cnt) $(c.expect) $good\n"
      break
    end
    if fl
      good += 1
    end
    if fl && max_good > -1
      if max_good < good
        @vprint :ClassGroup 2 "found enough $(I.cnt)\n"
        break
      end
    end
    if false && fl && good > 0 && (rank(c.M)- rk+1)/good < rat
      @vprint :ClassGroup 2 "rank too slow $(I.cnt) $good $(rank(c.M)) $rk\n"
      break
    end
  end
end

function class_group_via_lll{T}(c::ClassGrpCtx{T}, rat::Float64 = 0.2)
  O = order(c.FB.ideals[1])
  nb = nbits(abs(discriminant(O)))
  nb = div(nb, 2) + 30

  rt = time_ns()
  I = class_group_small_lll_elements_relation_start(c, O)
  single_env(c, I, nb, rat/10, -1)
  @vprint :ClassGroup 1 "search in order:  $((time_ns()-rt)*1e-9) rel mat:  $(c.M.bas_gens)\nin $(c.M)"

  @vtime :ClassGroup 1 h, piv = class_group_get_pivot_info(c)
  if h == 1 return c; end

#  for p = piv
#    I = class_group_small_lll_elements_relation_start(c, c.FB.ideals[p])
#    single_env(c, I, nb, rat, 1)
#  end
#
#  @vprint :ClassGroup 1 "search in ideals:  $((time_ns()-rt)*1e-9) rel mat:  $(c.M.bas_gens)\n"

#  @vtime :ClassGroup 1 h, piv = class_group_get_pivot_info(c)
#  if h == 1 return c; end

#  @vprint :ClassGroup 1 "Now with random...\n"
#  @vprint :ClassGroup 1 "length(piv) = $(length(piv)) and h = $h\n"
#  @vprint :ClassGroup 1 "$(piv)\n"

  class_group_new_relations_via_lll(c, rat, extra = -1)

  return c
end

function class_group_new_relations_via_lll{T}(c::ClassGrpCtx{T}, rat::Float64 = 0.2; extra::Int = 5, rand_exp::Int = 1)

  O = order(c.FB.ideals[1])
  nb = nbits(abs(discriminant(O)))
  nb = div(nb, 2) + 30

  st = c.rel_cnt

  @vtime :ClassGroup 1 h, piv = class_group_get_pivot_info(c)

  @vprint :ClassGroup 1 "Now with random...\n"
  @vprint :ClassGroup 1 "length(piv) = $(length(piv)) and h = $h\n"
  @vprint :ClassGroup 1 "$(piv)\n"
  if length(piv) == 0
    for i=1:5
      push!(piv, rand(1:length(c.FB.ideals)))
    end
    @vprint :ClassGroup 1 "piv was empty, supplemented it to\n"
    @vprint :ClassGroup 1 "$(piv)\n"
  end

  start = max(1, max(div(length(c.FB.ideals), 2)+1, length(c.FB.ideals)-10*(1+div(rand_exp, 3))))
  stop = length(c.FB.ideals)
  if isdefined(c, :randomClsEnv)
    rand_env = c.randomClsEnv
    random_extend(rand_env, 2.0)
    @vprint :ClassGroup 1 "re-using and extending random: $(nbits(norm(rand_env.rand))) and $(rand_env.exp)\n"
  else
    @vprint :ClassGroup 1 "want $(stop-start) primes for random. Try distinct rational primes...\n"
    JJ = [c.FB.ideals[stop]]
    start += length(c.FB.ideals) - stop
    stop  += length(c.FB.ideals) - stop
    i = stop
    while i>1 && length(JJ) < stop - start
      i -= 1
      if degree(c.FB.ideals[i]) > 1
        continue
      end
      if minimum(c.FB.ideals[i]) in [minimum(x) for x = JJ]
        continue
      else
        push!(JJ, c.FB.ideals[i])
      end
    end
    while length(JJ) < stop - start
      AA = rand(c.FB.ideals)
      if !(AA in JJ)
        push!(JJ, AA)
      end
    end
    rand_env = random_init(JJ, lb = root(abs(discriminant(O)), 2)^1, ub = abs(discriminant(O))^1, reduce = false)
    c.randomClsEnv = rand_env
  end

  if h > 0
    rand_exp += 1
    while gcd(h, rand_exp) > 1
      rand_exp += 1
    end
  end

  while true
    st = c.rel_cnt
    while (c.rel_cnt - st < 2)
      for p = piv
        @vprint :ClassGroup 1 "p: $p $rand_exp $(length(rand_env.base))\n"
        @vtime :ClassGroup 3 J = random_get(rand_env, reduce = false)
  #      @show nbits(norm(J)), rand_env.exp, rand_exp
        @vtime :ClassGroup 3 J *= c.FB.ideals[p]^rand_exp
        @vtime :ClassGroup 3 I = class_group_small_lll_elements_relation_start(c, J)
        @vtime :ClassGroup 3 single_env(c, I, nb, 0.8, -1)
        if h > 0 && c.rel_cnt - st > 2
          return
        else
          if h == 0 && c.rel_cnt - st > length(piv) + 2
            break
          end
        end
      end
    end

    @vprint :ClassGroup 1 "eval info\n"
    @vtime :ClassGroup 1 h, piv_new = class_group_get_pivot_info(c)
    @vprint :ClassGroup 1 "length(piv) = $(length(piv_new)) and h = $h\n"
    @vprint :ClassGroup 1 "$(piv_new)\n"

    if piv_new == piv
      if h > 0
        extra = 5
#        J = [rand(c.FB.ideals) for x=1:10]
#        println("extending rand")
#        random_extend(rand_env, J)
        random_extend(rand_env, root(abs(discriminant(O)), 2))
      end
      rand_exp += 1
      rand_exp = min(rand_exp, 13)
      if h>0
        while gcd(rand_exp, h) > 1
          rand_exp += 1
        end
      end

      if rand_exp % 3 == 0
        start = max(start -10, 1)
      end
    end
    piv = piv_new
    if h == 1 return end
  end
end


text:
##  Test continuous wavelet transform and wavelet power spectrum
##  with the example of the ElNino time series
##  again closely following Torrence & Compo
##
##  Author Bernd Blasius (bernd.blasius@gmail.com)

module W

using CSV, PyPlot

include("../src/cwt.jl")

function test_elnino()
  # read in the data
  df = CSV.read("../data/ElNino_data.csv"; comment="#")
  t   = df[:,1]  # time (year)
  sst = df[:,2]  # el nino signal (sea surface temperature SST)
  x   = df[:,3]  # el nino temperature anomaly

  figure("elnino time series")
  clf()
  subplot(2,1,1)
  plot(t,sst)
  xlabel("time")
  ylabel("SST")
  subplot(2,1,2)
  plot(t,x)
  xlabel("time")
  ylabel("anomaly")

  dt = mean(diff(t))   # actually we would need perfectly regular sampling intervals
  nvoices = 16  # total number of voices
  noctave = 5   # total number of octaves
  s0 = 1.0      # starting scale at 1 years

  wave, period, power, global_ws, variance, scale, coi =
    cwt(x,dt,noctave = noctave, nvoices= nvoices, s0=s0)

  # Significance levels: (variance=1 for the normalized SST)
  lag1 = 0.72    # lag-1 autocorrelation for red noise background
  signif,fft_theor = wave_signif(variance,dt,scale,0,lag1)
  sig95 = signif .* ones(length(x))'  # expand signif --> (J+1)x(N) array
  sig95 = power ./ sig95              # where ratio > 1, power is significant

  # significance level for global spectrum
  dof = length(x) .- scale  # the -scale corrects for padding at edges
  global_signif,fft_theor = wave_signif(variance,dt,scale,1,lag1,-1,dof)

  # plotting scale range
  lminp = log2(minimum(period))
  lmaxp = log2(maximum(period))
  lrange = (trunc(lminp) :  trunc(lmaxp))

  figure("El Nino wavelet power spectrum")
  clf()
  # either use pcolormesh ..
    pcolormesh(t, log2.(period),power, cmap=PyPlot.cm.jet, alpha = 0.5)
  # ... or use contours
  levels = range(0, stop=maximum(power), length=200)
  contourf(t, log2.(period),power,levels, cmap=PyPlot.cm.jet )
  #
  plot(t,log2.(coi),"k")           # cone-of-influence, anything "below" is dubious
  contour(t, log2.(period),sig95, [-99,1], colors="black")   # 95% significance value
  ylabel("period length")
  xlabel("time")
  ylim(lmaxp,lminp)
  yticks(lrange, 2 .^ lrange)

  figure("global spectrum")
  clf()
  plot(log2.(period),log2.(global_ws))
  plot(log2.(period),log2.(global_signif),"k.")  # 95% significance line
  xlabel("period")
  ylabel("log2(power)")
  xlim(lmaxp,lminp)
  xticks(lrange, 2 .^ lrange)
end

test_elnino()

end # module

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
abstract type Integrator{dType, tType} end

abstract type DeterministicIntegrator{dType, tType} <: Integrator{dType, tType} end

abstract type ODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
abstract type DAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
abstract type PODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
abstract type PDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end

abstract type IODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end
abstract type IDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end
abstract type HODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end
abstract type HDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end
abstract type LODEIntegrator{dType, tType} <: IODEIntegrator{dType, tType} end
abstract type LDAEIntegrator{dType, tType} <: IDAEIntegrator{dType, tType} end

parameters(integrator::Integrator) = error("parameters() not implemented for ", typeof(integrator))
equations(integrator::Integrator) = error("equations() not implemented for ", typeof(integrator))
equation(integrator::Integrator) = error("equation() not implemented for ", typeof(integrator))
equation(integrator::Integrator, i::Int) = error("equation() not implemented for ", typeof(integrator))
timestep(integrator::Integrator) = error("timestep() not implemented for ", typeof(integrator))
Base.ndims(integrator::Integrator) = error("ndims() not implemented for ", typeof(integrator))
GeometricBase.nconstraints(integrator::Integrator) = error("nconstraints() not implemented for ", typeof(integrator))
nstages(integrator::Integrator) = error("nstages() not implemented for ", typeof(integrator))

eachdim(integrator::Integrator) = 1:ndims(integrator)

"""
```julia
get_internal_variables(::Integrator) = NamedTuple()
```
Returns a `NamedTuple` containing all internal variables of an integrator that
shall be stored in an [`AtomicSolution`](@ref). If there is no method for a
specific integrator implemented an empty `NamedTuple()` is returned.
"""
get_internal_variables(::Integrator) = NamedTuple()
get_internal_variables(::Nothing) = NamedTuple()

# Create AtomicSolution for ODE.
function Solutions.AtomicSolution(solution::SolutionODE, integrator::Integrator)
    AtomicSolutionODE(get_initial_conditions(solution, 1)..., get_internal_variables(integrator))
end

# Create AtomicSolution for partitioned ODE.
function Solutions.AtomicSolution(solution::SolutionPODE, integrator::Integrator)
    AtomicSolutionPODE(get_initial_conditions(solution, 1)..., get_internal_variables(integrator))
end

# Create AtomicSolution for DAE.
function Solutions.AtomicSolution(solution::SolutionDAE, integrator::Integrator)
    AtomicSolutionDAE(get_initial_conditions(solution, 1)..., get_internal_variables(integrator))
end

# Create AtomicSolution for partitioned DAE.
function Solutions.AtomicSolution(solution::SolutionPDAE, integrator::Integrator)
    AtomicSolutionPDAE(get_initial_conditions(solution, 1)..., get_internal_variables(integrator))
end

abstract type Parameters{DT,TT} end

function_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT)
solution_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT)

initialize!(::Integrator, ::AtomicSolution) = nothing

"""
Performs one time step with a given integrator.

```julia
integrate_step!(integrator::Integrator, asol::AtomicSolution)
```

The function accepts two arguments: an integrator and an appropriate [`AtomicSolution`](@ref),
which contains the state of the system at the beginning and the end of the time step and possibly
additional information like solver output or the solution at internal stages of a Runge-Kutta
method.
"""
integrate_step!(integrator::Integrator, ::AtomicSolution) = error("integrate_step()! not implemented for ", typeof(integrator))


text:
## keepnames.jl  methods for NamedArray that keep the names (some checking may be done)

## (c) 2013, 2014, 2017, 2018 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

# Keep names for consistently named vectors, or drop them
function Base.hcat(N::NamedVecOrMat...)
    keepnames=true
    N1 = N[1]
    firstnames = names(N1, 1)
    for i in 2:length(N)
        keepnames &= names(N[i], 1) == firstnames
    end
    a = hcat(map(n -> n.array, N)...)
    if keepnames
        colnames = defaultnamesdict(size(a,2))
        NamedArray(a, (N1.dicts[1], colnames), (N1.dimnames[1], :hcat))
    else
        NamedArray(a)
    end
end

function Base.vcat(N::NamedMatrix...)
    keepnames=true
    N1=N[1]
    firstnames = names(N1, 2)
    for i=2:length(N)
        keepnames &= names(N[i], 2)==firstnames
    end
    a = vcat(map(a -> a.array, N)...)
    if keepnames
        rownames = defaultnamesdict(size(a,1))
        NamedArray(a, (rownames, N1.dicts[2]), (:vcat, N1.dimnames[2]))
    else
        NamedArray(a)
    end
end

function Base.vcat(N::NamedVector...)
    a = vcat(map(n -> n.array, N)...)
    anames = vcat(map(n -> names(n, 1), N)...)
    if length(unique(anames)) == length(a)
        return NamedArray(a, (anames,), (:vcat,))
    else
        return NamedArray(a, (defaultnamesdict(length(a)),), (:vcat,))
    end
end

## broadcast v1.0
Base.BroadcastStyle(::Type{A}) where {A <: NamedArray} = Broadcast.ArrayStyle{NamedArray}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{A}},
                      ::Type{T}) where {A <: NamedArray, T}
    namedarray = find_namedarray(bc)
    return similar(namedarray, T)
end
"`find_namedarray(As)` returns the first NamedArray among the arguments."
find_namedarray(bc::Base.Broadcast.Broadcasted) = find_namedarray(bc.args)
function find_namedarray(args::Tuple)
    find_namedarray(find_namedarray(args[1]), Base.tail(args))
end
find_namedarray(x::NamedArray) = x
find_namedarray(a::NamedArray, rest) = a
find_namedarray(::Any, rest) = find_namedarray(rest)
find_namedarray(x::Any) = ()

## reorder names
import Base: sort, sort!
function sort!(v::NamedVector; kws...)
    i = sortperm(v.array; kws...)
    newnames = names(v, 1)[i]
    empty!(v.dicts[1])
    for (ind, k) in enumerate(newnames)
        v.dicts[1][k] = ind
    end
    v.array = v.array[i]
    return v
end

sort(v::NamedVector; kws...) = sort!(copy(v); kws...)

## Note: I can't think of a sensible way to define sort!(a::NamedArray, dim>1)

## drop name of sorted dimension, as each index along that dimension is sorted individually
function sort(n::NamedArray; dims::Integer, kws...)
    nms = names(n)
    nms[dims] = [string(i) for i in 1:size(n, dims)]
    NamedArray(sort(n.array; dims=dims, kws...), tuple(nms...), n.dimnames)
end


text:
# bounds by quantiles, and standardizes and normalizes around median
function TransformStats(data, info)
    q01,q50,q99,iqr = info
    data = max.(data, q01')
    data = min.(data, q99')
    data = (data .- q50') ./ iqr'
    return data
end

# draw neural moments
function NeuralMoments(θ, reps, model::SNMmodel, nnmodel, nninfo)
    z = model.auxstat(θ, reps)
    mean([NeuralMoments(z[i], model, nnmodel, nninfo) for i = 1:reps])
end
# neural moments given statistic
function NeuralMoments(z, model::SNMmodel, nnmodel, nninfo)
    Float64.(nnmodel(TransformStats((z[:])', nninfo)'))
end

# estimate covariance
function EstimateΣ(θ, reps, model::SNMmodel, nnmodel, nninfo)
    z = model.auxstat(θ, reps)
    cov([min.(max.(NeuralMoments(z[i], model, nnmodel, nninfo), model.lb), model.ub) for i = 1:reps])
end

# moments and covariance
function mΣ(θ, reps, model::SNMmodel, nnmodel, nninfo)
    z = model.auxstat(θ, reps)
    Zs = [NeuralMoments(z[i], model, nnmodel, nninfo) for i = 1:reps]
    m = mean(Zs)
    c = cov([min.(max.(Zs[i], model.lb), model.ub) for i = 1:reps])
    m, c
end

# method with identity weight
function H(θ, m, reps, model::SNMmodel, nnmodel, nninfo)
    k = size(θ,1)
    invΣ = Matrix(1.0I, k, k)
    H(θ, m, reps, model, nnmodel, nninfo, invΣ)
end

# log likelihood (GMM-form) with fixed weight matrix
function H(θ, m, reps, model::SNMmodel, nnmodel, nninfo, invΣ::Matrix)
    x = m - NeuralMoments(θ, reps, model, nnmodel, nninfo)
    -0.5*dot(x,invΣ*x)
end

# this is for CUE version
function H(θ, m, reps, model::SNMmodel, nnmodel, nninfo, do_cue::Bool)
    mbar, Σ = mΣ(θ, reps, model, nnmodel, nninfo)
    x = m - mbar
    -0.5*dot(x,inv(Σ)*x)
end


text:
#Ost - "E6" - восток, Nord - "N8" - север, West - "W4" - запад, Sud - "S2" - юг
function main_function7!(r)
    arr=[]
    exit_to_the_lower_left_cell!(r, arr)
    moves!(r, Nord)
    counter_Nord=moves!(r, Sud)
    side=Ost
    mark_line!(r, side)
    for _ in 1:counter_Nord
        if !ismarker(r)
            move!(r, Nord)
            side=inverse(side)
            mark_line!(r, side)
        elseif ismarked(r)
            ismarked(r)
            move!(r, Nord)
            side=inverse(side)
            move!(r, side)
            mark_line!(r, side)
        end
    end
    moves!(r, West)
    moves!(r, Sud)
    return_to_starting_point!(r, arr)
    putmarker!(r)
end
function mark_line!(r, side)
    while !isborder(r, side)
        putmarker!(r)
        move!(r, side)
        if !isborder(r, side)
            move!(r, side)
        end
    end
end
function exit_to_the_lower_left_cell!(r, arr)
    while !(isborder(r, West)&&isborder(r, Sud))
        if !isborder(r, West)
            move!(r, West)
            push!(arr, West)
        end
        if !isborder(r, West)
            move!(r, West)
            push!(arr, West)
        end
    end
end
function return_to_starting_point!(r, arr)
    while length(arr)>0
        side=pop!(arr)
        move!(r, inverse(side))
    end
end
function inverse(side)
    HorizonSide(mod(Int(side)+2, 4))
end
function moves!(r, side)
    counter=0
    while !isborder(r, side)
        move!(r, side)
        counter +=1
    end
    return counter
end


text:
##### MODIFIED FROM https://github.com/ali-ramadhan/DocumenterCitations.jl #######
using Bibliography: xnames, xyear, xlink, xtitle, xin
using Documenter
using Documenter.Anchors
using Documenter.Builder
using Documenter.Documents
using Documenter.Selectors
using Documenter.Utilities
using Documenter.Expanders
using Markdown
using DataStructures: OrderedDict

struct CitationBibliography <: Documenter.Plugin
    bib::OrderedDict{String,<:Bibliography.AbstractEntry}
end

abstract type BibliographyBlock <: Expanders.ExpanderPipeline end

Selectors.order(::Type{BibliographyBlock}) = 12.0  # Expand bibliography last
Selectors.matcher(::Type{BibliographyBlock}, node, page, doc) = Expanders.iscode(node, r"^@bibliography")

linkify(text, link) = isempty(link) ? text : "$text"

function Selectors.runner(::Type{BibliographyBlock}, x, page, doc)
    @info "Expanding bibliography."
    raw_bib = "
"
    for (index,(id, entry)) in enumerate(doc.plugins[CitationBibliography].bib)
        # Add anchor that citations can link to from anywhere in the docs.
        Anchors.add!(doc.internal.headers, entry, entry.id, page.build)

        authors = xnames(entry)
        link = xlink(entry)
        title = xtitle(entry)
        published_in = xin(entry)

        raw_bib *= """

$authors, $(linkify(title, link)), $published_in
"""
    end
    raw_bib *= "\n
"

    page.mapping[x] = Documents.RawNode(:html, raw_bib)
end

abstract type Citations <: Builder.DocumentPipeline end

Selectors.order(::Type{Citations}) = 3.1  # After cross-references

function Selectors.runner(::Type{Citations}, doc::Documents.Document)
    @info "Citations: building citations."
    expand_citations(doc)
end

function expand_citations(doc::Documents.Document)
    for (src, page) in doc.blueprint.pages
        empty!(page.globals.meta)
        for element in page.elements
            expand_citation(page.mapping[element], page, doc)
        end
    end
end

function expand_citation(elem, page, doc)
    Documents.walk(page.globals.meta, elem) do link
        expand_citation(link, page.globals.meta, page, doc)
    end
end

function expand_citation(link::Markdown.Link, meta, page, doc)
    link.url !== "@cite" && return false

    if length(link.text) === 1 && isa(link.text[1], String)
        citation_name = link.text[1]

        if haskey(doc.plugins[CitationBibliography].bib, citation_name)
            entry = doc.plugins[CitationBibliography].bib[citation_name]
            headers = doc.internal.headers
            all_bibs= collect(keys(doc.plugins[CitationBibliography].bib))
            index = findfirst(k -> k==citation_name,all_bibs)
            if Anchors.exists(headers, entry.id)
                if Anchors.isunique(headers, entry.id)
                    # Replace the `@cite` url with a path to the referenced header.
                    anchor   = Anchors.anchor(headers, entry.id)
                    path     = relpath(anchor.file, dirname(page.build))
                    link.text = "[$index]"
                    link.url = string(path, Anchors.fragment(anchor))
                    return true
                else
                    push!(doc.internal.errors, :citations)
                    @warn "'$(entry.id)' is not unique in $(Utilities.locrepr(page.source))."
                end
            else
                push!(doc.internal.errors, :citations)
                @warn "reference for '$(entry.id)' could not be found in $(Utilities.locrepr(page.source))."
            end
        else
            error("Citation not found in bibliography: $(citation_name)")
        end
    else
        error("Invalid citation: $(link.text)")
    end
    return false
end

expand_citation(other, meta, page, doc) = true # Continue to `walk` through element `other`.


text:
include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))
include(joinpath(LOOPVECBENCHDIR, "loadsharedlibs.jl"))

using PrettyTables, BenchmarkTools
struct SizedResults{V <: AbstractVector} <: AbstractMatrix{String}
    results::Matrix{Float64}
    sizes::V
end
function Base.size(sr::SizedResults)
    M, N = size(sr.results)
    N, M+1
end
struct BenchmarkResult{V}
    tests::Vector{String}
    sizedresults::SizedResults{V}
end
function BenchmarkResult(tests, sizes)
    ntests = length(tests); nsizes = length(sizes)
    BenchmarkResult(
        append!(["Size"], tests),
        SizedResults(Matrix{Float64}(undef, ntests, nsizes), sizes)
    )
end
function Base.getindex(br::SizedResults, row, col)
    col == 1 ? string(br.sizes[row]) : string(br.results[col - 1, row])
end
Base.setindex!(br::BenchmarkResult, v, i...) = br.sizedresults.results[i...] = v

const HIGHLIGHT_BEST = Highlighter(
    (br,i,j) -> (j > 1 && maximum(@view(br.results[:, i])) == br.results[j-1,i]),
    foreground = :green
);
function Base.show(io::IO, br::BenchmarkResult)
    pretty_table(
        io, br.sizedresults, br.tests, crop = :none, highlighters = (HIGHLIGHT_BEST,)
    )
end

tothreetuple(i::Int) = (i,i,i)
tothreetuple(i::NTuple{3,Int}) = i
function benchmark_gemm(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "GFort-intrinsic", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        M, K, N = tothreetuple(s)
        C = Matrix{Float64}(undef, M, N)
        A = rand(M, K)
        B = rand(K, N)
        n_gflop = M*K*N*2e-9
        br[1,i] = n_gflop / @belapsed mul!($C, $A, $B)
        Cblas = copy(C)
        br[2,i] = n_gflop / @belapsed jgemm_nkm!($C, $A, $B)
        @assert C ≈ Cblas "Julia gemm wrong?"
        br[3,i] = n_gflop / @belapsed cgemm_nkm!($C, $A, $B)
        @assert C ≈ Cblas "Polly gemm wrong?"
        br[4,i] = n_gflop / @belapsed fgemm_nkm!($C, $A, $B)
        @assert C ≈ Cblas "Fort gemm wrong?"
        br[5,i] = n_gflop / @belapsed fgemm_builtin!($C, $A, $B)
        @assert C ≈ Cblas "Fort intrinsic gemm wrong?"
        br[6,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
        @assert C ≈ Cblas "LoopVec gemm wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end
function benchmark_AtmulB(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "GFort-intrinsic", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        M, K, N = tothreetuple(s)
        C = Matrix{Float64}(undef, M, N)
        At = rand(K, M)
        B = rand(K, N)
        n_gflop = M*K*N*2e-9
        br[1,i] = n_gflop / @belapsed mul!($C, $At', $B)
        Cblas = copy(C)
        br[2,i] = n_gflop / @belapsed jAtmulB!($C, $At, $B)
        @assert C ≈ Cblas "Julia gemm wrong?"
        br[3,i] = n_gflop / @belapsed cAtmulB!($C, $At, $B)
        @assert C ≈ Cblas "Polly gemm wrong?"
        br[4,i] = n_gflop / @belapsed fAtmulB!($C, $At, $B)
        @assert C ≈ Cblas "Fort gemm wrong?"
        br[5,i] = n_gflop / @belapsed fAtmulB_builtin!($C, $At, $B)
        @assert C ≈ Cblas "Fort intrinsic gemm wrong?"
        br[6,i] = n_gflop / @belapsed jAtmulBavx!($C, $At, $B)
        @assert C ≈ Cblas "LoopVec gemm wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end

function benchmark_dot(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        a = rand(s); b = rand(s);
        n_gflop = s * 2e-9
        br[1,i] = n_gflop / @belapsed dot($a, $b)
        dotblas = dot(a, b)
        br[2,i] = n_gflop / @belapsed jdot($a, $b)
        @assert jdot(a,b) ≈ dotblas "Julia dot wrong?"
        br[3,i] = n_gflop / @belapsed cdot($a, $b)
        @assert cdot(a,b) ≈ dotblas "Polly dot wrong?"
        br[4,i] = n_gflop / @belapsed fdot($a, $b)
        @assert fdot(a,b) ≈ dotblas "Fort dot wrong?"
        br[5,i] = n_gflop / @belapsed jdotavx($a, $b)
        @assert jdotavx(a,b) ≈ dotblas "LoopVec dot wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end
function benchmark_selfdot(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        a = rand(s);
        n_gflop = s * 2e-9
        br[1,i] = n_gflop / @belapsed dot($a, $a)
        dotblas = dot(a, a)
        br[2,i] = n_gflop / @belapsed jselfdot($a)
        @assert jselfdot(a) ≈ dotblas "Julia dot wrong?"
        br[3,i] = n_gflop / @belapsed cselfdot($a)
        @assert cselfdot(a) ≈ dotblas "Polly dot wrong?"
        br[4,i] = n_gflop / @belapsed fselfdot($a)
        @assert fselfdot(a) ≈ dotblas "Fort dot wrong?"
        br[5,i] = n_gflop / @belapsed jselfdotavx($a)
        @assert jselfdotavx(a) ≈ dotblas "LoopVec dot wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end
totwotuple(i::Int) = (i,i)
totwotuple(i::Tuple{Int,Int}) = i
function benchmark_gemv(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        M, N = totwotuple(s)
        x = Vector{Float64}(undef, M); A = rand(M, N); y = rand(N);
        n_gflop = M*N * 2e-9
        br[1,i] = n_gflop / @belapsed mul!($x, $A, $y)
        xblas = copy(x)
        br[2,i] = n_gflop / @belapsed jgemv!($x, $A, $y)
        @assert x ≈ xblas "Julia wrong?"
        br[3,i] = n_gflop / @belapsed cgemv!($x, $A, $y)
        @assert x ≈ xblas "Polly wrong?"
        br[4,i] = n_gflop / @belapsed fgemv!($x, $A, $y)
        @assert x ≈ xblas "Fort wrong?"
        br[5,i] = n_gflop / @belapsed jgemvavx!($x, $A, $y)
        @assert x ≈ xblas "LoopVec wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end
function benchmark_dot3(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        M, N = totwotuple(s)
        x = rand(M); A = rand(M, N); y = rand(N);
        n_gflop = M*N * 3e-9
        br[1,i] = n_gflop / @belapsed dot($x, $A, $y)
        dotblas = dot(x, A, y)
        br[2,i] = n_gflop / @belapsed jdot3($x, $A, $y)
        @assert jdot3(x, A, y) ≈ dotblas "Julia dot wrong?"
        br[3,i] = n_gflop / @belapsed cdot3($x, $A, $y)
        @assert cdot3(x, A, y) ≈ dotblas "Polly dot wrong?"
        br[4,i] = n_gflop / @belapsed fdot3($x, $A, $y)
        @assert fdot3(x, A, y) ≈ dotblas "Fort dot wrong?"
        br[5,i] = n_gflop / @belapsed jdot3avx($x, $A, $y)
        @assert jdot3avx(x, A, y) ≈ dotblas "LoopVec dot wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end
function sse!(Xβ, y, X, β)
    mul!(copyto!(Xβ, y), X, β, 1.0, -1.0)
    dot(Xβ, Xβ)
end
function benchmark_sse(sizes)
    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        N, P = totwotuple(s)
        y = rand(N); β = rand(P)
        X = randn(N, P)
        Xβ = similar(y)
        n_gflop = 2e-9*(P*N + 2N)
        br[1,i] = n_gflop / @belapsed sse!($Xβ, $y, $X, $β)
        lpblas = sse!(Xβ, y, X, β)
        br[2,i] = n_gflop / @belapsed jOLSlp($y, $X, $β)
        @assert jOLSlp(y, X, β) ≈ lpblas "Julia wrong?"
        br[3,i] = n_gflop / @belapsed cOLSlp($y, $X, $β)
        @assert cOLSlp(y, X, β) ≈ lpblas "Polly wrong?"
        br[4,i] = n_gflop / @belapsed fOLSlp($y, $X, $β)
        @assert fOLSlp(y, X, β) ≈ lpblas "Fort wrong?"
        br[5,i] = n_gflop / @belapsed jOLSlp_avx($y, $X, $β)
        @assert jOLSlp_avx(y, X, β) ≈ lpblas "LoopVec wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end

function benchmark_exp(sizes)
    tests = ["Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        a = rand(s); b = similar(a)
        n_gflop = 1e-9*s # not really gflops
        br[1,i] = n_gflop / @belapsed @. $b = exp($a)
        baseb = copy(b)
        br[2,i] = n_gflop / @belapsed cvexp!($b, $a)
        @assert b ≈ baseb "Clang wrong?"
        br[3,i] = n_gflop / @belapsed fvexp!($b, $a)
        @assert b ≈ baseb "Fort wrong?"
        br[4,i] = n_gflop / @belapsed @avx @. $b = exp($a)
        @assert b ≈ baseb "LoopVec wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end

function benchmark_aplusBc(sizes)
    tests = ["Julia", "Clang-Polly", "GFort-loops", "LoopVectorization"]
    br = BenchmarkResult(tests, sizes)
    for (i,s) ∈ enumerate(sizes)
        M, N = totwotuple(s)
        a = rand(M); B = rand(M,N); c = rand(N);
        c′ = c'; D = similar(B)
        n_gflop = 2e-9 * M*N
        br[1,i] = n_gflop / @belapsed @. $D = $a + $B * $c′
        Dcopy = copy(D)
        br[2,i] = n_gflop / @belapsed caplusBc!($D, $a, $B, $c)
        @assert D ≈ Dcopy "Polly wrong?"
        br[3,i] = n_gflop / @belapsed faplusBc!($D, $a, $B, $c)
        @assert D ≈ Dcopy "Fort wrong?"
        br[4,i] = n_gflop / @belapsed @avx @. $D = $a + $B * $c′
        @assert D ≈ Dcopy "LoopVec wrong?"
        # if i % 10 == 0
            # percent_complete = round(100i/ length(sizes), sigdigits = 4)
            # @show percent_complete
        # end
    end
    br
end


text:
# LP_tools.jl
# Contains helpful functions to debug Gurobi output
# Authored By Arthur J Delarue on 6/10/15

function print_iis_gurobi(m::Model)
	"""
	Taken from JuMP examples online. Given a LP instance, uses Gurobi to troubleshoot infeasibility.
	Outputs IIS to find "bad" constraints.
	"""
    grb = MathProgBase.getrawsolver(getInternalModel(m))
    Gurobi.computeIIS(grb)
    numconstr = Gurobi.num_constrs(grb)
    numvar = Gurobi.num_vars(grb)

    iisconstr = Gurobi.get_intattrarray(grb, "IISConstr", 1, numconstr)
    iislb = Gurobi.get_intattrarray(grb, "IISLB", 1, numvar)
    iisub = Gurobi.get_intattrarray(grb, "IISUB", 1, numvar)

    println("Irreducible Inconsistent Subsystem (IIS)")
    println("Variable bounds:")
    for i in 1:numvar
        v = Variable(m, i)
        if iislb[i] != 0 && iisub[i] != 0
            println(getLower(v), " <= ", getName(v), " <= ", getUpper(v))
        elseif iislb[i] != 0
            println(getName(v), " >= ", getLower(v))
        elseif iisub[i] != 0
            println(getName(v), " <= ", getUpper(v))
        end
    end

    println("Constraints:")
    for i in 1:numconstr
        if iisconstr[i] != 0
            println(m.linconstr[i])
        end
    end
    println("End of IIS")
end

function print_iis_gurobi(m::Model, im::Gurobi.GurobiMathProgModel)
    """
    Taken from JuMP examples online. Given a LP instance, uses Gurobi to troubleshoot infeasibility.
    Outputs IIS to find "bad" constraints.
    """
    grb = MathProgBase.getrawsolver(im)
    Gurobi.computeIIS(grb)
    numconstr = Gurobi.num_constrs(grb)
    numvar = Gurobi.num_vars(grb)

    iisconstr = Gurobi.get_intattrarray(grb, "IISConstr", 1, numconstr)
    iislb = Gurobi.get_intattrarray(grb, "IISLB", 1, numvar)
    iisub = Gurobi.get_intattrarray(grb, "IISUB", 1, numvar)

    println("Irreducible Inconsistent Subsystem (IIS)")
    println("Variable bounds:")
    lowerBounds = MathProgBase.getvarLB(im)
    upperBounds = MathProgBase.getvarUB(im)
    for i in 1:numvar
        v = Variable(m, i)
        if iislb[i] != 0 && iisub[i] != 0
            println(lowerBounds[i], " <= ", getName(v), " <= ", upperBounds[i])
        elseif iislb[i] != 0
            println(getName(v), " >= ", lowerBounds[i])
        elseif iisub[i] != 0
            println(getName(v), " <= ", upperBounds[i])
        end
    end

    cLowerBounds = MathProgBase.getconstrLB(im)
    cUpperBounds = MathProgBase.getconstrUB(im)
    println("Constraints:")
    for i in 1:numconstr
        if iisconstr[i] != 0
            println(cLowerBounds[i], " <= ", split(string(m.linconstr[i])," = ")[1], " <= ", cUpperBounds[i])
        end
    end
    println("End of IIS")
end

function writeDataToFile(directory::AbstractString,
    model_type::AbstractString="",
    max_rounds::Int=0,
    min_rides::Int=0,
    radius::Int=0,
    regularization::Bool=false,
    preprocess::Bool=false,
    num_clusters::Int=0,
    sample_size::Int=0,
    turnCost::Float64=0.,
    lambda::Float64=0.)
    """
    Writes relevant information for the LP to a file in the provided directory.
    """
    fileName = string(directory, "/info.txt")
    outputFile = open(fileName, "w")
    write(outputFile, string(fileName, "\n"))
    write(outputFile, string(now(), "\n"))
    write(outputFile, string("Model type    :\t", model_type, "\n"))
    write(outputFile, string("Max iterations:\t", max_rounds, "\n"))
    write(outputFile, string("Min # of rides:\t", min_rides, "\n"))
    write(outputFile, string("Node radius   :\t", radius, "m\n"))
    if regularization
        write(outputFile, "Regularization:\tON\n")
        write(outputFile, string("Lambda        :\t",lambda, "\n"))
    else
        write(outputFile, "Regularization:\tOFF\n")
        write(outputFile, string("Lambda        :\tN/A\n"))
    end
    if preprocess
        write(outputFile, "Preprocessing :\tON\n")
        write(outputFile, string("# of clusters :\t", num_clusters, "\n"))
        write(outputFile, string("Max samp. size:\t", sample_size, "\n"))
    else
        write(outputFile, "Preprocessing :\tOFF\n")
        write(outputFile, "# of clusters :\tN/A\n")
        write(outputFile, "Max samp. size:\tN/A\n")
    end
    write(outputFile, string("Left turn cost:\t", turnCost, "s\n"))
    close(outputFile)
end

function writeDataToFile(directory::AbstractString,
    model_type::AbstractString="",
    max_rounds::Int=0,
    min_rides::Int=0,
    radius::Int=0,
    preprocess::Bool=false,
    num_clusters::Int=0,
    sample_size::Int=0,
    turnCost::Float64=0.)
    """
    Writes relevant information for the LP to a file in the provided directory.
    """
    fileName = string(directory, "/info.txt")
    outputFile = open(fileName, "w")
    write(outputFile, string(fileName, "\n"))
    write(outputFile, string(now(), "\n"))
    write(outputFile, string("Model type    :\t", model_type, "\n"))
    write(outputFile, string("Max iterations:\t", max_rounds, "\n"))
    write(outputFile, string("Min # of rides:\t", min_rides, "\n"))
    write(outputFile, string("Node radius   :\t", radius, "m\n"))
    if preprocess
        write(outputFile, "Preprocessing :\tON\n")
        write(outputFile, string("# of clusters :\t", num_clusters, "\n"))
        write(outputFile, string("Max samp. size:\t", sample_size, "\n"))
    else
        write(outputFile, "Preprocessing :\tOFF\n")
        write(outputFile, "# of clusters :\tN/A\n")
        write(outputFile, "Max samp. size:\tN/A\n")
    end
    write(outputFile, string("Left turn cost:\t", turnCost, "s\n"))
    close(outputFile)
end

function appendDataToFile(directory::AbstractString,
    numBrokenConstraints::Int)
    """
    Adds the number of broken constraints through the relaxation parameter to info.txt
    """
    fileName = string(directory, "/info.txt")
    outputFile = open(fileName, "a")
    write(outputFile, string("Constraints violated by regularization: ", numBrokenConstraints, "\n"))
    close(outputFile)
end

function findWorstPathIndex(paths::Array{Array{Int}}, numExpensiveTurns::Array{Int}, turnCost::Float64, times::AbstractArray{Float64,2})
    """
    Given some paths (and extra turn costs) figure out which path is the longest
    """
    pathTime = zeros(length(numExpensiveTurns))
    for i = 1:length(pathTime)
        pathTime[i] = sum([times[paths[i][a],paths[i][a+1]] for a = 1:(length(paths[i])-1)]) + numExpensiveTurns[i] * turnCost
    end
    return indmax(pathTime)
end

function count_number_edges_used(totalPaths::Array{Any,1}, manhattan::Manhattan)
    """
    Given the array of paths being used, decide how many (and which) edges are being used
    """
    graph = manhattan.network
    edges_in = falses(nv(graph), nv(graph))
    for i = 1:length(totalPaths)
        for j = 1:length(totalPaths[i])
            for k = 1:(length(totalPaths[i][j]) - 1)
                edges_in[totalPaths[i][j][k], totalPaths[i][j][k+1]] = true
            end
        end
    end
    return sum(edges_in)
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing