text:
"""
This algorithm evaluate and solve the Unit Commitment optimization.

Author: Mauricio Hernandez
Date: 12/15/2020
References:
- https://jump.dev/JuMP.jl/0.17/refexpr.html

...
# Arguments
- `filePath::String`: the path if the CSV files with the generators information.
...

# Examples
```julia-repl
julia> solveUC("File name")
"Unit Commitment was solved correctly"
```
"""

#
using CSV, DelimitedFiles, DataFrames
using JuMP, CPLEX
using BenchmarkTools, Logging, Dates

#Enabling debugging code, use ENV["JULIA_DEBUG"] = "" to desable Debugging code
ENV["JULIA_DEBUG"] = "all"

const N_GEN = 145
const N_HRS = 168
const INITIAL_DAY = 1
const FINAL_DAY = 1

#curentTime = Dates.format(now(), "yyyy-mm-dd_HH-MM-SS");
# Logging file
io_log = open(string(".//outputs//log//log_", Dates.format(now(), "yyyy-mm-dd_HH-MM-SS") ,".txt"), "w")
logger = SimpleLogger(io_log)
flush(io_log)

global_logger(logger)

@info("Log file for UC model for $N_GEN generators")

#Loading input data
dfGenerator = CSV.read(".//inputs//data_generators.csv", DataFrame)
dfDemand = CSV.read(".//inputs//demand_reserves.csv", DataFrame)

@info("Data from generators and demand downloaded correcly")

typeFileAccess = "w+";

t1 = time_ns()
@time begin
for day in INITIAL_DAY:FINAL_DAY
    println("---------------------------")
    @info println("Solving UC for day: $day")
    println("---------------------------")

    dfDayDemand = dfDemand[(day-1)*N_HRS+1:(N_HRS*day)+1, [:Hour, :Demand, :Reserves]];

    @debug begin
      pathTempDemandFile = string(".//outputs//demandTEMP_day", day, ".csv");
      CSV.write(pathTempDemandFile, dfDayDemand)
    end

    model = Model(CPLEX.Optimizer)
    set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-5)

    @variable(model, genOnOff[1:N_GEN, 1:N_HRS], Bin)
    @variable(model, genOut[1:N_GEN, 1:N_HRS] >= 0) # Power bound constraint Equation 4
    @variable(model, genStartUp[1:N_GEN, 1:N_HRS], Bin)
    @variable(model, genShutDown[1:N_GEN, 1:N_HRS], Bin)

    #Setting objective function
    @objective(model, Min, sum(sum(dfGenerator.FixedCost[g]*genOnOff[g,t] + dfGenerator.VariableCost[g]*genOut[g,t] +
               dfGenerator.StartUpCost[g]*genStartUp[g,t] + dfGenerator.ShutdownCost[g]*genShutDown[g, t] for g in 1:N_GEN)
               for t in 1:N_HRS) )

    #Generators should have an initial and defined value at t=0
    @constraint(model, conInit[g=1:N_GEN], genOnOff[g,1] == dfGenerator.Uini[g])

    #Initial values of Generator shutdown and startup values should be defined
    if day == 1
      @constraint(model, conShutDn[g=1:N_GEN], genShutDown[g,1] == 0)
      @constraint(model, conShutUp[g=1:N_GEN], genStartUp[g,1] == 0)
    end

    #Logical Conditions
    @constraint(model, conLogical[t=1:N_HRS, g=1:N_GEN], genStartUp[g,t] + genShutDown[g, t] <= 1)

    #Start-up and Shut-down - Equation 3
    @constraint(model, conStartUpAndDn[t=2:N_HRS, g=1:N_GEN], (genStartUp[g,t] - genShutDown[g,t]) ==
               (genOnOff[g,t] - genOnOff[g,t-1]))

    #Power bounds - Equation 4
    @constraint(model, conMaxPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] <= dfGenerator.MaxPowerOut[g]*genOnOff[g,t] )
    @constraint(model, conMinPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] >= dfGenerator.MinPowerOut[g]*genOnOff[g,t] )

    #Ramping limits
    #Ramp Up - Equation 5
    @constraint(model, conRampUp[t=2:N_HRS, g=1:N_GEN], genOut[g,t] - genOut[g,t-1] <=
                dfGenerator.RampUpLimit[g]*genOnOff[g, t-1] + dfGenerator.RampStartUpLimit[g]*genStartUp[g,t])
    if day == 1
        @constraint(model, conIniRampUp[g=1:N_GEN], genOut[g,1] - dfGenerator.PowerIni[g] <=
                dfGenerator.RampUpLimit[g]*dfGenerator.Uini[g] + dfGenerator.RampStartUpLimit[g]*genStartUp[g,1])
    end

    #Ramp Down - Equation 6
    @constraint(model, conRampDown[t=2:N_HRS, g=1:N_GEN], (genOut[g,t-1] - genOut[g,t]) <=
               (dfGenerator.RampDownLimit[g]*genOnOff[g,t] + dfGenerator.RampShutDownLimit[g]*genShutDown[g,t]))

   if day == 1
      @constraint(model, conIniRampDown[g=1:N_GEN], dfGenerator.PowerIni[g] - genOut[g,1]  <=
              dfGenerator.RampDownLimit[g]*genOnOff[g,1] + dfGenerator.RampShutDownLimit[g]*genShutDown[g,1])
   end

    #Power balance - Equation 7
    @constraint(model, conDemand[t=1:N_HRS], sum(genOut[g, t] for g=1:N_GEN) == dfDayDemand.Demand[t])

    #Reserves - Equation 8
    @constraint(model, conReserves[t=1:N_HRS], sum(dfGenerator.MaxPowerOut[g]*genOnOff[g,t] for g=1:N_GEN) >= (dfDayDemand.Demand[t] + dfDayDemand.Reserves[t]))

    write(io_log, "Optimization Execution for day $day: \n")
    @info "Time elapsed (s):" @elapsed optimize!(model);

    println("------------------------------------")
    println("------- TERMINATION STATUS -------")
    if termination_status(model) == MOI.OPTIMAL
        #optimal_solution = value.(genOnOff)
        @info "Optimization objective: " optimal_objective = JuMP.objective_value(model)
        println("Optimal result")
    elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
        #suboptimal_solution = value.(genOnOff)
        suboptimal_objective = JuMP.objective_value(model)
        println("Time limit was exceeded")
        close(io_log)
    else
        error("The model was not solved correctly.")
        close(io_log)
    end

    #Printing general results
    println("------------------------------------")
    println("------- OBJECTIVE VALUE -------")
    println("Objective value for day $day: ", JuMP.objective_value(model))

    println("------------------------------------")
    println("------- PRIMAL STATUS -------")
    println(primal_status(model))

    println("------------------------------------")
    println("------- DUAL STATUS -------")
    println(JuMP.dual_status(model))

    println("Day: $day solved")
    println("---------------------------")

    #Storing decision variables results
    open(".//outputs//resultsObjective.csv", typeFileAccess) do io
               write(io,  string("Objective value for day: ", day, ",", JuMP.objective_value(model), "/n"))
               close(io)
        end;

    open(".//outputs//resultsGenOut.csv", typeFileAccess) do io
               writedlm(io,  transpose(JuMP.value.(genOut)), ',')
        end;

    open(".//outputs//resultsGenCommit.csv", typeFileAccess) do io
                   writedlm(io, transpose(JuMP.value.(genOnOff)), ',')
             end;

    open(".//outputs//resultsGenStartUp.csv", typeFileAccess) do io
                  writedlm(io, transpose(JuMP.value.(genStartUp)), ',')
            end;

    open(".//outputs//resultsGenShutDown.csv", typeFileAccess) do io
                  writedlm(io, transpose(JuMP.value.(genShutDown)), ',')
             end;

#=    arrGenOnOff = JuMP.value.(genOnOff)
    println("Type JuMP.value.(genOnOff)")
    println(typeof(arrGenOnOff))
    println("Number of columns")
    println(size(arrGenOnOff,2))
    println("DF Generator Size")
    println(size(JuMP.value.(genOnOff)))
    println("JuMP.value.(genOnOff) first column size")
    println(size(JuMP.value.(genOnOff)[:,size(arrGenOnOff,2)]))
=#
    dfGenerator.Uini = JuMP.value.(genOnOff)[:,size(JuMP.value.(genOnOff),2)]

  global typeFileAccess = "a+";
end
end #time
#close(ioGenOut)
#Releasing memory

t2 = time_ns()
elapsedTime = (t2 -t1)/1.0e9;
@info "UC Optimization solved in (s):" elapsedTime = (t2 -t1)/1.0e9;

#Cleaning garbage collector
GC.gc();
close(io_log); #

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
@testset "@lencheck" begin
    x, y = ones(10), ones(10)
    z = ones(9)

    @test (@lencheck 10 x) ≡ nothing

    @test (@lencheck 10 x y) ≡ nothing

    @test_throws DimensionError @lencheck 10 z

    @test_throws DimensionError @lencheck 10 x z
end

@testset "@contract" begin
    abstract type AbstractFoo end

    struct Foo <: AbstractFoo end

    func(f::AbstractFoo) = @contract AbstractFoo :func

    # testing that func(T) where T<:AbstractFoo throws until method T exists
    @test_throws ContractError func(Foo())

    func(f::Foo) = nothing

    @test func(Foo()) ≡ nothing
end

function test_restricted_type(type)

    @test_throws AccessError type.afield

    @test propertynames(type) ≡ ()

    return nothing
end

@testset "@restrict" begin
    struct Bar
        x::Any
    end

    Base.getproperty(b::Bar, s::Symbol) = @restrict Bar

    Base.propertynames(b::Bar) = ()

    bar = Bar(nothing)

    test_restricted_type(bar)
