text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/25/2014

# Lincoln Laboratory Airspace Encounter Model

module LLAEMImpl

export
    AddObserver,

    update,

    generateEncounter,
    getInitialState,
    getNextCommand,
    getTrajectory,

    LLAEM,
    LLAEMInitialState,
    LLAEMCommand

import Compat.ASCIIString

using AbstractEncounterModelImpl
using CommonInterfaces
using RLESUtils, Observers

import CommonInterfaces.addObserver
import CommonInterfaces.update
import AbstractEncounterModelInterfaces.generateEncounter
import AbstractEncounterModelInterfaces.getInitialState
import AbstractEncounterModelInterfaces.getNextCommand
import AbstractEncounterModelInterfaces.getTrajectory

type LLAEMInitialState

    t::Float64
    x::Float64
    y::Float64
    h::Float64
    v::Float64
    psi::Float64
    h_d::Float64
end

type LLAEMCommand

    t::Float64
    v_d::Float64
    h_d::Float64
    psi_d::Float64
end

type LLAEM <: AbstractEncounterModel

    initial_sample_filename::ASCIIString
    number_of_initial_samples::Int
    f_init::Union{IOStream, Void}

    transition_sample_filenames::Vector{ASCIIString}
    number_of_transition_samples::Vector{Int}
    f_tran::Vector{Union{IOStream, Void}}

    number_of_aircraft::Int

    A::Int
    L::Int
    C::Vector{Int}
    v_init::Vector{Float64}

    geometry_at_TCA::Vector{Float64}

    states::Array{Float64, 3}
    state_index::Vector{Int}

    dynamic_states::Array{Float64, 3}
    dn_state_index::Vector{Int}

    observer::Observer

    function LLAEM{T <: AbstractString}(
                    number_of_aircraft::Int,
                    initial_sample_filename::AbstractString,
                    number_of_initial_samples::Int,
                    transition_sample_filenames::Vector{T},
                    number_of_transition_samples::Vector{Int})

        obj = new()

        obj.initial_sample_filename = initial_sample_filename
        obj.number_of_initial_samples = number_of_initial_samples
        obj.f_init = nothing

        obj.transition_sample_filenames = transition_sample_filenames
        obj.number_of_transition_samples = number_of_transition_samples
        obj.f_tran = [nothing, nothing]

        obj.number_of_aircraft = number_of_aircraft

        obj.C = zeros(obj.number_of_aircraft)
        obj.v_init = zeros(obj.number_of_aircraft)

        obj.state_index = zeros(obj.number_of_aircraft)
        obj.states = zeros(obj.number_of_aircraft, maximum(obj.number_of_transition_samples), 4)

        obj.dn_state_index = zeros(obj.number_of_aircraft)
        obj.dynamic_states = zeros(obj.number_of_aircraft, maximum(obj.number_of_transition_samples) + 1, 6)

        obj.observer = Observer()

        return obj
    end
end

addObserver(aem::LLAEM, f::Function) = add_observer(aem.observer, f)
addObserver(aem::LLAEM, tag::AbstractString, f::Function) = add_observer(aem.observer, tag, f)

function generateEncounter(aem::LLAEM, sample_number::Int)

    if sample_number > 0
        reset_sample_from_file(aem)

        for i = 1:sample_number
            initial, transitions = read_sample_from_file(aem, aem.number_of_initial_samples, aem.number_of_transition_samples)

            if initial == nothing || transitions == nothing
                return -1
            end
        end
    else
        initial, transitions = read_sample_from_file(aem, aem.number_of_initial_samples, aem.number_of_transition_samples)

        if initial == nothing || transitions == nothing
            return -1
        end
    end

    # A, L, chi, beta, C1, C2, v1, v2, v1d, v2d, h1d, h2d, psi1d, psi2d, hmd, vmd, h
    convert_units_for_initial(initial)

    # t, h_d, psi_d, v_d, v, x, y, h, psi
    convert_units_for_transitions(transitions)

    A, L, chi, beta_, C1, C2, v1, v2, v1d, v2d, h1d, h2d, psi1d, psi2d, hmd, vmd, h = initial

    aem.A = A
    aem.L = L
    aem.C[1] = C1
    aem.C[2] = C2
    aem.v_init[1] = v1
    aem.v_init[2] = v2

    aem.geometry_at_TCA = [chi, beta_, hmd, vmd, h]

    for i = 1:aem.number_of_aircraft
        aem.states[i, :, 1] = transitions[i, 1:end-1, 1]  # time
        aem.states[i, :, 2] = transitions[i, 1:end-1, 4]  # v_d
        aem.states[i, :, 3] = transitions[i, 1:end-1, 2]  # h_d
        aem.states[i, :, 4] = transitions[i, 1:end-1, 3]  # psi_d
        aem.state_index[i] = 0
    end

    for i = 1:aem.number_of_aircraft
        aem.dn_state_index[i] = 0

        for j = 1:(maximum(aem.number_of_transition_samples) + 1)
            t, h_d, psi_d, v_d, v, x, y, h, psi = transitions[i, j, :]

            aem.dn_state_index[i] += 1
            aem.dynamic_states[i, j, :] = [t, x, y, h, v, psi]
        end
    end

    #for i = 1:(maximum(aem.number_of_transition_samples) + 1)
    #    print(reshape(aem.dynamic_states[1, i, :], 6)')
    #end

    #for i = 1:maximum(aem.number_of_transition_samples)
    #    print(reshape(aem.states[1, i, :], 4)')
    #end
