text:
"""
    BostonHousing(; as_df = true, dir = nothing)

The classical Boston Housing tabular dataset.

Sources:
   (a) Origin:  This dataset was taken from the StatLib library which is
                maintained at Carnegie Mellon University.
   (b) Creator:  Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the
                 demand for clean air', J. Environ. Economics & Management,
                 vol.5, 81-102, 1978.
   (c) Date: July 7, 1993

Number of Instances: 506

Number of Attributes: 13 continuous attributes (including target
                            attribute "MEDV"), 1 binary-valued attribute.

# Arguments

$ARGUMENTS_SUPERVISED_TABLE

# Fields

$FIELDS_SUPERVISED_TABLE

# Methods

$METHODS_SUPERVISED_TABLE

# Examples

```julia-repl
julia> using MLDatasets: BostonHousing

julia> dataset = BostonHousing()
BostonHousing:
  metadata => Dict{String, Any} with 5 entries
  features => 506×13 DataFrame
  targets => 506×1 DataFrame
  dataframe => 506×14 DataFrame

julia> dataset[1:5][1]
5×13 DataFrame
 Row │ CRIM     ZN       INDUS    CHAS   NOX      RM       AGE      DIS      RAD    TAX    PTRATIO  B        LSTAT
     │ Float64  Float64  Float64  Int64  Float64  Float64  Float64  Float64  Int64  Int64  Float64  Float64  Float64
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 0.00632     18.0     2.31      0    0.538    6.575     65.2   4.09        1    296     15.3   396.9      4.98
   2 │ 0.02731      0.0     7.07      0    0.469    6.421     78.9   4.9671      2    242     17.8   396.9      9.14
   3 │ 0.02729      0.0     7.07      0    0.469    7.185     61.1   4.9671      2    242     17.8   392.83     4.03
   4 │ 0.03237      0.0     2.18      0    0.458    6.998     45.8   6.0622      3    222     18.7   394.63     2.94
   5 │ 0.06905      0.0     2.18      0    0.458    7.147     54.2   6.0622      3    222     18.7   396.9      5.33

julia> dataset[1:5][2]
5×1 DataFrame
Row │ MEDV
    │ Float64
────┼─────────
  1 │    24.0
  2 │    21.6
  3 │    34.7
  4 │    33.4
  5 │    36.2

julia> X, y = BostonHousing(as_df=false)[]
([0.00632 0.02731 … 0.10959 0.04741; 18.0 0.0 … 0.0 0.0; … ; 396.9 396.9 … 393.45 396.9; 4.98 9.14 … 6.48 7.88], [24.0 21.6 … 22.0 11.9])
```
"""
struct BostonHousing <: SupervisedDataset
    metadata::Dict{String, Any}
    features
    targets
    dataframe
end

function BostonHousing(; as_df = true, dir = nothing)
    @assert dir === nothing "custom `dir` is not supported at the moment."
    path = joinpath(@__DIR__, "..", "..", "..", "data", "boston_housing.csv")
    df = read_csv(path)
    features = df[!, Not(:MEDV)]
    targets = df[!, [:MEDV]]

    metadata = Dict{String, Any}()
    metadata["path"] = path
    metadata["feature_names"] = names(features)
    metadata["target_names"] = names(targets)
    metadata["n_observations"] = size(targets, 1)
    metadata["description"] = BOSTONHOUSING_DESCR

    if !as_df
        features = df_to_matrix(features)
        targets = df_to_matrix(targets)
        df = nothing
    end

    return BostonHousing(metadata, features, targets, df)
end

const BOSTONHOUSING_DESCR = """
The Boston Housing Dataset.

Sources:
   (a) Origin:  This dataset was taken from the StatLib library which is
                maintained at Carnegie Mellon University.
   (b) Creator:  Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the
                 demand for clean air', J. Environ. Economics & Management,
                 vol.5, 81-102, 1978.
   (c) Date: July 7, 1993

Number of Instances: 506

Number of Attributes: 13 continuous attributes (including target
                            attribute "MEDV"), 1 binary-valued attribute.

Features:

    1. CRIM      per capita crime rate by town
    2. ZN        proportion of residential land zoned for lots over 25,000 sq.ft.
    3. INDUS     proportion of non-retail business acres per town
    4. CHAS      Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
    5. NOX       nitric oxides concentration (parts per 10 million)
    6. RM        average number of rooms per dwelling
    7. AGE       proportion of owner-occupied units built prior to 1940
    8. DIS       weighted distances to five Boston employment centres
    9. RAD       index of accessibility to radial highways
    10. TAX      full-value property-tax rate per 10,000 dollars
    11. PTRATIO  pupil-teacher ratio by town
    12. B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
    13. LSTAT    % lower status of the population

Target:

    14. MEDV     Median value of owner-occupied homes in 1000's of dollars

Note: Variable #14 seems to be censored at 50.00 (corresponding to a median price of \\\$50,000);
Censoring is suggested by the fact that the highest median price of exactly \\\$50,000 is reported in 16 cases,
while 15 cases have prices between \\\$40,000 and \\\$50,000, with prices rounded to the nearest hundred.
Harrison and Rubinfeld do not mention any censoring.

The data file stored in this repo is a copy of the UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/
"""

# Deprecated in v0.6,
function Base.getproperty(::Type{BostonHousing}, s::Symbol)
    if s == :features
        @warn "BostonHousing.features() is deprecated, use `BostonHousing().features` instead."
        return () -> BostonHousing(as_df=false).features
    elseif s == :targets
        @warn "BostonHousing.targets() is deprecated, use `BostonHousing().targets` instead."
        return () -> BostonHousing(as_df=false).targets
    elseif s == :feature_names
        @warn "BostonHousing.feature_names() is deprecated, use `BostonHousing().feature_names` instead."
        return () -> lowercase.(BostonHousing().metadata["feature_names"])
    else
        return getfield(BostonHousing, s)
    end
end

ext: jl


text:
using ViewerGL, LinearAlgebraicRepresentation
GL = ViewerGL
Lar = LinearAlgebraicRepresentation

V = [0.3307458 0.2735878 0.2527472 0.1955892 0.8735746 0.9604089 1.3600696 1.4469038 0.8952041 0.4087921 0.7697594 0.2833474 0.4465003 0.3177087 0.1293396 0.000548 1.1021824 0.6979647 1.0492281 0.6450104; 1.80716725 1.73722735 1.87091125 1.80097145 1.10328515 1.58712345 1.01597415 1.49981235 1.15705855 1.04345955 1.69419185 1.58059285 2.30620035 1.78225795 2.38416245 1.86022005 1.24012995 1.20346525 1.82393505 1.78727025]
EV = Array{Int64,1}[[1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8], [9, 10], [11, 12], [9, 11], [10, 12], [13, 14], [15, 16], [13, 15], [14, 16], [17, 18], [19, 20], [17, 19], [18, 20]]
VV = [[k] for k=1:size(V,2)]

GL.VIEW( GL.numbering(.05)((V,[VV, EV]),GL.COLORS[1]) );

# arrangement of input data
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV)
W, copEV, copFE = Lar.planar_arrangement(W, cop_EV)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# compute euler characteristic
χ = Lar.euler_characteristic(W, copEV, copFE)
println("χ = $χ ; bicon_comps = $(length(bicon_comps))");

# triangulation
triangulated_faces = Lar.triangulate2D(W, [copEV, copFE]);
FVs = convert(Array{Lar.Cells}, triangulated_faces);
V = convert(Lar.Points, W');
GL.VIEW( GL.GLExplode(V,FVs, 1.5,1.5,1.5,99,0.6) );

EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
GL.VIEW( GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1) );

# visualization of component graphs
# EW = Lar.cop2lar(copEV)
# w = convert(Lar.Points, W')
# comps = [ GL.GLLines(w,EW[comp],GL.COLORS[(k-1)%12+1]) for (k,comp) in enumerate(bicon_comps) ]
# GL.VIEW(comps);


text:
using Altro
using TrajectoryOptimization
using RobotDynamics
using Test
using BenchmarkTools
using SparseArrays
using LinearAlgebra
using RobotDynamics: KnotPoint, SampledTrajectory
using SparseArrays
using Altro.Cqdldl
const TO = TrajectoryOptimization
const RD = RobotDynamics

# Generate problem
N = 11
prob, opts = Problems.DubinsCar(:parallel_park, N=N)

opts.constraint_tolerance = opts.projected_newton_tolerance
al = Altro.ALSolver(prob, opts, show_summary=false)
solve!(al)
@test iterations(al) == 16

Zsol = get_trajectory(al)
max_violation(al)

pn = Altro.ProjectedNewtonSolver(prob, opts)
copyto!(pn.Z, Zsol)
copyto!(pn.Z̄data, pn.Zdata)
@test max_violation(pn) ≈ max_violation(al)

# Check trajectories
@test pn.Z ≈ Zsol
@test pn.Z̄ ≈ Zsol

# Check sizes
nx,nu = RD.dims(prob)
N = TO.horizonlength(prob)
Np = Altro.num_primals(pn)
@test Np == sum(nx) + sum(nu)

@test size(pn.hess[end]) == (nx[end]+nu[end], nx[end]+nu[end])
@test size(pn.grad[end]) == (nx[end]+nu[end],)

# Check Hessian and regularization
Altro.cost_hessian!(pn)
iz = pn.iz
@test pn.Atop[iz[end], iz[end]] ≈ cat(prob.obj[end].Q, prob.obj[end].R*0, dims=(1,2))
pn.Atop[iz[end],iz[end]]
Altro.primalregularization!(pn)
ρ_primal = opts.ρ_primal
@test pn.Atop[iz[end], iz[end]] ≈
    cat(prob.obj[end].Q, prob.obj[end].R*0, dims=(1,2)) + I * ρ_primal

# Compute the factorization
Altro.evaluate_constraints!(pn)
Altro.constraint_jacobians!(pn)
Altro.update_active_set!(pn)

A = Altro.getKKTMatrix(pn)
resize!(pn.qdldl, size(A,1))
Cqdldl.eliminationtree!(pn.qdldl, A)
Cqdldl.factor!(pn.qdldl)
F = Cqdldl.QDLDLFactorization(pn.qdldl)