end

ext: jl


text:
function get_sep(v_buf::Array{UInt8,1}, vi::Int8)
  z = zero(Int8)
  o = one(Int8)
  i = z
  while i < vi
    i += o
    y = getindex(v_buf, i)
    if y != 0x20
      return y
    end
  end
end

function geocsv_mkid(v_buf::Array{UInt8,1}, vi::Int8)
  # SID
  i = 0x00
  while i < vi
    i += 0x01
    if getindex(v_buf, i) == 0x5f
      setindex!(v_buf, 0x2e, i)
    end
  end
  return String(v_buf[0x01:vi])
end

function geocsv_assign!(C::SeisChannel,
                        k_buf::Array{UInt8,1},
                        v_buf::Array{UInt8,1},
                        ki::Int8,
                        vi::Int8)
  o = one(Int8)
  k = String(k_buf[o:ki])
  ptr = pointer(v_buf[o:vi])

  if k == "sample_rate_hz"
    setfield!(C, :fs, parse(Float64, unsafe_string(ptr)))
  elseif k == "latitude_deg"
    setfield!(C.loc, :lat, parse(Float64, unsafe_string(ptr)))
  elseif k == "longitude_deg"
    setfield!(C.loc, :lon, parse(Float64, unsafe_string(ptr)))
  elseif k == "elevation_m"
    setfield!(C.loc, :el, parse(Float64, unsafe_string(ptr)))
  elseif k == "azimuth_deg"
    setfield!(C.loc, :az, parse(Float64, unsafe_string(ptr)))
  elseif k == "dip_deg"
    setfield!(C.loc, :inc, 90.0-parse(Float64, unsafe_string(ptr)))
  elseif k == "depth_m"
    setfield!(C.loc, :dep, parse(Float64, unsafe_string(ptr)))
  elseif k == "scale_factor"
    setfield!(C, :gain, parse(Float64, unsafe_string(ptr)))
  elseif k == "scale_frequency_hz"
    C.misc[k] = parse(Float64, unsafe_string(ptr))
  elseif k == "scale_units"
    setfield!(C, :units, lowercase(unsafe_string(ptr)))
  else
    C.misc[k] = unsafe_string(ptr)
  end
  return nothing
end

function mkhdr(io::IO, c::UInt8, k_buf::Array{UInt8,1}, v_buf::Array{UInt8,1})
  k = true
  o = one(Int8)
  i = zero(Int8)
  j = zero(Int8)

  # skip space after a new line
  while c == 0x20; c = fastread(io); end

  while c != 0x0a
    if c == 0x23
      c = fastread(io)
      while c == 0x20
        c = fastread(io)
      end
    # transition at 0x3a
    elseif c == 0x3a && k == true
      k = false
      c = fastread(io)
      while c == 0x20
        c = fastread(io)
      end
    elseif k
      i += o
      setindex!(k_buf, c, i)
      c = fastread(io)
    else
      j += o
      setindex!(v_buf, c, j)
      c = fastread(io)
    end
  end
  return i,j
end

function read_geocsv_slist!(S::GphysData, io::IO)
  o = one(Int16)
  oo = one(Int64)
  z = zero(Int16)

  c = 0x00
  j = z
  sep = 0x2c
  t = oo
  t_old = oo
  t_exp = 5
  Δ = oo
  Δ_gap = oo
  i = oo
  x = zero(Float32)
  k = ""
  v = ""
  nx = zero(UInt64)

  C = SeisChannel()
  C.loc = GeoLoc()
  X = Array{Float32,1}(undef, 0)

  k_buf = Array{UInt8,1}(undef, 80)
  v_buf = Array{UInt8,1}(undef, 80)
  t_buf = Array{UInt8,1}(undef, 32)
  t_ptr = pointer(t_buf)

  reading_data = false
  is_float = false

  # Time structure
  tm = TmStruct()

  while !fasteof(io)
    c = fastread(io)

    # new line ----------------------------------
    (c == 0x0a) && continue

    # parse header ------------------------------
    if c == 0x23
      if reading_data == true

      # '#' after newline
        C.t = vcat(C.t, [i-1 zero(Int64)])
        append!(C.x, X)
        push!(S, C)
        C = SeisChannel()
        C.loc = GeoLoc()
        i = oo
        reading_data = false
      end
      (ki, vi) = mkhdr(io, c, k_buf, v_buf)

      if ki == Int8(9)
        sep = get_sep(v_buf, vi)
      elseif ki == Int8(12)
        if k_buf[1] == 0x73 && k_buf[2] == 0x61
          nx = buf_to_uint(v_buf, vi)
        else
          geocsv_assign!(C, k_buf, v_buf, ki, vi)
        end
      elseif ki == Int8(3)
        id = geocsv_mkid(v_buf, vi)
        ii = findid(S, id)
        if ii > 0
          C = pull(S, ii)
        else
          C.id = id
        end
        Nt = size(C.t, 1)
        if Nt > 0
          t_old = endtime(C.t, round(Int64, sμ/C.fs))
          i = C.t[Nt, 1] + 1
          if C.t[Nt, 2] == 0
            C.t = C.t[1:Nt-1,:]
          end
        else
          t_old = oo
        end
      else
        geocsv_assign!(C, k_buf, v_buf, ki, vi)
      end

    # any other character after newline
    else
      if reading_data == false
        # transition from header to data
        t_str = get(C.misc, "start_time", "1970-01-01T000000.000000Z")
        t = String.(split(t_str, ('T', '.')))
        ts = (Date(t[1]).instant.periods.value)*86400000000 +
              div(Time(t[2]).instant.value, 1000) -
              dtconst

        #= weirdly, *much* more efficient than *either*
        Array{Float32,1}(undef, 0) with push! or
        Array{Float32,1}(undef, x) s setindex!
        =#
        X = Float32[]; sizehint!(X, nx)
        Δ = round(Int64, sμ/getfield(C, :fs))
        if isempty(C.t)
          C.t = [1 ts]
        else
          C.t = vcat(C.t, [i ts-t_old-Δ])
        end
        reading_data = true
        while is_u8_digit(c) == false
          c = fastread(io)
        end
      end
      x = stream_float(io, c)
      push!(X, x)
      i += 1
    end
  end
  append!(C.x, X)
  push!(S, C)
  return nothing
end

function read_geocsv_tspair!(S::GphysData, io::IO)
  o = one(Int16)
  oo = one(Int64)
  z = zero(Int16)

  c = 0x00
  j = z
  sep = 0x2c
  t = oo
  t_old = oo
  t_exp = 5
  Δ = oo
  Δ_gap = oo
  i = oo
  x = zero(Float32)
  k = ""
  v = ""
  nx = zero(UInt64)

  C = SeisChannel()
  C.loc = GeoLoc()
  X = Array{Float32,1}(undef, 0)

  k_buf = Array{UInt8,1}(undef, 80)
  v_buf = Array{UInt8,1}(undef, 80)
  t_buf = Array{UInt8,1}(undef, 32)
  t_ptr = pointer(t_buf)

  read_state = 0x00
  reading_data = false
  is_float = false

  # Time structure
  tm = TmStruct()

  # read_state:
  # 0x00    new line
  # 0x01    hdr (subroutine)
  # 0x02    time
  # 0x03    fractional-second
  # 0x04    data

  while !fasteof(io)
    c = fastread(io)

    # new line ----------------------------------
    if c == 0x0a
      read_state = 0x00

      # No parsing of c
      continue
    end

    # determine next read state -----------------
    if read_state == 0x00

      # '#' after newline
      if c == 0x23
        # transition from data to header
        if reading_data == true

          # finish current SeisChannel
          C.t = vcat(C.t, [i-1 zero(Int64)])
          append!(C.x, X)
          push!(S, C)
          C = SeisChannel()
          C.loc = GeoLoc()
          i = oo
          reading_data = false
        end
        (ki, vi) = mkhdr(io, c, k_buf, v_buf)
        if ki == Int8(9)
          sep = get_sep(v_buf, vi)
        elseif ki == Int8(12)
          if k_buf[1] == 0x73 && k_buf[2] == 0x61
            nx = buf_to_uint(v_buf, vi)
          else
            geocsv_assign!(C, k_buf, v_buf, ki, vi)
          end
        elseif ki == Int8(3)
          id = geocsv_mkid(v_buf, vi)
          ii = findid(S, id)
          if ii > 0
            C = pull(S, ii)
          else
            C.id = id
          end
          Nt = size(C.t, 1)
          if Nt > 0
            t_old = endtime(C.t, round(Int64, sμ/C.fs))
            i = C.t[Nt, 1] + 1
            if C.t[Nt, 2] == 0
              C.t = C.t[1:Nt-1,:]
            end
          else
            t_old = oo
          end
        else
          geocsv_assign!(C, k_buf, v_buf, ki, vi)
        end
        read_state = 0x00

        continue

      # any other character after newline
      else
        read_state = 0x02
        fill!(t_buf, 0x20)
        j = z

        # Requires transition from header to time
        if reading_data == false
          #= weirdly, *much* more efficient than *either*
          Array{Float32,1}(undef, 0) with push! or
          Array{Float32,1}(undef, x) s setindex!
          =#
          X = Float32[]; sizehint!(X, nx)
          Δ = round(Int64, sμ/getfield(C, :fs))
          Δ_gap = div(3*Δ,2)

          # Flag that we're now reading data
          reading_data = true
        end
      end
    end

    # 0x02 = time state -------------------------
    # elseif read_state == 0x02
    if read_state == 0x02

      # '.'
      if c == 0x2e
        # only happens at a decimal in a time block
        ccall(:strptime, Cstring, (Cstring, Cstring, Ref{TmStruct}), t_ptr, "%FT%T", tm)
        t = (Sys.iswindows() ? ccall(:_mkgmtime, Int64, (Ref{TmStruct},), tm) : ccall(:timegm, Int64, (Ref{TmStruct},), tm))* 1000000
        read_state = 0x03
        continue

      else
        j += o
        setindex!(t_buf, c, j)
      end

    # 0x03 = fractional-second state ------------
    elseif read_state == 0x03

      # ',', typically
      if c == sep
        t_exp = 5
        j = z

        if t-t_old > Δ_gap
          if t_old == oo
            C.t = vcat(C.t, [i t])
          else
            C.t = vcat(C.t, [i t-t_old-Δ])
          end
        end
        t_old = t

        read_state = 0x04
        continue

      # intent: ignore the time zone, users can fix manually
      # zero out t_exp when we reach timezone crap
      elseif c == 0x02b || c == 0x02d || c == 0x05a
        t_exp = 0
      elseif t_exp > oo
        t += Int64(c-0x30)*10^t_exp
        t_exp -= oo
      else
        continue
      end

    # 0x04 = data state -------------------------
    elseif read_state == 0x04
      x = stream_float(io, c)
      push!(X, x)
      read_state = 0x00
      i += 1
    else
      error("indeterminate read state")
    end

  end
  C.t = vcat(C.t, [i-1 zero(Int64)])
  append!(C.x, X)
  push!(S, C)
  return nothing
