text:
using JuMP, Ipopt, Plots, LinearAlgebra, SparseArrays

#A JuMP model is initialized
mod = Model(optimizer_with_attributes( Ipopt.Optimizer ,  "max_iter" => 100,
            "mumps_mem_percent" => 500))

#The parameters of the 1-D rod
L = 0.1                                 #Length of rod
λ = 45.0                                # Conductivity
c = 460.0                               # Specific heat capacitivity
ρ = 7800.0                              # Material density
α = λ / (c*ρ)                           # Diffusivity

#Number of spatial discretization points and length of each
N = 11                                 #Discretization points / nodes
Δx = L / (N-1)                         # Δx = x[i+1] - x[i]

#Integration from t0 = 0 to tf
tf = 10000.0;                          #Final time
Δt = 10^(-1);                          #Sampling period
k = round(Int,tf / Δt);                #k = number of steps

#System's matrix: θ' = Aθ + Bu
A = spdiagm(-1 => ones( N-1), 0 => -2*ones(N), 1 => ones(N-1));
A[2,1] = 2.0;
A[end - 1,end] = 2.0;
A = (α/( Δx^2)) * A;

#input matrix
B = spzeros(N);
B[1] = 1;
B = (2/ (c * ρ * Δx)) * B;

#output matrix
C = spzeros(1,N)
C[1,end] = 1.0;

#Initial and reference temperatures in kelvin
θinit = 273.0
θref = 500.0
#bounded input values[u_min,u_max]
u_min = 0.0
u_max = 15000.0