end

generateEncounter(aem::LLAEM) = generateEncounter(aem, 0)

function getInitialState(aem::LLAEM, aircraft_number::Int)

    state = aem.dynamic_states[aircraft_number, 1, :]   # t, x, y, h, v, psi
    command = aem.states[aircraft_number, 1, :]         # t, v_d, h_d, psi_d

    return LLAEMInitialState(state[1], state[2], state[3], state[4], state[5], state[6], command[3])
end

getNextCommand(aem::LLAEM, aircraft_number::Int) = get_next_state(aem, aircraft_number)

function update(aem::LLAEM, aircraft_number::Int)

    command = getNextCommand(aem, aircraft_number)

    if command != nothing
        return LLAEMCommand(command...)
    else
        return nothing
    end
end

function getTrajectory(aem::LLAEM, aircraft_number::Int)

    return reshape(aem.dynamic_states[aircraft_number, :, 1:4], aem.dn_state_index[aircraft_number], 4)
end

function convert_units_for_initial(values)

    # A, L, chi, beta, C1, C2, v1, v2, v1d, v2d, h1d, h2d, psi1d, psi2d, hmd, vmd, h
    values[7] *= 1.68780
    values[8] *= 1.68780
    values[9] *= 1.68780
    values[10] *= 1.68780
    values[11] /= 60
    values[12] /= 60
    values[15] *= 6076.12

    return values
end

function convert_units_for_transitions(values)

    # t, h_d, psi_d, v_d, v, x, y, h, psi
    values[:, :, 2] /= 60
    values[:, :, 4] *= 1.68780
    values[:, :, 5] *= 1.68780

    return values
end

function get_next_state(aem, aircraft_number)

    aem.state_index[aircraft_number] += 1

    if aem.state_index[aircraft_number] > aem.number_of_transition_samples[aircraft_number]
        return nothing
    else
        return reshape(aem.states[aircraft_number, aem.state_index[aircraft_number], :], 4)
    end
end

# TODO support to read number of samples
# returns a sample at a time, currently

function read_sample_from_file(aem, number_of_initial_samples, number_of_transition_samples)

    # 1       id
    # 2       A = [1:4]
    # 3       L = [1:5]
    # 4       chi = [90, 270]   : deg
    # 5       beta = [0:30:360] : deg
    # 6       C1 = [1, 2]
    # 7       C2 = [1, 2]
    # 8,9     v1, v2 = [50, 100, 200, 300, 400, 500, 600]   : kt
    # 10, 11  v1d, v2d = [-5, -2, -0.25, 0.25, 2, 5]        : kt/s
    # 12, 13  h1d, h2d = [-5000, -3000, -2000, -1000, -400, 400, 1000, 2000, 3000, 5000]    : ft/min
    # 14, 15  psi1d, psi2d = [-8, -6, -3.5, -1, -0.25, 0.25, 1, 3.5, 6, 8]  : deg/s
    # 16      hmd = [0, 0.0822896, 0.5, 1, 3]   : nm
    # 17      vmd = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 6000]  : ft
    # 18      h = altitude  : ft

    if aem.f_init == nothing
        aem.f_init = open(aem.initial_sample_filename, "r")
    end

    line = readline(aem.f_init)

    if line == ""
        close(aem.f_init)
        aem.f_init = nothing

        initial = nothing
    else
        initial = split(chomp(line))

        A_ = copy(initial[2])

        if A_ == "b"
            A = 1
        elseif A_ == "c"
            A = 2
        elseif A_ == "d"
            A = 3
        elseif A_ == "o"
            A = 4
        end

        initial = [A, float(initial[3:end])]
    end

    # 1     id
    # 2     time step
    # 3     h1d         : ft/min
    # 4     psi_d       : deg/s
    # 5     v_d         : kt/s
    # 6     v           : kt
    # 7     x           : ft
    # 8     y           : ft
    # 9     h           : ft
    # 10    psi         : deg

    transitions = nothing

    for i = 1:aem.number_of_aircraft
        if aem.f_tran[i] == nothing
            aem.f_tran[i] = open(aem.transition_sample_filenames[i], "r")
        end

        line = readline(aem.f_tran[i])

        if line == ""
            close(aem.f_tran[i])
            aem.f_tran[i] = nothing
        else
            values = float(split(chomp(line)))

            if transitions == nothing
                transitions = Array(Float64, aem.number_of_aircraft, maximum(aem.number_of_transition_samples) + 1, length(values) - 1)
            end

            transitions[i, 1, :] = values[2:end]

            for j = 2:(aem.number_of_transition_samples[i] + 1)
                values = float(split(chomp(readline(aem.f_tran[i]))))

                transitions[i, j, :] = values[2:end]
            end
        end
    end

    return initial, transitions