# Try a step
b = Altro.update_b!(pn)
Np = Altro.num_primals(pn)
@test b[1:Np] ≈ zeros(Np)
@test b[Np+1:end] ≈ -pn.d[pn.active[Np+1:end]]
dY = F \ b
pn.Z̄data .= pn.Zdata .+ dY[1:Np]
Altro.evaluate_constraints!(pn)
viol = max_violation(pn, nothing)
@test viol < max_violation(al)
pn.Zdata .= pn.Z̄data  # accept the step

# Take a second step
b = Altro.update_b!(pn)
dY2 = F \ b
pn.Z̄data .= pn.Zdata .+ dY2[1:Np]
Altro.evaluate_constraints!(pn)
viol2 = max_violation(pn, nothing)
@test viol2 < viol

## Solve the entire thing
N = 11
prob, opts = Problems.DubinsCar(:parallel_park, N=N)

opts.constraint_tolerance = opts.projected_newton_tolerance
al = Altro.ALSolver(prob, opts, show_summary=false)
solve!(al)

Zsol = get_trajectory(al)
pn = Altro.ProjectedNewtonSolver(prob, opts)
copyto!(pn.Z, Zsol)
pn.Z[end].dt = Inf
Z0data = copy(pn.Zdata)

opts.constraint_tolerance = 1e-10
viol = solve!(pn)
@test viol < 1e-10

# Test allocations
bpn = @benchmark let pn = $pn
    copyto!(pn.Zdata, $Z0data)
    solve!(pn)
end samples=1 evals=1
@test bpn.allocs == 0

#################################
## Basic tests
#################################

# Create a trajectory with views into a single array
nx,nu = RD.dims(prob)
N = horizonlength(prob)
h = 0.1
NN = N*nx[1] + (N-1)*nu[1]
Zdata = zeros(NN + nu[1])
n,m = nx[1], nu[1]
ix = [(1:n) .+ (k-1)*(n+m) for k = 1:N]
iu = [n .+ (1:m) .+ (k-1)*(n+m) for k = 1:N-1]
iz = [(1:n+m) .+ (k-1)*(n+m) for k = 1:N]
Z = SampledTrajectory([KnotPoint{n,m}(view(Zdata, iz[k]), (k-1)*h, h) for k = 1:N])
Z[end].dt = 0.0
RD.state(Z[2]) .= 1
RD.control(Z[3]) .= 2.1
@test Zdata[ix[2]] ≈ fill(1,n)
@test Zdata[iu[3]] ≈ fill(2.1,m)
Nc = sum(TO.num_constraints(prob))  # number of stage constraints
Np = N*n + (N-1)*m  # number of primals
Nd = N*n + Nc       # number of duals (no constraints)

## Create PNConstraintSet
A = spzeros(Np + Nd, Np + Nd)
b = zeros(Np +  Nd)
a = trues(Nd)

# Add Hessian blocks
blocks = Altro.SparseBlocks(size(A)...)
for k = 1:N-1
    isdiag = TO.is_diag(prob.obj[k])
    Altro.addblock!(blocks, iz[k], iz[k], isdiag)
end
let k = N
    isdiag = TO.is_diag(prob.obj[k])
    Altro.addblock!(blocks, ix[k], ix[k], isdiag)
end

con = prob.constraints[1]
Altro.getinputinds(con, n, m)

pnconset = Altro.PNConstraintSet(prob.constraints, Z, opts, A, b, a, blocks)
alconset = Altro.ALConstraintSet{Float64}()
Altro.initialize!(alconset, prob.constraints, prob.Z, opts, zeros(N))
@test length(pnconset) == length(alconset)

copyto!(prob.Z, Z)
prob.Z ≈ Z
Altro.evaluate_constraints!(pnconset, Z)
Altro.evaluate_constraints!(alconset, prob.Z)
for i = 1:length(alconset)
    @test pnconset[i].vals ≈ alconset[i].vals
end

Altro.constraint_jacobians!(pnconset, Z)
Altro.constraint_jacobians!(alconset, prob.Z)
for i = 1:length(alconset) - 1
    @test pnconset[i].jac ≈ alconset[i].jac
end

# Create a PN Solver
al = Altro.ALSolver(prob)
ilqr = al.ilqr
pn = Altro.ProjectedNewtonSolver(prob)
Np = Altro.num_primals(pn)

pnconset = pn.conset
@test pnconset.cinds[1][1] == n .+ (1:2m)
@test pnconset.cinds[end][1] == 1:n
@test pnconset.cinds[end][2] == n + 2m .+ (1:n)
@test pnconset.cinds[1][2] == n + 2m + n .+ (1:2m)
@test pnconset.cinds[2][1] == n + 2m + n + 2m .+ (1:4)
@test pnconset.cinds[end][3] == n + 2m + n + 2m + 4 .+ (1:n)

# Cost expansion
copyto!(pn.Z, ilqr.Z)
copyto!(pn.Z̄data, pn.Zdata)
Altro.cost_gradient!(pn)
Altro.cost_hessian!(pn)
Altro.cost_expansion!(ilqr.obj.obj, ilqr.Efull, ilqr.Z)
for k = 1:N
    @test ilqr.Efull[k].hess ≈ pn.hess[k]
    @test ilqr.Efull[k].grad ≈ pn.grad[k]
end

# Dynamics Jacobian
Altro.constraint_jacobians!(pn)
Altro.dynamics_expansion!(ilqr)
for k = 1:N-1
    @test ilqr.D[k].∇f ≈ pn.∇f[k]
end

# Dynamics error
Altro.dynamics_error!(pn)
for k = 1:N-1
    @test pn.e[k] ≈ RD.discrete_dynamics(prob.model[1], prob.Z[k]) - RD.state(prob.Z[k+1])
end

# Update everything
let pn = pn
    Altro.evaluate_constraints!(pn)
    Altro.constraint_jacobians!(pn)
    Altro.cost_gradient!(pn)
    Altro.cost_hessian!(pn)
    Altro.update_active_set!(pn)
end


text:
# fom Integer to Posit
Posit8(x::Int64) = ccall((:i64_to_p8, SoftPositPath), Posit8, (Int64,),x)
Posit8(x::Int32) = ccall((:i32_to_p8, SoftPositPath), Posit8, (Int32,),x)

Posit16(x::Int64) = ccall((:i64_to_p16, SoftPositPath), Posit16, (Int64,),x)
Posit16(x::Int32) = ccall((:i32_to_p16, SoftPositPath), Posit16, (Int32,),x)

Posit32(x::Int64) = ccall((:i64_to_p32, SoftPositPath), Posit32, (Int64,),x)
Posit32(x::Int32) = ccall((:i32_to_p32, SoftPositPath), Posit32, (Int32,),x)

# from Int64 to Posit_1
Posit8_1(x::Int64) = ccall((:i64_to_pX1, SoftPositPath), Posit8_1, (Int64,Int64),x,8)
Posit16_1(x::Int64) = ccall((:i64_to_pX1, SoftPositPath), Posit16_1, (Int64,Int64),x,16)
Posit24_1(x::Int64) = ccall((:i64_to_pX1, SoftPositPath), Posit24_1, (Int64,Int64),x,24)

# from Int64 to Posit_2
Posit8_2(x::Int64) = ccall((:i64_to_pX2, SoftPositPath), Posit8_2, (Int64,Int64),x,8)
Posit16_2(x::Int64) = ccall((:i64_to_pX2, SoftPositPath), Posit16_2, (Int64,Int64),x,16)
Posit24_2(x::Int64) = ccall((:i64_to_pX2, SoftPositPath), Posit24_2, (Int64,Int64),x,24)

# from Int32 to Posit_1
Posit8_1(x::Int32) = ccall((:i32_to_pX1, SoftPositPath), Posit8_1, (Int32,Int64),x,8)
Posit16_1(x::Int32) = ccall((:i32_to_pX1, SoftPositPath), Posit16_1, (Int32,Int64),x,16)
Posit24_1(x::Int32) = ccall((:i32_to_pX1, SoftPositPath), Posit24_1, (Int32,Int64),x,24)

# from Int32 to Posit_2
Posit8_2(x::Int32) = ccall((:i32_to_pX2, SoftPositPath), Posit8_2, (Int32,Int64),x,8)
Posit16_2(x::Int32) = ccall((:i32_to_pX2, SoftPositPath), Posit16_2, (Int32,Int64),x,16)
Posit24_2(x::Int32) = ccall((:i32_to_pX2, SoftPositPath), Posit24_2, (Int32,Int64),x,24)

# promotion
promote_rule(::Type{Int64},::Type{Posit32}) = Posit32
promote_rule(::Type{Int32},::Type{Posit32}) = Posit32
promote_rule(::Type{Int64},::Type{Posit16}) = Posit16
promote_rule(::Type{Int32},::Type{Posit16}) = Posit16
promote_rule(::Type{Int64},::Type{Posit8}) = Posit8
promote_rule(::Type{Int32},::Type{Posit8}) = Posit8

promote_rule(::Type{Int64},::Type{Posit24_1}) = Posit24_2
promote_rule(::Type{Int32},::Type{Posit24_1}) = Posit24_2
promote_rule(::Type{Int64},::Type{Posit16_1}) = Posit16_1
promote_rule(::Type{Int32},::Type{Posit16_1}) = Posit16_1
promote_rule(::Type{Int64},::Type{Posit8_1}) = Posit8_1
promote_rule(::Type{Int32},::Type{Posit8_1}) = Posit8_1

promote_rule(::Type{Int64},::Type{Posit24_2}) = Posit24_2
promote_rule(::Type{Int32},::Type{Posit24_2}) = Posit24_2
promote_rule(::Type{Int64},::Type{Posit16_2}) = Posit16_2
promote_rule(::Type{Int32},::Type{Posit16_2}) = Posit16_2
promote_rule(::Type{Int64},::Type{Posit8_2}) = Posit8_2
promote_rule(::Type{Int32},::Type{Posit8_2}) = Posit8_2


text:
using HDF5
using DataFrames
using Dates