#Decision variables of the optimization problem
@variable(mod,θ[1:N,1:k]);                              # temperature at x[1:N] at time step k
@variable(mod,u_min <= u[1:k - 1] <= u_max);            # bounded constrained Input u
@variable(mod, y[1:k]);                                 # Output y
#Initial values
@constraint(mod, θ[:,1] .== θinit)
@constraint(mod, y[1] .== θinit)
#System dynamics
for j in 1:k - 1,i in 1:N
    @constraint(mod,θ[i,j + 1] - (θ[i,j] + Δt * ((A[:,i]' * θ[:,j]) + B[i] * u[j])).== 0.0)
end

for j in 1:k - 1
    @constraint(mod,y[j + 1] .== (C * θ[:,j + 1]))
end
#Weighing matrices
Q = C' * C * 1000.0;
R = 0.0001;

#err = sum(e' Q e)
err= @NLexpression(mod,(sum(Q[i,i] * (θref - θ[i,j]) ^ 2
                        for j in 1:k-1,i in 1:N)))
#in_err = sum(u'Ru)
in_err = @NLexpression(mod,(sum(R * u[j] ^ 2 for j in 1:k-1)))

J = @NLexpression(mod,0.5 * Δt * (err + in_err))
#defining the cost function
@NLobjective(mod, Min, J)
#optimizes the model
optimize!(mod)

#plotting values
states_θstart = zeros(k)
states_θend = zeros(k)
input_u = zeros(k - 1)
for i = 1 : k-1
    states_θstart[i] = JuMP.value(θ[1,i])
    states_θend[i] = JuMP.value(θ[end,i])
    input_u[i] = JuMP.value(u[i])
end

states_θend[end] = JuMP.value(θ[end,end])
states_θstart[end] = JuMP.value(θ[1,end])

plot(states_θstart,label = "θ1")
plot!(states_θend,label="y = Cθ")
savefig("lq-optimal_1D_Q1000_R0001_output-bounded.png")
p3 = plot(input_u,label="input,Q = 1000,R = 0.001")
savefig("lq-optimal_1D_Q1000_R0001-input-bounded.png")


text:
using Millboard
using Base.Test

@test """
+==================+=======+
| Base.WorkerState | value |
+==================+=======+
|        W_CREATED |     0 |
+------------------+-------+
|      W_CONNECTED |     1 |
+------------------+-------+
|    W_TERMINATING |     2 |
+------------------+-------+
|     W_TERMINATED |     3 |
+------------------+-------+""" == table(Base.WorkerState) |> string

@test """
+==================+=======+
| Base.WorkerState | value |
+==================+=======+
|        W_CREATED |     0 |
+------------------+-------+
|      W_CONNECTED |     1 |
+------------------+-------+
|    W_TERMINATING |     2 |
+------------------+-------+
|     W_TERMINATED |     3 |
+------------------+-------+""" == table(Base.W_CREATED) |> string

ext: jl


text:
using Base.Test

using Distances
using Clustering

srand(34568)

d = 3
n = 200
k = 10

X = rand(d, n)
costs = pairwise(SqEuclidean(), X)
@assert size(costs) == (n, n)

R = kmedoids(costs, k)
@test isa(R, KmedoidsResult)
@test length(R.medoids) == length(unique(R.medoids))
@test all(R.assignments .>= 1) && all(R.assignments .<= k)
@test R.assignments[R.medoids] == collect(1:k) # Every medoid should belong to its own cluster
@test sum(R.counts) == n
@test R.acosts == costs[sub2ind((n, n), R.medoids[R.assignments], 1:n)]
@test_approx_eq sum(R.acosts) R.totalcost
@test R.converged

# this data set has three obvious groups:
# group 1: [1, 3, 4], values: [1, 2, 3]
# group 2: [2, 5, 7], values: [6, 7, 8]
# group 3: [6, 8, 9], values: [21, 20, 22]
#

X = reshape(@compat(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22])), 1, 9)
costs = pairwise(SqEuclidean(), X)

R = kmedoids!(costs, [1, 2, 6])
@test isa(R, KmedoidsResult)
@test R.medoids == [3, 5, 6]
@test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
@test R.counts == [3, 3, 3]
@test_approx_eq R.acosts [1, 1, 0, 1, 0, 0, 1, 1, 1]
@test_approx_eq R.totalcost 6.0
@test R.converged


text:
module PicoSAT

export solve, itersolve

import Base: convert, iterate, IteratorSize

@static if Sys.isunix()
    const libpicosat = joinpath(dirname(@__FILE__), "..", "deps", "libpicosat.so")
end
@static if Sys.iswindows()
    throw(ErrorException("PicoSAT.jl does not currently work on Windows"))
end

const UNKNOWN = 0
const SATISFIABLE = 10
const UNSATISFIABLE = 20

struct PicoPtr
    ptr::Ptr{Cvoid}
end
convert(::Type{Ptr{Cvoid}}, p::PicoPtr) = p.ptr

version()   = bytestring(ccall((:picosat_version, libpicosat), Ptr{Cchar}, ()))
config()    = bytestring(ccall((:picosat_config,  libpicosat), Ptr{Cchar}, ()))
copyright() = bytestring(ccall((:picosat_copyright, libpicosat), Ptr{Cchar}, ()))

# constructor
picosat_init() = ccall((:picosat_init, libpicosat), PicoPtr, ())

# destructor
picosat_reset(p::PicoPtr) = ccall((:picosat_reset, libpicosat), Cvoid, (PicoPtr,), p)

"""
Measure all time spent in all calls in the solver.  By default only the
time spent in 'picosat_sat' is measured.  Enabling this function might
for instance triple the time needed to add large CNFs, since every call
to 'picosat_add' will trigger a call to 'getrusage'.
"""
picosat_measure_all_calls(p::PicoPtr) =
    ccall((:picosat_measure_all_calls, libpicosat), Cvoid, (PicoPtr,), p)

"""
Set the prefix used for printing verbose messages and statistics.
(Default is "c ")
"""
picosat_set_prefix(p::PicoPtr, str::String) =
    ccall((:picosat_set_prefix, libpicosat), Cvoid, (PicoPtr,Ptr{Cchar}), p, str)

"""
Set verbosity level
A verbosity level >= 1 prints more and more detailed progress reports to the output file.
Verbose messages are prefixed with the string set by 'picosat_set_prefix'
"""
picosat_set_verbosity(p::PicoPtr, level::Integer) =
    ccall((:picosat_set_verbosity, libpicosat), Cvoid, (PicoPtr,Cint), p, level)

"""
Disable (set_plain == true) / Enable all prerocessing.
"""
picosat_set_plain(p::PicoPtr, v::Bool) =
    ccall((:picosat_set_plain, libpicosat), Cvoid, (PicoPtr,Cint), p, v ? 1 : 0)

"""
If you know a good estimate on how many variables you are going to use
then calling this function before adding literals will result in less
resizing of the variable table.  But this is just a minor optimization.
Beside exactly allocating enough variables it has the same effect as
calling 'picosat_inc_max_var'.
"""
picosat_adjust(p::PicoPtr, max_idx::Integer) =
    ccall((:picosat_adjust, libpicosat), Cvoid, (PicoPtr,Cint), p, max_idx)

"""
As alternative to a decision limit you can use the number of propagations as limit.
This is more linearly related to execution time. This has to
be called after 'picosat_init' and before 'picosat_sat'.
"""
picosat_set_propagation_limit(p::PicoPtr, limit::Integer) =
    ccall((:picosat_set_propagation_limit, libpicosat), Cvoid, (PicoPtr,Culonglong), p, limit)

"""
Add a literal of the next clause.  A zero terminates the clause.  The
solver is incremental.  Adding a new literal will reset the previous
assignment.   The return value is the original clause index to which
this literal respectively the trailing zero belong starting at 0.
"""
picosat_add(p::PicoPtr, lit::Integer) =
    ccall((:picosat_add, libpicosat), Cint, (PicoPtr,Cint), p, lit)

"""
Call the main SAT solver.
A negative decision limit sets no limit on the number of decisions.
"""
picosat_sat(p::PicoPtr, limit::Integer) =
    ccall((:picosat_sat, libpicosat), Cint, (PicoPtr,Cint), p, limit)

"""p cnf <m> n"""
picosat_variables(p::PicoPtr) =
    ccall((:picosat_variables, libpicosat), Cint, (PicoPtr,), p)

"""
After 'picosat_sat' was called and returned 'PICOSAT_SATISFIABLE', then
the satisfying assignment can be obtained by 'dereferencing' literals.
The value of the literal is return as '1' for 'true',  '-1' for 'false'
and '0' for an unknown value.
"""
picosat_deref(p::PicoPtr, lit::Integer) =
    ccall((:picosat_deref, libpicosat), Cint, (PicoPtr,Cint), p, lit)

function add_clause(p::PicoPtr, clause)
    for lit in clause
        v = convert(Cint, lit)
        v == 0 && throw(ErrorException("PicoSAT Error: non zero integer expected"))
        picosat_add(p, v)
    end
    picosat_add(p, 0)
    return
end

function add_clauses(p::PicoPtr, clauses)
    for item in clauses
        add_clause(p, item)
    end
    return
end

function picosat_setup(clauses, vars::Integer, verbose::Integer, proplimit::Integer)
    p = picosat_init()
    if p.ptr === C_NULL
        picosat_reset(p)
        throw(ErrorException("PicoSAT Error: failed to initialize PicoSAT library"))
    end
    if verbose < 0
        picosat_reset(p)
        throw(ArgumentError("verbose must be ≥ 0"))
    end
    picosat_set_verbosity(p, verbose)
    if vars != -1
        picosat_adjust(p, vars)
    end
    if proplimit < 0
        picosat_reset(p)
        throw(ArgumentError("proplimit must be ≥ 0"))
    end
    if proplimit > 0
        picosat_set_propagation_limit(p, proplimit)
    end
    try
        add_clauses(p, clauses)
    catch ex
        picosat_reset(p)
        rethrow(ex)
    end
    return p
end

function get_solution(p::PicoPtr)
    nvar = picosat_variables(p)
    if nvar < 0
        picosat_reset(p)
        throw(ErrorException("number of solution variables < 0"))
    end
    sol = zeros(Int, nvar)
    for i = 1:nvar
        v = picosat_deref(p, i)
        @assert(v == 1 || v == -1)
        sol[i] = v * i
    end
    return sol
end
"""
`solve(clauses; vars::Integer=-1, verbose::Integer=0, proplimit::Integer=0)`
   - `vars` - the number of variables
   - `verbose` - prints solver logs to `STDOUT` when `verbose > 0` with increasing detail.
   - `proplimit` - helps to bound the execution time.  The number of propagations and the solution time are roughly linearly related.  A value of 0 (default) allows for an unbounded number of propagations.

Returns a solution if the problem is satisfiable.
Satisfiable solutions are represented as a vector of signed integers.
If the problem is not satisfiable the method returns an `:unsatisfiable` symbol.
If a solution cannot be found within the defined propagation limit, an `:unknown` symbol is returned.
```julia
julia> import PicoSAT
julia> cnf = Any[[1, -5, 4], [-1, 5, 3, 4], [-3, -4]];
julia> PicoSAT.solve(cnf)
5-element Array{Int64,1}:
   1
  -2
  -3
  -4
   5
```
"""
function solve(clauses;
               vars::Integer=-1,
               verbose::Integer=0,
               proplimit::Integer=0)
    p = picosat_setup(clauses, vars, verbose, proplimit)
    res = picosat_sat(p, -1)
    if res == SATISFIABLE
        result = get_solution(p)
    elseif res == UNSATISFIABLE
        result = :unsatisfiable
    elseif res == UNKNOWN
        result = :unknown
    else
        picosat_reset(p)
        throw(ErrorException("PicoSAT Error: return value $res"))
    end
    picosat_reset(p)
    return result
end

mutable struct PicoSolIterator
    ptr::PicoPtr
    vars::Vector{Int}

    function PicoSolIterator(p::PicoPtr)
        @assert p.ptr !== C_NULL
        iter = new(p, Int[])
        finalizer(i -> picosat_reset(i.ptr), iter)
        return iter
    end
end
"""
`itersolve(clauses; vars::Integer=-1, verbose::Integer=0, proplimit::Integer=0)`
   - `vars` - the number of variables
   - `verbose` - prints solver logs to `STDOUT` when `verbose > 0` with increasing detail.
   - `proplimit` - helps to bound the execution time.  The number of propagations and the solution time are roughly linearly related.  A value of 0 (default) allows for an unbounded number of propagations.

Returns an iterable object over all solutions.
When a user-defined propagation limit is specified, the iterator may not produce all feasible solutions.

```julia
julia> import PicoSAT
julia> cnf = Any[[1, -5, 4], [-1, 5, 3, 4], [-3, -4]];
julia> PicoSAT.itersolve(cnf)
julia> for sol in PicoSAT.itersolve(cnf)
           println(sol)
       end
[1,-2,-3,-4,5]
[1,-2,-3,4,-5]
[1,-2,-3,4,5]
[1,-2,3,-4,-5]
...
```
"""
function itersolve(clauses;
                   vars::Integer=-1,
                   verbose::Integer=0,
                   proplimit::Integer=0)
    p = picosat_setup(clauses, vars, verbose, proplimit)
    return PicoSolIterator(p)
end

# Add inverse of current solution to the clauses
function blocksol(it::PicoSolIterator)
    nvar = picosat_variables(it.ptr)
    if nvar < 0
        throw(ErrorException("PicoSAT Error: number of solution variables < 0"))
    end
    if length(it.vars) < nvar
        resize!(it.vars, nvar)
    end
    for i = 1:nvar
        it.vars[i] = picosat_deref(it.ptr, i) > 0 ? 1 : -1
    end
    for i = 1:nvar
        picosat_add(it.ptr, it.vars[i] < 0 ? i : -i)
    end
    picosat_add(it.ptr, 0)
    return
end

function next_solution(it::PicoSolIterator)
    res = picosat_sat(it.ptr, -1)
    if res == SATISFIABLE
        result = get_solution(it.ptr)
        # add inverse sol for next iter
        blocksol(it)
    elseif res == UNSATISFIABLE
        result = :unsatisfiable
    elseif res == UNKNOWN
        result = :unknown
    else
        throw(ErrorException("PicoSAT Errror: return value $res"))
    end
    return result
end

satisfiable(sol) = sol !== :unknown && sol !== :unsatisfiable

function Base.iterate(it::PicoSolIterator, state=nothing)
    sol = next_solution(it)
    if satisfiable(sol)
        return (sol, nothing)
    else
        return nothing
    end
end

IteratorSize(it::PicoSolIterator) = Base.SizeUnknown()
end # module


text:
#########
# TRACE #
#########

function update!(tr::OptimizationTrace{T,O}, state::S,
                 iteration::Integer, val::T, dt::Dict,
                 store_trace::Bool, show_trace::Bool, show_every::Int = 1,
                 callback = nothing) where {S<:AbstractOptimizerState, T, O}
    os = OptimizationTraceRecord{T,O}(iteration, val, dt)
    store_trace && push!(tr, os)
    if show_trace
        if iteration % show_every == 0
            show(os)
            print("\n")
            flush(stdout)
        end
    end
    if callback != nothing && (iteration % show_every == 0)
        if store_trace
            stopped = callback(tr)
        else
            stopped = callback(os)
        end
    else
        stopped = false
    end
    stopped
end

function trace!(tr, iteration, objfun, state, population, method, options, curr_time=time())
    dt = Dict{String,Any}()
    dt["time"] = curr_time
    # set additional trace value
    trace!(dt, objfun, state, population, method, options)
    update!(tr,
            state,
            iteration,
            value(objfun),
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

"""
    trace!(record::Dict{String,Any}, objfun, state, population, method, options)

Update the trace `record`. This function allows to supplement an additional information into the optimization algorithm trace by modifing a trace `record`. It can be overwiden by specifing particular parameter types.
"""
trace!(record::Dict{String,Any}, objfun, state, population, method, options) = ()

#########
# STATE #
#########

abschange(objfun::O, state::S) where {O<:AbstractObjective, S<:AbstractOptimizerState} =
    abschange(value(objfun), value(state))
abschange(curr::T, prev) where T = abs(curr - prev)
relchange(objfun::O, state::S) where {O<:AbstractObjective, S<:AbstractOptimizerState} =
    relchange(value(objfun), value(state))
relchange(curr::T, prev) where T = abs(curr - prev)/abs(curr)

function assess_convergence(objfun, state::AbstractOptimizerState, method, options::Options)
    converged = false

    if abschange(objfun, state) ≤ options.abstol
        converged = true
    end
    if relchange(objfun, state) ≤ options.reltol
        converged = true
    end

    return converged
end

##############
# POPULATION #
##############

"""
    initial_population(method, individual::AbstractVector)

Initialize population by replicating the `inividual` vector.
"""
initial_population(method::M, individual::I) where {M<:AbstractOptimizer, I<:AbstractVector} =
    return [copy(individual) for i in 1:population_size(method)]

"""
    initial_population(method, individual::Function)

Initialize population from the `inividual` function which returns an individual object.
"""
initial_population(method::M, individualFunc::Function) where {M<:AbstractOptimizer} =
    [individualFunc() for i in 1:population_size(method)]

"""
    initial_population(method, individuals::AbstractMatrix)

Initialize population from the `inividuals` matrix where each individual is a column.
"""
function initial_population(method::M, individuals::I) where {M<:AbstractOptimizer, I<:AbstractMatrix}
    n = population_size(method)
    @assert size(individuals,2) >= n "Size of initial population must be no smaller then $n"
    return [individuals[:,i] for i in 1:n]
end

"""
    initial_population(method, bounds::ConstraintBounds)

Initialize a random population within the individual `bounds`.
"""
function initial_population(method::M, bounds::ConstraintBounds) where {M<:AbstractOptimizer}
    n = population_size(method)
    cn = nconstraints_x(bounds)
    indv = rand(cn, n)
    if length(bounds.eqx) > 0
        indv[bounds.eqx,:] .= bounds.valx
    end
    T = eltype(bounds)
    rngs = Dict(i=>zeros(T,2) for i in bounds.ineqx)
    for (j,v,s) in zip(bounds.ineqx, bounds.bx, bounds.σx)
        rngs[j][1] -= s*v
        if s > 0
            rngs[j][2] += v
        end
    end
    for i in 1:n
        for (j,(r,l)) in rngs
            indv[j,i] *= r
            indv[j,i] += l
        end
    end
    initial_population(method, indv)
end


text:
"""
    plot_spectrum(eeg::SSR, chan::Int; kwargs...)
    plot_spectrum(eeg::SSR, chan::AbstractString; kwargs...)

Plot the spectrum of a steady state response measurement.
"""
function plot_spectrum(eeg::SSR, chan::Int; targetFreq::Number = modulationrate(eeg))

    channel_name = channelnames(eeg)[1]

    # Check through the processing to see if we have done a statistical test at target frequency
    signal = nothing
    result_idx = find_keys_containing(eeg.processing, "statistics")

    for r = 1:length(result_idx)
        result = get(eeg.processing, collect(keys(eeg.processing))[result_idx[r]], 0)
        if result[!, :AnalysisFrequency][1] == targetFreq

            result_snr = result[!, :SNRdB][chan]
            signal = result[!, :SignalAmplitude][chan]^2
            noise = result[!, :NoiseAmplitude][chan]^2
            title = "Channel $(channel_name). SNR = $(round(result_snr, sigdigits=4)) dB"
        end
    end

    if signal == nothing
        title = "Channel $(channel_name)"
        noise = 0
        signal = 0
    end

    title = replace(title, "_" => " ")

    avg_sweep = dropdims(Statistics.mean(eeg.processing["sweeps"], dims = 2), dims = 2)
    avg_sweep = avg_sweep[:, chan]
    avg_sweep = convert(Array{Float64}, vec(avg_sweep))

    p = plot_spectrum(
        avg_sweep,
        eeg.header["sampRate"][1];
        titletext = title,
        targetFreq = targetFreq,
        noise_level = noise,
        signal_level = signal,
    )

    return p
end

function plot_spectrum(
    eeg::SSR,
    chan::AbstractString;
    targetFreq::Number = modulationrate(eeg),
)

    return plot_spectrum(
        eeg,
        something(findfirst(isequal(chan), channelnames(eeg)), 0),
        targetFreq = targetFreq,
    )
end


text:
function checkWaveSpec!(inSpec::waveSpec)
if isdefined(inSpec,:f)
  nf=length(inSpec.f)
  if any(isnan(inSpec.f))
    #println("checkWaveSpec: NaNs found in f")
  end
  if !isdefined(inSpec,:df)
    inSpec.df=makeWaveSpecFreq(inSpec.f)
  end
else
  #println("checkWaveSpec: f not defined")
end

if isdefined(inSpec,:th)
  nth=length(inSpec.th)
  if any(isnan(inSpec.th))
    #println("checkWaveSpec: NaNs found in th")
  end
  if !isdefined(inSpec,:dth)
    inSpec.dth=makeWaveSpecDirn(inSpec.th)
  end
else
  #println("checkWaveSpec: th not defined")
end

if isdefined(inSpec,:Sf)
  inSpec.Sf=vec(inSpec.Sf)
  if isdefined(inSpec,:f)
    if length(inSpec.Sf)!=nf
       #println("checkWaveSpec: Sf is different lengths to f")
       return
    end
  end
end

if isdefined(inSpec,:Sth)
  inSpec.th=vec(inSpec.th)
  if isdefined(inSpec,:th)
    if length(inSpec.th)!=nth
       #println("checkWaveSpec: Sth is different lengths to th")
       return
    end
  end
end

if isdefined(inSpec,:f) & isdefined(inSpec,:th)
   fm=repmat(inSpec.f[:],1,nth)
   dfm=repmat(inSpec.df[:],1,nth)
   thm=repmat(pi/180*inSpec.th[:]',nf,1)
   dthm=repmat(pi/180*inSpec.dth[:]',nf,1)
end

if isdefined(inSpec,:S)
  nS=collect(size(inSpec.S))
  if isdefined(inSpec,:f) & isdefined(inSpec,:th)
    if any(nS.==nf) & any(nS.==nth)
      if nS[1]==nth
        inSpec.S=inSpec.S'
      end
    else
       #println("checkWaveSpec: S is different size to f and/or th")
       return
    end
  end
  if !isdefined(inSpec,:Sf)
    if isdefined(inSpec,:th)
      # Loop over freq
      inSpec.Sf=zeros(size(inSpec.S,1))
      for i=1:size(inSpec.S,1)
        sumS=0.
	# Sum all directions
        for j=1:size(inSpec.S,2)
          sumS+=dthm[j]*inSpec.S[i,j]
        end
	inSpec.Sf[i]=sumS
      end
      inSpec.Sf=vec(sum(dthm.*inSpec.S,2))
    end
  end
  if !isdefined(inSpec,:Sth)
    if isdefined(inSpec,:f)
      inSpec.Sth=zeros(size(inSpec.S,2))
      # Loop over dirns
      for i=1:size(inSpec.S,2)
        sumS=0.
	# Sum all frequencies
        for j=1:size(inSpec.S,1)
          sumS+=dfm[j]*inSpec.S[j,i]
        end
	inSpec.Sth[i]=sumS
      end
      inSpec.Sth=vec(sum(dfm.*inSpec.S,1))
    end
  end
end
if !isdefined(inSpec,:S) & isdefined(inSpec,:Sf)
  inSpec.S=inSpec.Sf[:,:]
end

return inSpec

end

"""
  Extends Base copy to copy MetoceanTools.waveSpec
"""
function Base.copy(inSpec::waveSpec)::waveSpec
  oSpec=waveSpec()
  isdefined(inSpec,:f) ? oSpec.f=copy(inSpec.f) : 0
  isdefined(inSpec,:df) ? oSpec.df=copy(inSpec.df) : 0
  isdefined(inSpec,:th) ? oSpec.th=copy(inSpec.th) : 0
  isdefined(inSpec,:dth) ? oSpec.dth=copy(inSpec.dth) : 0
  isdefined(inSpec,:Sf) ? oSpec.Sf=copy(inSpec.Sf) : 0
  isdefined(inSpec,:S) ? oSpec.S=copy(inSpec.S) : 0
  return oSpec
end

"""
  Function to split a full direction spectrum into two direction spectra
    by cutting along direction space specified by direction boundaries

  INPUT:
    inSpec - type waveSpec containing fields at least th and S both
               defining vector of directions (th) and direction spectra
         either full directon/frequency or direction spectrum
         Will copy over f,df and will split dth if is defined
    dSplit - defines the directions at which to split the spectrm [stDir,edDir]
  OUTPUT:
    o1,o2  - two partitions of type waveSpec

"""
function wvSpec_dirSplit(inSpec::waveSpec,dSplit)
    th=rem(inSpec.th+720,360)
    dSplitr=rem(dSplit+720,360)
    c1=indmin(abs(th-dSplitr[1]))
    c2=indmin(abs(th-dSplitr[2]))
    if c1==c2
       error("No split required as split directions found are identical at indices $(c1) and $(c2)")
    end
    o1=waveSpec()
    o2=waveSpec()
    o1.f=copy(inSpec.f)
    o2.f=copy(inSpec.f)
    o1.df=copy(inSpec.df)
    o2.df=copy(inSpec.df)
    nth=length(th)
    if c1>c2
      t=c2
      c2=c1
      c1=t
    end
    o1.th=inSpec.th[c1:(c2-1)]
    o2.th=[inSpec.th[c2:end];inSpec.th[1:(c1-1)]]
    o1.th=vec(o1.th)
    o2.th=vec(o2.th)
    if isdefined(inSpec,:df)
      o1.dth=inSpec.dth[c1:(c2-1)]
      o2.dth=[inSpec.dth[c2:end];inSpec.dth[1:(c1-1)]]
      o1.dth=vec(o1.dth)
      o2.dth=vec(o2.dth)
    end
    o1.S=inSpec.S[:,c1:(c2-1)]
    o2.S=[inSpec.S[:,c2:end] inSpec.S[:,1:(c1-1)]]
    checkWaveSpec!(o1)
    checkWaveSpec!(o2)
    return o1,o2
end
function wvSpec_dirSplit(inSpec::Array{waveSpec,1},dSplit)
  nSpec=length(inSpec)
  o1=Array{waveSpec,1}(nSpec)
  o2=Array{waveSpec,1}(nSpec)
  for i=1:nSpec
    o1[i],o2[i]=wvSpec_dirSplit(inSpec[i],dSplit)
  end
  return o1,o2
end

function makeWaveSpecDelDirn(th::Array{Float64,1})
  dth=similar(th)
  nth=length(th)
  for i=2:nth-1
    th1=rem(th[i-1]+720,360)
    th2=rem(th[i]+720,360)
    th3=rem(th[i+1]+720,360)
    dth1=abs(gmAngDiff(th2,th1))
    dth2=abs(gmAngDiff(th3,th2))
    dth[i]=0.5*(dth1+dth2)
  end
  dth[1]=dth[2]
  dth[nth]=dth[nth-1]
  return dth
end

function makeWaveSpecFreq()
  f=[1/x for x=30:-0.05:2]
  return f
end
function makeWaveSpecDelFreq(f::Array{Float64,1})
  df=copy(f)
  nf=length(f)
  df[1]=f[2]-f[1]
  df[2:(nf-1)]=( f[3:nf]-f[2:(nf-1)]     )/2 +
  ( f[2:(nf-1)]-f[1:(nf-2)] )/2
  df[nf]=f[nf]-f[nf-1]
  return df
end

#######################################################################
# Make a new wave spectrum
#######################################################################
#function makeSpecJONSWAP(Hs,Tp,Gamma,SigmaA,SigmaB)::waveSpec
function makeSpecJONSWAP(inParm::Array{Float64,1})::waveSpec
  # 0. Map input parameters
  (Hs,Tp,Gamma,SigmaA,SigmaB)=inParm
  # A. Make a new spectrum object
  s=waveSpec()
  # B. Create JONSWAP
  g::Float64=9.81
  fp::Float64=1./Tp
  nf=length(s.f)
  s.Sf=zeros(nf)
  for i=1:nf
    if s.f[i]<fp
      sigma=SigmaA
    else
      sigma=SigmaB
    end
    s.Sf[i]= g^2*(2*pi)^-4* s.f[i]^-5* exp(-5/4*(s.f[i]/fp)^-4)*
          Gamma^(exp(-1 * (((s.f[i]-fp)^2) /
	  (2*sigma^2*fp^2)) ))
  end
  sumS=0.
  for i=1:nf
    sumS+=s.Sf[i]*s.df[i]
  end
  alpha=Hs.^2./(16.*sumS)
  for i=1:nf
    s.Sf[i]=alpha.*s.Sf[i]
  end
  return s
end

#######################################################################
# Calculate wave spectrum integrated parameters
#######################################################################
function calcSpecPar(mySpec::waveSpec)

  # 0. Remove the zero frequency if present
  fZero=find(mySpec.f.==0)
  if !isempty(fZero)
    shift!(mySpec.f)
    shift!(mySpec.df)
    mySpec.S=mySpec.S[2:end,:]
    mySpec.S=mySpec.S[:,:]
  end

  # A. Detemine frequency cut-off
  fcut::Float64=mySpec.f[end]+mySpec.df[end]/2
  ecut::Float64=mySpec.S[end]

  # B. Calculate moment contributions from the tail
  m0t::Float64=1/4*fcut.*ecut
  m1t::Float64=1/3*fcut^2.*ecut
  m2t::Float64=1/2*fcut^3.*ecut

  # C. Calculate moments of the spectrum
  preCalc0::Array{Float64,1}=mySpec.df
  preCalc1::Array{Float64,1}=preCalc0.*mySpec.S[:,1]
  m0::Float64  = sum(sum(         preCalc1) + m0t)
  m1::Float64  = sum(sum(mySpec.f    .* preCalc1) + m1t)
  m2::Float64  = sum(sum(mySpec.f.^2 .* preCalc1) + m2t)
  mm1::Float64 = sum(sum(1./mySpec.f .* preCalc1))
  fe4::Float64 = sum(sum(mySpec.f .* preCalc0.*mySpec.S.^4))
  e4::Float64  = sum(sum(      preCalc0.*mySpec.S.^4))

  # D. Calculate integrated parameters
  Data=zeros(6)
  Name=["Hs","Tp","T01","T02","Tm01","T4"]
  Data[1]=4.*sqrt(m0);
  iTp::Int64=indmax(mySpec.S)
  Data[2]=1./(mySpec.f[iTp])
  Data[3]=m0./m1        #[s] Mean Period
  Data[4]=sqrt(m0./m2)  #[s] Zero-crossing Period
  Data[5]=mm1./m0      #[s] Energy Period
  # Calculate T4
  #  Weighted Peak by recommend by Young, Ocean Eng. Vol22 No.7 pp 669 to 686.
  #  The Determination of confidence limits associated with estimates of the
  #    spectral peak frequency
  Data[6]=1./(fe4./e4)

  return Data,Name
end

"""
  Part of:  fitSpecJONSWAPsingle

  Returns the RMS error of fit between single JONSWAP and input discrete spectrum

  Written by: Jason McConochie
  Revision 1, 30Dec2016, Jason McConochie
"""
function errorfitSpecJONSWAPsingle(wParm::Array{Float64,1},Sf::Array{Float64,1},fParm,fMask)
  maskReplace!(wParm,fParm,fMask)
  if any(wParm.<0.0)
    return 1000.
  end
  if wParm[3]<0.8
    return 100./(wParm[3])
  end
  solverSpec=makeSpecJONSWAP(wParm)
  sumS=0.
  for i=1:length(Sf)
    sumS+=(Sf[i]-solverSpec.Sf[i])^2
  end
  errRMS=sqrt(sumS)
  #errRMS=sqrt( sum((Sf.-solverSpec.Sf).^2) )
  return errRMS
end
"""
  fitSpecJONSWAPsingle

  Makes a single fit of the JONSWAP single spectrum to an input discrete spectrum
    interpolating the input spectrum to frequencies defined by waveSpec() initialisation

  Input:
     f,S - input spectrum
     fParm - array of [Hs,Tp,Gamma,SigmaA,SigmaB]
     fMask - Boolean array of mask indicating which of the fParm should be fixed in the fitting process
  Output:
     An array of parameters like fParm containing the best fit parameters

  Written by: Jason McConochie
  Revision 1, 30Dec2016, Jason McConochie
"""
function fitSpecJONSWAPsingle(f::Array{Float64,1},Sf::Array{Float64,1},fParm::Array{Float64,1},fMask::Array{Bool,1})
  # A. Interpolate to the standard waveSpec frequencies
  nS=waveSpec()
  fi=nS.f
  Si=interp1(f,Sf,fi)
  Si[Si.<0.]=0.
  r = optimize(b -> errorfitSpecJONSWAPsingle(b,Si,fParm,fMask), fParm)
  r1 = r.minimizer
  return r1
end
"""
  fitSpecJONSWAPsingle

  Makes a fit of the JONSWAP single spectrum to an input discrete spectrum
    interpolating the input spectrum to frequencies defined by waveSpec() initialisation

  This function takes an array of spectra defintions and returns a tsParm structure

  Input:
     f,S - input spectrum
     fParm - array of [Hs,Tp,Gamma,SigmaA,SigmaB]
     fMask - Boolean array of mask indicating which of the fParm should be fixed in the fitting process
  Output:
     A tsParm type of parameters
        sDate - dates starting from 1/1/2000 and then every hour
        Data - array containing Hs,Tp,Gamma,SigmaA,SigmaB on each row
        Name - Names of parameters

     This is suitable to be sent to tsQuickPlot

  Written by: Jason McConochie
  Revision 1, 30Dec2016, Jason McConochie
"""
function fitSpecJONSWAPsingle(inSpec::Array{waveSpec,1},fParm::Array{Float64,1},fMask::Array{Bool,1})
  o=tsParm()
  o.Data=zeros(Float64,(length(inSpec),5))
  o.Name=["Hs","Tp","Gamma","SigmaA","SigmaB"]
  o.sDate=Array{DateTime}(length(inSpec))
  for i=1:length(inSpec)
    o.Data[i,:]=fitSpecJONSWAPsingle(inSpec[i].f,inSpec[i].Sf,fParm,fMask)
    o.sDate[i]=DateTime(2000,1,1,0,0,0)+Dates.Day(i)
    solvedSpec=makeSpecJONSWAP(o.Data[i,:])
    #PyPlot.clf()
    #plot(inSpec[i].f,vec(inSpec[i].Sf),marker="o")
    #plot(solvedSpec.f,solvedSpec.Sf)
    #savefig(*("Jason",@sprintf("%015d",i)))
  end
  return o
end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
using Test, FileIO, ModelConstructors
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))
rp = joinpath(dirname(@__FILE__), "../reference/SLM") # reference path

# Set up inputs
input_inc = load(joinpath(rp, "find_solution_input.jld2"))
input_dec = load(joinpath(rp, "find_solution_input_decrease.jld2"))
input_sin = load(joinpath(rp, "find_solution_input_sine.jld2"))
in_inc = load(joinpath(rp, "design.jld2"))
in_dec = load(joinpath(rp, "design_decrease.jld2"))
in_sin = load(joinpath(rp, "design_sine.jld2"))
d = Dict()
default_slm_kwargs!(d)
d[:left_value] = input_inc["p0_norun"]
d[:right_value] = input_inc["p1_norun"]
d[:min_value] = input_inc["p0_norun"]
d[:max_value] = input_inc["p1_norun"]

# Scale problem
@testset "Scaling" begin
    d_false = deepcopy(d)
    d_false[:scaling] = false
    yscale_inc = load(joinpath(rp, "scaleproblem.jld2"))
    yscale_dec = load(joinpath(rp, "scaleproblem_decrease.jld2"))
    yscale_sin = load(joinpath(rp, "scaleproblem_sine.jld2"))
    leftval_inc = load(joinpath(rp, "left_value.jld2"))
    rightval_inc = load(joinpath(rp, "right_value.jld2"))
    for (i, in_data, yscale_data) in zip(1:3, [input_inc, input_dec, input_sin], [yscale_inc, yscale_dec, yscale_sin])
        if i == 1
            ŷ = scale_problem!(vec(in_data["w"]), vec(in_data["p_sol"]), d)
        elseif i == 2
            ŷ = scale_problem!(vec(in_data["w"]), vec(in_data["rev_p_sol"]), d)
        else
            ŷ = scale_problem!(vec(in_data["x"]), vec(in_data["y"]), d)
        end

        @test d[:y_scale] == yscale_data["YScale"]
        @test d[:y_shift] == yscale_data["YShift"]
        @test ŷ           == vec(yscale_data["yhat"])

        if i == 1
            @test d[:left_value] == leftval_inc["left_value"]
            @test d[:right_value] == rightval_inc["right_value"]
            @test d[:min_value] == in_data["p0_norun"] * d[:y_scale] + d[:y_shift]
            @test d[:max_value] == in_data["p1_norun"] * d[:y_scale] + d[:y_shift]
        end
    end
    @test vec(in_inc["y"]) == scale_problem!(vec(in_inc["x"]), vec(in_inc["y"]), d_false)
    @test d_false[:left_value] == input_inc["p0_norun"]
    @test d_false[:right_value] == input_inc["p1_norun"]
    @test d_false[:y_shift] == 0.
    @test d_false[:y_scale] == 1.
end

# Design matrix
Mdes = Dict()
rhs = Dict()
Mineq = Dict()
rhsineq = Dict()
Meq = Dict()
rhseq = Dict()
for (name, in_data) in zip([:inc, :dec, :sin], [in_inc, in_dec, in_sin])
    Mineq[name] = zeros(0, Int(in_data["nc"]))
    Meq[name] = zeros(0, Int(in_data["nc"]))
    rhsineq[name] = Vector{Float64}(undef, 0)
    rhseq[name] = Vector{Float64}(undef, 0)
    Mdes[name] = construct_design_matrix(vec(in_data["x"]), vec(in_data["knots"]),
                                         vec(in_data["dx"]), Int.(vec(in_data["xbin"])),
                                         Int(in_data["nx"]), Int(in_data["nk"]), Int(in_data["nc"]))
    rhs[name] = in_data["y"]
end
@testset "Design matrix" begin
    for (name, in_data) in zip([:inc, :dec, :sin], [in_inc, in_dec, in_sin])
        @test @test_matrix_approx_eq Mdes[name] in_data["Mdes"]
        @test rhs[name] == in_data["rhs"]
    end
end

# Regularizer matrix
Mreg = Dict()
for (name, in_data) in zip([:inc, :dec, :sin], [in_inc, in_dec, in_sin])
    Mreg[name] = construct_regularizer(vec(in_data["dx"]), Int(in_data["nk"]))
end
reg_inc = load(joinpath(rp, "regularizer.jld2"))
reg_dec = load(joinpath(rp, "regularizer_decrease.jld2"))
reg_sin = load(joinpath(rp, "regularizer_sine.jld2"))
@testset "Regularizer" begin
    for (name, reg_data) in zip([:inc, :dec, :sin], [reg_inc, reg_dec, reg_sin])
        @test @test_matrix_approx_eq Mreg[name] reg_data["Mreg"]
        @test all(reg_data["rhsreg"] .== 0.)
    end
end

# Test C2
C2_inc = load(joinpath(rp, "C2.jld2"))
C2_dec = load(joinpath(rp, "C2_decrease.jld2"))
C2_sin = load(joinpath(rp, "C2_sine.jld2"))
for (name, in_data) in zip([:inc, :dec, :sin], [in_inc, in_dec, in_sin])
    Meq[name], rhseq[name] = C2_matrix(Int(in_data["nk"]), Int(in_data["nc"]), vec(in_data["dx"]), Meq[name], rhseq[name])
end
@testset "C2" begin
    for (name, C2_data) in zip([:inc, :dec, :sin], [C2_inc, C2_dec, C2_sin])
        @test @test_matrix_approx_eq Meq[name] C2_data["Meq"]
        @test @test_matrix_approx_eq rhseq[name] C2_data["rhseq"]
    end
end

# Left and right values
left_inc = load(joinpath(rp, "left_value.jld2"))
left_dec = load(joinpath(rp, "left_value.jld2")) # did not create a "decreasing" version, so we just use the same file as increasing
left_sin = load(joinpath(rp, "left_value_sine.jld2")) # and adjust accordingly, noting that the left and right values are the reverse ones
right_inc = load(joinpath(rp, "right_value.jld2"))
right_dec = load(joinpath(rp, "right_value.jld2"))
right_sin = load(joinpath(rp, "right_value_sine.jld2"))
for (i, name, side_data, in_data) in zip(1:3, [:inc, :dec, :sin], [left_inc, right_dec, left_sin], [in_inc, in_dec, in_sin])
    if i == 2
        Meq[name], rhseq[name] = set_right_value(side_data["right_value"], Int(in_data["nc"]), Int(in_data["nk"]), Meq[name], rhseq[name])
    else
        Meq[name], rhseq[name] = set_left_value(side_data["left_value"], Int(in_data["nc"]), Meq[name], rhseq[name])
    end
end
for (i, name, side_data, in_data) in zip(1:3, [:inc, :dec, :sin], [right_inc, left_dec, right_sin], [in_inc, in_dec, in_sin])
    if i == 2
        Meq[name], rhseq[name] = set_left_value(side_data["left_value"], Int(in_data["nc"]), Meq[name], rhseq[name])
    else
        Meq[name], rhseq[name] = set_right_value(side_data["right_value"], Int(in_data["nc"]), Int(in_data["nk"]), Meq[name], rhseq[name])
    end
end
@testset "Left and right value" begin
    for (name, in_data) in zip([:inc, :sin], [right_inc, right_sin])
        @test @test_matrix_approx_eq Meq[name] in_data["Meq"]
        @test @test_matrix_approx_eq rhseq[name] in_data["rhseq"]
    end
end

# Global minimum and maximum
minmax_inc = load(joinpath(rp, "min_max.jld2"))
minmax_dec = load(joinpath(rp, "min_max_decrease.jld2"))
minmax_sin = load(joinpath(rp, "min_max_sine.jld2"))
for (name, in_data, minmax_data) in zip([:inc, :dec], [in_inc, in_dec], [minmax_inc, minmax_dec])
    Mineq[name], rhsineq[name] = set_min_value(minmax_data["min_value"], Int(in_data["nk"]), Int(in_data["nc"]),
                                               vec(in_data["dx"]), Mineq[name], rhsineq[name])
    Mineq[name], rhsineq[name] = set_max_value(minmax_data["max_value"], Int(in_data["nk"]), Int(in_data["nc"]),
                                               vec(in_data["dx"]), Mineq[name], rhsineq[name])
end
@testset "Global minimum and maximum" begin
    for (name, minmax_data) in zip([:inc, :dec], [minmax_inc, minmax_dec])
        @test @test_matrix_approx_eq Mineq[name]  minmax_data["Mineq"]
        @test @test_matrix_approx_eq rhsineq[name] vec(minmax_data["rhsineq"])
    end
end

# Monotonicity
mono_inc = load(joinpath(rp, "monotone.jld2"))
mono_dec = load(joinpath(rp, "monotone_decrease.jld2"))
mono_sin = load(joinpath(rp, "monotone_sine.jld2"))
inc_int = [0 π/2; 3*π/2 5*π/2; 7*π/2 4*pi]
dec_int = [π/2 3*π/2; 5*π/2 7*π/2]
for (i, name, in_data, mono_data) in zip(1:3, [:inc, :dec, :sin], [in_inc, in_dec, in_sin], [mono_inc, mono_dec, mono_sin])
    monotone_settings = Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}}(undef, 0)
    total_monotone_intervals = 0
    if i == 1
        total_monotone_intervals += monotone_increasing!(monotone_settings, Int(in_data["nk"]))
    elseif i == 2
        total_monotone_intervals += monotone_decreasing!(monotone_settings, Int(in_data["nk"]))
    elseif i == 3
        total_monotone_intervals += increasing_intervals_info!(monotone_settings, vec(in_data["knots"]), inc_int, Int(in_data["nk"]))
        total_monotone_intervals += decreasing_intervals_info!(monotone_settings, vec(in_data["knots"]), dec_int, Int(in_data["nk"]))
    end
    Mineq[name], rhsineq[name] = construct_monotonicity_matrix(monotone_settings, Int(in_data["nc"]), Int(in_data["nk"]),
                                                              vec(in_data["dx"]), total_monotone_intervals, Mineq[name], rhsineq[name])
end
@testset "Monotonicity" begin
    for (name, mono_data) in zip([:inc, :dec, :sin], [mono_inc, mono_dec, mono_sin])
        @test @test_matrix_approx_eq Mineq[name]  mono_data["Mineq"]
        @test @test_matrix_approx_eq rhsineq[name] mono_data["rhsineq"]
    end
end

# Curvature
curv_inc = load(joinpath(rp, "curvature.jld2"))
curv_dec = load(joinpath(rp, "curvature_decrease.jld2"))
curv_sin = load(joinpath(rp, "curvature_sine.jld2"))
cu_int  = [π 2*π; 3*π 4*π]
cd_int  = [0 π; 2*π 3*π]
for (i, name, in_data) in zip(1:3, [:inc, :dec, :sin], [in_inc, in_dec, in_sin])
    curvature_settings = Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{Float64}}}}(undef, 0)

    if i == 1
        concave_down_info!(curvature_settings)
    elseif i == 2
        concave_up_info!(curvature_settings)
    elseif i == 3
        concave_up_intervals_info!(curvature_settings, cu_int)
        concave_down_intervals_info!(curvature_settings, cd_int)
    end
    Mineq[name], rhsineq[name] = construct_curvature_matrix(curvature_settings, Int(in_data["nc"]), Int(in_data["nk"]), vec(in_data["knots"]),
                                                              vec(in_data["dx"]), Mineq[name], rhsineq[name])