end

function reset_sample_from_file(aem)

    if aem.f_init != nothing
        close(aem.f_init)
        aem.f_init = nothing
    end

    for i = 1:aem.number_of_aircraft
        if aem.f_tran[i] != nothing
            close(aem.f_tran[i])
            aem.f_tran[i] = nothing
        end
    end
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
# # Custom ReLU layer

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/custom-relu.jl)

# We demonstrate how DiffOpt can be used to generate a simple neural network
# unit - the ReLU layer. A neural network is created using Flux.jl and
# trained on the MNIST dataset.

# This tutorial uses the following packages

using JuMP
import DiffOpt
import Ipopt
import ChainRulesCore
import Flux
import Statistics
import Base.Iterators: repeated
using LinearAlgebra

# ## The ReLU and its derivative

# Define a relu through an optimization problem solved by a quadratic solver.
# Return the solution of the problem.
function matrix_relu(
    y::AbstractArray{T};
    model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
) where T
    _x = zeros(size(y))
    N = length(y[:, 1])
    empty!(model)
    set_silent(model)
    @variable(model, x[1:N] >= 0)
    for i in 1:size(y, 2)
        @objective(
            model,
            Min,
            dot(x, x) -2dot(y[:, i], x)
        )
        optimize!(model)
        _x[:, i] = value.(x)
    end
    return _x
end

# Define the backward differentiation rule, for the function we defined above.
function ChainRulesCore.rrule(
    ::typeof(matrix_relu),
    y::AbstractArray{T};
    model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
) where T
    pv = matrix_relu(y, model = model)
    function pullback_matrix_relu(dl_dx)
        ## some value from the backpropagation (e.g., loss) is denoted by `l`
        ## so `dl_dy` is the derivative of `l` wrt `y`
        x = model[:x] ## load decision variable `x` into scope
        dl_dy = zeros(T, size(dl_dx))
        dl_dq = zeros(T, size(dl_dx)) ## for step-by-step explanation
        for i in 1:size(y, 2)
            ## set sensitivities
            MOI.set.(
                model,
                DiffOpt.BackwardInVariablePrimal(),
                x,
                dl_dx[:, i]
            )
            ## compute grad
            DiffOpt.backward(model)
            ## return gradient wrt objective function parameters
            obj_exp = MOI.get(
                model,
                DiffOpt.BackwardOutObjective()
            )
            dl_dq[:, i] = JuMP.coefficient.(obj_exp, x) ## coeff of `x` in q'x = -2y'x
            dq_dy = -2 ## dq/dy = -2
            dl_dy[:, i] = dl_dq[:, i] * dq_dy
        end
        return (ChainRulesCore.NoTangent(), dl_dy,)
    end
    return pv, pullback_matrix_relu
end

# For more details about backpropagation, visit [Introduction, ChainRulesCore.jl](https://juliadiff.org/ChainRulesCore.jl/dev/).
# ## prepare data
import MLDatasets
N = 1000
imgs = MLDatasets.MNIST.traintensor(1:N)
labels = MLDatasets.MNIST.trainlabels(1:N);

# Preprocessing
train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), N)) ## stack all the images
train_Y = Flux.onehotbatch(labels, 0:9);

test_imgs = MLDatasets.MNIST.testtensor(1:N)
test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), N))
test_Y = Flux.onehotbatch(MLDatasets.MNIST.testlabels(1:N), 0:9);

# ## Define the Network

# Network structure

inner = 15

m = Flux.Chain(
    Flux.Dense(784, inner), #784 being image linear dimension (28 x 28)
    matrix_relu,
    Flux.Dense(inner, 10), # 10 beinf the number of outcomes (0 to 9)
    Flux.softmax,
)