end

function read_geocsv_file!(S::GphysData, fname::String, tspair::Bool, memmap::Bool)
  io = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")
  if tspair == true
    read_geocsv_tspair!(S, io)
  else
    read_geocsv_slist!(S, io)
  end
  close(io)
  return nothing
end


text:
@userplot Calibration_Plot
@userplot Phylogeny_Calibration

@recipe function f(c::Calibration_Plot)
    legend := false
    grid --> false
    seriescolor --> "#A9A9A9"
    markerstrokecolor --> "#A9A9A9"
    xguide --> "observed event percentage"
    yguide --> "bin midpoint"
    formatter := x -> string(Int(round(x * 100))) * "%"

    (c, c.args...)
end

@recipe function f(
    ::Calibration_Plot, result::HLAModelResult
)

    N = result.sf.data["N"]
    rs = result.sf.data["rs"]
    idx = result.sf.data["idx"]
    sf = stanfit(result)
    d = sf.data
    posterior = extract(sf)

    theta_pred = Vector{Float64}(undef, N)
    y = result.sf.data["y"]

    @threads for i in 1:N
        theta = mean(theta_i(sf, posterior, i))
        theta_pred[i] = theta
    end

    @series begin
        seriestype := :path
        linecolor := "black"
        linestyle := :dash
        [0, 1], [0, 1]
    end

    @series begin
        seriestype := :calibration
        theta_pred, y
    end
end

@recipe function f(c::Phylogeny_Calibration)
    legend := false
    grid --> false
    seriescolor --> "#B2001D"
    markerstrokecolor --> "#B2001D"
    xguide --> "observed event percentage"
    yguide --> "bin midpoint"
    formatter := x -> string(Int(round(x * 100))) * "%"

    (c, c.args...)
end

@recipe function f(
    ::Phylogeny_Calibration, result::HLAModelResult
)

    if !("phy" ∈ keys(result.sf.data))
        error("HLAModelResult does not contain phylogeny information.")
    end

    N = result.sf.data["N"]
    rs = result.sf.data["rs"]
    idx = result.sf.data["idx"]

    theta = [result.sf.data["phy"][rs[i], idx[i]] for i in 1:N]
    y = result.sf.data["y"]

    @series begin
        seriestype := :calibration
        theta, y
    end

    @series begin
        seriestype := :path
        linecolor := "black"
        linestyle := :dash
        [0, 1], [0, 1]
    end
end

@recipe function f(::Type{Val{:calibration}}, x, y, z)
    bins = plotattributes[:bins] == :auto ? 10 : plotattributes[:bins]
    df = binned_intervals(x, Bool.(y), bins = bins)

    # error bars
    for row in eachrow(df)
        @series begin
            seriestype := :path
            linecolor := plotattributes[:linecolor]
            x := [row[:expected], row[:expected]]
            y := [row[:lower], row[:upper]]
            ()
        end

    end

    # scatter points
    seriestype := :scatter
    x := df[!, :expected]
    y := df[!, :observed]
end

function binned_intervals(
    theta, y;
    bins::Int = 10, lower::Real = 0.025, upper = 0.975
)
    check_calibration_arguments(theta, y)

    sort_idx = sortperm(theta)
    theta = theta[sort_idx]
    y = y[sort_idx]

    df = DataFrame(
        expected = Float64[], observed = Float64[],
        lower = Float64[], upper = Float64[]
    )

    cutpoints = range(0, 1, length = bins)

    for (lo, hi) in zip(cutpoints[1:end-1], cutpoints[2:end])
        idx = findall(x -> lo <= x <= hi, theta)
        y_p = y[idx]

        push!(df, [mean(y_p), lo + (hi - lo) / 2, lo, hi])
    end

    return df
end

function check_calibration_arguments(x, y)
    if length(x) != length(y)
        error("x ($(length(x))) and y ($(length(y))) must be of same length.")
    elseif !all(0 .<= x .<= 1)
        error("elements of x must be between 0 and 1.")
    end
end

function binned_intervals(result::HLAModelResult; bins = 20, samples = 200)
    stan_input = Escape.stan_input(result)

    theta = posterior_predictions(result.sf, stan_input, samples = samples)
    binned_indices = Escape.binned_indices(theta, bins = bins)
    predictions = rand.(Bernoulli.(theta))

    df = DataFrame(
        expected = Float64[], observed = Float64[],
        lower = Float64[], upper = Float64[]
    )

    for indices in binned_indices
        length(indices) > 0 || continue

        expected = map(x -> mean(predictions[x, indices]), 1:samples)
        observed = mean(stan_input["y"][indices])
        lower = minimum(expected)
        upper = maximum(expected)

        push!(df, [mean(expected), observed, lower, upper])
    end

    return df
end

function binned_indices(theta::AbstractMatrix; bins = 20)
    cutpoints = range(0, 1, length = bins)
    bins = [Int[] for i in 1:bins]

    for i in 1:size(theta, 2)
        m = median(theta[:, i])

        for (j, (low, high)) in enumerate(zip(cutpoints[1:end-1], cutpoints[2:end]))
            if low <= m <= high
                push!(bins[j], i)
            end
        end
    end

    return bins
end

function posterior_predictions(sf::StanInterface.Stanfit, stan_input::Dict; samples = 200)
    posterior = StanInterface.extract(sf)
    theta = Matrix{Float64}(undef, samples, stan_input["N"])

    p = Progress(size(theta, 2), desc = "theta_i ")

    @threads for i in 1:size(theta, 2)
        theta[:, i] .= Escape.theta_i(sf, posterior, i)[1:samples]
        next!(p)
    end

    return theta
end


text: include("tripartite/tripartite.jl")


text:
module Portfolio

import ..Clarus

function cash(;params...)
  return Clarus.api_request("Portfolio","Cash";params...)
end

function cashbydate(;params...)
  return Clarus.api_request("Portfolio","CashByDate";params...)
end

function fixings(;params...)
  return Clarus.api_request("Portfolio","Fixings";params...)
end

function mtm(;params...)
  return Clarus.api_request("Portfolio","MTM";params...)
end

function notional(;params...)
  return Clarus.api_request("Portfolio","Notional";params...)
end

function summary(;params...)
  return Clarus.api_request("Portfolio","Summary";params...)
end

function trades(;params...)
  return Clarus.api_request("Portfolio","Trades";params...)
end

end #MODULE-END


text:
struct Data
    name::String

    capacity::Int64
    vehicles::Int64

    vertices::Vector{Vertex}
    edges::Vector{Edge}
    requireds::Vector{Edge}

    lb::Int64
    ub::Int64
end