end
@testset "Curvature" begin
    for (name, curv_data) in zip([:inc, :dec, :sin], [curv_inc, curv_dec, curv_sin])
        @test @test_matrix_approx_eq Mineq[name]  curv_data["Mineq"]
        @test @test_matrix_approx_eq rhsineq[name] vec(curv_data["rhsineq"])
    end
end


text:
function run_binary_BH_simulation(;

        #Time variables
        tspan::Array{Float64 , 1}   =   [0. , 0.]   ,
        dt::Float64                 =   0.05        ,
        out_every::Float64          =   1.0         ,
        max_run_time::String        =   "24:00:00"  ,

        #Spatial variables
        Nnodes:: Int64              =   128     ,
        Rcavity::Float64            =   40.     ,

        #Black hole parameters
        Mbh::Float64                =   1.0     ,
        Rorbit::Float64             =   6.0     ,
        alpha::Float64              =   10.0    ,

        #Initial configuration
        Gaussian_pulse:: Array{Float64, 1} = [ 3.5 , 3.7 , 20.  , 0.1 , 2 ],

        output_folder              = ""
        )

    if(tspan[2] - tspan[1] < 0)
       println("ERROR: t_final must be larger than t_initial")
        return
    end

    #Define variables
    Ωorbit = sqrt(2*Mbh/Rorbit^3)

    #filename
    folder_name = output_folder

    if(output_folder == "")
        folder_name = string("/", Dates.today() , "/")
    end

    filename = string("binary_bh_M_", Mbh, "_alpha_", alpha , "_Rorbit_", Rorbit , "_Rcavity_",Rcavity , "_N_" , Nnodes ,"_ti_" , tspan[1] , "_tf_" , tspan[2] )

    #Setup masses
    p0 = Param(

    #Temporal variables
    t_sim_init      = tspan[1]          ,
    t_sim_final     = tspan[2]          ,
    out_every_t     = out_every      ,
    deltat          = dt             ,

    max_runtime     = max_run_time  ,
    #Spacial variables
    xmin           =    -Rcavity    ,
    xmax           =    Rcavity     ,
    xnodes         =    Nnodes      ,

    #Spacial variables
    ymin           =    -Rcavity    ,
    ymax           =    Rcavity     ,
    ynodes         =    Nnodes      ,

    #Initial Conditions
    A0              =   Gaussian_pulse[1] , #3.5         ,
    σ               =   Gaussian_pulse[2] , #4.5         ,
    r0              =   Gaussian_pulse[3] , #40.0        ,
    ω               =   Gaussian_pulse[4] , #0.5         ,
    m               =   Gaussian_pulse[5] , #2           ,

    #Klein Gordon Mass
    μ               =   0.0         ,

    #Dissipation terms
    dissipation     =  :true        ,

    R_orbit         = Rorbit        ,
    R1              = 2*Mbh         ,
    R2              = 2*Mbh         ,

    Ω               = Ωorbit        ,
    alpha1          = alpha         ,
    alpha2          = alpha         ,

    #Periodic BC?
    Boundaries      = :radial       , # :square :periodic

    #Data variables
    folder          = folder_name   ,
    fname           = filename

    )

    solve_wave_equation_2D(p0)