# Define input data
# The original data is repeated `epochs` times because `Flux.train!` only
# loops through the data set once

epochs = 5

dataset = repeated((train_X, train_Y), epochs);

# Parameters for the network training

# training loss function, Flux optimizer
custom_loss(x, y) = Flux.crossentropy(m(x), y)
opt = Flux.ADAM()
evalcb = () -> @show(custom_loss(train_X, train_Y))

# Train to optimize network parameters

@time Flux.train!(custom_loss, Flux.params(m), dataset, opt, cb = Flux.throttle(evalcb, 5));

# Although our custom implementation takes time, it is able to reach similar
# accuracy as the usual ReLU function implementation.

# Average of correct guesses
accuracy(x, y) = Statistics.mean(Flux.onecold(m(x)) .== Flux.onecold(y));

# Training accuracy

accuracy(train_X, train_Y)

# Test accuracy

accuracy(test_X, test_Y)

# Note that the accuracy is low due to simplified training.
# It is possible to increase the number of samples `N`,
# the number of epochs `epoch` and the connectivity `inner`.

ext: jl


text:
"""
    PerformanceOpt

Abstract supertype for performance optimisation flags.
"""
abstract type PerformanceOpt end

"""
    SparseOpt <: PerformanceOpt

    SparseOpt()

An optimisation flag that ignores all zero values in the grid.

For low-density simulations performance may improve by
orders of magnitude, as only used cells are run.

This is complicated for optimising neighborhoods - they
must run if they contain just one non-zero cell.

Specifiy with:

```julia
ruleset = Ruleset(rule; opt=SparseOpt())
# or
output = sim!(output, rule; opt=SparseOpt())
```

`SparseOpt` is best demonstrated with this simulation, where the grey areas do not
run except where the neighborhood partially hangs over an area that is not grey:

![SparseOpt demonstration](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/complexlife_spareseopt.gif)
"""
struct SparseOpt{F<:Function} <: PerformanceOpt
    f::F
end
SparseOpt() = SparseOpt(==(0))

@inline can_skip(opt::SparseOpt{<:Function}, val) = opt.f(val)

"""
    NoOpt <: PerformanceOpt

    NoOpt()

Flag to run a simulation without performance optimisations besides basic high performance
programming. Still fast, but not intelligent about the work that it does: all cells are run
for all rules.

`NoOpt` is the default `opt` method.
"""
struct NoOpt <: PerformanceOpt end

"""
    Processor

Abstract supertype for selecting a hardware processor, such as ia CPU or GPU.
"""
abstract type Processor end

"""
    CPU <: Processor

Abstract supertype for CPU processors.
"""
abstract type CPU <: Processor end

"""
    SingleCPU <: CPU

    SingleCPU()

[`Processor`](@ref) flag that specifies to use a single thread on a single CPU.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=SingleCPU())
# or
output = sim!(output, rule; proc=SingleCPU())
```
"""
struct SingleCPU <: CPU end

"""
    ThreadedCPU <: CPU

    ThreadedCPU()

[`Processor`](@ref) flag that specifies to use a `Threads.nthreads()` CPUs.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=ThreadedCPU())
# or
output = sim!(output, rule; proc=ThreadedCPU())
```
"""
struct ThreadedCPU{L} <: CPU
    spinlock::L
end
ThreadedCPU() = ThreadedCPU(Base.Threads.SpinLock())
Base.Threads.lock(opt::ThreadedCPU) = lock(opt.spinlock)
Base.Threads.unlock(opt::ThreadedCPU) = unlock(opt.spinlock)

"""
    BoundaryCondition

Abstract supertype for flags that specify the boundary conditions used in the simulation,
used in [`inbounds`](@ref) and to update [`NeighborhoodRule`](@ref) grid padding.
These determine what happens when a neighborhood or jump extends outside of the grid.
"""
abstract type BoundaryCondition end

"""
    Wrap <: BoundaryCondition

    Wrap()

[`BoundaryCondition`](@ref) flag to wrap cordinates that boundary boundaries back to the
opposite side of the grid.

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Wrap())
# or
output = sim!(output, rule; boundary=Wrap())
```
"""
struct Wrap <: BoundaryCondition end

"""
    Remove <: BoundaryCondition

    Remove()

[`BoundaryCondition`](@ref) flag that specifies to assign `padval` to cells that overflow
grid boundaries. `padval` defaults to `zero(eltype(grid))` but can be assigned as a keyword
argument to an [`Output`](@ref).

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Remove())
# or
output = sim!(output, rule; boundary=Remove())
```
"""
struct Remove <: BoundaryCondition end