text:
"""
   out = adaptive_pg(theta0, prior, state, n; kwargs...)

Particle Gibbs (PG) with adaptive proposal on the parameters.

# Arguments
- `theta0`: Initial parameter vector
- `prior`: Function returning prior log density values for parameters
- `state`: `SMCState` data structure
- `n`: Number of iterations

# Optional keyword arguments
- `b`: Burn-in length; default `0.1n`
- `thin`: Thinning interval; default `1` (no thinning)
- `save_paths`: Whether to save (samples of) latent trajectories; default `false`
- `backward_sampling`: Whether to perform backward sampling; default `true`
- `show_progress`: Show progress every `show_progress` seconds (default: `false`).
- `algorithm`: The adaptation algorithm; either Robust Adaptive Metropolis
   `:ram` (default) or Adaptive Scaling within Adaptive Metropolis `:aswam`.

The hidden Markov model and SMC state is defined in `state`; see `?SMCState`.

The output `out` is a `NamedTuple` containing the simulation output: `out.acc`
is the mean acceptance rate, `out.Theta` contains the simulated parameter
values (each column is a parameter vector). If requested (by `save_paths=true`),
 `out.X[i][k]` contains the simulated state corresponding to `out.Theta[:,i]`
 at time `k`.
"""
function adaptive_pg(theta0::ParamT, prior::Function, state, n::Int;
    b::Int=Int(ceil(0.1n)), thin::Int=1,
    save_paths::Bool=false, acc_opt::FT=FT(0.234), backward_sampling::Bool=true,
    show_progress::Real=false, algorithm=:ram) where {FT <: AbstractFloat, ParamT<:AbstractVector{FT}}

    if isnothing(state.lM)
        error("SMC state `state` must supply log-transition density!")
    end

    # Run the particle filter for HMM w/ param theta0:
    _set_model_param!(state, theta0)
    _run_smc!(state)
    # ... and initialise reference trajectory to a sample from it:
    _pick_particle!(state)

    # Initialise random walk Metropolis state
    r = RWMState(theta0)
    # Initialise adaptation
    if algorithm == :ram
        s = RobustAdaptiveMetropolis(theta0, acc_opt)
    elseif algorithm == :aswam
        s = AdaptiveScalingWithinAdaptiveMetropolis(theta0, acc_opt)
    else
        error("Unknown algorithm: ", algorithm)
    end
    # Allocate space for simulated parameters
    nsim = Int(floor((n-b)/thin))
    Theta = zeros(FT, length(theta0), nsim)
    acc = zero(FT)
    if save_paths
        X = _init_path_storage(state, nsim)
    else
        X = Missing
    end
    if show_progress > 0
        progress = Progress(n, Float64(show_progress))
    else
        progress = NullProgress()
    end
    for k = 1:n
        # Propose parameter update:
        ###########################
        draw!(r, s); pr_ = prior(r.y)
        if pr_ > -Inf
            p = prior(r.x) + _reference_log_likelihood(state)
            _set_model_param!(state, r.y)
            p_ = pr_ + _reference_log_likelihood(state)
            #alpha = min(one(FT), exp(p_ - p))
            alpha = _accept_prob(p, p_, FT)
            if rand() <= alpha
                accept!(r)
                acc += 1
            else
                _set_model_param!(state, r.x)
            end
        else # Early rejection if prior fails:
            alpha = zero(FT)
        end
        adapt!(s, r, alpha, k)

        # State update with CPF (& optional backward sampling)
        ############################
        _run_csmc!(state, backward_sampling)

        @inbounds if k>b && rem(k-b, thin)==0
            i = Int((k-b)/thin)
            Theta[:,i] = r.x
            if save_paths
                _save_reference!(state)
                _copy_reference!(X[i], state)
            end
        end
        next!(progress)
    end
    (Theta=Theta, acc=acc/n, X=X, adapt_state=s, theta0=theta0)
end


text:
@testset "OpenStreetMap space" begin
    space = OpenStreetMapSpace(TEST_MAP)
    @test length(space.s) == 1799
    @test sprint(show, space) ==
          "OpenStreetMapSpace with 1456 roadways and 1799 intersections"

    Random.seed!(689)
    model = ABM(Agent10, space)

    @test osm_random_road_position(model) != osm_random_road_position(model)
    intersection = random_position(model)
    @test intersection[1] == intersection[2]
    @test intersection[3] == 0.0
    ll = osm_latlon(intersection, model)
    @test intersection == osm_intersection(ll, model)

    start_latlon = (39.534773980413505, -119.78937575923226)
    start_i = osm_intersection(start_latlon, model)
    i_diff = sum(abs.(osm_latlon(start_i, model) .- start_latlon))
    start_r = osm_road(start_latlon, model)
    r_diff = sum(abs.(osm_latlon(start_r, model) .- start_latlon))
    @test i_diff > r_diff

    finish_latlon = (39.52530416953533, -119.76949287425508)
    finish_i = osm_intersection(finish_latlon, model)
    finish_r = osm_road(finish_latlon, model)

    route = osm_plan_route(start_r, finish_r, model)
    @test length(route) == 20
    add_agent!(start_r, model, route, finish_r)
    add_agent!(finish_i, model, [], finish_i)

    @test osm_latlon(model[2], model) == osm_latlon(finish_i[1], model)
    np = nearby_positions(model[2], model)
    @test length(np) == 3
    @test osm_latlon(np[1], model) == (39.5259267, -119.76681090000002)

    rand_road = osm_random_road_position(model)
    rand_intersection = random_position(model)
    rand_route = osm_plan_route(rand_road, rand_intersection, model; return_trip = true)
    @test length(rand_route) == 18

    @test osm_road_length(model[1].pos, model) ≈ 106.49016546202557
    @test osm_road_length(finish_r[1], finish_r[2], model) ≈ 395.8120895006937

    @test length(osm_plan_route(start_r[2], finish_r[1], model)) == 22
    @test length(osm_plan_route(start_r, finish_r[1], model)) == 21
    @test length(osm_plan_route(start_r[2], finish_r, model)) == 21
    @test osm_plan_route(start_r, finish_r[1], model) !=
          osm_plan_route(start_r[2], finish_r, model)
    @test length(osm_plan_route(start_r, finish_r, model)) == 20

    @test osm_map_coordinates(model[1], model) == (-2904.337825035879, 1444.5549163962387)
    @test osm_map_coordinates(model[2], model) == (-1475.8022717621393, 461.34058791100955)

    @test !osm_is_stationary(model[1])
    move_agent!(model[1], model, 403.2)
    @test length(model[1].route) == 16
    move_agent!(model[1], model, 1500)
    @test isempty(model[1].route)
    @test !osm_is_stationary(model[1])
    move_agent!(model[1], model, 200)
    @test !osm_is_stationary(model[1])
    move_agent!(model[1], model, 200)
    @test !osm_is_stationary(model[1])
    move_agent!(model[1], model, 200)
    @test osm_is_stationary(model[1])

    @test osm_is_stationary(model[2])
    @test move_agent!(model[2], model, 50) == nothing

    add_agent!(start_r, model, route, finish_r)
    osm_random_route!(model[3], model)
    @test model[3].destination != finish_r
    move_agent!(model[3], model, 100.0)

    for i in 1:5
        s = (start_r[1:2]..., start_r[3] - i)
        route = osm_plan_route(s, finish_r, model)
        add_agent!(s, model, route, finish_r)
    end

    @test sort!(nearby_ids(model[6], model, 2)) == [4, 5, 7, 8]
    @test sort!(nearby_ids(model[6], model, 800.0)) == [3, 4, 5, 7, 8]
end


text:
function h2vtable(m, postag::String)
        ih = findfirst(isequal(postag), m.tags)
        if isnothing(ih) error("$(postag) not found in ctb") end
        vs = collect(m.h2v[ih])
        sort!(vs; by=last, rev=true)
        nv = length(vs)
        mat = Matrix{Any}(undef, (nv, 3))
        word_ids, probs = map(first, vs), exp.(map(last, vs))
        mat[:,1] = m.words[word_ids]
        mat[:,2] = map(p -> trunc(p, digits=6), probs)
        mat[:,3] = map(i -> (i > length(m.words) - m.usr_words ? "user.dict" : "CTB"), word_ids)
        return UselessTable(mat; cnames=["word", "prob.", "source"], foots=["$(postag) has $(length(vs)) words", "Unknown words porb. = $(1. - sum(probs))"])
end

function v2htable(m, word::String)
        word_id = findfirst(isequal(word), m.words)
        if isnothing(word_id) error("$(word) not found in m") end
        np = length(m.tags)
        hs = [(m.tags[first(p)], haskey(last(p), word_id) ? exp(last(p)[word_id]) : 0.) for p in enumerate(m.h2v)]
        sort!(hs; by=last, rev=true)
        mat = Matrix{Any}(undef, (np, 1))
        mat[:,1] = map(d -> trunc(d, digits=6), map(last, hs))
        return UselessTable(mat; cnames=["Prob."], rnames=map(first, hs), topleft="From\\To")
end

function h2htable(m)
        np = length(m.tags)
        mat = map(d->Int(round(exp(d)*100)), m.h2h)
        return UselessTable(mat; cnames=m.tags, rnames=m.tags, topleft="From\\To(%)")
end

function hprtable(m)
        np = length(m.tags)
        mat = Matrix{Any}(undef, (np, 1))
        mat[:,1] = map(d->trunc(exp(d), digits=6), m.hpr)
        return UselessTable(mat; cnames=["prob."], rnames=m.tags)
end

function postable(m)
        tsv = readdlm(joinpath(pathof(KongYiji), "..", "..", "data", "postable.tsv"), '\t', String)
        return UselessTable(tsv[2:end,:]; cnames=tsv[1,:], heads=["CTB postable"])
end

function normalize!(h2v::Vector{Dict{Int,Tv}}, INF::Vector{Tv})
        for (ih, vs) in enumerate(h2v)
                tot = sum(values(vs)) + EPS * (length(vs) + 1)
                for (k, v) in vs
                        vs[k] = (v + EPS) / tot
                end
                INF[ih] = EPS / tot
        end
end

function normalize!(hpr::Vector{Tv})
        hpr .+= EPS
        hpr ./= sum(hpr)
end

function normalize!(h2h::Matrix{Tv})
        h2h .+= EPS
        h2h ./= sum(h2h; dims=2)
end

log!(hpr::Vector{Tv}) = hpr .= log.(hpr)
log!(h2h::Matrix{Tv}) = h2h .= log.(h2h)
function log!(h2v::Vector{Dict{Int,Tv}})
        for (ih, vs) in enumerate(h2v)
                for (k, v) in vs
                        vs[k] = log(v)
                end
        end
        h2v
end


text:
################################################################################
#
#       NfOrd/Ideal.jl : Ideals in orders of absolute number fields
#
# This file is part of Hecke.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann
#  Copyright (C) 2015, 2016, 2017 Claus Fieker
#
################################################################################

include("Ideal/Ideal.jl")
include("Ideal/TwoElement.jl")
include("Ideal/Arithmetic.jl")
include("Ideal/Prime.jl")
include("Ideal/FacElem.jl")
include("Ideal/Enum.jl")
include("Ideal/Relative.jl")


text:
"""
  Penalty for BR2
"""
struct InviscidPenalty{Tsol, Tres} <: AbstractDiffusionPenalty

  function InviscidPenalty{Tsol, Tres}(mesh, sbp, opts) where {Tsol, Tres}

    return new()
  end