end


text:
using MAT, LaTeXStrings, Plots, ProgressBars, Statistics, StatsBase, StatsPlots, Ranges
cd(dirname(@__FILE__))
file = matopen("i140703-001_lfp-spikes.mat")
vars = matread("i140703-001_lfp-spikes.mat")

#   block->segments-> anaglogsignals (LFB signals_Q8) cells:(1:96)
#                  -> spiketrains
#Part A
neurons = vars["block"]["segments"][1,1]["spiketrains"]
events = vars["block"]["segments"][1,1]["events"][1,1]
all_neurons_spiketimes=[[] for i in 1:length(neurons)]
for index in 1:length(neurons)
    all_neurons_spiketimes[index]= neurons[1, index]["times"][1, :] ./ 30000
end

scatter()
for i in 1:length(neurons)
    scatter!(all_neurons_spiketimes[i], [i for j in 1:length(all_neurons_spiketimes[i])],
    markershape=:vline, size = (4000,2500), color=:black, tickfontsize=30)
end
plot!(legend=false, dpi =400, ylabel=L"trial id", xlabel="time(ms)")
TS_ON = events["times"][findall(x -> x == "65296", events["labels"])] ./ 3000
WS_ON = events["times"][findall(x -> x == "65360", events["labels"])] ./ 3000
CUE_ON_PG = events["times"][findall(x -> x == "65365", events["labels"])]  ./ 3000
#CUE_OFF = events["times"][findall(x -> x == "65360", events["labels"])]  ./ 3000
GO_ON_LF = events["times"][findall(x -> x == "65369", events["labels"])]  ./ 3000
RW_ON_CONF_PG = events["times"][findall(x -> x == "65509", events["labels"])]  ./ 3000
plot!(TS_ON)