"""
    ParameterSource

Abstract supertypes for parameter source wrappers. These allow
parameters to be retreived from auxilliary data or from other grids.
"""
abstract type ParameterSource end

"""
    Aux <: ParameterSource

    Aux{K}()
    Aux(K::Symbol)

Use auxilary array with key `K` as a parameter source.

Implemented in rules with:

```julia
get(data, rule.myparam, index...)
```

When an `Aux` param is specified at rule construction with:

```julia
rule = SomeRule(; myparam=Aux{:myaux})
output = ArrayOutput(init; aux=(myaux=myauxarray,))
```

If the array is a DimensionalData.jl `DimArray` with a `Ti` (time)
dimension, the correct interval will be selected automatically,
precalculated for each timestep so it has no significant overhead.
"""
struct Aux{K} end
Aux(key::Symbol) = Aux{key}()

_unwrap(::Aux{X}) where X = X
_unwrap(::Type{<:Aux{X}}) where X = X

"""
    Grid <: ParameterSource

    Grid{K}()
    Grid(K::Symbol)

Use grid with key `K` as a parameter source.

Implemented in rules with:

```julia
get(data, rule.myparam, index...)
```

And specified at rule construction with:

```julia
SomeRule(; myparam=Grid{:somegrid})
```
"""
struct Grid{K} end
Grid(key::Symbol) = Grid{key}()

_unwrap(::Grid{X}) where X = X
_unwrap(::Type{<:Grid{X}}) where X = X


text:
__precompile__()

module Colors

using FixedPointNumbers, ColorTypes, Reexport, Compat
@reexport using ColorTypes

typealias AbstractGray{T} Color{T,1}
typealias Color3{T} Color{T,3}

import Base: ==, +, -, *, /
import Base: convert, eltype, hex, isless, linspace, show, typemin, typemax

# Additional exports, not exported by ColorTypes
export weighted_color_mean,
       hex, @colorant_str,
       protanopic, deuteranopic, tritanopic,
       distinguishable_colors, whitebalance,
       colordiff, DE_2000, DE_94, DE_JPC79, DE_CMC, DE_BFD, DE_AB, DE_DIN99, DE_DIN99d, DE_DIN99o,
       MSC, sequential_palette, diverging_palette, colormap,
       colormatch, CIE1931_CMF, CIE1964_CMF, CIE1931J_CMF, CIE1931JV_CMF

# Early utilities
include("utilities.jl")

# Include other module components
include("conversions.jl")
include("algorithms.jl")
include("parse.jl")
include("differences.jl")
include("colormaps.jl")
include("display.jl")
include("colormatch.jl")

@deprecate rgba RGBA
@deprecate hsva HSVA
@deprecate hsla HSLA
@deprecate xyza XYZA
@deprecate xyYa xyYA
@deprecate laba LabA
@deprecate luva LuvA
@deprecate lchaba LCHabA
@deprecate lchuva LCHuvA
@deprecate din99a DIN99A
@deprecate din99da DIN99dA
@deprecate din99oa DIN99oA
@deprecate lmsa LMSA
@deprecate argb32 ARGB32

end # module


text:
# ## Loading modules
# using Distributed
# using DataFrames
# using CSV
#
# addprocs(6)
#
# ## Load modules in parallel
# @everywhere include("RP_Cournot/modules/MyModules.jl")
#
# @everywhere begin
#     using .MyModules
#     using Random
#     using SharedArrays
#     ## Setting seed
#     Random.seed!(66636)
# end
#
# ## Settin up directories
#
# dir = @__DIR__ ;
# datadir=dir*"/data/oildata.csv" ;
# resultsdir=dir*"/results/"
#
# ## Loading data
#
# pt, qt = load_data_select(datadir)

## Initializing array for results
table=SharedArray{Float64,2}(zeros(36,3))

## Testing for each year for the

@time MyModules.shared_fill_table!(table, pt, qt)

## Saving results to csv

annual=DataFrame(Year= collect(1973:2008),
                 TS=table[:,1],
                 pval=table[:,2],
                 stars=mystars.(table[:,2]))
# CSV.write(resultsdir*"annual.csv", annual)
CSV.write(resultsdir*"annual_sub_$set.csv", annual)

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
testing_units = Int == Int64

include("colortypes.jl")
include("core.jl")
include("map.jl")
include("overlays.jl")
include("io.jl")
if testing_units
    include("readnrrd.jl")
end
include("readremote.jl")
include("readstream.jl")
@osx_only include("readOSX.jl")
include("algorithms.jl")
include("edge.jl")
include("parallel.jl")