end

function applyPenalty(penalty::InviscidPenalty{Tsol, Tres}, sbp, params::AbstractParamType,
                      sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix{Tsol}, theta::AbstractMatrix{Tres},
                      qL_el::AbstractMatrix, qR_el::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      alphas::AbstractVector,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L::AbstractMatrix, res1R::AbstractMatrix,
                      res2L::AbstractMatrix, res2R::AbstractMatrix
                     ) where {Tsol, Tres}

  fill!(res1L, 0); fill!(res1R, 0)
  fill!(res2L, 0); fill!(res2R, 0)

  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = length(jacL)
  dim = size(nrm_face, 1)

  for i=1:numNodesPerFace
    fac = calcLength(params, ro_sview(nrm_face, :, i))*sbpface.wface[i]
    for j=1:numDofPerNode
      res1L[j, i] =  fac*delta_w[j]
      res1R[j, i] = -fac*delta_w[j]
    end
  end

  return nothing
end

function applyPenalty_diff(penalty::InviscidPenalty{Tsol, Tres}, sbp,
                      params::AbstractParamType, sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix{Tsol},
                      delta_w_dotL::Abstract4DArray, delta_w_dotR::Abstract4DArray,
                      theta::AbstractMatrix{Tres},
                      theta_dotL::Abstract4DArray, theta_dotR::Abstract4DArray,
                      qL_el::AbstractMatrix, qR_el::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      alphas::AbstractVector,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L_dotL::Abstract4DArray, res1L_dotR::Abstract4DArray,
                      res1R_dotL::Abstract4DArray, res1R_dotR::Abstract4DArray,
                      res2L::AbstractMatrix, res2R::AbstractMatrix,
                      res2L_dotL::Abstract4DArray, res2L_dotR::Abstract4DArray,
                      res2R_dotL::Abstract4DArray, res2R_dotR::Abstract4DArray,
                     ) where {Tsol, Tres}

  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = length(jacL)
  dim = size(nrm_face, 1)

  for i=1:numNodesPerFace
    fac = calcLength(params, ro_sview(nrm_face, :, i))
    for j=1:numDofPerNode
      #res1L[j, i] =  fac*delta_w[j]
      #res1R[j, i] = -fac*delta_w[j]
    end
  end

  for q=1:numNodesPerElement
    for p=1:numNodesPerFace
      fac = calcLength(params, ro_sview(nrm_face, :, p))*sbpface.wface[p]
      for j=1:numDofPerNode
        for i=1:numDofPerNode
          res1L_dotL[i, j, p, q] =  fac*delta_w_dotL[i, j, p, q]
          res1L_dotR[i, j, p, q] =  fac*delta_w_dotR[i, j, p, q]
          res1R_dotL[i, j, p, q] = -fac*delta_w_dotL[i, j, p, q]
          res1R_dotL[i, j, p, q] = -fac*delta_w_dotR[i, j, p, q]
        end
      end
    end
  end

  fill!(res2L_dotL, 0); fill!(res2L_dotR, 0)
  fill!(res2R_dotL, 0); fill!(res2R_dotR, 0)

  return false
end


text:
# clustering julia

# define the accuracy function
function acc(pred::AbstractVector, truth::AbstractVector)
    return sum(pred .== truth) / length(pred);
end

# load the libraries
using GLMNet, RDatasets, MLBase, Plots, DecisionTree, Distances, NearestNeighbors, Random, LinearAlgebra, LIBSVM

# load iris data
iris = dataset("datasets", "iris");
describe(iris)

X = Matrix{Float64}(iris[!, 1:end-1]);
labels = iris[!, end];
labelsmap = labelmap(labels);
Y = labelencode(labelsmap, labels);

# using the first two classes to perform feature reduction
# similar to 04. dim reduction
# PCA
using ScikitLearn, MultivariateStats, StatsBase, Statistics
x = X[1:100, :];
y = Y[1:100];