savefig("../../Fundamentals-of-Neuroscience_Fall1400/Term_Project/Figs/rasterplot.pdf")


text:
# GenerateAllPartitionArray.jl
module GenerateAllPartitionArray

    using Combinatorics, DelimitedFiles, Statistics

    export
        # preprocessing functions
        isInputMatrixProper,
        findNumOfPartitions,
        requestUserInputToSelectHeuristics,
        printInputMatrixSummary,
        getUserInput,
        computeListOfAllCompositions,
        computeLexicalOrderCompostitions,
        groupCompositionsByPartition,
        computeListOfAllCombinationsOfOverlaps,
        findOverlapRowIndices,
        findPartitionsToReserveForLast,

        # core functions
        generateAllPartitionArray,
        searchMatrixForCoveringOfCandidates,
        makeInputMatrixSuitableForSearch,
        findUnUsedCompositions,
        findUnUsedPartitions,
        findSumOfUnUsedParts,
        findCandidateCompositionsAtPosition,
        isCompositionPossible,
        isRegionPossible,
        findMissingRegionPcs,
        findDuplicateRegionPcs,
        findPossibleOverlapsForRegion,
        findCandidateOverlapsForRegion,
        selectHeuristicsAndCompute,
        judgeQualityOfCandidate,
        sortCandidatesByHeuristicScore,
        success,
        backtrack,
        printSolution

        # files containing preprocessing and core functions
        include("./Preprocess.jl")
        include("./Core.jl")