function getDailyDF(f::String) :: DataFrame
    h5f  = h5open(f,"r")
    dd0  = h5f["daily/block0_values"] # Integer64 (datetime)
    dd1  = h5f["daily/block1_values"] # Floats
    dd2  = h5f["daily/block2_values"] # Integers (count and daylight indicator)
    dd0n = h5f["daily/block0_items"]
    dd1n = h5f["daily/block1_items"]
    dd2n = h5f["daily/block2_items"]
    dfd0 = DataFrame(read(dd0)')
    dfd1 = DataFrame(read(dd1)')
    dfd2 = DataFrame(read(dd2)')
    ddn = vcat(read(dd0n), read(dd1n), read(dd2n))

    #+
    # Concatenate dataframes, add column labels and convert datetime
    #-
    dfd = hcat(dfd0,dfd1,dfd2)
    names!(dfd,Symbol.(ddn))
    dfd[:date] = unix2datetime.(div.(dfd[:date],1000000000))

    #+
    # Convert the integer daytime indicator to a categorical, descriptive
    # variable.
    #-
    dfd[:daytime] = CategoricalArray(map(x-> if x==1 return "Day" else return "Night" end , dfd[:daytime]))
    close(h5f)
    return dfd
end

function getHourlyDF(f::String) :: DataFrame
    h5f  = h5open(f,"r")
    dh0 =  h5f["hourly/block0_values"] # Floats dataset
    dh1  = h5f["hourly/block1_values"] # integer dataset (count)
    dh2  = h5f["hourly/block2_values"] # integer64 dataset (time)
    dh0n = h5f["hourly/block0_items"]
    dh1n = h5f["hourly/block1_items"]
    dh2n = h5f["hourly/block2_items"]

    #+
    # Create dataframes from the dataset contents and prepare an
    # array with the field labels
    #
    dfh0 = DataFrame(read(dh0)')
    dfh1 = DataFrame(read(dh1)')
    dfh2 = DataFrame(read(dh2)')
    dfn = vcat(read(dh0n), read(dh1n), read(dh2n))

    #+
    # Concatenate the datasets, assign column names and convert int64 to datetime
    #-
    dfh = hcat(dfh0,dfh1,dfh2)
    names!(dfh,Symbol.(dfn))
    dfh[:time] = unix2datetime.(div.(dfh[:time],1000000000))
    close(h5f)
    return dfh
end

function getRawDF(f::String) :: DataFrame
    h5f  = h5open(f,"r")

    dr0  = h5f["rawdata/block0_values"] # Floats dataset
    dr1  = h5f["rawdata/block1_values"] # integer64 (datetime) dataset
    dr0n = h5f["rawdata/block0_items"]
    dr1n = h5f["rawdata/block1_items"]

    #+
    # Create DataFrames for the block0 (floats) and block1 (integer64 containing
    # times in nanosecods). Note the tick used to transpose the read matrices.
    # Build an array with the field labels from both datasets.
    #-
    dfr0 = DataFrame(read(dr0)')
    dfr1 = DataFrame(read(dr1)')
    drn  = vcat(read(dr1n),read(dr0n))

    #+
    # Concatenate both DataFrames and assign the column names.
    #-
    dfr = hcat(dfr1,dfr0)
    names!(dfr,Symbol.(drn))

    #+
    # Convert the interger64 values to proper times. Remember those
    # times are UTC, but the Julia datatime is plain (has no TZ info)
    #-
    dfr[:time] = unix2datetime.(div.(dfr[:time],1000000000))
    close(h5f)
    return dfr
end

max_stars_repo_name: jguillaumes/MeteoClient

max_issues_repo_name: jguillaumes/MeteoClient

max_forks_repo_name: jguillaumes/MeteoClient


text:
type MarketModelPathwiseCashRebate <: MarketModelPathwiseMultiProduct
  evolution::EvolutionDescription
  paymentTimes::Vector{Float64}
  amounts::Matrix{Float64}
  numberOfProducts::Int
  currentIndex::Int
end

function MarketModelPathwiseCashRebate(evolution::EvolutionDescription,
                                      paymentTimes::Vector{Float64},
                                      amounts::Matrix{Float64},
                                      numberOfProducts::Int)

  check_increasing_times(paymentTimes)
  amounts_rows, amounts_cols = size(amounts)

  amounts_rows == numberOfProducts || error("the number of rows in the matrix must equal the number of products")
  amounts_cols == length(paymentTimes) || error("the number of cols in the matrix must equal the number of payment times")
  length(evolution.evolutionTimes) == length(paymentTimes) || error("the number of evolution times must equal the number of payment times")

  return MarketModelPathwiseCashRebate(evolution, paymentTimes, amounts, numberOfProducts, 1)
end

possible_cash_flow_times(mm::MarketModelPathwiseCashRebate) = mm.paymentTimes

max_number_of_cashflows_per_product_per_step(::MarketModelPathwiseCashRebate) = 1

reset!(mm::MarketModelPathwiseCashRebate) = mm.currentIndex = 1

function next_time_step!(mm::MarketModelPathwiseCashRebate, ::CurveState, numberCashFlowsThisStep::Vector{Int}, cashFlowsGenerated::Vector{Vector{MarketModelPathWiseCashFlow}})
  for i = 1:mm.numberOfProducts
    numberCashFlowsThisStep[i] = 1
    cashFlowsGenerated[i][1].timeIndex = mm.currentIndex
    cashFlowsGenerated[i][1].amount[1] = mm.amounts[1, mm.currentIndex]

    for k = 2:mm.evolution.numberOfRates + 1
      cashFlowsGenerated[i][1].amount[k] = 0.0
    end
  end

  mm.currentIndex += 1
  return true
end

# Clone #
clone(mm::MarketModelPathwiseCashRebate) = MarketModelPathwiseCashRebate(clone(mm.evolution), copy(mm.paymentTimes), copy(mm.amounts), mm.numberOfProducts, mm.currentIndex)


text:
import SpecialFunctions
import SpecialFunctions: besselj

@testset "Bessel Roots" begin
    @test_throws DomainError approx_besselroots(0.0, -1)

    msg = "`besselroots` has been renamed to `approx_besselroots`, and `besselroots` will be removed in the next major release."
    @test_logs (:warn, msg) besselroots(0.2,3)

    # Check if besselj(ν, approx_besselroots(ν, n) ) is small:
    tol = 1e-11

    for ν = 0.:0.1:5.
        n = 10
        @test norm( besselj.(ν, approx_besselroots(ν, n) ), Inf ) < tol
    end
end


text:
# julia
using FASTX
# main
print("print the full header line from fasta file\nusage: julia fasta_to_header.jl input_fasta output_txt")
reader = open(FASTA.Reader,ARGS[1])
open(ARGS[2],"a") do io
   for record in reader
      println(io,FASTA.identifier(record),"\t",FASTA.description(record))
   end
end
close(reader)


text:
# # A vector of SVectors, denoting the edge ID each face is connected to.
# function face_edge_connectivity(mesh::QuadrilateralMesh{Dim,T,U}) where {Dim,T,U}
#     if length(mesh.edges) === 0
#         @error "Mesh does not have edges!"
#     end
#     # A vector of MVectors of zeros for each face
#     # Each MVector is the length of the number of edges
#     face_edge = [MVector{4, U}(0, 0, 0, 0) for _ in eachindex(mesh.faces)]
#     # For each face in the mesh, generate the edges.
#     # Search for the index of the edge in the mesh.edges vector
#     # Insert the index of the edge into the face_edge connectivity vector
#     for i in eachindex(mesh.faces)
#         for (j, edge) in enumerate(edges(mesh.faces[i]))
#             face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
#         end
#         if any(x->x === U(0), face_edge[i])
#             @error "Could not determine the face/edge connectivity of face $i"
#         end
#     end
#     return [SVector(sort!(conn).data) for conn in face_edge]
# end
#
# # A vector of SVectors, denoting the edge ID each face is connected to.
# function face_edge_connectivity(mesh::TriangleMesh{Dim,T,U}) where {Dim,T,U}
#     if length(mesh.edges) === 0
#         @error "Mesh does not have edges!"
#     end
#     # A vector of MVectors of zeros for each face
#     # Each MVector is the length of the number of edges
#     face_edge = [MVector{3, U}(0, 0, 0) for _ in eachindex(mesh.faces)]
#     # For each face in the mesh, generate the edges.
#     # Search for the index of the edge in the mesh.edges vector
#     # Insert the index of the edge into the face_edge connectivity vector
#     for i in eachindex(mesh.faces)
#         for (j, edge) in enumerate(edges(mesh.faces[i]))
#             face_edge[i][j] = searchsortedfirst(mesh.edges, SVector(edge.data))
#         end
#         if any(x->x === U(0), face_edge[i])
#             @error "Could not determine the face/edge connectivity of face $i"
#         end
#     end
#     return [SVector(sort!(conn).data) for conn in face_edge]
# end
#
# #
# # # Return a mesh with boundary edges
# # function add_boundary_edges(mesh::M; boundary_shape="Unknown"
# #     ) where {M <: UnstructuredMesh_2D}
# #     if 0 === length(mesh.edge_face_connectivity)
# #         mesh = add_connectivity(mesh)
# #     end
# #     return M(name = mesh.name,
# #              points = mesh.points,
# #              edges = mesh.edges,
# #              materialized_edges = mesh.materialized_edges,
# #              faces = mesh.faces,
# #              materialized_faces = mesh.materialized_faces,
# #              edge_face_connectivity = mesh.edge_face_connectivity,
# #              face_edge_connectivity = mesh.face_edge_connectivity,
# #              boundary_edges = boundary_edges(mesh, boundary_shape),
# #              face_sets = mesh.face_sets
# #             )
# # end
# #
#
# # # Return a vector of the faces adjacent to the face of ID face
# # function adjacent_faces(face::UInt32, mesh::UnstructuredMesh_2D)
# #     edges = mesh.face_edge_connectivity[face]
# #     the_adjacent_faces = UInt32[]
# #     for edge in edges
# #         faces = mesh.edge_face_connectivity[edge]
# #         for face_id in faces
# #             if face_id != face && face_id != 0
# #                 push!(the_adjacent_faces, face_id)
# #             end
# #         end
# #     end
# #     return the_adjacent_faces
# # end
# #
# # # Return a vector containing vectors of the edges in each side of the mesh's bounding shape, e.g.
# # # For a rectangular bounding shape the sides are North, East, South, West. Then the output would
# # # be [ [e1, e2, e3, ...], [e17, e18, e18, ...], ..., [e100, e101, ...]]
# # function boundary_edges(mesh::UnstructuredMesh_2D, boundary_shape::String)
# #     # edges which have face 0 in their edge_face connectivity are boundary edges
# #     the_boundary_edges = UInt32.(findall(x->x[1] === 0x00000000, mesh.edge_face_connectivity))
# #     if boundary_shape == "Rectangle"
# #         # Sort edges into NESW
# #         bb = boundingbox(mesh.points)
# #         y_north = bb.ymax
# #         x_east  = bb.xmax
# #         y_south = bb.ymin
# #         x_west  = bb.xmin
# #         p_NW = Point_2D(x_west, y_north)
# #         p_NE = bb.tr
# #         p_SE = Point_2D(x_east, y_south)
# #         p_SW = bb.bl
# #         edges_north = UInt32[]
# #         edges_east = UInt32[]
# #         edges_south = UInt32[]
# #         edges_west = UInt32[]
# #         # Insert edges so that indices move from NW -> NE -> SE -> SW -> NW
# #         for edge ∈  the_boundary_edges
# #             epoints = edgepoints(mesh.edges[edge], mesh.points)
# #             if all(x->abs(x[2] - y_north) < 1e-4, epoints)
# #                 insert_boundary_edge!(edge, p_NW, edges_north, mesh)
# #             elseif all(x->abs(x[1] - x_east) < 1e-4, epoints)
# #                 insert_boundary_edge!(edge, p_NE, edges_east, mesh)
# #             elseif all(x->abs(x[2] - y_south) < 1e-4, epoints)
# #                 insert_boundary_edge!(edge, p_SE, edges_south, mesh)
# #             elseif all(x->abs(x[1] - x_west) < 1e-4, epoints)
# #                 insert_boundary_edge!(edge, p_SW, edges_west, mesh)
# #             else
# #                 @error "Edge $edge could not be classified as NSEW"
# #             end
# #         end
# #         return [ edges_north, edges_east, edges_south, edges_west ]
# #     else
# #         return [ convert(Vector{UInt32}, the_boundary_edges) ]
# #     end
# # end
# #
#
#
# # # Find the faces which share the vertex of ID v.
# # function faces_sharing_vertex(v::Integer, mesh::UnstructuredMesh_2D)
# #     shared_faces = UInt32[]
# #     for i ∈ 1:length(mesh.faces)
# #         if v ∈  mesh.faces[i]
# #             push!(shared_faces, UInt32(i))
# #         end
# #     end
# #     return shared_faces
# # end
# #
# #
# # # Insert the boundary edge into the correct place in the vector of edge indices, based on
# # # the distance from some reference point
# # function insert_boundary_edge!(edge_index::UInt32, p_ref::Point_2D, edge_indices::Vector{UInt32},
# #                                mesh::UnstructuredMesh_2D)
# #     # Compute the minimum distance from the edge to be inserted to the reference point
# #     insertion_distance = minimum(distance.(Ref(p_ref),
# #                                            edgepoints(mesh.edges[edge_index], mesh.points)))
# #     # Loop through the edge indices until an edge with greater distance from the reference point
# #     # is found, then insert
# #     nindices = length(edge_indices)
# #     for i ∈ 1:nindices
# #         epoints = edgepoints(mesh.edges[edge_indices[i]], mesh.points)
# #         edge_distance = minimum(distance.(Ref(p_ref), epoints))
# #         if insertion_distance < edge_distance
# #             insert!(edge_indices, i, edge_index)
# #             return nothing
# #         end
# #     end
# #     insert!(edge_indices, nindices+1, edge_index)
# #     return nothing
# # end
# #
# #
#
# # # If a point is a vertex
# # function isvertex(p::Point_2D, mesh::UnstructuredMesh_2D)
# #     for point in mesh.points
# #         if p ≈ point
# #             return true
# #         end
# #     end
# #     return false
# # end
# #
#
# #
# # function remap_points_to_hilbert(points::Vector{Point_2D})
# #     bb = boundingbox(points)
# #     npoints = length(points)
# #     # Generate a Hilbert curve
# #     hilbert_points = hilbert_curve(bb, npoints)
# #     nhilbert_points = length(hilbert_points)
# #     # For each point, get the index of the hilbert points that is closest
# #     point_indices = Vector{Int64}(undef, npoints)
# #     for i ∈ 1:npoints
# #         min_distance = 1.0e10
# #         for j ∈ 1:nhilbert_points
# #             pdistance = distance²(points[i], hilbert_points[j])
# #             if pdistance < min_distance
# #                 min_distance = pdistance
# #                 point_indices[i] = j
# #             end
# #         end
# #     end
# #     return sortperm(point_indices)
# # end
# #
# # function reorder_points_to_hilbert!(mesh::UnstructuredMesh_2D)
# #     # Points
# #     # point_map     maps  new_points[i] == mesh.points[point_map[i]]
# #     # point_map_inv maps mesh.points[i] == new_points[point_map_inv[i]]
# #     point_map  = remap_points_to_hilbert(mesh.points)
# #     point_map_inv = UInt32.(sortperm(point_map))
# #     # reordered to resemble a hilbert curve
# #     permute!(mesh.points, point_map)
# #
# #     # Adjust face indices
# #     # Point IDs have changed, so we need to change the point IDs referenced by the faces
# #     new_faces_vec = [ point_map_inv[face] for face in mesh.faces]
# #     for i in 1:length(mesh.faces)
# #         mesh.faces[i] = SVector(point_map_inv[mesh.faces[i]])
# #     end
# #     return nothing
# # end
# #
# # function reorder_faces_to_hilbert!(mesh::UnstructuredMesh_2D)
# #     face_map  = UInt32.(remap_points_to_hilbert(centroid.(materialize_faces(mesh))))
# #     permute!(mesh.faces, face_map)
# #     return nothing
# # end
# #
# # function reorder_to_hilbert!(mesh::UnstructuredMesh_2D)
# #     reorder_points_to_hilbert!(mesh)
# #     reorder_faces_to_hilbert!(mesh)
# #     return nothing
# # end
# #
# # # Return the ID of the edge shared by two adjacent faces
# # function shared_edge(face1::UInt32, face2::UInt32, mesh::UnstructuredMesh_2D)
# #     for edge1 in mesh.face_edge_connectivity[face1]
# #         for edge2 in mesh.face_edge_connectivity[face2]
# #             if edge1 == edge2
# #                 return edge1
# #             end
# #         end
# #     end
# #     return 0x00000000
# # end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
import Rotations.RotMatrix

"""
    degree_rotation(deg)

Create a rotation matrix for the angle given in degrees.
"""
function degree_rotation(deg::Real)
    RotMatrix(pi * (deg / 180))
end

"""
    frac_rotation(frac)

Create a rotation matrix for the angle given as a fraction,
i.e. 1.0 represents a full rotation and 0.5 a half rotation.
"""
function frac_rotation(frac::Real)
    RotMatrix(2 * pi * frac)
end

"""
    radian_rotation(rad)

Create a rotation matrix for the angle given in radians.
"""
function radian_rotation(rad::Real)
    RotMatrix(rad)
end

export degree_rotation, frac_rotation, radian_rotation


text:
# X-ray Window

abstract type AbstractWindow end

Base.show(io::IO, wnd::AbstractWindow) = print(io, wnd.name)

struct NoWindow <: AbstractWindow
    name::String
    NoWindow() = new("No window")
end

NeXLCore.transmission(wnd::NoWindow, energy::Float64, angle::Float64 = π / 2) = 1.0

struct LayerWindow <: AbstractWindow
    name::String
    layers::Vector{Film}
    support::Film
    openfraction::Float64
end

NeXLCore.transmission(wnd::LayerWindow, energy::Float64, angle::Float64 = π / 2) =
    (
        length(wnd.layers) > 0 ?
        mapreduce(lyr -> NeXLCore.transmission(lyr, energy, angle), *, wnd.layers) : 1.0
    ) * (
        wnd.openfraction +
        (1.0 - wnd.openfraction) * NeXLCore.transmission(wnd.support, energy, angle)
    )

"""
    AP33Model()

Construct a modeled window for the Moxtek AP3.3 window.
"""

function AP33Model()
    support, openarea = Film(pure(n"Si"), 0.038), 0.77
    paralene = Film(parse(Material, "C10H8O4N", density = 1.39), 3.0e-5)
    aluminum = Film(pure(n"Al"), 4.0e-6)
    return LayerWindow("Moxtek AP3.3 model", [aluminum, paralene], support, openarea)
end

function AP5Model()
    support, openarea = Film(pure(n"C"), 0.0265), 0.78
    paralene = Film(parse(Material, "C10H8O4N", density = 1.39), 3.0e-5)
    aluminum = Film(pure(n"Al"), 4.0e-6)
    return LayerWindow("Moxtek AP5 model", [aluminum, paralene], support, openarea)
end

"""
    Beryllium(thickness=5.0e-4)

Construct a beryllium window.
"""
function Beryllium(thickness = 5.0e-4)
    support, openarea = Film(pure(n"Be"), thickness), 0.0
    return LayerWindow("$(thickness*1.0e4) μm Be window", [], support, openarea)
end

"""
    function AmptekC1()
    function AmptekC2()

Create modeled windows for the Amptek C1 or C2 Si₃N₄ windows according to the specifications here:
https://www.amptek.com/products/accessories-for-xrf-eds/c-series-low-energy-x-ray-windows#Specifications
"""

function AmptekC1()
    si3n4 = Film(parse(Material, "Si3N4", density = 3.44), 150.0e-7)
    al = Film(pure(n"Al"), 250.0e-7)
    support, openarea = Film(pure(n"Si"), 1.50e-4), 0.80
    return LayerWindow("Amptek C1", [si3n4, al], support, openarea)
end

function AmptekC2()
    si3n4 = Film(parse(Material, "Si3N4", density = 3.44), 40.0e-7)
    al = Film(pure(n"Al"), 30.0e-7)
    support, openarea = Film(pure(n"Si"), 1.50e-4), 0.80
    return LayerWindow("Amptek C2", [si3n4, al], support, openarea)
end

struct TabulatedWindow <: AbstractWindow
    name::String
    interpolation::AbstractInterpolation
    extrapolation::LayerWindow
    match::Float64
end

function NeXLCore.transmission(
    wnd::TabulatedWindow,
    energy::Float64,
    angle::Float64 = π / 2,
)
    bnds = bounds(parent(wnd.interpolation))[1]
    return (energy >= bnds[1]) && (energy < bnds[2]) ? wnd.interpolation(energy) : #
           wnd.match * transmission(wnd.extrapolation, energy, angle)
end

"""
    AP33Tabulation()
    AP5Tabulation()

Construct tabulated window models for the Moxtek AP3.3 and AP5 windows.
"""
function AP33Tabulation()
    data = CSV.read(
        joinpath(dirname(pathof(@__MODULE__)), "AP3_3_mod.csv"),
        DataFrame,
        comment = "//"
    )
    inter = LinearInterpolation(data[:, 1], data[:, 2])
    extra = AP33Model()
    return TabulatedWindow("Moxtek AP3.3", inter, extra, 1.01582935271)
end

function AP5Tabulation()
    data = CSV.read(joinpath(dirname(pathof(@__MODULE__)), "AP5.csv"), DataFrame)
    inter = LinearInterpolation(data[:, 1], data[:, 2])
    extra = AP5Model()
    return TabulatedWindow("Moxtek AP5", inter, extra, 1.0)
end
""

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
abstract PixelType

macro wraptype(T)
  quote
    type $T <: PixelType end
    export $T
  end
end

suffix(T) = symbol(split(string(T),".")[2])

macro dlsym(func, lib)
  z, zlocal = gensym(string(func)), gensym()
  eval(current_module(),:(global $z = C_NULL))
  z = esc(z)
  quote
    let $zlocal::Ptr{Void} = $z::Ptr{Void}
      if $zlocal == C_NULL
        $zlocal = Libdl.dlsym($(esc(lib))::Ptr{Void}, $(esc(func)))
        global $z = $zlocal
      end
      $zlocal
    end
  end
end

macro itkcall(fname, T, rettype, argtypes, args...)
  quote
    ccall(@dlsym(symbol($fname,suffix($T)),
            Libdl.dlopen("libcitkffi")),
          $rettype, $argtypes, $(args...))
  end
end

type Image{T <: PixelType}
  ptr::Ptr{Void}
  function Image(p::Ptr{Void})
    self = new(p)
    finalizer(self, function(x)
      if x.ptr != 0
        @itkcall("DeleteImage", T, Void, (Ptr{Void},), x.ptr)
        x.ptr = 0
      end
    end)
    self
  end
end

newimg{T}(::Type{T}) = begin
  ptr = @itkcall("", T, Ptr{Void}, ())
  Image{T}(ptr)
end

newimg{T}(image::Image{T}) = begin
  ptr = @itkcall("", T, Ptr{Void}, ())
  Image{T}(ptr)
end

pixeltype{T}(image::Image{T}) = begin
  T
end

readimg{T}(image::Image{T}, filename::AbstractString) = begin
  @itkcall("ReadImage", T, Void, (Ptr{Void},Ptr{UInt8}), image.ptr, filename)
end

readimg{T}(::Type{T}, filename::AbstractString) = begin
  ptr = @itkcall("", T, Ptr{Void}, ())
  @itkcall("ReadImage", T, Void, (Ptr{Void},Ptr{UInt8}), ptr, filename)
  Image{T}(ptr)
end

writeimg{T}(image::Image{T}, filename::AbstractString) = begin
  @itkcall("WriteImage", T, Void, (Ptr{Void},Ptr{UInt8}), image.ptr, filename)
end

export newimg
export pixeltype
export readimg
export writeimg

# image
# createanother
# pixeltype
# delete
# allocate
# dim
# setdim
# spacing
# setspacing
# origin
# setorigin
# direction
# setdirectoin
# bufferedsize
# setbufferedsize
# size
# setsize
# read
# write
# toarr
# fromarr

max_stars_repo_name: orobix/itkffi

max_issues_repo_name: orobix/itkffi

max_forks_repo_name: orobix/itkffi


text:
mutable struct PDataMA97{T} <: PDataFact{T}
    PLD::Ma97          # H = (PL)D(PL)'  or  LDL' = P'HP
    Δ::Array{T,1}  # diagonal, eigenvalues of D
    D::Any
    H::Any
    pivot::Any
    #Q::Array{Float64,2}  # orthogonal matrix, eigenvectors of D:  should be sparse
    # QΔQ'  =  D
    Q::SparseMatrixCSC{T,Int64}
    g̃::Array{T,1}  # transformed gradient
    g::Array{T,1}  # untransformed gradient
    λ::T
    success::Bool        # previous iteration was successfull
    OK::Bool             # preprocess success

    PDataMA97() = new{Nothing}()
    PDataMA97(PLD, Δ, D, H, pivot, Q, g̃, g, l, success, OK) =
        new{eltype(g̃)}(PLD, Δ, D, H, pivot, Q, g̃, g, l, success, OK)
end

function preprocess(PData::PDataMA97, H, g, n1, n2)
    PLD = Ma97
    pivot = Array{Int32,1}
    vD1 = Array{Float64,2}
    vD2 = Array{Float64,2}

    H97 = convert(SparseMatrixCSC{Float64,Int}, H)
    try
        PLD = Ma97(H97, print_level = -1)
        ma97_factorize!(PLD)
    catch
        println("*******   Problem in MA97")
        res = PDataMA97()
        res.OK = false
        return res
    end
    try
        ret = ma97_inquire(PLD)
        vD = ret[2]
        pivot = ret[1]
        # Take care of zero valued pivots
        z(a) = a == 0
        pivot[findall(z, pivot)] = collect(1:length(pivot))[findall(z, pivot)]

        vD1 = copy(vD[1:1, :])'[:, 1:1]
        vD2 = copy(vD[2:2, 1:end-1])'[:, 1:1]
    catch
        println("*******   Problem after MA97")
        println(" Cond(H) = $(cond(full(H)))")
        res = PDataMA97()
        res.OK = false
        return res
    end

    vD1 = vec(vD1)
    vD2 = vec(vD2)
    D = SparseArrays.spdiagm(0 => vD1, -1 => vD2, 1 => vD2)

    #################  Future object BlockDiag operator?
    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
    vQ2 = zeros(length(vD2))      #
    vQ2 = zeros(length(vD2))      #
    vQ2m = zeros(length(vD2))     #
    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
    # if D diagonal, nothing more will be computed

    i = 1
    while i < length(vD1)
        if vD2[i] == 0.0
            i += 1
        else
            mA = [vD1[i] vD2[i]; vD2[i] vD1[i+1]]  #  2X2 submatrix
            # DiagmA, Qma = eig(mA)                 #  spectral decomposition of mA
            X = eigen(mA)
            DiagmA = X.values
            Qma = X.vectors
            veig[i] = DiagmA[1]
            vQ1[i] = Qma[1, 1]
            vQ2[i] = Qma[1, 2]
            vQ2m[i] = Qma[2, 1]
            vQ1[i+1] = Qma[2, 2]
            veig[i+1] = DiagmA[2]
            i += 2
        end
    end

    Q = sparse(SparseArrays.spdiagm(0 => vQ1, -1 => vQ2m, 1 => vQ2)) # sparse representation of Q

    Q = Q[abs.(pivot), abs.(pivot)]
    veig = veig[abs.(pivot)]

    Δ = 1.0 ./ veig #'[:,1]             # don't forget inquire returns D^-1
    l_m, = findmin(Δ)

    ĝ = copy(g)
    ma97_solve!(PLD, ĝ, job = :PL)

    g̃ = Q' * ĝ

    λ = max(-l_m, 0.0)
    return PDataMA97(PLD, Δ, D, H, pivot, Q, g̃, g, λ, true, true)
end

function AInv(X::PDataMA97, d̃::Array{T,1}) where {T}
    d̂ = X.Q * d̃

    ma97_solve!(X.PLD, d̂, job = :LPS)

    return d̂
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

import Compat.ASCIIString
using CSVFeatures
using DataFrameFeatures
using RLESUtils, MathUtils, LookupCallbacks, FileUtils, StringUtils
using DataFrames

const FEATURE_MAP = LookupCallback[
    LookupCallback("ra_detailed.ra_active", x -> Bool(x)),
    LookupCallback("ra_detailed.ownInput.dz", x -> Float64(x)),
    LookupCallback(["ra_detailed.ownInput.z", "ra_detailed.intruderInput[1].z"], (z1, z2) -> Float64(z2 - z1)), #intruderInput.z is quantized!!!
    LookupCallback("ra_detailed.ownInput.psi", x -> Float64(x)),
    LookupCallback("ra_detailed.intruderInput[1].sr", x -> Float64(x)),
    LookupCallback("ra_detailed.intruderInput[1].chi", x -> Float64(x)),
    LookupCallback("ra_detailed.intruderInput[1].vrc", x -> x == 0), #split categorical to 1-hot
    LookupCallback("ra_detailed.intruderInput[1].vrc", x -> x == 1), #split categorical to 1-hot
    LookupCallback("ra_detailed.intruderInput[1].vrc", x -> x == 2), #split categorical to 1-hot
    LookupCallback("ra_detailed.ownOutput.cc", x -> bin(Int(x), 3)[1] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.cc", x -> bin(Int(x), 3)[2] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.cc", x -> bin(Int(x), 3)[3] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.vc", x -> bin(Int(x), 3)[1] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.vc", x -> bin(Int(x), 3)[2] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.vc", x -> bin(Int(x), 3)[3] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.ua", x -> bin(Int(x), 3)[1] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.ua", x -> bin(Int(x), 3)[2] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.ua", x -> bin(Int(x), 3)[3] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.da", x -> bin(Int(x), 3)[1] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.da", x -> bin(Int(x), 3)[2] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.da", x -> bin(Int(x), 3)[3] == '1'), #categorical 3-bit
    LookupCallback("ra_detailed.ownOutput.target_rate", x -> Float64(x)),
    LookupCallback("ra_detailed.ownOutput.crossing", x -> Bool(x)),
    LookupCallback("ra_detailed.ownOutput.alarm", x -> Bool(x)),
    LookupCallback("ra_detailed.ownOutput.alert", x -> Bool(x)),
    LookupCallback("ra_detailed.intruderOutput[1].vrc", x -> x == 0), #split categorical to 1-hot
    LookupCallback("ra_detailed.intruderOutput[1].vrc", x -> x == 1), #split categorical to 1-hot
    LookupCallback("ra_detailed.intruderOutput[1].vrc", x -> x == 2), #split categorical to 1-hot
  LookupCallback("ra_detailed.intruderOutput[1].tds", x -> Float64(x)),
    LookupCallback("response.state", x -> x == "none"), #split categorical to 1-hot
    LookupCallback("response.state", x -> x == "stay"), #split categorical to 1-hot
    LookupCallback("response.state", x -> x == "follow"), #split categorical to 1-hot
    LookupCallback("response.timer", x -> Float64(x)),
    LookupCallback("response.h_d", x -> Float64(x)),
    LookupCallback("response.psi_d", x -> Float64(x)),
    LookupCallback("adm.v", x -> Float64(x)),
    LookupCallback("adm.h", x -> Float64(x)) #same as z
    ]

const FEATURE_NAMES = ASCIIString[
    "RA",
    "vert_rate",
    "alt_diff",
    "psi",
    "intr_sr",
    "intr_chi",
    "intr_vrc0",
    "intr_vrc1",
    "intr_vrc2",
    "cc0",
    "cc1",
    "cc2",
    "vc0",
    "vc1",
    "vc2",
    "ua0",
    "ua1",
    "ua2",
    "da0",
    "da1",
    "da2",
    "target_rate",
    "crossing",
    "alarm",
    "alert",
    "intr_out_vrc0",
    "intr_out_vrc1",
    "intr_out_vrc2",
    "intr_out_tds",
    "response_none",
    "response_stay",
    "response_follow",
    "response_timer",
    "response_h_d",
    "response_psi_d",
    "v",
    "h"
    ]

const COC_STAYS_MAP = LookupCallback[
    LookupCallback(["RA_1", "response_none_1", "response_stay_1"],
                    (ra, none, stay) -> (ra, ra ? none : true, ra ? stay : false)),
    LookupCallback(["RA_2", "response_none_2", "response_stay_2"],
                    (ra, none, stay) -> (ra, ra ? none : true, ra ? stay : false))
    ]

function is_converging(psi1::Float64, chi1::Float64, psi2::Float64, chi2::Float64)
    #println("psi1=$psi1, chi1=$chi1, psi2=$psi2, chi2=$chi2")
    if abs(chi1) > pi/2 && abs(chi2) > pi/2 #flying away from each other
        return false
    end
    z1 = to_plusminus_pi(psi2 - psi1)
    z2 = to_plusminus_pi(psi1 - psi2)
    isconverge = z1 * chi1 <= 0 && z2 * chi2 <= 0
    isconverge
end

const ADD_FEATURE_MAP = LookupCallback[
    LookupCallback(["psi_1", "intr_chi_1", "psi_2", "intr_chi_2"], is_converging),
    LookupCallback(["h_1", "h_2"], (h1, h2) -> abs(h2 - h1)),
    LookupCallback(["intr_sr_1", "abs_alt_diff"], (sr, altdiff) -> sqrt(sr^2 - altdiff^2))
    ]

const ADD_FEATURE_NAMES = ASCIIString[
    "converging",
    "abs_alt_diff",
    "horizontal_range"
    ]

function csvs2dataframes(in_dir::AbstractString, out_dir::AbstractString)
    if !isdir(out_dir) #create output dir if it doesn't exist
        mkpath(out_dir)
    end
    csvfiles = readdir_ext("csv", in_dir)
    df_files = csv_to_dataframe(csvfiles, FEATURE_MAP, FEATURE_NAMES, outdir=out_dir)
    add_features!(df_files, ADD_FEATURE_MAP, ADD_FEATURE_NAMES, overwrite=true)
    transform(df_files, (:psi_1, rad2deg), (:psi_2, rad2deg), (:intr_chi_1, rad2deg), (:intr_chi_2, rad2deg),
                overwrite=true)
end

#Note: This can be improved using DataFramesMeta
function correct_coc_stays!(dir::AbstractString)
    df_files = readdir_ext(".csv.gz", dir)
    transform(df_files, COC_STAYS_MAP, overwrite=true)
end

#TODO: remove hardcodings
function get_id(filename::AbstractString)
    s = basename(filename)
    s = split(s, '.')[1]
    s = replace(s, "trajSaveMCTS_ACASX_EvE_", "")
    s = replace(s, "trajSaveMCTS_ACASX_Multi_", "")
    s = replace(s, "trajSaveMCTS_ACASX_GM_", "")
    s = replace(s, "_dataframe", "")
    s = replace(s, "_filled", "")
    encounter_id = parse(Int64, s) #encounter number
    encounter_id
end

function merge_dataframes(dir::AbstractString, outfile::AbstractString)
    Ds = get_dataframes(dir)
    D = vcat(Ds...)
    writetable(outfile, D)
end

function add_encounter_info!(dir::AbstractString)
    files = readdir_ext(".csv.gz", dir)
    for (i, f) in enumerate(files)
        fileroot = split(basename(f), '.')[1]
        df = readtable(f)
        df[:t] = Int64[1:nrow(df);]
        df[:encounter_id] = fill(get_id(fileroot), nrow(df))
        writetable(f, df)
    end
    files
end

function make_dataset(dfdir::AbstractString, jsondir::AbstractString,
    outdir::AbstractString)

    mkpath(outdir)

    dffiles = readdir(dfdir)
    sort!(dffiles, by=get_id) #FIXME: ordering is important to sync with encounter_meta
    map!(f->joinpath(dfdir, f), dffiles)

    jsonfiles = readdirGZs(jsondir)
    sort!(jsonfiles, by=get_id) #this order is important

    meta = metadf(length(jsonfiles))
    meta[:encounter_id] = map(get_id, jsonfiles)
    meta[:nmac] = map(is_nmac, jsonfiles)

    #make sure things match up
    @assert all(map(get_id, dffiles) .== map(get_id, jsonfiles))

    Ds = DFSet(meta, map(readtable, dffiles))
    save_csvs(outdir, Ds)
    Ds
end


text:
@inline hasstacks(interp::AbstractPushInterpreter, i::PushInstruction) =
    all(st -> hasstack(interp, st), neededstacks(i))

@inline inputtypes(i::PushInstruction) = DataType[]
@inline outputtypes(i::PushInstruction) = DataType[]
@inline neededstacks(i::PushInstruction) = unique(vcat(inputtypes(i), outputtypes(i)))
@inline arity(i::PushInstruction) = length(inputtypes(i))
@inline numoutputs(i::PushInstruction) = length(outputtypes(i))

name(i::PushInstruction) = name(typeof(i))
desc(i::PushInstruction) = name(i)
opname(i::PushInstruction) = opname(typeof(i))
oldopname(i::PushInstruction) = oldopname(typeof(i))
oldname(i::PushInstruction) = oldname(typeof(i))
names(i::PushInstruction) = names(typeof(i))

oldopname(::Type{I}) where {I<:PushInstruction} = uppercase(opname(I)) # Default is that the old Push name is same as current op name. Override for exception like add/+ etc.
names(::Type{I}) where {I<:PushInstruction} = unique([name(I), oldname(I)])

@inline function name(::Type{I}) where {I<:PushInstruction}
    oname, stname = split_type_name_in_parts(I)
    (stname === nothing) ? opname(I) : (lowercase(stname) * "_" * opname(I))
end
@inline function oldname(::Type{I}) where {I<:PushInstruction}
    oname, stname = split_type_name_in_parts(I)
    (stname === nothing) ? oldopname(I) : (uppercase(stname) * "." * oldopname(I))
end
@inline opname(::Type{I}) where {I<:PushInstruction} =
    lowercase(first(split_type_name_in_parts(I)))

const TypeNameRegexp = r"([^{]+){([^,]+)(.*)}"

function split_type_name_in_parts(::Type{I}) where {I<:PushInstruction}
    m = match(TypeNameRegexp, string(I))
    if m === nothing
        return string(I), nothing, nothing
    elseif length(m.captures) == 2 || m[3] == ""
        return m[1], m[2], nothing
    elseif length(m.captures) == 3
        return m[1], m[2], map(strip, filter(s -> length(s) > 0, split(m[3], ",")))
    end
end

@inline execute!(i::PushInstruction, interp::AbstractPushInterpreter) = nothing # No-Op is default

# We use a naming convention where types are indicated A, B, C and so.
# So an instruction taking 2 input arguments of type A and producing one output of the same
# type A would be of type PushInstructionA2tA1.
# An instruction taking one input of type A, one input of type B and producing an output
# of type C would be of type: PushInstructionA1B1tC1

abstract type PushInstructionA1tA1{T} <: PushInstruction end
inputtypes(i::PushInstructionA1tA1{T}) where {T} = DataType[T]
outputtypes(i::PushInstructionA1tA1{T}) where {T} = DataType[T]

abstract type PushInstructionA1tB1{I,O} <: PushInstruction end
inputtypes(i::PushInstructionA1tB1{I,O}) where {I,O} = DataType[I]
outputtypes(i::PushInstructionA1tB1{I,O}) where {I,O} = DataType[O]

abstract type PushUnaryOp{T} <: PushInstructionA1tA1{T} end

abstract type PushInstructionA2tB1{I,O} <: PushInstruction end
inputtypes(i::PushInstructionA2tB1{I,O}) where {I,O} = DataType[I, I] # Note that both are given so arity is right
outputtypes(i::PushInstructionA2tB1{I,O}) where {I,O} = [O]
neededstacks(i::PushInstructionA2tB1{I,O}) where {I,O} = DataType[I, O]

@inline function execute!(i::PushInstructionA2tB1{I,O}, interp::AbstractPushInterpreter) where {I,O}
    stackhas(interp, I, 2) || return(nothing)
    a, b = pop2!(interp, I)
    res = execute_op(i, a, b) # Should produce a result of type O (but we do not check it)
    push!(interp, res)
end

# Binary Op is an A2tB1 where input type (A) and output type (B) are the same.
abstract type PushBinaryOp{T} <: PushInstructionA2tB1{T,T} end
neededstacks(i::PushBinaryOp{T}) where {T} = DataType[T]
name(i::PushBinaryOp{T}) where {T} = lowercase(string(T)) * "_" * opname(typeof(i))

abstract type PushInstructionA0tA1{T} <: PushInstruction end
inputtypes(i::PushInstructionA0tA1{T}) where {T} = DataType[]
outputtypes(i::PushInstructionA0tA1{T}) where {T} = DataType[T]

# Since Vector are used as Push programs we might need to explicitly
# indicate when they are literals. We can do this with the Lit instruction.
struct Lit{T} <: PushInstructionA0tA1{T}
    value::T
    Lit(x::T) where {T} = new{T}(x)
end
@inline execute!(l::Lit, interp::AbstractPushInterpreter) = push!(interp, l.value)
@inline desc(i::Lit) = string(i.value)

# Load instructions for specific areas:
include("instructions/stack_instructions.jl")
include("instructions/number_instructions.jl")
include("instructions/bool_instructions.jl")
include("instructions/conversion_instructions.jl")


text:
using LinearAlgebra, CUDA

import Base: *, transpose, adjoint, permutedims, size, axes, IndexStyle
# import Base: A_mul_B!
# import Base: A_mul_Bt, A_mul_Bt!, A_mul_Bc, A_mul_Bc!
# import Base: At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!
# import Base: At_mul_Bt, At_mul_Bt!, Ac_mul_Bc, Ac_mul_Bc!
import LinearAlgebra.BLAS: gemm!, scal!
import LinearAlgebra: rmul!, lmul!, axpy!, norm
# import Base.LinAlg: scale! `scale!(a::Number, B::AbstractArray)` is deprecated, use `lmul!(a, B)` instead.
# export axpy!

# AutoGrad defines: @primitive1 *(x1,x2),dy  (dy*x2')  (x1'*dy)
# We specialize it below to avoid transposes
# Full-scale lazy transpose requires a lot more things to work with Adjoint(::KnetArray)
(*)(A::KnetMatrix{T},B::KnetMatrix{T}) where {T} = gemm!('N','N',one(T),A,B,zero(T),similar(A,(size(A,1),size(B,2))))
A_mul_Bt(A::KnetMatrix{T}, B::KnetMatrix{T}) where {T} = gemm!('N','T',one(T),A,B,zero(T),similar(A,size(A,1),size(B,1)))
At_mul_B(A::KnetMatrix{T}, B::KnetMatrix{T}) where {T} = gemm!('T','N',one(T),A,B,zero(T),similar(A,size(A,2),size(B,2)))
At_mul_Bt(A::KnetMatrix{T}, B::KnetMatrix{T}) where {T} = gemm!('T','T',one(T),A,B,zero(T),similar(A,size(A,2),size(B,1)))
@primitive1 *(x1::KnetMatrix,x2::KnetMatrix),dy  A_mul_Bt(dy,x2)  At_mul_B(x1,dy)
@primitive1 Knet.A_mul_Bt(x1::KnetMatrix,x2::KnetMatrix),dy  (dy*x2)  At_mul_B(dy,x1)
@primitive1 Knet.At_mul_B(x1::KnetMatrix,x2::KnetMatrix),dy  A_mul_Bt(x2,dy)  (x1*dy)
@primitive1 Knet.At_mul_Bt(x1::KnetMatrix,x2::KnetMatrix),dy  At_mul_Bt(x2,dy)  At_mul_Bt(dy,x1)

# Allow 1-D vectors as (N,1) in matmul:
(*)(A::KnetVector{T},B::KnetMatrix{T}) where {T} = reshape(A,:,1) * B
(*)(A::KnetMatrix{T},B::KnetVector{T}) where {T} = (C = A * reshape(B,:,1); size(A,1) == 1 ? C[1] : vec(C))

# deprecated:
# A_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','N',one(T),A,B,zero(T),C)
# A_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','T',one(T),A,B,zero(T),C)
# A_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(similar(A,(size(A,1),size(B,1))),A,B)
# A_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(C,A,B)
# A_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt(A,B)

# At_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','N',one(T),A,B,zero(T),C)
# At_mul_B{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(similar(A,(size(A,2),size(B,2))),A,B)
# Ac_mul_B!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(C,A,B)
# Ac_mul_B{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B(A,B)

# At_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','T',one(T),A,B,zero(T),C)
# At_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(similar(A,(size(A,2),size(B,1))),A,B)
# Ac_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(C,A,B)
# Ac_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt(A,B)

function gemm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T}) where {T}
    cublasop(c::Char)=CUBLAS.cublasOperation_t(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    size2(x,i)=(if ndims(x)<=2; size(x,i); elseif i==1; div(length(x),size(x,ndims(x))); elseif i==2; size(x,ndims(x)); else 1; end)
    if transA == 'N'
        m=size2(A,1); k=size2(A,2)
    else
        m=size2(A,2); k=size2(A,1)
    end
    if transB == 'N'
        n=size2(B,2); k==size2(B,1) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    else
        n=size2(B,1); k==size2(B,2) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    (m == size2(C,1) && n == size(C,2)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    lda = size2(A,1); ldb = size2(B,1); ldc = size2(C,1)
    if T<:Float64
        # @cublas(cublasDgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
        CUBLAS.cublasDgemm_v2(CUBLAS.handle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    elseif T<:Float32
        # @cublas(cublasSgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
        CUBLAS.cublasSgemm_v2(CUBLAS.handle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    # elseif T<:Float16
    #     @cublas(cublasHgemm, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

function axpy!(n::Integer, alpha::Number, x::KnetArray{T}, incx::Integer, y::KnetArray{T}, incy::Integer) where {T}
    length(x) == length(y) || throw(DimensionMismatch("$(map(size,(x,y)))"))
    alpha = T[alpha]
    if T<:Float32
        # @cublas(cublasSaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx, y, incy)
        CUBLAS.cublasSaxpy_v2(CUBLAS.handle(), n, alpha, x, incx, y, incy)
    elseif T<:Float64
        # @cublas(cublasDaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx, y, incy)
        CUBLAS.cublasDaxpy_v2(CUBLAS.handle(), n, alpha, x, incx, y, incy)
    else
        error("$T not supported")
    end
    return y
end

axpy!(alpha::Number, x::KnetArray{T}, y::KnetArray{T}) where {T} = axpy!(length(x),alpha,x,1,y,1)

function scal!(n::Integer, alpha::Number, x::KnetArray{T}, incx::Integer) where {T}
    alpha = T[alpha]
    if T<:Float32
        # @cublas(cublasSscal_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx)
        CUBLAS.cublasSscal_v2(CUBLAS.handle(), n, alpha, x, incx)
    elseif T<:Float64
        # @cublas(cublasDscal_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx)
        CUBLAS.cublasDscal_v2(CUBLAS.handle(), n, alpha, x, incx)
    else
        error("$T not supported")
    end
    return x
end

lmul!(alpha::Number, x::KnetArray{T}) where {T} = scal!(length(x),alpha,x,1)
rmul!(x::KnetArray{T}, alpha::Number) where {T} = scal!(length(x),alpha,x,1)

transpose(x::KnetVecOrMat)=_transpose(x)
adjoint(x::KnetVecOrMat)=_transpose(x)
_transpose(x::KnetVector) = copy(reshape(x,1,:))
_transpose(x::KnetMatrix) = _transpose!(similar(x,(size(x,2),size(x,1))),x)

function _transpose!(y::KnetMatrix{T}, x::KnetMatrix{T}) where {T}
    if T<:Float32
        # @cublas(cublasSgeam, (Cptr,UInt32,UInt32,Cint,Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Cint),
        #         cublashandle(),1,1,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
        op = CUBLAS.cublasOperation_t(1)
        CUBLAS.cublasSgeam(CUBLAS.handle(),op,op,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
    elseif T<:Float64
        # @cublas(cublasDgeam, (Cptr,UInt32,UInt32,Cint,Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Cint),
        #         cublashandle(),1,1,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
        op = CUBLAS.cublasOperation_t(1)
        CUBLAS.cublasDgeam(CUBLAS.handle(),op,op,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
    else
        error("CUBLAS does not support $T")
    end
    return y
end

#= TODO: use the lazy transpose:
using LinearAlgebra: Adjoint, Transpose, AdjOrTrans
transpose(x::KnetArray)=Transpose(x)
adjoint(x::KnetArray)=Adjoint(x)
const AdjointKnetVec{T} = Adjoint{T,<:KnetVector}
const TransposeKnetVec{T} = Transpose{T,<:KnetVector}
const AdjOrTransKnetVec{T} = AdjOrTrans{T,<:KnetVector}
const AdjOrTransKnetMat{T} = AdjOrTrans{T,<:KnetMatrix}
size(v::AdjOrTransKnetVec) = (1, length(v.parent))
size(A::AdjOrTransKnetMat) = reverse(size(A.parent))
axes(v::AdjOrTransKnetVec) = (Base.OneTo(1), axes(v.parent)...)
axes(A::AdjOrTransKnetMat) = reverse(axes(A.parent))
IndexStyle(::Type{<:AdjOrTransKnetVec}) = IndexLinear()
IndexStyle(::Type{<:AdjOrTransKnetMat}) = IndexCartesian()
=#

"""

    mat(x; dims = ndims(x) - 1)

Reshape `x` into a two-dimensional matrix by joining the first dims dimensions, i.e.
`reshape(x, prod(size(x,i) for i in 1:dims), :)`

`dims=ndims(x)-1` (default) is typically used when turning the output of a 4-D convolution
result into a 2-D input for a fully connected layer.

`dims=1` is typically used when turning the 3-D output of an RNN layer into a 2-D input for
a fully connected layer.

`dims=0` will turn the input into a row vector, `dims=ndims(x)` will turn it into a column
vector.

"""
mat(x; dims::Int=ndims(x)-1)=reshape(x, (dims > 0 ? prod(size(x,i) for i in 1:dims) : 1), :)

# conv: reshape(x, (:,xn)): rowdims=ndims-1
# rnns: reshape(x, (x1,:)): rowdims=1
# general: reshape(x, (x1*x2..xi, x[i+1]*...*xn))
# specify the first rowdims are joined, the remaining are joined
# 1-D input can be turned into a rowvec (rowdims=0), or colvec (rowdims=1).
# default dims=ndims(x)-1 will turn vec into a rowvec but dims=1 will not work for conv.

# Low level gemm! call with pointers: CPU conv4 uses this. Based on julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1105

using LinearAlgebra
using LinearAlgebra.BLAS: libblas, BlasInt
using LinearAlgebra.BLAS: @blasfunc

# C := alpha*op(A)*op(B) + beta*C, where:
# op(X) is one of op(X) = X, or op(X) = XT, or op(X) = XH,
# alpha and beta are scalars,
# A, B and C are matrices:
# op(A) is an m-by-k matrix,
# op(B) is a k-by-n matrix,
# C is an m-by-n matrix.

for (gemm, elty) in ((:dgemm_,:Float64), (:sgemm_,:Float32))
    @eval begin
        function gemm!(transA::AbstractChar, transB::AbstractChar, M::Int, N::Int, K::Int, alpha::($elty), A::Ptr{$elty}, B::Ptr{$elty}, beta::($elty), C::Ptr{$elty})
            if transA=='N'; lda=M; else; lda=K; end
            if transB=='N'; ldb=K; else; ldb=N; end
            ldc = M;
            # ccall((@blasfunc($gemm), libblas), Nothing,
            #       (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt},
            #        Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
            #        Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
            #        Ptr{BlasInt}),
            #       &transA, &transB, &M, &N, &K,
            #       &alpha, A, &lda, B, &ldb, &beta, C, &ldc)

            ccall((@blasfunc($gemm), libblas), Cvoid,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                   Ref{BlasInt}),
                  transA, transB, M, N, K,
                  alpha, A, lda,
                  B, ldb,
                  beta, C, ldc)
        end
    end
end

function norm(x::KnetArray{T}, p::Real=2) where {T}
    if length(x) == 0
        zero(T)
    elseif p == 2
        sqrt(sum(abs2,x))
    elseif p == 1
        sum(abs,x)
    elseif p == Inf
        maximum(abs,x)
    elseif p == 0
        countnz(x)
    elseif p == -Inf
        minimum(abs,x)
    else
        sum(abs.(x).^p)^(1/p)
    end
end


text:
# Convenience constructor for creating StateSpace objects

"""`sys = ss(A, B, C, D, Ts=0) -> sys`

Create a state-space model `sys::StateSpace{T, MT<:AbstractMatrix{T}}`
where `MT` is the type of matrixes `A,B,C,D` and `T` the element type.

This is a continuous-time model if Ts is omitted or set to 0.
Otherwise, this is a discrete-time model with sampling period Ts.
Set Ts=-1 for a discrete-time model with unspecified sampling period.

`sys = ss(D[, Ts, ...])` specifies a static gain matrix D."""
function ss(A::AbstractArray, B::AbstractArray, C::AbstractArray, D::AbstractArray, Ts::Real=0)
    # Check the kwargs for metadata
    nu = size(B, 2)
    ny, nx = size(C, 1), size(C, 2)
    return StateSpace(A, B, C, D, Ts)
end

# Function for accepting scalars
function ss(A::Union{Number,AbstractArray}, B::Union{Number,AbstractArray}, C::Union{Number,AbstractArray}, D::Union{Number,AbstractArray}, Ts::Real=0)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    A = to_matrix(T, A)
    B = to_matrix(T, B)
    C = to_matrix(T, C)
    if D == 0
        D = fill(zero(T), size(C,1), size(B,2))
    else
        D = to_matrix(T, D)
    end
    ss(A, B, C, D, Ts)
end

# Function for creation of static gain
function ss(D::AbstractArray{T}, Ts::Real=0) where {T<:Number}
    ny, nu = size(D, 1), size(D, 2)
    A = fill(zero(T), 0, 0)
    B = fill(zero(T), 0, nu)
    C = fill(zero(T), ny, 0)

    return ss(A, B, C, D, Ts)
end
ss(d::Number, Ts::Real=0; kwargs...) = ss([d], Ts)

# ss(sys) converts to StateSpace
ss(sys::LTISystem) = convert(StateSpace, sys)


text:
# Weighted histogram analysis method.

# Kumar, S., Rosenberg, J. M., Bouzida, D., Swendsen, R. H., & Kollman, P. A.
# (1992). The weighted histogram analysis method for free‐energy calculations
# on biomolecules. I. The method. Journal of Computational Chemistry, 13(8),
# 1011-1021.

function unbias_histograms(T::Float64, bias_force_constant::Float64,
                           bin_size::Float64,
                           trajectories::Dict{Float64,Vector{Float64}},
                           num_iter::Int; aux_output::IO=stderr)
    # Most extreme positions of restrained coordinate.
    min_xi = Inf
    max_xi = -Inf

    for trajectory in values(trajectories)
        min_xi = min(min_xi, minimum(trajectory))
        max_xi = max(max_xi, maximum(trajectory))
    end

    # Round to one decimal place (0.1 nm).
    min_xi = round(min_xi, RoundDown; digits=1)
    max_xi = round(max_xi, RoundUp; digits=1)

    bins = min_xi:bin_size:max_xi
    num_edges = length(bins)+1
    edges = range(min_xi - bin_size / 2, max_xi + bin_size / 2;
                  length=num_edges)

    total = zeros(Float64, size(bins))
    windows = Float64[]
    sizes = Float64[]

    for (window, trajectory) in sort(collect(trajectories))
        _, autocorrelation_steps = stderr_bin(1.0, trajectory)
        inefficiency = 1 + 2 * autocorrelation_steps
        hist = fit(Histogram, trajectory, edges).weights / inefficiency
        total .+= hist
        push!(windows, window)
        push!(sizes, sum(hist))
    end

    offsets = bins .- windows'
    exppot = exp.(-beta(T) .* (0.5 .* bias_force_constant .* offsets.^2))

    if size(exppot) != (length(total), length(sizes))
        error("Mismatched dimensions.")
    end

    # Apply WHAM equations iteratively.
    prob = nothing
    weights = ones(Float64, size(sizes))

    meter = Progress(num_iter, output=aux_output)

    for _ in ProgressWrapper(1:num_iter, meter)
        factors = dropdims(sum(sizes' .* exppot ./ weights'; dims=2); dims=2)
        prob = total ./ factors
        # Normalize.
        prob ./= sum(prob)
        weights = dropdims(sum(prob .* exppot; dims=1); dims=1)
    end

    pmf = -log.(prob) ./ beta(T)

    bins, pmf
end


text:
using Gnuplot
x = LinRange(0,1,100)
@gp x x.^2 "w filledcu y = 0 lc 'red' t 'area fill'"
@gp :- "set style fill transparent solid 0.2"
save(term ="pngcairo size 600,400", output ="areafill.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
const col_types = Dict(
    :date => String,
    :name => String,
    :role => String,
    :arena_place => Union{Missing, Float64},
    :arena_victories => Union{Missing, Float64},
    :arena_defences => Union{Missing, Float64},
    :total_power => Union{Missing, Float64},
    :maxed_heroes => Union{Missing, Float64},
    :clan_tokens => Union{Missing, Float64},
    :raid_points => Union{Missing, Float64},
    :portal_stones => Union{Missing, Float64},
    :Artus => Union{Missing, Float64},
    :Trorin => Union{Missing, Float64}
    )

function append_stats(clan::AbstractDict)
    dir = joinpath(stats_path,clan[:name])
    files = joinpath.(dir,readdir(dir))
    filter!(t -> !occursin.("DS_Store",t),files)
    full = DataFrame()
    for f in files
        df = read(f;types = col_types, normalizenames=true) |> DataFrame
        df.date = [Date(x,"dd-mm-yyyy") for x in df.date]
        if isempty(full)
            full = df
        else
            append!(full,df)
        end
    end
    filter!(r -> !in(r[:name],clan[:past_members]),full)
    if :Artus in names(full)
        select!(full,Not([:Artus,:Trorin]))
    end
    return full
end

## add missing values for dates were new members have no data
function complete_stats(full)
    sort!(full,(:date,:name))
    df = combine(groupby(full,:name)) do dd
        days = union(full.date)
        for d in days
            if !in(d,dd.date)
                println("missing data for $(dd.name[1]) day $(dd.date[1])")
                push!(full,[d,dd.name[1],dd.role[end], missings(ncol(full)-3)...])
            end
        end
    end
    sort!(full,(:date,:name))
end

function get_stats(clan::AbstractDict)
    df = append_stats(clan)
    complete_stats(df)
end


text:
#############################################
########### Output Functions ################
#############################################
function generate_output_data(CRS, savedir)
    unique_name = savedir*"/"*string(hash(CRS.parameters))*Random.randstring()*".csv"
    parameter_df = DataFrame()
    for key in keys(CRS.parameters)
        parameter_df[key] = CRS.parameters[key]
    end
    parameter_df[:save_name] = unique_name
    return unique_name, parameter_df
end

function InitializeOutput(CRS)
    outDF = DataFrame()
    ### Unclear if this always alines correctly. Could use join(..., by= ) to ensure that
    outDF[Symbol("mID")] = 1:length(CRS.molecule_list)
    outDF[Symbol("molecule")] = CRS.molecule_list
    return outDF
end

function df2json(df::DataFrame)
  len = length(df[:,1])
  indices = names(df)
  jsonarray = [Dict([string(index) => df[index][i] for index in indices])
               for i in 1:len]
  return JSON.json(jsonarray)
end

function writejson(path::String,df::DataFrame)
  open(path,"w") do f
    write(f,df2json(df))
  end
end

max_stars_repo_name: croningp/clusterautocatalyticsets

max_issues_repo_name: croningp/clusterautocatalyticsets

max_forks_repo_name: croningp/clusterautocatalyticsets