# standardize the data and use PCA
data = (x .- mean(x, dims=1)) ./ std(x, dims=1);
p = fit(PCA, data', maxoutdim=2);
P = projection(p);
pca_reduction = P' * data';
scatter(pca_reduction[1, 1:50], pca_reduction[2, 1:50], lab="setosa");
scatter!(pca_reduction[1, 51:end], pca_reduction[2, 51:end], lab="versicolor")
@show norm(data' - reconstruct(p, pca_reduction))

# using tSNE
@sk_import manifold : TSNE;
tsnepy = TSNE(n_components=2 ,perplexity=20.0,early_exaggeration=50);
@time tsne_reduction = tsnepy.fit_transform(data);
scatter(tsne_reduction[1:50, 1], tsne_reduction[1:50, 2], lab="setosa");
scatter!(tsne_reduction[51:100, 1], tsne_reduction[51:100, 2], lab="versicolor")

# using UMAP
using UMAP
debye = pairwise(Euclidean(), data, data, dims=1);
@time umap_reduction = umap(debye.^2, 2);
scatter(umap_reduction[1, 1:50], umap_reduction[2, 1:50], lab="setosa");
scatter!(umap_reduction[1, 51:end], umap_reduction[2, 51:end], lab="cersicolor")

corr = cor(data')
@time umap_reduction = umap(10*corr, 2);  # note that corr.^2 does not work
scatter(umap_reduction[1, 1:50], umap_reduction[2, 1:50], lab="setosa");
scatter!(umap_reduction[1, 51:end], umap_reduction[2, 51:end], lab="cersicolor")

# back to clustering
X;
Y;
labeldecode(labelsmap, Y);

# perclass splitting
# test-train split
function perclass_split(y::AbstractVector, at::T) where T<:Real
    uids = unique(y);
    keepids = [];
    for ui in uids
        curids = findall(y .== ui);
        rowids = randsubseq(curids, at);
        push!(keepids, rowids...);
    end
    return keepids;
end

trainids = perclass_split(y, 0.7);
testids = setdiff(collect(1:length(y)), trainids);

function assign_class(pred::T) where T<:Real
    return argmin(abs.(collect(1:3) .- pred));
end

# lasso and elastic net
path = glmnet(X[trainids,:], y[trainids]);
cv = glmnetcv(X[trainids,:], y[trainids]); # choosing the best lammbda (reg const)
lam = path.lambda[argmin(cv.meanloss)];

path = glmnet(X[trainids, :], Y[trainids], lambda = [lam]);
pred_lasso = GLMNet.predict(path, X[testids, :]);
@show acc(assign_class.(pred_lasso[:,1]), y[testids]) # acc = 1.0

# Ridge
path = glmnet(X[trainids,:], y[trainids], alpha=0.0);
cv = glmnetcv(X[trainids,:], y[trainids], alpha=0.0); # choosing the best lammbda (reg const)
lam = path.lambda[argmin(cv.meanloss)];

path = glmnet(X[trainids, :], Y[trainids], lambda = [lam], alpha=0.0);
pred_ridge = GLMNet.predict(path, X[testids, :]);
@show acc(assign_class.(pred_ridge[:,1]), y[testids])

# Elastic Net
#combination of lasso and ridge
path = glmnet(X[trainids,:], y[trainids], alpha=0.5);
cv = glmnetcv(X[trainids,:], y[trainids], alpha=0.5); # choosing the best lammbda (reg const)
lam = path.lambda[argmin(cv.meanloss)];
path = glmnet(X[trainids, :], Y[trainids], lambda = [lam], alpha=0.5);
pred_ridge = GLMNet.predict(path, X[testids, :]);
@show acc(assign_class.(pred_ridge[:,1]), y[testids])

# Tree Based models
# decision tree
dtclassifier = DecisionTreeClassifier(max_depth=2);
@time DecisionTree.fit!(dtclassifier, X[trainids, :], Y[trainids]);
pred_dt = DecisionTree.predict(dtclassifier, X[testids, :]);
@show acc(pred_dt, Y[testids])

# random forest
rfclassifier = RandomForestClassifier(n_trees=20);
@time DecisionTree.fit!(rfclassifier, X[trainids, :], Y[trainids]);
pred_rf = DecisionTree.predict(rfclassifier, X[testids, :]);
@show acc(pred_rf, Y[testids]);

# xgboost
using XGBoost
@time xgbclassifier = xgboost(X[trainids, :], 100, eta=0.2, max_depth=2, silent=false, label=Y[trainids]);
pred_xgb = XGBoost.predict(xgbclassifier, X[testids, :]);
@show acc(assign_class.(pred_xgb), Y[testids])

# SVM
svm = svmtrain(X[trainids, :]', Y[trainids]);
pred_svm, __ = svmpredict(svm, X[testids, :]');
@show acc(pred_svm, Y[testids])

max_stars_repo_path: 06. Classification.jl

max_issues_repo_path: 06. Classification.jl

max_forks_repo_path: 06. Classification.jl


text:
@testset "Test Maxwell solvers" begin
    include(joinpath("MaxwellSolver", "test_poisson_fourier_1d.jl"))
    include(joinpath("MaxwellSolver", "test_fdtd_1d.jl"))
end


text:
# Computes the LDLinv, along with the adjaceny matrix, both need to be passed
# around to create a solver for a parallelized thread
function computeLDLinv(sddm)
    a,d = adj(sddm)
    a1 = extendMatrix(a,d)
    ldli = approxchol_lap_pc(a1)
    la = lap(a1)
    return ldli, la
end

# Wrap the LDLinv object as a preallocated linear operator
function approxcholOperator(ldli::LDLinv{Tind,Tval},buff::Vector{Tval}) where {Tind,Tval}
    prod = @closure rhs -> LDLsolver!(buff,ldli,rhs)
    return PreallocatedLinearOperator{Tval}(length(ldli.d), length(ldli.d), true, true, prod, nothing, nothing)
end

# Compute a solver for a grounded system (SDDM matrix) with a PreallocatedLinearOperator, and an adjacency matrix
function approxcholSolver(P::PreallocatedLinearOperator, la::AbstractArray; tol::Real=1e-6, maxits=300, verbose=false)

    tol_=tol
    maxits_=maxits
    verbose_=verbose

    f = function(b;tol=tol_, maxits=maxits_, verbose=verbose_)
        xaug = Krylov.cg(la,[b; -sum(b)] .- mean([b; -sum(b)]), M=P, rtol = tol, itmax=maxits, verbose=verbose)[1]
        xaug .= xaug .- xaug[end]
        return xaug[1:end-1]
    end

    return f

end

# Compute a solver for a grounded system (SDDM matrix) with a LDLinv object, and an adjacency matrix
function approxcholSolver(ldli::LDLinv, la::AbstractArray; tol::Real=1e-6, maxits=300, verbose=false)

    buff = zeros(length(ldli.d))
    P = approxcholOperator(ldli,buff)

    tol_=tol
    maxits_=maxits
    verbose_=verbose

    f = function(b;tol=tol_, maxits=maxits_, verbose=verbose_)
        xaug = Krylov.cg(la,[b; -sum(b)] .- mean([b; -sum(b)]), M=P, rtol = tol, itmax=maxits, verbose=verbose)[1]
        xaug .= xaug .- xaug[end]
        return xaug[1:end-1]
    end

    return f

end

# Compute a solver for a grounded system (SDDM matrix)
function approxcholSolver(sddm::AbstractArray; tol::Real=1e-6, maxits=300, verbose=false)

    # Compute LDLinv object and adjacency matrix
    ldli,la = computeLDLinv(sddm)

    # Wrap the LDLinv object as a preallocated linear operator
    buff = zeros(length(ldli.d))
    P = approxcholOperator(ldli,buff)

    tol_=tol
    maxits_=maxits
    verbose_=verbose

    f = function(b;tol=tol_, maxits=maxits_, verbose=verbose_)
        xaug = Krylov.cg(la,[b; -sum(b)] .- mean([b; -sum(b)]), M=P, rtol = tol, itmax=maxits, verbose=verbose)[1]
        xaug .= xaug .- xaug[end]
        return xaug[1:end-1]
    end

    return f

end


text:
# This code generates the runtests.jl file automatically for some lucky exercises that work with this script.

# NOTE: This is a very experimental version of this idea. The code is not all that well written,
# and may not work on most exercises. If it works, it might save you hours of time writing that exercise.
# If it doesn't it will throw an error and quit.

using JSON, HTTP

# Define 404 Not found Exception so test_JSONtoTests.jl can handle non existent slugs.
struct NotFoundException <: Exception
    var::String
end

"""
    array_to_str(array::AbstractArray)

Convert an array to a array in string form.

## Examples

Array with different types of elements:
```julia-repl
julia> array_to_str([1, 2, "hey"])
"[1, 2, \"hey\"]"
```

Empty array:
```julia-repl
julia> array_to_str([])
"[]"
```
"""
function array_to_str(array::AbstractArray)
    # Start array with opening bracket
    str = ["["]
    # If array is empty, return an empty array
    if length(array) == 0
        return "[]"
    end

    # If list contains items, for each item in array, add it to array
    for item in array
        # If item is a string, enclose it in quotations and add it to our array
        if typeof(item) == String
            push!(str, "\"$item\", ")
        # If item is an array, call function recursively to add it to our array
        elseif typeof(item) == AbstractArray
            push!(str, array_to_str(item))
        # Else, if item is something else, just write it to array
        else
            push!(str, "$item, ")
        end
    end

    # Make it a string and remove the last extra comma and space from array
    truncated_str = join(str)[1:end-2]

    # Return final array string by closing with bracket
    return string(truncated_str, "]")
end

function write_test_code(testcase_array)
    # Initialize array to add code so far to
    lines_written_so_far = []

    for testcases in testcase_array

        # Get the function name that we want to test
        function_to_test = get(testcases, "property", nothing)

        # Write the testcase name, and description
        push!(lines_written_so_far, "@testset \"" * string(get(testcases, "description", nothing)) * "\" begin\n")

        # Get the parameters that need to be input into the function for testing
        input_params = []

        for parameter in keys(get(testcases, "input", nothing))

            # If the parameters are arrays, format them as arrays in String using our function array_to_str()
            # e.g. [a, b, c, d]
            if isa(get(get(testcases, "input", nothing), parameter, nothing), AbstractArray)
                # Start list string form
                parameter_temp = array_to_str(get(get(testcases, "input", nothing), parameter, nothing))
                # Add the parameter with its variable name
                push!(input_params, "$parameter = " *  parameter_temp * ", ")
            else
                # If parameters are not lists, get the parameters seperated by commas.

                # If parameters are strings, enclose them in quotations
                if isa(get(get(testcases, "input", nothing), parameter, nothing), String)
                    push!(input_params, "$parameter = " * "\"" * string(get(get(testcases, "input", nothing), parameter, nothing)) * "\"" * ", ")

                # Else, just put them in as they were
                else
                    push!(input_params, "$parameter = " * string(get(get(testcases, "input", nothing), parameter, nothing)) * ", ")
                end
            end
        end

        # Finalise parameters to input into function
        input_temp = string(join(input_params))[1:lastindex(join(input_params))-2]

        # Create expected output to compare with
        expected_output = []

        # If expected output is an array, write it as an array like before
        if isa(get(testcases, "expected", nothing), AbstractArray)
            expected_output = array_to_str(get(testcases, "expected", nothing))
        # If expected_output is a Dict, write it as a Dict
        elseif isa(get(testcases, "expected", nothing), Dict)
            expected_output = string(get(testcases, "expected", nothing))
        # If expected output is a Bool, write it without quotes
        elseif isa(get(testcases, "expected", nothing), Bool)
            expected_output = string(get(testcases, "expected", nothing))
        else
        # If expected output is not any of the handled cases, best choice is to print it out in quotes.
            expected_output = "\"" * string(get(testcases, "expected", nothing)) * "\""
        end

        # Write the actual code that tests our function.
        push!(lines_written_so_far, "    @test $function_to_test(" * input_temp * ") == " * strip(expected_output))
        push!(lines_written_so_far, "\nend\n\n")
    end
    return lines_written_so_far
end

function main(exercise_slug)
    # Get the exercise the user is working on from command line argument
    current_exercise = exercise_slug

    # Download the test data for the exercise the user is working on and parse the JSON.
    try
        get_json = HTTP.get("https://raw.githubusercontent.com/exercism/problem-specifications/master/exercises/$current_exercise/canonical-data.json")
        json_code = String(get_json.body)
        data = JSON.parse(json_code)
    catch error
        if isa(error, HTTP.ExceptionRequest.StatusError)
            throw(NotFoundException("Sorry, there was a problem processing the data.\nMaybe the exercise is not in the exercism problem-specifications list?"))
        end

        throw("Something went wrong.")
    end

    # If exercise exists in problem-specifications, go ahead with the conversion.
    get_json = HTTP.get("https://raw.githubusercontent.com/exercism/problem-specifications/master/exercises/$current_exercise/canonical-data.json")
    json_code = String(get_json.body)
    data = JSON.parse(json_code)

    # Create an array to push our julia code text to.
    lines = String[]

    # Write the initial two lines for runtets.jl
    push!(lines, "# canonical data version: " * string(get(data, "version", nothing)) * "\n\nusing Test\n\n")
    push!(lines, "include(\"$current_exercise.jl\")\n\n")

    # Some exercises have a slightly different json structure: the testcases are stored in data["cases"][i]["cases"].
    # I don't know why this structure is the way it is, but this generator needs to handle it.
    # If the json has the weird double case nested structure, handle that
    if (get(data["cases"][1], "cases", nothing) !== nothing)
        for i in data["cases"]
            testcase_array = i["cases"]
            # For every testcase in the nested case, write the test code
            push!(lines, join(write_test_code(testcase_array)))
        end
    else
        # If the structure is not the nested on we discussed, handle cases as usual
        testcase_array = get(data, "cases", nothing)
        # For every testcase in the canonical data json file, write the test code
        push!(lines, join(write_test_code(testcase_array)))
    end

    # Finalise the runtests code, and return it.
    final_code = join(lines)
    return final_code
end

# If this code is run on its own, show the code
if length(ARGS) > 0
    test_code = main(ARGS[1])
    println(test_code)
else
    # Else, if this code is run by the coverage tester, do nothing. The coverage tester will detect errors on its own.
    # println("Please mention the slug of the exercise as an argument to this script")
end

# For now, print code. Later when code works on 100% of exercises, we will save the code to its folder.

#=
Write final code to folder
RUN ONLY IF CODE WORKS ON EXERCISE

open("../exercises/$(ARGS[1])/runtests.jl", "w") do io
    write(io, final_code)
end;
=#

max_stars_repo_name: soumitradev/julia

max_issues_repo_name: soumitradev/julia

max_forks_repo_name: soumitradev/julia


text:
module TestGroupGenerics

using Test
using RLEVectors

@testset begin

# compare group
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test (x .< 3) == RLEVector([true,false],[4,6])

# math group
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test x .+ 5 == 5 .+ x
@test x .+ 4 == RLEVector([5,5,6,6,8,8])
@test div.(vec,2) == collect(div.(x,2))
@test median(RLEVector([1,2,3,2,1,5,4])) == median([1,2,3,2,1,5,4])
@test median(RLEVector([1,2,3,2,1,5])) == median([1,2,3,2,1,5])
@test median(RLEVector([3,2,1])) == 2.0
@test median(RLEVector([2, 2, 2, 3, 3, 3])) == 2.5 # median is average of end of run and next value
@test sum(RLEVector([4, 4, 5, 5, 6, 6])) == 30
@test mean(RLEVector([4, 4, 5, 5, 6, 6])) == 5.0
@test x .^ 2 == RLEVector( [1, 1, 4, 4, 16, 16] )
@test x .^ 3 == RLEVector( [1, 1, 8, 8, 64, 64] )
@test x + x == RLEVector( [2, 2, 4, 4, 8, 8] )
@test x - x == RLEVector( [0, 0, 0, 0, 0, 0] )
rle = RLEVector( [4, 4, 9, 9, 16, 16] )
@test sqrt.(rle) == RLEVector( [2, 2, 3, 3, 4, 4] )

# math on bools
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test x .+ true == x .+ 1
@test x .+ false == x
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test x .+ true == x .+ 1
@test x .+ false == x

# broadcast with different sizes
x = RLEVector([1,1,3])
@test x .+ 4 == RLEVector([5,5,7])
@test x .+ [4 5 6] == [[5 6 7] ; [5 6 7] ; [7 8 9]]
@test x .+ [4, 5, 6] == RLEVector([5, 6, 9])
@test x .+ [[2 2 2] ; [4 4 4]]' == [[3 3 5] ; [ 5 5 7]]'

# findmax, findmin
@test findmin(RLEVector([1,2,3,4,1,1])) == findmin([1,2,3,4,1,1])
@test findmax(RLEVector([1,2,3,4,1,1])) == findmax([1,2,3,4,1,1])

# indexin
foo = IntegerRle( collect(1:1000), collect(5:5:5000))
x = RLEVector([2,2,4,4,3,3])
y = RLEVector([0,0,0,3,3,3,4,4])
@test indexin(x,y) == indexin(collect(x), collect(y))
@test indexin(x,collect(3:11)) == indexin(collect(x),collect(3:11))
@test indexin([200,200,1,1,5,5],foo) == indexin([200,200,1,1,5,5],foo)

@test findall(in(RLEVector(collect(3:10))), RLEVector([1,1,2,2,3,3])) == collect(5:6)
@test findall(in(3:10), RLEVector([1,1,2,2,3,3])) == collect(5:6)
@test findall(in(collect(3:10)), RLEVector([1,1,2,2,3,3])) == collect(5:6)
@test findall(in(RLEVector(collect(1:4))), [3,4,5]) == [1,2]

# in
@test in(3, RLEVector( [ 1,2,2,3 ] )) == true
@test in(4, RLEVector( [ 1,2,2,3 ] )) == false
@test setdiff( Set([1,2,3,4,5]), RLEVector([1,1,2,2,4,4,5,5]) ) == Set([3])

end # testset

end # module


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

include("model_CloudError.jl")
include("model_ErrorAdditionalInfo.jl")
include("model_ErrorResponse.jl")
include("model_Identity.jl")
include("model_ParameterValues.jl")
include("model_ParameterValuesValue.jl")
include("model_PolicyAssignment.jl")
include("model_PolicyAssignmentListResult.jl")
include("model_PolicyAssignmentProperties.jl")
include("model_PolicySku.jl")


text:
module ClimatePlots

using Dates
using Reexport
@reexport using ClimateBase
using PyCall
using PyPlot
using Statistics

const basemap = PyNULL()
const mpl = PyNULL()
const cmocean = PyNULL()
# const scipy = PyNULL()

function __init__()
    copy!(mpl, pyimport_conda("matplotlib", "matplotlib", "conda-forge"))
    copy!(basemap, pyimport_conda("mpl_toolkits.basemap", "basemap", "conda-forge"))
    copy!(cmocean, pyimport_conda("cmocean", "cmocean", "conda-forge"))
  # copy!(scipy, pyimport_conda("scipy.interpolate", "scipy"))
end

include("mapping.jl")

export plot, mapclimgrid, hist

end # module


text:
export samplingFreq, roundFreq, numSamplesPerPeriod, phaseShift

function samplingFreq(rp::RedPitaya,dec::Integer)
  baseFreq = 125e6
  freq = div(baseFreq,dec)
  return freq
end

function roundFreq(rp::RedPitaya,dec::Integer,freq)
  numPeriods = round(Int64,samplingFreq(rp,dec) / freq)
  return div(samplingFreq(rp,dec), numPeriods)
end

function numSamplesPerPeriod(rp::RedPitaya,dec::Integer,freqR)
  return Int(div(samplingFreq(rp,dec),freqR))
end

"""
Return the number of time points the signal u has the be shifted
in order to have phase pi/2 (i.e. be a sine wave)
"""
function phaseShift(u, numPeriods)
  numSampPerPeriod = div(length(u),numPeriods)
  complexSine = exp(-2*pi*im*numPeriods*(0:length(u)-1)/length(u))
  fourierCoeff = dot(complexSine,u)
  return round(Int,numSampPerPeriod*(angle(fourierCoeff) -pi/2) / (2pi))
end

function test2(rp)
  a=zeros(30)
  for l=1:30
    wpos = parse(Int32,query(rp,"ACQ:WPOS?"))
    u = receiveAnalogSignalLowLevel(rp,1,0,16384,"STA",true,true)
    a[l] = wpos
  end
  return a
end


text:
const GOOGLE_URL = "http://finance.google.com/finance/historical?"  # for querying google finance's servers

"""
Download stock price data from Google Finance into a TS object.

`google(symb::String; from::String="2000-01-01", thru::String=string(Dates.today()))::TS`

# Arguments
- `symb` ticker symbol of the stock
- `from` starting date of the historical data request (string formatted as yyyy-mm-dd)
- `thru` ending date of the historical data request (string formatted as yyyy-mm-dd)

# Example

```
julia> google("IBM", from="2010-06-09", thru=string(Dates.today()))
1756x5 Temporal.TS{Float64,Date}: 2010-06-09 to 2017-05-30
Index       Open    High    Low     Close   Volume
2010-06-09  124.83  125.84  123.58  123.9   7.800309e6
2010-06-10  125.99  128.22  125.8   127.68  7.47961e6
2010-06-11  126.73  128.8   126.44  128.45  5.827093e6
2010-06-14  128.5   129.97  128.49  128.5   6.753113e6
2010-06-15  128.93  129.95  128.37  129.79  6.652612e6
⋮
2017-05-23  152.57  153.68  151.92  152.03  2.564503e6
2017-05-24  152.21  152.76  151.23  152.51  3.732399e6
2017-05-25  153.25  153.73  152.95  153.2   2.582815e6
2017-05-26  152.85  153.0   152.06  152.49  2.443507e6
2017-05-30  151.95  152.67  151.59  151.73  3.666032e6
```
"""
function google(symb::String;
                from::String="2000-01-01",
                thru::String=string(Dates.today()))::TS
    from_date = parse(Date, from, Dates.DateFormat("yyyy-mm-dd"))
    thru_date = parse(Date, thru, Dates.DateFormat("yyyy-mm-dd"))
    url = string("$(GOOGLE_URL)q=$(symb)",
                 "&startdate=$(Dates.monthabbr(Dates.month(from_date)))",
                 "+$(@sprintf("%.2d",Dates.dayofmonth(from_date)))",
                 "+$(Dates.year(from_date))",
                 "&enddate=$(Dates.monthabbr(thru_date))",
                 "+$(@sprintf("%.2d",Dates.dayofmonth(thru_date)))",
                 "+$(Dates.year(thru_date))&output=csv")
    response = HTTP.get(url)
    indata = Temporal.csvresp(response)
    return TS(indata[1], indata[2], indata[3][2:end])
end


text:
using ServerMetrics
using Dates

statsd_exporter = ServerMetrics.StatsdExporter(
    statsd_backend=ServerMetrics.UDPBackend("127.0.0.1", 8125),
    send_interval=Dates.Second(1)
)

Base.@kwdef struct UptimeMetrics <: AbstractMetricCollection
    server_uptime_seconds = Gauge()
    server_heartbeats_total = Counter()
end

@info "Starting statsd exporter"
ServerMetrics.start_statsd_exporter!(statsd_exporter)

metrics = UptimeMetrics()
ServerMetrics.publish_metrics_from(metrics; overwrite=true)
startup_timestamp = Dates.datetime2unix(Dates.now())
ServerMetrics.@spawn_sticky_periodic_task "UptimeTracker" Dates.Second(1) begin
    time_delta = Dates.datetime2unix(Dates.now()) - startup_timestamp
    set!(metrics.server_uptime_seconds, time_delta)
    inc!(metrics.server_heartbeats_total)
end

sleep(10)
@info "Stopping statsd exporter"
ServerMetrics.stop_statsd_exporter!(statsd_exporter)


text:
using OrderedCollections

appendSymbol(path::Nothing, name::Symbol) = name
appendSymbol(path         , name::Symbol) = :( $path.$name )

"""
    buildModia3D(model)

Generate the following type of code:

```
model | Model(_id = rand(Int),
              start = Map(_qdd = zeros(1)),
              mbs_equations = :[
                 zeros(1) = multibodyResiduals(_id, _leq_mode, instantiatedModel, _qdd, rev.variables)
                 rev.qdd = _qdd[1]])
```
"""
function buildModia3D(model)    #(model; path="")
    jointInfo = []
    getJointInfo!(model, jointInfo, nothing)

    if length(jointInfo) == 0
        ndofTotal = 0    # what to do???
        args = nothing
        error("No joints present ->  so no states")
    else
        ndofTotal = sum(j.ndof for j in jointInfo)
        args = [appendSymbol(e.path, :variables) for e in jointInfo]
    end

    code = Expr[]
    i = 1
    for e in jointInfo
        if e.ndof == 1
            push!(code, :( $(e.path).qdd = _qdd[$i] ))
        else
            iend = i + e.ndof - 1
            push!(code, :( $(e.path).qdd = _qdd[$i:$iend] ))
        end
        i += e.ndof
    end

    mbsCode = Model(_id = rand(Int),
                    _qdd = Var(start = zeros(ndofTotal)),
                    mbs_equations = :[0 = multibodyResiduals(_id, _leq_mode, instantiatedModel, time, _qdd, $(args...))
                                      $(code...)])

                    #mbs_equations = :[zeros($ndofTotal) = multibodyResiduals(_id, _leq_mode, instantiatedModel, time, _qdd, $(args...))
                    #                  $(code...)])

    return mbsCode | model
end

"""
    getJointInfo!(model, jointInfo, path)

Recursively traverse the OrderedDict `model` and return a vector of tuples.
A tuple contains the path to a joint (where the path is defined as
quoted expression, for example :(a.b.c)), and the length of the `qdd` vector.
"""
function getJointInfo!(model, jointInfo, path)::Nothing
    if haskey(model, :_constructor)
        constructor = model[:_constructor]
        if typeof(constructor) <: OrderedDict && haskey(constructor, :ndof)
            ndof = constructor[:ndof]
            push!(jointInfo, (path=path, ndof=ndof))
            return
        end
    end

    for (key,value) in model # zip(keys(model), model)
        if typeof(value) <: OrderedDict
            getJointInfo!(value, jointInfo, appendSymbol(path, key))
        end
    end
    return nothing
end


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #
function _lsqfit_local_regression_model(X::Array{Float64,1},Y::Array{Float64,1})::PSResult

    # setup the model -
    @. model(x, p) = p[1]+p[2]*x+p[3]*x^2

    # setup the fit -
    p0 = [-1.0, 3.0, -2.0]

    # run the fit -
    fit_bounds = curve_fit(model, X, Y, p0)

    # Wrap -
    a = fit_bounds.param

    local_model = LocalExpectationRegressionModel(a[1], a[2], a[3], 0.0)

    # grab the result -
    return PSResult(local_model)
end

function _evaluate_local_regression_model(model::LocalExpectationRegressionModel,X::Array{Float64,1})::PSResult

    # initialize -
    f_value_array = Array{Float64,1}()

    # get the parameters for the model -
    a0 = model.a0
    a1 = model.a1
    a2 = model.a2
    a3 = model.a3

    # compute -
    for value in X
        term_1 = a0
        term_2 = a1*value
        term_3 = a2*(value)^2
        term_4 = a3*(value)^3
        f_value = term_1+term_2+term_3+term_4
        push!(f_value_array,f_value)
    end

    # return -
    return PSResult{Array{Float64,1}}(f_value_array)
end

function _calculate_intrinsic_value_trade_legs(contractSet::Set{PSAbstractAsset},underlyingPrice::Float64)::PSResult

    # initialize -
    tmp_array = Float64[]

    # main loop -
    for contract in contractSet
        result = intrinsic_value(contract,underlyingPrice)
        if (isa(result.value,Exception) == true)
            return result
        end
        iv_value = result.value.intrinsic_value
        push!(tmp_array,iv_value)
    end

    # total intrinsic value -
    total_intrinsic_value = sum(tmp_array)

    # return -
    return PSResult{Float64}(total_intrinsic_value)
end

function _calculate_options_cost_table(contractSet::Set{PSAbstractAsset}, underlying_price_table::Array{Float64,2};
    riskFreeRate::Float64 = 0.015, timeMultiplier=1.0, earlyExercise::Bool = false)::PSResult

    # now that we have the underlying price table, let's operate on this table to estimate values -
    # what is the size of the underlying price table?
    (number_of_paths,number_of_time_steps) = size(underlying_price_table)

    # initialize an empty option_cost_table -
    option_excercise_reward_table = zeros(number_of_paths,number_of_time_steps)

    # initialize an empty profit_loss_table -
    intrinsic_value_table = zeros(number_of_paths,number_of_time_steps)

    # ok, so before we do anything about excerising the option contract, lets complete the P/L table -
    for path_index = 1:number_of_paths
        for time_index = 1:number_of_time_steps

            # grab the potential price -
            underlying_price_value = underlying_price_table[path_index, time_index]

            # process each leg of the trade -
            result = _calculate_intrinsic_value_trade_legs(contractSet, underlying_price_value)
            if (isa(result.value,Exception) == true)
                return result
            end
            total_intrinsic_value = result.value

            # capture this value -
            intrinsic_value_table[path_index,time_index] = total_intrinsic_value
        end
    end

    # ok depending upon whether this is an American (USA,USA!) or an European option
    # we do different things = earlyExercise = false is a European option
    if (earlyExercise == false)

        # what is the deltaT?
        ΔT = (number_of_time_steps*timeMultiplier)*(1.0/364.5)
        d = exp(-1*riskFreeRate*ΔT)

        # compute the value -
        for path_index = 1:number_of_paths
            value = d*(intrinsic_value_table[path_index,end])
            option_excercise_reward_table[path_index,end] = value
        end
    else

        # initialize the last col of the cost table -
        for path_index = 1:number_of_paths
            value = 1.0*(intrinsic_value_table[path_index,end])
            option_excercise_reward_table[path_index,end] = value
        end

        # so the option cost table
        backward_index_collection = collect(range(number_of_time_steps,step=-1,stop=2))
        for time_index in backward_index_collection

             # What is my d?
            ΔT = (1.0*timeMultiplier)
            d = exp(-1*riskFreeRate*ΔT)

            # get all the Y's for the regression -
            Y = intrinsic_value_table[:,time_index]

            # get all of the X's (we'll filter out OTM rows later) -
            X = underlying_price_table[:,time_index-1]

            # which of the X's are ITM?
            itm_index_array = findall(x->x>0, intrinsic_value_table[:,time_index-1])
            Xdata = X[itm_index_array]
            Ydata = Y[itm_index_array].*d

            # compute the local model -
            result = _lsqfit_local_regression_model(Xdata,Ydata)
            if (isa(result.value,Exception) == true)
                return result
            end
            local_model = result.value

            # which paths will have early an early excercise event?
            # lets compare what we would get if we excercised now, versus waiting -
            result = _evaluate_local_regression_model(local_model,Xdata)
            if (isa(result.value,Exception) == true)
                return result
            end
            Ycontinuation = result.value

            # ok, so let's compare the excercise value versus the Ycontinuation -
            for (index,itm_index) in enumerate(itm_index_array)

                excercise_value = intrinsic_value_table[itm_index,time_index-1]
                continuation_value = Ycontinuation[index]

                if (excercise_value>=continuation_value)
                    option_excercise_reward_table[itm_index,time_index-1] = excercise_value
                    #option_excercise_reward_table[itm_index,time_index] = 0.0

                    # we excercised - so all future times are equal to zero
                    for local_time_index = time_index:number_of_time_steps
                        option_excercise_reward_table[itm_index,local_time_index] = 0.0
                    end
                else
                    option_excercise_reward_table[itm_index,time_index-1] = 0.0
                end
            end
        end
    end

    # compute the option price -
    (number_of_paths,number_of_time_steps) = size(option_excercise_reward_table)
    final_option_cost_array = Array{Float64,1}()
    d = ones(number_of_time_steps)
    for time_index = 1:number_of_time_steps
        ΔT = (time_index)
        d[time_index] = exp(-riskFreeRate*ΔT)
    end

    # compute the discounted payout -
    discounted_payout_array = option_excercise_reward_table*d

    # setup names tuple -
    result_tuple = (option_payout_table=option_excercise_reward_table,discount=d)

    # return - what is going on?
    return PSResult(result_tuple)
end
# ----------------------------------------------------------------------------------------------------------- #

# --- PUBLIC METHODS ---------------------------------------------------------------------------------------- #
# function longstaff_option_contract_price(contractSet::Set{PSAbstractAsset}, model::PSGeometricBrownianMotionModelParameters, initialCondition::Float64, tspan::Tuple{Float64,Float64}, timeStep::Float64;
#     number_of_trials::Int64=10000, return_time_step::Float64 = 1.0, earlyExercise::Bool = false, riskFreeRate::Float64 = 0.015, timeMultiplier=1.0)::PSResult

#     # TODO: checks - all ok with args?

#     # simulate the underlying asset price dynamics -
#     result = evaluate(model,initialCondition,tspan,timeStep;number_of_trials=number_of_trials,return_time_step=return_time_step)
#     if (isa(result.value,Exception) == true)
#         return result
#     end

#     # we have data: this call returns a named tuple w/T,X,μ,σ
#     X = result.value.X

#     # setup the price table -
#     underlying_price_table = transpose(X)

#     # call the valuation method -
#     return longstaff_option_contract_price(contractSet, underlying_price_table;
#         earlyExercise=earlyExercise, riskFreeRate=riskFreeRate, timeMultiplier=timeMultiplier)
# end

"""
    longstaff_option_contract_price
"""
function longstaff_option_contract_price(contractSet::Set{PSAbstractAsset}, underlyingPriceTable::Array{Float64,2};
    earlyExercise::Bool = false, riskFreeRate::Float64 = 0.015, timeMultiplier=1.0)::PSResult

    # update the option cost table -
    result = _calculate_options_cost_table(contractSet, underlyingPriceTable;
        earlyExercise = earlyExercise, riskFreeRate=riskFreeRate,timeMultiplier=timeMultiplier)
    if (isa(result.value,Exception) == true)
        return result
    end
    option_price_table = result.value

    # return -
    return PSResult(option_price_table)
end
# ----------------------------------------------------------------------------------------------------------- #