#= Read in input matrix from file =#
const inputMatrix = readdlm("Babbitt_6_96.txt", '\t', Int8, '\n')   # as found e.g., in Babbitt's Arie da Capo
#const inputMatrix = readdlm("Smalley_6_96.txt", '\t', Int8, '\n')  # as found e.g., in Babbitt's Sheer Pluck (constructed by David Smalley)
#const inputMatrix = readdlm("Babbitt_4_96.txt", '\t', Int8, '\n')  # as found e.g., in Babbitt's My Ends are My Beginings

#= Define primary types =#
"""
    Immutable composite type called Candidate

A left-to-right sequence of K Candidates are needed to form a cover (i.e., a solution and complete all-partition array).
"""
struct Candidate
    partition::Int8
    composition::Array{Int8,1}
    overlaps::BitArray{2}
end

"""
    Mutable composite type called Region
"""
mutable struct Region
    pcs::Array{Int8,1}
    inPushPcs::Array{Int8,1} # indexed by overlapsAllowedHere
    outPushPcs::Array{Int8,1} # indexed by overlapsAllowedHere
    missingPcs::Array{Int8,1} #
    duplicatePcs::Array{Int8,1} #
    # inner constructor of default values
    # Region() = new(zeros(Int8,12), zeros(Int8,6), zeros(Int8,6), Array{Int8,1}(), Array{Int8,1}()) # remove?
end

"""
    Immutable composite type called APArray
"""
struct APArray
    numRows::Int8
    numCols::Int8
    numElements::Int16
    numOverlaps::Int16
    numPartitions::Int8
    aggregate::Array{Int8,1}
end

"""
    Immutable composite type called Heuristic
"""
struct Heuristic
    offset::Int8
    heuristics::Array{Int8,1}
    reserves::Int8
end

"""
    Mutable composite type called Selected
"""
mutable struct Selected
    candidates::Array{Int32,1}
    overlaps::Array{Int16,1}
end

#= Run main algorithm =#
selected, allCandidateLists, cnt, position = generateAllPartitionArray(inputMatrix)

end
#= End of module =#


text:
!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSPRemezWrap, Compat, Compat.Test

#
# Length 151 LPF (Low Pass Filter).
#
@testset "remez_151_lpf" begin
    h = remez_wrap(151, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_151_lpf.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 152 LPF. Non-default "weight" input.
#
#    from scipy.signal import remez
#    lpf = remez(152, [0, 0.475, 0.5, 1.0], [1.0, 0.0], weight=[1,2], Hz=2.0)
#    lpf.tofile('remez_152_lpf.txt', sep='\n')
#
@testset "remez_152_lpf" begin
    h = remez_wrap(152, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; weight=[1.0,2.0], Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_152_lpf.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 51 HPF (High Pass Filter).
#
#    from scipy.signal import remez
#    hpf = remez(51, [0, 0.75, 0.8, 1.0], [0.0, 1.0], Hz=2.0)
#    hpf.tofile('remez_51_hpf.txt', sep='\n')
#
@testset "remez_51_hpf" begin
    h = remez_wrap(51, [0, 0.75, 0.8, 1.0], [0.0, 1.0]; Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_51_hpf.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 180 BPF (Band Pass Filter).
#
#    from scipy.signal import remez
#    bpf = remez(180, [0, 0.375, 0.4, 0.5, 0.525, 1.0], [0.0, 1.0, 0.0], Hz=2.0)
#    bpf.tofile('remez_180_bpf.txt', sep='\n')
#
@testset "remez_180_bpf" begin
    h = remez_wrap(180, [0, 0.375, 0.4, 0.5, 0.525, 1.0], [0.0, 1.0, 0.0]; Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_180_bpf.txt"),'\t')
    @test h ≈ h_scipy
end

#
#  Odd-symmetric filters - hilbert and differentiators type.
#  Even length - much better approximation since it is not constrained to 0 at
#  the nyquist frequency
#
# Length 20 hilbert
#
#    from scipy.signal import remez
#    h = remez(20, [0.1, 0.95], [1], type="hilbert", Hz=2.0)
#    h.tofile('remez_20_hilbert.txt', sep='\n')
#
@testset "remez_20_hilbert" begin
    h = remez_wrap(20, [0.1, 0.95],[1.0]; filter_type=filter_type_hilbert, Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_20_hilbert.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 21 hilbert
#
#    from scipy.signal import remez
#    h = remez(21, [0.1, 0.95], [1], type="hilbert", Hz=2.0)
#    h.tofile('remez_21_hilbert.txt', sep='\n')
#
@testset "remez_21_hilbert" begin
    h = remez_wrap(21, [0.1, 0.95],[1.0]; filter_type=filter_type_hilbert, Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_21_hilbert.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 200 differentiator
#
#    from scipy.signal import remez
#    h = remez(200,[0.01, 0.99],[1],type="differentiator" Hz=2.0)
#    h.tofile('remez_200_differentiator.txt', sep='\n')
#
@testset "remez_200_differentiator" begin
    h = remez_wrap(200, [0.01, 0.99],[1.0]; filter_type=filter_type_differentiator, Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_200_differentiator.txt"),'\t')
    @test h ≈ h_scipy
end

#
# Length 201 differentiator
#
#    from scipy.signal import remez
#    h = remez(201,[0.05, 0.95],[1],type="differentiator" Hz=2.0)
#    h.tofile('remez_201_differentiator.txt', sep='\n')
#
@testset "remez_201_differentiator" begin
    h = remez_wrap(201, [0.05, 0.95],[1.0]; filter_type=filter_type_differentiator, Hz=2.0);
    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_201_differentiator.txt"),'\t')
    @test h ≈ h_scipy
end


text:
using Symplectic
using LinearAlgebra
using ForwardDiff

const bs = beamSplitter
const ϕ = π/4
θ = π / 4
invBS = inv(beamSplitter(π/4)) ⊕ I(2*2)
S = bs(ϕ, 1, 2, 4)*bs(θ, 2, 3, 4)*bs(ϕ, 3, 4, 4)*bs(ϕ, 1, 4, 4)*bs(ϕ, 1, 2, 4)*bs(θ, 2, 3, 4)*bs(ϕ, 3, 4, 4)*bs(ϕ, 1, 4, 4)
S = *(generify(S; randomize = randomPhaseShifiting )...)

dseq2 = decoupleSequence([S for i in 1:4]...);
dseq1 = decoupleSequence([i==4 ? S*invBS : S for i in 1:4]...);
S2 = *(dseq2...);
S1 = *(dseq1...);
S2 = S2[1:2, 1:2] ⊕ *(generify( S2[3:end, 3:end], randomize = randomLocalSymplecticMatrix)...)

genericity(S2[3:end, 3:end])
round.(S2[3:end, 3:end]; digits = 3)

seq = ibs(S; T = BS);
round.(*(seq...); digits = 3)
round.(seq[5]; digits=3)

tmp = [ 1 1 1 1;
        -1 1 -1 1;
        -1  -1 1 1;
        1 -1 -1 1] / 2

nonSymplecticity(tmp)

transpose(tmp) * Ω * tmp

a, b, c = preIwasawaFactorization(tmp)
a
b
c

lstf(tmp[:,1], tmp[1,:])

S * lstf(S[:, 1], S[1, :]) * S

opnorm(round.(transpose(S) * S - I))

function func(θs::Vector)
        T = teleportation(phaseShifting([0, 0, θs[1:4]...])*Matrix(S)*phaseShifting([θs[1:4]..., 0, 0]), [1, 2], [1, 2])
        return tr((T - I(4)) * transpose(T - I(4)))
end

g = x -> ForwardDiff.gradient(func, x)

g([0,0,0,0])


text:
using GPAR: GPARegressor, FiniteGPARegressor, is_closed_downwards, replaces, imputes

@testset "gpar" begin

    @testset "GPARegressor" begin
        f1, f2 = GP(EQ(), GPC()), GP(EQ(), GPC())
        gpar = GPARegressor([f1, f2])

        @test replaces(gpar) == false
        @test imputes(gpar) == false
    end

    @testset "indexing" begin
        rng, N = MersenneTwister(123456), 10
        x = randn(rng, N)
        f = GP(EQ(), GPC())
        gpar = GPARegressor([f, f]; replaces=true, imputes=true)
        finite_gpar = gpar(x)

        @test finite_gpar isa FiniteGPARegressor
    end

    @testset "rand" begin
        rng, N, S = MersenneTwister(123456), 10, 3
        x = randn(rng, N)
        f1 = GP(ZeroMean(), EQ(), GPC())
        f2 = GP(OneMean(), EQ(), GPC())
        gpar = GPARegressor([f1, f2])(x)

        @test size(rand(rng, gpar, 1)) == (N, 2, 1)
        @test size(rand(rng, gpar, S)) == (N, 2, S)
    end

    @testset "logpdf" begin

    end
end


text:
export mTFI

function mTFI(θ,m=20)
    σw = θ[1]
    σε = θ[2]
    λ   = exp(θ[3])
    μ0 = θ[4]
    σ0 = θ[5]

    πs = [-ω(j,λ,5/6) for j in 1:m]

    F = [πs'; Matrix{Float64}(I,m-1,m-1) zeros(m-1)]
    Θ = [1.0; zeros(m-1,1)]
    G = zeros(1,1)

    H = [1.0 zeros(1,m-1)]
    Γ = zeros(1,1)

    Q = fill(σw^2,1,1)
    R = fill(σε^2,1,1)

    μ = fill(μ0,m)
    Σ0 = σ0^2*Matrix{Float64}(I,m,m)

    LinearStateSpaceModel(F,Q,G,H,R,Γ,Θ,μ,Σ0)
end

# function _gradient(θ,m=20)
#     σw = θ[1]
#     σv = θ[2]
#     σε = θ[3]
#     λ   = exp(θ[4])

#     dQdσw = [2*σw 0; 0 0]
#     dQdσv = [0 0; 0 2*σv]

#     dRdσε = fill(2*σε,1,1)

#     πs = λ*[-dω(j,λ,5/6) for j in 1:m]

#     dFdλ = [zeros(1,m+2); zeros(1,m+2); 0 0 πs'; zeros(m-1,m+2)]

#     dF = cat(zeros(m+2,m+2,3),dFdλ,dims=[3])
#     dQ = cat(dQdσw,dQdσv,zeros(2,2,2),dims=[3])
#     dH = zeros(1,m+2,4)
#     dR = cat(zeros(1,1,2),dRdσε,zeros(1,1),dims=[3])

#     dΘ = zeros(m+2,2,4)

#     dμ = zeros(m+2,1,4)
#     dΣ0 = zeros(m+2,m+2,4)

#     LinearStateSpaceGradient(dF,dQ,dH,dR,dΘ,dμ,dΣ0)
# end