text:
module TerraDG
using WriteVTK
using Printf
using Logging
using LinearAlgebra
import YAML

include("configuration.jl")
include("basis.jl")
include("equations.jl")
include("grid.jl")
include("kernels/surface.jl")
include("kernels/volume.jl")
include("kernels/filtering.jl")
include("kernels/time.jl")
include("plotters.jl")
include("error_writer.jl")
include("global_matrices.jl")
include("kernels/limiter.jl")

"""
    evaluate_rhs(eq, scenario, filter, globals, du, dofs, grid)

Evalutes the right-hand-side of the equation `eq` for
scenario `scenario`, with filter `filter`,
collection of global matrices `globals`, update
`du`, degrees of freedom `dofs` and grid `grid`.

Updates `du` in place.
"""
function evaluate_rhs(eq, scenario, filter, globals, du, dofs, grid, config)
    buffers_face = BuffersFaceIntegral(grid.basis, get_ndofs(eq))
    buffers_volume = BuffersVolume(grid.basis, get_ndofs(eq))
    linearindex = vec(collect(Iterators.product(1:grid.basis.order, 1:grid.basis.order))) #linear indexing
    var = get_ndofs(eq)
    reference_massmatrix = massmatrix(grid.basis, grid.basis.dimensions)
    du .= 0.0
    maxeigenval = -Inf

    ∇ = globals.reference_derivative_matrix
    for i in eachindex(grid.cells)
        @views cell = grid.cells[i]
        @views data = dofs[:,:, cell.dataidx]
        @views flux = grid.flux[:,:,cell.dataidx]

        #evaluate_flux(eq, data, flux)

        #ATTEMPT 1
        """
        if config.source == "without_source"
            evaluate_flux(eq, data, flux)
        else
            evaluate_flux_with_source(eq, data, flux, cell, grid.basis)
        end
        """

        x, y = globalposition(cell, (0.5, 0.5))

        if config.scenario_name == "stiff_inclusion"
            evaluate_flux_with_stiffinclusion(eq, scenario, data, flux, (x,y))
        else
            evaluate_flux(eq, data, flux)
        end

        # Volume matrix is zero for FV/order=1
        if length(grid.basis.quadpoints) > 1
            @views evaluate_volume(globals, buffers_volume, flux, grid.basis, inverse_jacobian(cell), volume(cell), du[:,:,cell.dataidx])
            #@info "after VI"
            #@info du[:,:,cell.dataidx]
        end
    end
    for i in eachindex(grid.cells)
        @views cell = grid.cells[i]
        @views data = dofs[:,:, cell.dataidx]
        @views flux = grid.flux[:,:,cell.dataidx]
        elem_massmatrix = volume(cell) * reference_massmatrix
        inv_massmatrix = inv(elem_massmatrix)

        # Here we also need to compute the maximum eigenvalue of each cell
        # and store it for each cell (needed for timestep restriction later!)
        faces = [left, top, right, bottom]
        for (i, neigh) in enumerate(cell.neighbors)
            @views dofsneigh = dofs[:,:,neigh.dataidx]
            @views fluxneigh = grid.flux[:,:,neigh.dataidx]

            # @info data
            # Project dofs and flux of own cell to face
            project_to_faces(globals, data, flux, buffers_face.dofsface, buffers_face.fluxface, faces[i])
            facetypeneigh = cell.facetypes[i]

            #@info data

            if facetypeneigh == regular
                # Project neighbors to faces
                # Neighbor needs to project to opposite face
                faceneigh = globals.oppositefaces[faces[i]]
                project_to_faces(globals, dofsneigh, fluxneigh, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh, faceneigh)
                #@info "from if"
                #@info data
            else
                #@info "from else"
                #@info data
                @assert(facetypeneigh == boundary)
                normalidx = globals.normalidxs[faces[i]]

                # For boundary cells, we operate directly on the dofsface
                evaluate_boundary(eq, scenario, faces[i], normalidx, buffers_face.dofsface, buffers_face.dofsfaceneigh)
                # Evaluate flux on face directly
                # Note: When extrapolating, this is not exact!
                # The error is given by the commutation error of face projection and flux!

                #evaluate_flux(eq, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh)

                # ATTEMPT 1
                """
                if config.source == "without_source"
                    evaluate_flux(eq, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh)
                else
                    evaluate_flux_with_source(eq, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh, neigh, grid.basis)
                end
                """
                x, y = globalposition(cell, (0.5, 0.5))
                if config.scenario_name == "stiff_inclusion"
                    evaluate_flux_with_stiffinclusion(eq, scenario, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh, (x,y))
                else
                    evaluate_flux(eq, buffers_face.dofsfaceneigh, buffers_face.fluxfaceneigh)
                end

                #@info data
            end
            @views cureigenval = evaluate_face_integral(eq, globals, buffers_face, cell, faces[i], du[:,:,cell.dataidx], config.riemann)
            maxeigenval = max(maxeigenval, cureigenval)
        end

        @views du[:,:,cell.dataidx] = inv_massmatrix * @views du[:,:,cell.dataidx] #
    end
    grid.maxeigenval = maxeigenval
end

"""
    main(configfile::String)

Runs a DG-simulation with configuration from `configfile`.
"""
function main(configfile::String)
    config = Configuration(configfile)
    filter = make_filter(config)
    eq = make_equation(config)
    scenario = make_scenario(config)
    grid = make_grid(config, eq, scenario)
    integrator = make_timeintegrator(config, grid)

    #@info size(grid.dofs)
    #@info grid.dofs[:,:,40]
    #@info grid.cells[5].size
    #@info grid.cells[5].facetypes

    @info "Initialising global matrices"
    globals = GlobalMatrices(grid.basis, filter, grid.basis.dimensions)
    @info "Initialised global matrices"

    filename = "output/plot"

    # Init everything
    for cell in grid.cells
        @views interpolate_initial_dofs(eq, scenario, grid.dofs[:,:,cell.dataidx],cell,grid.basis)
    end
    plotter = VTKPlotter(eq, scenario, grid, filename)

    grid.time = 0
    timestep = 0
    next_plotted = config.plot_start

    while grid.time < config.end_time
        if timestep > 0

            time_start = time()
            if config.equation_type == "non_linear"
                minmod_limiter(eq, scenario, grid.dofs, grid, grid.basis, get_ndofs(eq), 2)
            end
            dt = 1/(config.order^2+1) * config.cellsize[1] * config.courant * 1/grid.maxeigenval
            # Only step up to either end or next plotting
            dt = min(dt, next_plotted-grid.time, config.end_time - grid.time)
            @assert dt > 0

            @info "Running timestep" timestep dt grid.time
            step(integrator, grid, dt) do du, dofs, time
                # slope limiting before each timestep
                if config.equation_type == "non_linear"
                    minmod_limiter(eq, scenario, grid.dofs, grid, grid.basis, get_ndofs(eq), 2)
                end
                evaluate_rhs(eq, scenario, filter, globals, du, dofs, grid, config)
                #Here is the end of the PROBLEM A of Fractional Step method for a problem with source

                #Now we need to solve the PROBLEM B of Fractional Step Method for a problem with source
                #We will use Forward Euler for this (order=1)
                if config.source == "with_source"

                    #@info "from outside"
                    #@info sp[1](1,2,3)
                    for i in eachindex(grid.cells)
                        @views cell = grid.cells[i]
                        @views data = dofs[:,:, cell.dataidx]
                        sp = Array{Function,1}
                        globalx,globaly = globalposition(cell, (0.5, 0.5))
                        sp = evaluate_source(eq, grid.time, (globalx,globaly))
                        # Finding global Coordinates of Cell
#Calculation is specifically for order=1 and ndof=5

                        #@info globalx
                        #@info globaly
                        #globalz = 0.0 # as problem is 2 Dimensional
                        #@info size(data)
                        #@info data

#Calculation is specifically for order=1 and ndof=5
                        data[1,1] = data[1,1] + ( dt* (sp[1](globalx,globaly,grid.time)[1][1]) )
                        data[1,2] = data[1,2] + ( dt* (sp[2](globalx,globaly,grid.time)[1][1]) )
                        data[1,3] = data[1,3] + ( dt* (sp[3](globalx,globaly,grid.time)[1][1]) )
                        data[1,4] = data[1,4] + ( dt* (sp[4](globalx,globaly,grid.time)[1][1]) )
                        data[1,5] = data[1,5] + ( dt* (sp[5](globalx,globaly,grid.time)[1][1]) )

                        #@info data
                    end

                end
            end
            grid.time += dt
            time_end = time()
            time_elapsed = time_end - time_start
            @info "Timestep took" time_elapsed
        else
            # Compute initial eigenvalue (needed for dt)
            grid.maxeigenval = -1
            for cell in grid.cells
                @views celldata = grid.dofs[:,:,cell.dataidx]
                for normalidx=1:2
                    cureigenval = max_eigenval(eq, celldata, normalidx)
                    grid.maxeigenval = max(grid.maxeigenval, cureigenval)
                end
            end
        end

        if abs(grid.time - next_plotted) < 1e-10
            # slope limiting before plotting
            @info "Writing output" grid.time
            if config.equation_type == "non_linear"
                minmod_limiter(eq, scenario, grid.dofs, grid, grid.basis, get_ndofs(eq), 2)
            end
            plot(plotter)
            next_plotted = grid.time + config.plot_step
        end
        #break # For sibson to run only for t=0
        timestep += 1
    end
    save(plotter)
    if is_analytical_solution(eq, scenario)
        evaluate_error(eq, scenario, grid, grid.time)
    end
end

end


text:
nprocs()==1 && addprocs()

@everywhere logis(x) = 1 ./ (1 .+ exp(-x))
@everywhere dlogis(x) = logis(x) .* (1 .- logis(x))

@everywhere function initParams(x,y,n_hidden,max_iter)
    w1 = randn(n_hidden,size(x,1))
    w2 = randn(n_hidden,n_hidden)
    w3 = randn(size(y,1),n_hidden)
    B1 = rand(n_hidden,1)
    B2 = rand(n_hidden,1)
    errors = zeros(max_iter)
    outputs = zeros(size(y,2),max_iter)
    return (w1,w2,w3,errors,outputs,B1,B2)
end

@everywhere function ffward(x,y,w1,w2,w3)
    a1 = w1 * x
    h1 = logis(a1)
    a2 = w2 * h1
    h2 = logis(a2)
    ay = w3 * h2
    yhat = logis(ay)
    err = yhat - y
    return(a1,h1,a2,h2,ay,yhat,err)
end

@everywhere function updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err;lr=1)
    w1 .-= δa1 * x' .* lr
    w2 .-= δa2 * h1' .* lr
    w3 .-= err * h2' .* lr
end

@everywhere function bp(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,_,_)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa2 = (w3' * err) .* dlogis(a2)
        δa1 = (w2' * δa2) .* dlogis(a1)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

@everywhere function fa(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,B1,B2)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa2 = (B2.*err) .* dlogis(a2)
        δa1 = (B1.*δa2) .* dlogis(a1)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

@everywhere function dfa(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,B1,B2)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa2 = (B2.*err) .* dlogis(a2)
        δa1 = (B1.*err) .* dlogis(a1)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

@everywhere function ifa(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,B1,_)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa1 = (B1.*err) .* dlogis(a1)
        δa2 = (w2 * δa1) .* dlogis(a2)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

# Example: XOR problem

x = [1 1 ; 0 1 ; 1 0 ; 0 0]'
y = [0 1 1 0]
n_hidden = 15; # Number of hidden units
max_iter = 500; #Number of learning steps
trials = 30;

dfa_err=@parallel (hcat) for _=1:trials
    (_,dfa_errors) = dfa(x,y,n_hidden,max_iter)
    dfa_errors
end
(dfa_yhat,_) = dfa(x,y,n_hidden,max_iter)

bp_err = @parallel (hcat) for _=1:trials
    (_,bp_errors) = bp(x,y,n_hidden,max_iter)
    bp_errors
end
(bp_yhat,_) = bp(x,y,n_hidden,max_iter)

ifa_err=@parallel (hcat) for _=1:trials
    (_,ifa_errors) = ifa(x,y,n_hidden,max_iter)
    ifa_errors
end
(ifa_yhat,_) = ifa(x,y,n_hidden,max_iter);

fa_err=@parallel (hcat) for _=1:trials
    (_,fa_errors) = fa(x,y,n_hidden,max_iter)
    fa_errors
end
(fa_yhat,_) = fa(x,y,n_hidden,max_iter);

using Plots
l = @layout [
    a b
    c d
    e f
]
gr()
p1=plot(dfa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="DFA")
p2=plot(ifa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="IFA")
p3=plot(fa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="FA")
p4=plot(bp_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="BP")
p5=plot([dfa_err[:,1] ifa_err[:,1] fa_err[:,1] bp_err[:,1]],label=["DFA" "IFA" "FA" "BP"],title="First Trial Error")
p6=plot([mean(dfa_err,2) mean(ifa_err,2) mean(fa_err,2) mean(bp_err,2)],label=["DFA" "IFA" "FA" "BP"],title="Mean Trials Error")
plot(p1,p2,p3,p4,p5,p6,size=(900,900),layout=l)

ext: jl


text:
# A few headers displays
function print_header(Title)

    @printf(" \n\n %-15s\n", Title)
    @printf(
        "\n%-16s  %5s  %9s  %7s %7s  %5s  %5s  %5s  %6s  %s   %s\n",
        "Name",
        "nvar",
        "f",
        "‖∇f‖∞",
        "‖∇f‖₂",
        "#obj",
        "#grad",
        " H  ",
        "#hprod",
        "status",
        "time"
    )
end

function display_header_problems()
    @printf(
        "%-15s  %5s  %9s  %7s %7s  %5s  %5s  %6s  %s\n",
        "Name",
        "nvar",
        "f",
        "‖∇f‖∞",
        "‖∇f‖₂",
        "#obj",
        "#grad",
        "#hprod",
        "status"
    )
end

function display_header_iterations()
    @printf(" Iter    f          ||g||       λ        α       status \n",)
end

# A few displays for verbose iterations
function display_failure(iter, fnext, λ, α)
    @printf("%4d  %10.3e             %7.1e %7.1e    unsuccessful\n", iter, fnext, λ, α)
end

function display_v_success(iter, f, norm_g, λ, α)
    @printf("%4d  %10.3e %9.2e   %7.1e %7.1e Very successful\n", iter, f, norm_g, λ, α)
end

function display_success(iter, f, norm_g, λ, α)
    @printf("%4d  %10.3e %9.2e   %7.1e %7.1e      successful\n", iter, f, norm_g, λ, α)
end

function print_stats(prob, dim, f, gNorm, gnorm2, calls, status, timt)
    @printf(
        "%-16s  %5d  %9.2e  %7.1e  %7.1e  %5d  %5d %5d %6d  %s  %8.3f\n",
        prob,
        dim,
        f,
        gNorm,
        gnorm2,
        calls[1],
        calls[2],
        calls[3],
        calls[4],
        status,
        timt
    )
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    TableDistance(; normalize=true, weights=nothing)

Distance between rows of Tables.jl tables.

## Options

* `normalize` - whether or not to normalize the columns
                of the tables before computing distances
                (default to `true`)
* `weights`   - dictionary with weights for each column
                (default to uniform weights)

## Example

```julia
julia> pairwise(TableDistance(), table₁, table₂)
```
"""
struct TableDistance{W}
  normalize::Bool
  weights::W
end

TableDistance(; normalize=true, weights=nothing) =
  TableDistance{typeof(weights)}(normalize, weights)

function pairwise(td::TableDistance, table₁, table₂)
  ts, ds, ws = preprocess(td, table₁, table₂)

  function f((c, d))
    x = Tables.getcolumn(ts[1], c)
    y = Tables.getcolumn(ts[2], c)
    ws[c] * pairwise(d, x, y)
  end

  mapreduce(f, +, ds)
end

function pairwise(td::TableDistance, table)
  ts, ds, ws = preprocess(td, table)

  function f((c, d))
    x = Tables.getcolumn(ts[1], c)
    ws[c] * pairwise(d, x)
  end

  mapreduce(f, +, ds)
end

function preprocess(td, tables...)
  # retrieve parameters
  distances = [default_distances(table) for table in tables]
  weights   = isnothing(td.weights) ? default_weights(first(tables)) : td.weights

  # sanity checks
  @assert length(unique(distances)) == 1 "incompatible schema"
  @assert keys(weights) == keys(first(distances)) "invalid columns in weights"
  @assert all(values(weights) .> 0) "weights must be positive"

  # normalize weights
  c, w = keys(weights), values(weights)
  weights = Dict(c .=> w ./ sum(w))

  # normalize columns
  stdtables = td.normalize ? normalize(tables...) : tables

  stdtables, first(distances), weights
end


text:
function kRPCConnect(client_name::String)
    return kRPCConnect("localhost", 50000, 50001, client_name)
end

function kRPCConnect(port::Int64, client_name::String)
    return kRPCConnect("localhost", port, port+1, client_name)
end

function connect_or_error(conn, cr :: krpc.schema.ConnectionRequest)
    SendRawProto(conn, cr)
    resp = readproto(RecvRawProto(conn), krpc.schema.ConnectionResponse())
    if resp.status == krpc.schema.ConnectionResponse_Status.MALFORMED_MESSAGE
        throw("Malfomed Initial Message, Check Version")
    elseif resp.status == krpc.schema.ConnectionResponse_Status.TIMEOUT
        throw("Timeout, check network")
    elseif resp.status == krpc.schema.ConnectionResponse_Status.WRONG_TYPE
        throw("Wrong type, yell at kRPC.jl maintainer")
    end
    return resp
end

function generateStubs(conn)
	cachepath = joinpath(@__DIR__, "..", "deps")
	genpath = joinpath(cachepath, "gen.jl")
	if isfile(genpath)
		kRPC.eval(:(include($genpath)))
	else
		mkpath(cachepath)
	    info = GetServices(conn)
	    gennex = Codegen.generateHelpers(info, conn)
        open(genpath, "w") do f
            println(f, gennex)
        end
        include_dependency(genpath)
	    kRPC.eval(gennex)
	end
end

function kRPCConnect(client_name::String, host::String="localhost", port::Int64=50000, stream_port::Int64=50001)
    if lastindex(client_name) > 31
        throw(ArgumentError("client name too long"))
    end

    conn = connect(host, port)
    str_conn = connect(host, stream_port)

    resp = connect_or_error(conn, krpc.schema.ConnectionRequest(client_name=client_name,_type=krpc.schema.ConnectionRequest_Type.RPC))
    connect_or_error(str_conn, krpc.schema.ConnectionRequest(client_identifier=resp.client_identifier,_type=krpc.schema.ConnectionRequest_Type.STREAM))

    stream_running = Bool[true]
    conn = kRPCConnection(conn, str_conn, resp.client_identifier, (@async stream_listener(str_conn, stream_running)), stream_running)
    if !isdefined(kRPC, :Remote)
        generateStubs(conn)
    end

    kRPC.eval(:(Remote.change_connection($conn))) # whyyyyy

    return conn
end

function kRPCConnect(f::Function, client_name::String, host::String="localhost", port::Int64=50000, stream_port::Int64=50001)
    kc = kRPCConnect(client_name, host, port, stream_port)
    try
        f(kc)
    finally
        close(kc)
    end
end

function SendMessage(conn::kRPCConnection, pcall::kPC)
    res = SendBiMessage(conn, krpc.schema.Request(calls=[pcall.pcall]))
    handle_potential_error(res)
    if pcall.rty != Nothing
        return getJuliaValue(res.results[1].value,pcall.rty)
    end
end

function SendMessage(conn::kRPCConnection, pcalls::AbstractArray{kPC})
    res = SendBiMessage(conn, krpc.schema.Request(calls=map(p->p.pcall, pcalls)))
    handle_potential_error(res)
    return [if pcalls[i].rty != Nothing getJuliaValue(res.results[i].value, pcalls[i].rty) else Nothing() end for i=1:length(pcalls)]
end

function close(conn::kRPCConnection)
    conn.running[1] = false
    Base.close(conn.conn)
end

function handle_potential_error(result)
    if ProtoBuf.has_field(result, :error)
        throw(kRPC.make_error(result.error))
    end
end

function make_error(err::krpc.schema.Error)
    return "A remote RPC error occured.\nThe error occured in $(err.service).$(err.name):\n$(err.description)\n$(err.stack_trace)"
end


text:
"Problem 6 in the Hock-Schittkowski suite"
function hs6()

  nlp = Model()

  @variable(nlp, x[1:2])
  setvalue(x[1], -1.2)
  setvalue(x[2],  1.0)

  @NLobjective(
    nlp,
    Min,
    (1 - x[1])^2
  )

  @NLconstraint(
    nlp,
    10 * (x[2] - x[1]^2) == 0
  )

  return nlp
end


text:
permittivity(w) = @error "permittivity function is not set"
permeability(w) = @error "permeability function is not set"

function refractive_index(w)
    eps = permittivity(abs(w))
    mu = permeability(abs(w))
    n = sqrt(eps * mu + 0im)
    return n
end

function beta_func(w)
    n = refractive_index(w)
    if w < 0
        n = conj(n)
    end
    return n * w / C0
end

function k_func(w)
    beta = beta_func(w)
    return real(beta)
end

function k1_func(w)
    func(w) = k_func(w)
    return derivative(func, w, 1)
end

function k2_func(w)
    func(w) = k_func(w)
    return derivative(func, w, 2)
end

function k3_func(w)
    func(w) = k_func(w)
    return derivative(func, w, 3)
end

function absorption_coefficient(w)
    beta = beta_func(w)
    return imag(beta)
end

function phase_velocity(w)
    n = refractive_index(w)
    return C0 / real(n)
end

function group_velocity(w)
    k1 = k1_func(w)
    return 1 / k1
end

function diffraction_length(w, a0)
    k = k_func(w)
    return k * a0^2
end

function dispersion_length(w, t0)
    k2 = k2_func(w)
    if k2 == 0
        Ldisp = Inf
    else
        Ldisp = t0^2 / abs(k2)
    end
    return Ldisp
end

function dispersion_length3(w, t0)
    k3 = k3_func(w)
    if k3 == 0
        Ldisp3 = Inf
    else
        Ldisp3 = t0^3 / abs(k3)
    end
    return Ldisp3
end

function absorption_length(w)
    ga = absorption_coefficient(w)
    if ga == 0
        La = Inf
    else
        La = ga / 2
    end
    return La
end

function chi1_func(w)
    eps = permittivity(abs(w))
    return eps - 1
end

function chi3_func(w, n2)
    n = refractive_index(w)
    return 4 / 3 * real(n)^2 * EPS0 * C0 * n2
end

function critical_power(w, n2)
    Rcr = 3.79
    lam = 2 * pi * C0 / w
    n = refractive_index(w)
    return Rcr * lam^2 / (8 * pi * abs(real(n)) * abs(real(n2)))
end

function nonlinearity_length(w, I0, n2)
    if n2 == 0
        Lnl = Inf
    else
        Lnl = 1 / (abs(real(n2)) * I0 * w / C0)
    end
    return Lnl
end

"""Self-focusing distance by the Marburger formula (P in watts)."""
function selffocusing_length(w, a0, P, n2)
    Ld = diffraction_length(w, a0)
    PPcr = P / critical_power(w, n2)
    if PPcr > 1
        zf = 0.367 * Ld / sqrt((sqrt(PPcr) - 0.852)^2 - 0.0219)
    else
        zf = Inf
    end
    return zf
end

"""
N-th derivative of a function f at a point x.

The derivative is found using five-point stencil:
    http://en.wikipedia.org/wiki/Five-point_stencil
Additional info:
    http://en.wikipedia.org/wiki/Finite_difference_coefficients
"""
function derivative(f, x, n)
    if x == 0
        h = 0.01
    else
        h = 0.001 * x
    end

    if n == 1
        res = (f(x - 2 * h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2 * h)) /
              (12 * h)
    elseif n == 2
        res = (- f(x - 2 * h) + 16 * f(x - h) - 30 * f(x) + 16 * f(x + h) -
                 f(x + 2 * h)) / (12 * h^2)
    elseif n == 3
        res = (- f(x - 2 * h) + 2 * f(x - h) - 2 * f(x + h) +
                 f(x + 2 * h)) / (2 * h^3)
    elseif n == 4
        res = (f(x - 2 * h) - 4 * f(x - h) + 6 * f(x) - 4 * f(x + h) +
               f(x + 2 * h)) / (h^4)
    else
        error("Wrong derivative order.")
    end
    return res
end


text:
# Functions for changing the values of an associative array

#=
Reduce all value to logical, checking if that cell is empty
=#

using SparseArrays,LinearAlgebra

logical(A::Assoc) = Assoc(copy(A.row),copy(A.col),promote([1.0],A.val)[1],LinearAlgebra.fillstored!(dropzeros!(copy(A.A)),1))

function parsehelper(type, v)
    parsed = tryparse(type, v)
    if parsed===nothing
        return 1
    end
    return parsed
end

function str2num(A::Assoc, type::Type=Int)
    if ~(type<:Number)
        println("Type not number")
    return nothing
    end
    r,c,v = find(A)
    v = parsehelper.(type, v)
    # TODO do mixed types of values exist? If so, use them
    A = Assoc(r,c,v)
end

function str2num(type::Type, A::Assoc)
    return str2num(A, type)
end

# Converts numeric values to strings
function num2str(A::Assoc)::Assoc
	r, c, v = find(A)
	v1 = string.(v)
	return Assoc(r, c, v1)
end

#=
convert all entries to a certain type
=#

function convertvals(type, A::Assoc)
    r,c,v = find(A)
    v = convert.(type, v)
    A = Assoc(r,c,v)
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using IterativeSolvers
struct NewtonResults
    z::Vector # z = [x1;u1;...;xN]
    λ::Vector # constraint multipliers
    ν::Vector # dynamics constraints multipliers
    s::Vector # slack variables

    ∇²J::SparseMatrixCSC # Hessian of cost function (including penalty terms)
    ∇J::Vector # gradient of cost function (including penalty terms)

    ∇c::SparseMatrixCSC # Jacobian of constraints
    c::Vector # constraints

    ∇d::SparseMatrixCSC # Jacobian of dynamics constraints
    d::Vector # dynamics constraints

    A::SparseMatrixCSC # KKT Hessian matrix
    b::Vector # KKT gradient matrix

    active_set::Vector
    active_set_ineq::Vector

    x_eval::Vector

    δ::Vector

    ρ::Vector

    Nz::Int
    Np::Int
    Nx::Int
    NpI::Int
end

function NewtonResults(Nz::Int,Np::Int,Nx::Int,NpI::Int)
    z = zeros(Nz)
    λ = zeros(Np)
    ν = zeros(Nx)
    s = zeros(Np)

    ∇²J = spzeros(Nz,Nz)
    ∇J = zeros(Nz)

    ∇c = spzeros(Np,Nz)
    c = zeros(Np)

    ∇d = spzeros(Nx,Nz)
    d = zeros(Nx)

    A = spzeros(Nz+Np+Nx+Np,Nz+Np+Nx+Np)
    b = spzeros(Nz+Np+Nx+Np)

    active_set = zeros(Bool,Np)
    active_set_ineq = zeros(Bool,Np)
    x_eval = zeros(Nx)

    δ = zeros(Nz+Np+Nx+Np)

    ρ = zeros(1)

    Nz = Nz
    Np = Np
    Nx = Nx
    NpI = NpI

    NewtonResults(z,λ,ν,s,∇²J,∇J,∇c,c,∇d,d,A,b,active_set,active_set_ineq,x_eval,δ,ρ,Nz,Np,Nx,NpI)
end

function copy(r::NewtonResults)
    NewtonResults(copy(r.z),copy(r.λ),copy(r.ν),copy(r.s),copy(r.∇²J),copy(r.∇J),copy(r.∇c),copy(r.c),copy(r.∇d),copy(r.d),copy(r.A),copy(r.b),copy(r.active_set),copy(r.active_set_ineq),copy(r.x_eval),copy(r.δ),copy(r.ρ),copy(r.Nz),copy(r.Np),copy(r.Nx),copy(r.NpI))
end

function NewtonResults(solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    nm = nn + mm
    Nz = nn*N + mm*(N-1)
    Np = p*(N-1) + p_N
    Nx = N*nn
    NpI = pI*(N-1) + pI_N

    NewtonResults(Nz,Np,Nx,NpI)
end

function get_batch_sizes(solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    nm = nn + mm
    Nz = nn*N + mm*(N-1)
    Np = p*(N-1) + p_N
    Nx = N*n
    Nu = mm*(N-1) # number of control decision variables u
    NpI = pI*(N-1) + pI_N

    return Nz,Np,Nx,Nu,nm,NpI
end

function newton_active_set!(newton_results::NewtonResults,results::SolverIterResults,solver::Solver,tolerance::Float64=1.0e-3,slacks::Bool=true)
    # get problem dimensions
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
    Nz,Np,Nx,Nu,nm,NpI = get_batch_sizes(solver)

    active_set = newton_results.active_set
    active_set_ineq = newton_results.active_set_ineq

    s = newton_results.s

    active_set .= copy(vcat(results.active_set...))
    active_set_ineq .= 0 # reset

    for k = 1:N
        k != N ? p_idx = p : p_idx = p_N
        k != N ? pI_idx = pI : pI_idx = pI_N

        for i = 1:p_idx
            idx = (k-1)*p + i

            if i <= pI_idx
                if results.C[k][i] > -tolerance
                    active_set[idx] = 1
                    active_set_ineq[idx] = 1
                    s[idx] = sqrt(2.0*max(0.0,-results.C[k][i]))
                end
            # else
            #     active_set[idx] = 1
            end
        end
    end
end

function update_newton_results!(newton_results::NewtonResults,results::SolverIterResults,solver::Solver)
    # get problem dimensions
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
    Nz,Np,Nx,Nu,nm,NpI = get_batch_sizes(solver)

    # pull out results for convenience
    X = results.X
    U = results.U
    Iμ = results.Iμ
    C = results.C
    Cx = results.Cx
    Cu = results.Cu
    fdx = results.fdx
    fdu = results.fdu
    x0 = solver.obj.x0

    # pull out newton results for convenience
    z = newton_results.z
    λ = newton_results.λ
    ν = newton_results.ν
    s = newton_results.s

    ∇²J = newton_results.∇²J
    ∇J = newton_results.∇J

    ∇c = newton_results.∇c
    c = newton_results.c

    ∇d = newton_results.∇d
    d = newton_results.d

    active_set = newton_results.active_set
    active_set_ineq = newton_results.active_set_ineq

    x_eval = newton_results.x_eval

    Inn = sparse(Matrix(I,nn,nn))

    newton_results.ρ[1] = maximum(vcat(results.μ...))
    ρ = newton_results.ρ[1]

    # update batch matrices
    for k = 1:N
        x = results.X[k]

        if k != N
            u = results.U[k]
            Q,R,H,q,r = taylor_expansion(solver.obj.cost,x,u) .* solver.dt
        else
            Qf,qf = taylor_expansion(solver.obj.cost,x)
        end

        # Indices
        k != N ? idx = (((k-1)*nm + 1):k*nm) : idx = (((k-1)*nm + 1):Nz) # index over x and u
        k != N ? idx2 = (((k-1)*p + 1):k*p) : idx2 = (((k-1)*p + 1):Np) # index over p
        idx3 = ((k-1)*nm + 1):((k-1)*nm + nn) # index over x only
        idx4 = ((k-1)*nm + nn + 1):k*nm # index over u only
        idx5 = ((k-1)*mm + 1):k*mm # index through stacked u vector
        k != N ? idx6 = (((k-1)*p + 1):((k-1)*p + pI)) : idx6 = (((k-1)*p + 1):((k-1)*p + pI_N))# index over p only inequality indices
        idx7 = ((k-1)*nn + 1):(k*nn)
        idx8 = ((k-2)*nm + 1):((k-1)*nm + nn)
        idx9 = ((k-2)*nm + 1):((k-2)*nm + nn)
        idx10 = ((k-2)*nm + nn + 1):((k-1)*nm)

        k != N ? pI_idx = pI : pI_idx = pI_N

        as = active_set[idx2]

        # Assemble ∇²J, ∇J, ∇c, c, z
        if k != N
            ∇²J[idx3,idx3] = Q + ρ*Cx[k][as,:]'*Cx[k][as,:]
            ∇²J[idx3,idx4] = H' + ρ*Cx[k][as,:]'*Cu[k][as,:]
            ∇²J[idx4,idx3] = H + ρ*Cu[k][as,:]'*Cx[k][as,:]
            ∇²J[idx4,idx4] = R + ρ*Cu[k][as,:]'*Cu[k][as,:]

            ∇J[idx3] = q #+ ρ*Cx[k][as,:]'*C[k][as]
            ∇J[idx4] = r #+ ρ*Cu[k][as,:]'*C[k][as]

            ∇c[idx2,idx3] = Cx[k]
            ∇c[idx2,idx4] = Cu[k]

            c[idx2] = C[k]

            z[idx] = [x;u]
        else
            ∇²J[idx,idx] = Qf #+ ρ*Cx[N][as,:]'*Cx[N][as,:]
            ∇J[idx] = qf #+ ρ*Cx[N][as,:]'C[N][as,:]

            ∇c[idx2,idx] = results.Cx[N]
            c[idx2] = results.C[N]

            z[idx] = x
        end

        # assemble ∇d, d
        if k == 1
            ∇d[1:nn,1:nn] = Inn
            d[1:nn] = X[1] - solver.obj.x0
        else
            ∇d[idx7,idx9] = -fdx[k-1]
            ∇d[idx7,idx10] = -fdu[k-1]
            ∇d[idx7,idx3] = Inn

            solver.fd(view(x_eval,idx7),X[k-1][1:n],U[k-1][1:m])
            d[idx7] = X[k] - x_eval[idx7]
        end

        # assemble λ, ν
        λ[idx2] = results.λ[k]
        ν[idx7] = results.s[k]

    end

    return nothing
end

function update_results_from_newton_results!(results::SolverIterResults,newton_results::NewtonResults,solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    nm = nn + mm
    Nz = nn*N + mm*(N-1)
    Np = p*(N-1) + p_N
    Nu = mm*(N-1) # number of control decision variables u

    z = newton_results.z
    λ = newton_results.λ

    # update results with stack vector
    for k = 1:N
        idx = ((k-1)*nm + 1):k*nm # index over x and u
        k != N ? idx2 = (((k-1)*p + 1):k*p) : idx2 = (((k-1)*p + 1):Np) # index over p
        k != N ? idx3 = (((k-1)*nm + 1):((k-1)*nm + nn)) : idx3 = (((k-1)*nm + 1):Nz)# index over x only
        idx4 = ((k-1)*nm + nn + 1):k*nm # index over u only
        idx5 = ((k-1)*mm + 1):k*mm # index through stacked u vector
        idx6 = ((k-1)*p + 1):((k-1)*p + pI) # index over p only inequality indices

        results.X[k] = z[idx3]
        k != N ? results.U[k] = z[idx4] : nothing
        results.λ[k] = λ[idx2]
    end

    update_constraints!(results,solver)
    update_jacobians!(results,solver)

    return nothing
end

function solve_kkt!(newton_results::NewtonResults)
    z = newton_results.z
    λ = newton_results.λ
    ν = newton_results.ν
    s = newton_results.s

    ∇²J = newton_results.∇²J
    ∇J = newton_results.∇J

    ∇c = newton_results.∇c
    c = newton_results.c

    ∇d = newton_results.∇d
    d = newton_results.d

    A = newton_results.A
    b = newton_results.b

    active_set = newton_results.active_set
    active_set_ineq = newton_results.active_set_ineq

    δ = newton_results.δ
    ρ = newton_results.ρ[1]

    # get batch problem sizes
    Nz = newton_results.Nz
    Np = newton_results.Np
    Nx = newton_results.Nx

    Np_as = sum(active_set)

    # assemble KKT matrix,vector
    _idx1 = Array(1:Nz)
    _idx2 = Array((1:Np) .+ Nz)[active_set]
    _idx22 = Array((1:Np) .+ Nz)[active_set_ineq]
    _idx3 = Array((1:Nx) .+ (Nz + Np))
    _idx4 = Array((1:Np) .+ (Nz + Np + Nx))[active_set_ineq]

    _idx5 = [_idx1;_idx2;_idx3;_idx4]

    A[_idx1,_idx1] = ∇²J
    A[_idx1,_idx2] = ∇c[active_set,:]'
    A[_idx1,_idx3] = ∇d'
    A[_idx2,_idx1] = ∇c[active_set,:]
    A[_idx3,_idx1] = ∇d

    A[_idx2,_idx4] = Diagonal(s)[active_set_ineq,active_set]'
    A[_idx4,_idx2] = Diagonal(s)[active_set_ineq,active_set]
    A[_idx4,_idx4] = Diagonal(λ[active_set_ineq]) + ρ*Diagonal(c[active_set_ineq]) + 3/2*ρ*Diagonal(s[active_set_ineq].^2)

    A[_idx1,_idx4] = ρ*∇c[active_set_ineq,:]'*Diagonal(s[active_set_ineq])
    A[_idx4,_idx1] = ρ*Diagonal(s[active_set_ineq])*∇c[active_set_ineq,:]

    b[_idx1] = -(∇J + ∇c[active_set,:]'*λ[active_set] + ∇d'*ν + ρ*∇c[active_set_ineq,:]'*0.5*s[active_set_ineq].^2 + ρ*∇d'*d)
    b[_idx2] = -c[active_set]
    b[_idx22] -= 0.5*s[active_set_ineq].^2
    b[_idx3] = -d
    b[_idx4] = -(λ[active_set_ineq] .* s[active_set_ineq] + ρ*c[active_set_ineq] .* s[active_set_ineq] + 0.5*ρ*s[active_set_ineq].^3)

    ## indexing
    println("A:")
    println("rank: $(rank(Array(A)))")
    println("size: $(size(A))")
    println("cond: $(cond(Array(A)))")
    println("Active set A:")
    println("rank: $(rank(Array(A[_idx5,_idx5])))")
    println("size: $(size(A[_idx5,_idx5]))")
    println("cond: $(cond(Array(A[_idx5,_idx5])))")
    println("Nz: $Nz")
    println("Np: $Np")
    println("Nx: $Nx")
    println("active set: $(sum(active_set))")
    println("active set (ineq.): $(sum(active_set_ineq))")
    # println(rank(Array(A[[_idx1;_idx2;_idx3],[_idx1;_idx2;_idx3]])))
    # println(length([_idx1;_idx2;_idx3]))
    # println("No slack A:")
    # println(rank(Array(A[_idx_tmp,_idx_tmp])))
    # println(size(A[_idx_tmp,_idx_tmp]))
    # println("slacks")
    # println(sum(active_set))
    # println("reg.")
    # println(rank(Array(A)+Matrix(I,Nz+Np+Nx+Np,Nz+Np+Nx+Np)))
    # println("KKT size: $(length(_idx5))")

    # δ[_idx5] = A[_idx5,_idx5]\b[_idx5]
    δ[_idx5] = minres(A[_idx5,_idx5],b[_idx5])

    # z .+= alpha*δ[_idx1]
    # λ[active_set] += alpha*δ[_idx2]
    # ν .+= alpha*δ[_idx3]
    # s[active_set_ineq] += alpha*δ[_idx4]

    return nothing
end

function newton_step!(results::SolverIterResults,newton_results::NewtonResults,solver::Solver,alpha::Float64=1.0)
    update_newton_results!(newton_results,results,solver)
    solve_kkt!(newton_results)

    ###
    _idx1 = Array(1:newton_results.Nz)
    _idx2 = Array((1:newton_results.Np) .+ newton_results.Nz)[newton_results.active_set]
    _idx3 = Array((1:newton_results.Nx) .+ (newton_results.Nz + newton_results.Np))
    _idx4 = Array((1:newton_results.Np) .+ (newton_results.Nz + newton_results.Np + newton_results.Nx))[newton_results.active_set_ineq]

    newton_results.z .+= alpha*newton_results.δ[_idx1]
    newton_results.λ[newton_results.active_set] += alpha*newton_results.δ[_idx2]
    newton_results.λ[newton_results.active_set_ineq] = max.(0,newton_results.λ[newton_results.active_set_ineq])
    newton_results.ν .+= alpha*newton_results.δ[_idx3]
    newton_results.s[newton_results.active_set_ineq] += alpha*newton_results.δ[_idx4]
    ###

    update_results_from_newton_results!(results,newton_results,solver)
    newton_active_set!(newton_results,results,solver)
    return nothing
end

function newton_solve!(results::SolverIterResults,solver::Solver)
    # copy current results
    results = copy(results)
    results_new = copy(results)

    # instantiate Newton Results
    newton_results = NewtonResults(solver)

    # get initial cost and max constraint violation
    newton_active_set!(newton_results,results_new,solver)
    update_newton_results!(newton_results,results_new,solver)
    # J_prev = cost(solver,results_new)
    J_prev = newton_cost(results_new,newton_results,solver)

    c_max = max_violation(results_new)

    println("Newton initialization")
    println("Cost: $J_prev")
    println("c_max: $c_max")
    println("\n")

    α = 1.0
    max_iter = 100
    max_c = 1e-8
    iter = 1
    ls_param = 0.1
    J = Inf

    while (c_max > max_c) && iter <= max_iter
        newton_step!(results_new,newton_results,solver,α)
        # backwardpass!(results_new,solver)
        # rollout!(results_new,solver,1.0)
        # results_new.X .= deepcopy(results_new.X_)
        # results_new.U .= deepcopy(results_new.U_)
        # update_constraints!(results_new,solver)
        # J = cost(solver,results_new)
        J = newton_cost(results_new,newton_results,solver)
        c_max = max_violation(results_new)

        println("*Newton step: $iter")
        println("α = $α")
        println("Cost: $J \n    prev. Cost: $J_prev")
        println("ΔJ = $(J-J_prev)")
        println("c_max: $c_max")

        # iter == 1 ? J_prev = J : nothing

        if J <= J_prev + ls_param*newton_results.b'*α*newton_results.δ
            results = copy(results_new)
            J_prev = copy(J)
            println("improved")
            α = 1.0
        else
            println("α = $α increased cost")
            results_new = copy(results)
            newton_results = NewtonResults(solver)
            newton_active_set!(newton_results,results_new,solver)
            update_newton_results!(newton_results,results_new,solver)
            α = 0.5*α
        end
        println("-----\n")
        iter += 1
    end
    println("Newton solve complete")
    println("J: $J")
    println("c_max: $c_max")
end

function newton_cost(results::SolverIterResults,newton_results::NewtonResults,solver::Solver)
    ρ = newton_results.ρ[1]
    J = cost(solver,to_array(results.X),to_array(results.U))
    J += newton_results.λ'*(vcat(results.C...)+0.5*newton_results.s.^2)
    J += ρ*(vcat(results.C...) + 0.5*newton_results.s.^2)'*(vcat(results.C...) + 0.5*newton_results.s.^2)

    #########

    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
    Nz,Np,Nx,Nu,nm,NpI = get_batch_sizes(solver)

    # pull out results for convenience
    X = results.X
    U = results.U
    Iμ = results.Iμ
    C = results.C
    Cx = results.Cx
    Cu = results.Cu
    fdx = results.fdx
    fdu = results.fdu
    x0 = solver.obj.x0

    ∇d = newton_results.∇d
    d = newton_results.d

    for k = 1:solver.N

        k != N ? idx = (((k-1)*nm + 1):k*nm) : idx = (((k-1)*nm + 1):Nz) # index over x and u
        k != N ? idx2 = (((k-1)*p + 1):k*p) : idx2 = (((k-1)*p + 1):Np) # index over p
        idx3 = ((k-1)*nm + 1):((k-1)*nm + nn) # index over x only
        idx4 = ((k-1)*nm + nn + 1):k*nm # index over u only
        idx5 = ((k-1)*mm + 1):k*mm # index through stacked u vector
        k != N ? idx6 = (((k-1)*p + 1):((k-1)*p + pI)) : idx6 = (((k-1)*p + 1):((k-1)*p + pI_N))# index over p only inequality indices
        idx7 = ((k-1)*nn + 1):(k*nn)
        idx8 = ((k-2)*nm + 1):((k-1)*nm + nn)
        idx9 = ((k-2)*nm + 1):((k-2)*nm + nn)
        idx10 = ((k-2)*nm + nn + 1):((k-1)*nm)

        k != N ? pI_idx = pI : pI_idx = pI_N

        if k == 1
            # ∇d[1:nn,1:nn] = Inn
            d[1:nn] = X[1] - solver.obj.x0
        else
            # ∇d[idx7,idx9] = -fdx[k-1]
            # ∇d[idx7,idx10] = -fdu[k-1]
            # ∇d[idx7,idx3] = Inn

            solver.fd(view(newton_results.x_eval,idx7),X[k-1][1:n],U[k-1][1:m])
            d[idx7] = X[k] - newton_results.x_eval[idx7]
        end
    end
    ####

    J += newton_results.ν'*newton_results.d
    J += 0.5*ρ*newton_results.d'*newton_results.d
    return J
end

# """
# $(SIGNATURES)
#     Time Varying Discrete Linear Quadratic Regulator (TVLQR)
# """
# function lqr(A::Array, B::Array, Q::AbstractMatrix, R::AbstractMatrix, Qf::AbstractMatrix)::Array
#     n,m = size(B[1])
#     N_ = length(B)
#     K = [zeros(m,n) for k = 1:N_]
#     S = zeros(n,n)
#
#     # Boundary condition
#     S .= Qf
#
#     # Riccati
#     for k = 1:N_
#         # Compute control gains
#         K[k] = -(R + B[k]'*S*B[k])\(B[k]'*S*A[k])
#
#         # Calculate cost-to-go for backward propagation
#         S .= Q + A[k]'*S*A[k] - A[k]'*S*B[k]*K[k]
#     end
#     return K
# end
#
# function lqr(results::SolverResults, solver::Solver)::Array
#     Q, R, Qf = get_cost_matrices(solver)
#     A, B = results.fdx, results.fdu
#     return lqr(A,B,Q,R,Qf)
# end


text:
using FHist, StatsBase, Statistics
using Test

@testset "The basics" begin
    a = rand(10^3)
    sth1 = fit(Histogram, a)
    h1sqrt = Hist1D(sth1)
    h1p = Hist1D(sth1)
    @test h1sqrt.hist == sth1
    @test h1p.hist == sth1

    h1 = Hist1D(a)
    sth1 = fit(Histogram, a)
    @test h1.hist.weights ≈ sth1.weights
end

@testset "Fast route" begin
    a = rand(10^5)
    r1 = 0:0.1:1
    r2 = collect(r1)
    for r in (r1,r2)
        h1 = Hist1D(a, r)
        sth1 = fit(Histogram, a, r)
        @test all(h1.hist.weights .≈ sth1.weights)
    end
end

@testset "Slow route" begin
    a = rand(10^5)
    r1 = [0, 0.3, 0.4, 0.8, 1]
    h1 = Hist1D(a, r1)
    sth1 = fit(Histogram, a, r1)
    @test all(h1.hist.weights .≈ sth1.weights)
end

@testset "push! route" begin
    a = rand(10^5)
    r1 = [0, 0.3, 0.4, 0.8, 1]
    h1 = Hist1D(a, r1)
    h2 = Hist1D(Int; bins = r1)
    sth1 = fit(Histogram, a, r1)
    for ele in a
        push!(h2, ele)
    end

    @test all(h1.hist.weights .≈ sth1.weights)
    @test all(h1.hist.weights .≈ h2.hist.weights)
end

@testset "Weighted" begin
    a = rand(10^5)
    wgts1 = ones(10^5) |> weights
    wgts = rand(10^5) |> weights
    r1 = 0:0.1:1
    r2 = collect(r1)

    h1 = Hist1D(a)
    h2 = Hist1D(a, wgts1)
    @test h1==h2
    for r in (r1,r2)
        h1 = Hist1D(a, weights(wgts), r)
        sth1 = fit(Histogram, a, weights(wgts), r)
        @test all(h1.hist.weights .≈ sth1.weights)
    end
end

@testset "Hist2D" begin
    x = rand(10)
    y = rand(10)
    h = Hist2D((x,y))
    @test integral(h) == 10

    rx, ry = 0:0.1:1, 0:0.2:1
    wgts = weights(2*ones(length(x)))
    h = Hist2D((x,y), wgts, (rx, ry))
    @test integral(h) == sum(wgts)
    @test nbins(h) == (length(rx)-1, length(ry)-1)

    @test bincenters(Hist2D((x,y), (0:1,0:1))) == ([0.5], [0.5])
    @test bincenters(Hist2D((x,y), wgts, (0:1,0:1))) == ([0.5], [0.5])
    @test nbins(Hist2D((x,y), ([0,0.5,1],[0,0.5,1]))) == (2,2)
    @test nbins(Hist2D((x,y), ([0,0.3,1],[0,0.3,1]))) == (2,2)
    @test nbins(Hist2D((x,y), wgts, ([0,0.5,1],[0,0.5,1]))) == (2,2)
    @test nbins(Hist2D((x,y), wgts, ([0,0.3,1],[0,0.3,1]))) == (2,2)

    @test integral(Hist2D((x,y), wgts, nbins=(5,5))) == sum(wgts)
    @test integral(Hist2D((x,y), nbins=(5,5))) == length(x)

    h1 = Hist2D((x,y), wgts, (rx, ry))
    h2 = Hist2D(Float64; bins=(rx, ry))
    h3 = Hist2D(Float64; bins=(rx, ry))
    push!.(h2, x, y, wgts)
    atomic_push!.(h3, x, y, wgts)
    @test h1 == h2
    @test h1 == h3
end

@testset "Special bins" begin
    # integer values and integer binning
    a = floor.(Int,abs.(randn(10^6)))
    h1 = Hist1D(a, 0:5)
    h2 = Hist1D(fit(Histogram, a, 0:5))
    @test h1 == h2

    # integer values and integer binning
    a = floor.(Int,abs.(randn(10^6)))
    h1 = Hist1D(a, -6:2:6)
    h2 = Hist1D(fit(Histogram, a, -6:2:6))
    @test h1 == h2

    # test floating point rounding when
    # coercing explicit bins into StepRange
    bins = collect(-3:0.1:3.1)
    h = Hist1D(rand(100), bins)
    @test binedges(h) isa AbstractRange
    @test collect(binedges(h)) == bins
end

@testset "Normalize" begin
    a = rand(10^5)
    wgts1 = 2 .* ones(10^5) |> weights
    h1 = Hist1D(a, wgts1)
    @test integral(h1) ≈ sum(wgts1) atol=1e-8
    @test integral(normalize(h1)) ≈ 1 atol=1e-8

    h1 = Hist2D((a,a), wgts1, (0:0.1:1,0:0.1:1))
    @test integral(h1) ≈ sum(wgts1) atol=1e-8
    @test integral(normalize(h1)) ≈ 1 atol=1e-8
end

@testset "Cumulative" begin
    a = rand(10^5)
    wgts1 = 2 .* ones(10^5) |> weights
    h1 = Hist1D(a, wgts1)
    @test argmax(bincounts(cumulative(h1; forward=true))) == nbins(h1)
    @test argmax(bincounts(cumulative(h1; forward=false))) == 1
    h1 = Hist1D([1,2,2,3,3,3],1:4)
    @test bincounts(cumulative(h1, forward=true)) == [1, 3, 6]
    @test bincounts(cumulative(h1, forward=true)) == cumulative(h1, forward=true).sumw2
    @test bincounts(cumulative(h1, forward=false)) == [6, 5, 3]
    @test bincounts(cumulative(h1, forward=false)) == cumulative(h1, forward=false).sumw2
end

@testset "Bin errors" begin
    h1 = Hist1D(randn(100), -3:3)
    @test h1.sumw2 == binerrors(identity, h1)
    @test sqrt.(h1.sumw2) == binerrors(h1)
    h1 = Hist2D((randn(100), randn(100)), (-3:3,-3:3))
    @test h1.sumw2 == binerrors(identity, h1)
    @test sqrt.(h1.sumw2) == binerrors(h1)

    h2 = Hist1D(Float64; bins=0:1)
    push!(h2, 0.5, 0.5)
    push!(h2, 0.5, 0.5)

    @test h2.sumw2 == [0.5]
    @test binerrors(h2) == [sqrt(0.5)]
end

@testset "Statistics" begin
    r = -3:0.1:3
    N = 10^5
    h1 = Hist1D(randn(N), r)
    @test mean(h1) ≈ 0 atol=0.05
    @test std(h1) ≈ 1 atol=0.1
    @test nbins(h1) == length(r)-1
    @test median(h1) == quantile(h1, 0.5)
    @test quantile(h1, 0.0) == first(bincenters(h1))
    @test quantile(h1, 1.0) == last(bincenters(h1))

end

@testset "Lookup" begin
    h1 = Hist1D(randn(100))
    @test lookup.(Ref(h1), bincenters(h1)) == bincounts(h1)
    @test ismissing(lookup(h1, last(binedges(h1)) + 0.1))
    @test ismissing(lookup(h1, first(binedges(h1)) - 0.1))

    h1 = Hist2D((randn(100), randn(100)), (-3:3,-3:3))
    cx, cy = bincenters(h1)
    # (x,y) tuples of bin centers with same shape as counts
    tbc = ((x,y)->(x,y)).(cx, cy')
    f = x->lookup(h1,x...)
    @test f.(tbc) == bincounts(h1)
    @test ismissing(lookup(h1, 10, 10))
end

@testset "Sample" begin
    @test mean(FHist.sample(Hist1D(rand(10^5), 0:0.1:1), n=10^5)) ≈ 0.5 atol=0.1
end

@testset "Empty" begin
    a = rand(10^5)
    r = 0:0.1:1
    h1 = Hist1D(a, r)
    empty!(h1)
    @test maximum(h1.hist.weights) == 0
    @test maximum(h1.sumw2) == 0
    @test h1.hist.edges[1] == r

    h1 = Hist2D((randn(10),randn(10)), (-3:3,-3:3))
    empty!(h1)
    @test maximum(h1.hist.weights) == 0
    @test maximum(h1.sumw2) == 0
end

@testset "Atomic push" begin
    if get(ENV, "CI", "false") == "true"
        @test Threads.nthreads() > 1
    end
    a = randn(10^5)
    h1 = Hist1D(a, -3:0.2:3)
    h2 = Hist1D(Int; bins=-3:0.2:3)
    h3 = Hist1D(Int; bins=-3:0.2:3)
    h4 = Hist1D(Int; bins=-3:0.2:3)
    for i in a
        push!(h2, i)
        atomic_push!(h3, i)
    end
    Threads.@threads for i in a
        atomic_push!(h4, i)
    end
    @test h1 == h2
    @test h1 == h3
    @test h1 == h4
end

@testset "Broadcast" begin
    lookup(h::Hist1D, value) = h.hist.weights[FHist._edge_binindex(h.hist.edges[1], value)]

    h1 = Hist1D(Int; bins=0:3)

    push!.(h1, [0.5, 1.5])
    @test lookup.(h1, [0.5,1.5,2.5]) == [1, 1, 0]

    empty!(h1)
    push!.(h1, [0.5, 1.5], 2)
    @test lookup.(h1, [0.5,1.5,2.5]) == [2, 2, 0]
end

@testset "Arithmetic" begin
    @testset "Unweighted regular binning" begin
        h1 = Hist1D([0.5,1.5,1.5,2.5], 0:3)
        h2 = Hist1D([0.5,1.5,2.5,2.5], 0:3)

        h = h1 + h2
        @test h.hist.weights == [2.0, 3.0, 3.0]
        @test h.sumw2 == [2.0, 3.0, 3.0]

        h = h1 - h2
        @test h.hist.weights == [0.0, 1.0, -1.0]
        @test h.sumw2 == [2.0, 3.0, 3.0]

        h = h1 / h2
        @test h.hist.weights == [1.0, 2.0, 0.5]
        @test h.sumw2 == [2.0, 6.0, 0.375]

        h = h1 * 2
        @test h.hist.weights == [2.0, 4.0, 2.0]
        @test h.sumw2 == [4.0, 8.0, 4.0]

        h = h1/(h1+h2*2)
        @test h.hist.weights ≈ [0.333333, 0.5, 0.2] atol=1e-6
        @test h.sumw2 ≈ [0.17284, 0.21875, 0.0544] atol=1e-6
    end

    @testset "Weighted regular binning" begin
        h1 = Hist1D([0.5,1.5,1.5,2.5], weights([3,3,2,1]), 0:3)
        h2 = Hist1D([0.5,1.5,2.5,2.5], weights([3,3,2,1]), 0:3)

        h = h1/(h1+h2*2)
        @test h.hist.weights ≈ [0.333333, 0.454545, 0.142857] atol=1e-6
        @test h.sumw2 ≈ [0.17284, 0.191107, 0.029155] atol=1e-6
    end

    @testset "Weighted irregular binning" begin
        h1 = Hist1D([0.5,1.5,1.5,2.5], weights([3,3,2,1]), [0,1,2,4])
        h2 = Hist1D([0.5,1.5,2.5,2.5], weights([3,3,2,1]), [0,1,2,4])

        h = h1/(h1+h2*2)
        @test h.hist.weights ≈ [0.333333, 0.454545, 0.142857] atol=1e-6
        @test h.sumw2 ≈ [0.17284, 0.191107, 0.029155] atol=1e-6
    end

    @testset "Hist2D" begin
        h1 = Hist2D(([0.5,1.5,1.5,2.5], [0.5,0.5,0.5,0.5]), (0:3,0:1))
        h2 = Hist2D(([0.5,1.5,2.5,2.5], [0.5,0.5,0.5,0.5]), (0:3,0:1))
        h = h1/(h1+h2*2)
        @test vec(h.hist.weights) ≈ [0.333333, 0.5, 0.2] atol=1e-6
        @test vec(h.sumw2) ≈ [0.17284, 0.21875, 0.0544] atol=1e-6
    end

end

@testset "Merging" begin
    h1 = Hist1D(randn(100), -3:3)
    h2 = Hist1D(randn(100), -3:3)
    @test merge(h1,h2) == h1+h2

    h1 = Hist2D((randn(10),randn(10)), (-3:3,-3:3))
    h2 = Hist2D((randn(10),randn(10)), (-3:3,-3:3))
    @test merge(h1,h2) == h1+h2
end

@testset "Rebinning" begin
    h1 = Hist1D(rand(10^2), 0:0.1:1)
    @test h1 == rebin(h1, 1)
    @test integral(h1) == integral(rebin(h1, 5))
    @test sum(h1.sumw2) == sum(rebin(h1, 5).sumw2)
    @test binedges(rebin(h1, 5)) == [0, 0.5, 1.0]

    h2 = Hist1D(rand(10^2), [0.0, 0.1, 0.7, 0.9, 1.0])
    @test h2 == rebin(h2, 1)
    @test integral(h2) == integral(rebin(h2, 2))
    @test sum(h2.sumw2) == sum(rebin(h2, 2).sumw2)
    @test binedges(rebin(h2, 2)) == [0, 0.7, 1.0]

    @test rebin(h1, 2) == (h1 |> rebin(2))

    h1 = Hist2D((rand(10^2),rand(10^2)), (0:0.1:1,0:0.1:1))
    @test h1 == rebin(h1, 1, 1)
    @test integral(h1) == integral(rebin(h1, 5))
    @test sum(h1.sumw2) == sum(rebin(h1, 5).sumw2)
    @test binedges(rebin(h1, 5)) == ([0, 0.5, 1.0], [0, 0.5, 1.0])

    bins = [0.0, 0.1, 0.7, 0.9, 1.0]
    h2 = Hist2D((rand(10^2),rand(10^2)), (bins,bins))
    @test h2 == rebin(h2, 1) == (h2 |> rebin(1, 1))
    @test integral(h2) == integral(rebin(h2, 2))
    @test sum(h2.sumw2) == sum(rebin(h2, 2).sumw2)
    @test binedges(rebin(h2, 2)) == ([0, 0.7, 1.0], [0, 0.7, 1.0])

    h2 = Hist2D((rand(10^2),rand(10^2)), (0:0.1:1,0:0.5:1))
    @test nbins(rebin(h2, 10, 2)) == (1, 1)
    @test_throws AssertionError rebin(h2, 2, 10)
end

@testset "Profile" begin
    xy = collect(hcat([[-2.0, 1.5], [-2.0, -3.5], [-2.0, 1.5], [0.0, -2.0], [0.0, -2.0], [0.0, 0.0], [0.0, 2.0], [0.0, 4.0], [2.0, 1.5]]...)')
    h = Hist2D((xy[:,1],xy[:,2]), (-5:2:5,-5:2:5))

    hx = profile(h, :x)
    hy = profile(h, :y)

    @test hx == (h |> profile(:x))
    @test hx.sumw2 == [0.0, 2.6666666666666665, 1.088, 0.0, 0.0]
    @test hy.sumw2 == [0.0, 0.0, 0.0, 0.6875, 0.0]
    @test bincounts(hx) == [0.0, 0.0, 0.4, 2.0, 0.0]
    @test bincounts(hy) == [-2.0, 0.0, 0.0, -0.5, 0.0]
    @test binedges(hx) == -5:2:5
    @test binedges(hy) == -5:2:5
end

@testset "Projection" begin
    xs = rand(10)
    ys = rand(10)
    r = 0:0.1:1
    h1x = Hist1D(xs, r)
    h1y = Hist1D(ys, r)
    h2 = Hist2D((xs,ys), (r,r))
    @test project(h2, :x) == (h2 |> project(:x))
    @test h1x == project(h2, :x)
    @test h1y == project(h2, :y)
end

@testset "Transpose" begin
    h1 = Hist2D((randn(10),randn(10)), (-3:3,-3:3))
    t = FHist.transpose
    @test (t∘t)(h1) == h1
end

@testset "Repr" begin
    for h1 in (Hist1D(randn(100), -3:3),
               Hist2D((randn(10),randn(10)), (-3:3,-3:3)),
               profile(Hist2D((randn(10^5),randn(10^5)), (-3:0.1:3, -5:0.1:5)), :x))
        @test all(occursin.(["edges:", "total count:", "bin counts:"], repr(h1)))
        @test !occursin("<svg", repr(h1))
        @test all(occursin.(["edges:", "total count:", "bin counts:", "<svg"], repr("text/html", h1)))
        @test all(occursin.(["edges=", "integral=", "Hist"], repr(h1, context=:compact=>true)))
    end
end

@testset "Restrict" begin
    h = Hist1D(randn(500), -5:0.2:5)
    hleft = restrict(h, -Inf, 0.0)
    hright = restrict(h, 0.0, Inf)

    @test h == restrict(h)
    @test restrict(h, -1, 1) == (h |> restrict(-1,1))
    @test integral(hleft) + integral(hright) == integral(h)
    @test nbins(hleft) + nbins(hright) == nbins(h)
    @test sum(hleft.sumw2) + sum(hright.sumw2) == sum(h.sumw2)

    @test all(-1 .<= bincenters(restrict(h,-1,1)) .<= 1)
    @test_throws AssertionError restrict(h, 10, Inf)
end

@testset "Overflow" begin
    a = randn(10^5)
    b = randn(10^5)
    ϵ = 1e-6
    edges = -3:0.5:3
    aclip = clamp.(a,nextfloat(first(edges)),prevfloat(last(edges)))
    bclip = clamp.(b,nextfloat(first(edges)),prevfloat(last(edges)))
    sh1 = fit(Histogram, aclip, edges)
    h1 = Hist1D(a, edges, overflow=true)
    @test h1.hist == sh1
    @test h1.sumw2 == bincounts(h1)

    h1 = Hist1D(a, overflow=true)
    before = last(bincounts(h1))
    push!(h1, nextfloat(last(binedges(h1))))
    after = last(bincounts(h1))
    @test after == before + 1

    sh2 = fit(Histogram, (aclip,bclip), (edges,edges))
    h2 = Hist2D((a,b), (edges,edges), overflow=true)
    @test h2.hist == sh2
    @test h2.sumw2 == bincounts(h2)
    @test integral(project(h2, :x)) == length(aclip)
    @test rebin(h2, 2).overflow == true

end

@testset "Utils" begin
    @test all(FHist.pearson_err(10) .≈ (3.7015621187164243, 2.7015621187164243))
    @test FHist.sqrt_err(9) == (3,3)
    @test FHist._is_uniform_bins([1,2,3,4,5]) == true
    @test FHist._is_uniform_bins([1,2,3,4,5+1e-8]) == false
    @test FHist._is_uniform_bins(1:10) == true
    h = Hist1D(randn(10^3))
    @test FHist.hists_to_bars([h]) == (binedges(h)[1:end-1], bincounts(h), ones(nbins(h)))
end


text:
using ShapeFactory

app = dispatch("CATIA.Application")
prt = app.ActiveDocument.Part

# running rectangle example
#
function ex_a(num::Integer, rad::Real)

    pln = prt.FindObjectByName("xy plane")
    gst = prt.FindObjectByName("Geometrical Set.1")

    sk1 = gst.HybridSketches.Add(pln)
    f2d = sk1.OpenEdition()
    con = sk1.Constraints

    ang0 = prt.Parameters.CreateDimension("master_angle", "ANGLE", 0)
    ang0.ValuateFromString("7deg")

    pt0 = f2d.CreatePoint(0, 0)
    con.AddMonoEltCst.(0, pt0) # fixed

    ax0 = f2d.CreateLine(0, 0, 0, rad * 3)
    ax0.StartPoint = pt0
    ax0.Construction = true
    con.AddMonoEltCst(13, ax0) # vertical

    cor = [-1 1; 1 1; 1 -1; -1 -1] .* rad
    edg = []

    for n in 1:num

        pts = f2d.CreatePoint.(cor[:,1], cor[:,2])
        if n == 1
            con.AddMonoEltCst.(0, pts) # fixed
        else
            con.AddBiEltCst.(2, pts, edg) # on line
        end

        run = [pts circshift(pts, -1)]
        edg = f2d.CreateLine.(1:4, 0, 0, 0)

        setproperty!.(edg, :StartPoint, run[:, 1])
        setproperty!.(edg, :EndPoint, run[:, 2])

        con.AddBiEltCst.(11, edg[2:end], edg[1:end-1]) # perp

        if n != 1
            ax = f2d.CreateLine(0, 0, 0, 100)
            ax.StartPoint = pt0
            ax.Construction = true
            ang = con.AddBiEltCst(6, ax, ax0) # angle
            ang.AngleSector = 0
            ax0 = ax
            prt.Relations.CreateFormula("", "", ang.Dimension, "master_angle")
            con.AddBiEltCst(11, edg[1], ax) # perpendicular
        end

    end

    sk1.CloseEdition()
    prt.Update()

end

ex_a(10, 10)

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
using Gnuplot, Random
Random.seed!(123)
test = rand(50,50)
@gp "set multiplot layout 1,2"
# on top
@gp :- 1 "set auto fix" "set size square" :-
@gp :- test "w image notit"  """set cblabel "CBTitle \\n (my unit)" """ :-
@gp :- "set cblabel offset -6.5, 10 font ',8' textcolor lt 3 rotate by 0" :-
@gp :- "set tmargin at screen 0.85" :-
# on the side
@gp :- 2 "set auto fix" "set size square" :-
@gp :- test "w image notit"  """set cblabel "CBTitle (my unit)" """
@gp :- "set cblabel offset 0, 0 font ',8' textcolor lt 4 rotate by 90" :-
@gp :- palette(:viridis) :-
@gp :- "set rmargin at screen 0.9"
save(term = "pngcairo size 800,400", output = "colorbar_labels_layout.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using MPI

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Grids: halo_size

#####
##### Distributed incompressible model constructor
#####

function DistributedIncompressibleModel(; architecture, grid, model_kwargs...)
    i, j, k = architecture.local_index
    Rx, Ry, Rz = architecture.ranks
    my_connectivity = architecture.connectivity

    Nx, Ny, Nz = size(grid)
    Lx, Ly, Lz = length(grid)

    # Pull out endpoints for full grid.
    xL, xR = grid.xF[1], grid.xF[Nx+1]
    yL, yR = grid.yF[1], grid.yF[Ny+1]
    zL, zR = grid.zF[1], grid.zF[Nz+1]

    # Make sure we can put an integer number of grid points in each rank.
    # Will generalize in the future.
    # TODO: Check that we have enough grid points on each rank to fit the halos!
    @assert isinteger(Nx / Rx)
    @assert isinteger(Ny / Ry)
    @assert isinteger(Nz / Rz)

    nx, ny, nz = Nx÷Rx, Ny÷Ry, Nz÷Rz
    lx, ly, lz = Lx/Rx, Ly/Ry, Lz/Rz

    x₁, x₂ = xL + (i-1)*lx, xL + i*lx
    y₁, y₂ = yL + (j-1)*ly, yL + j*ly
    z₁, z₂ = zL + (k-1)*lz, zL + k*lz

    # FIXME? local grid might have different topology!
    my_grid = RegularRectilinearGrid(topology=topology(grid), size=(nx, ny, nz), x=(x₁, x₂), y=(y₁, y₂), z=(z₁, z₂), halo=halo_size(grid))

    ## Construct local model

    pressure_solver = haskey(model_kwargs, :pressure_solver) ? Dict(model_kwargs)[:pressure_solver] :
                                                               DistributedFFTBasedPoissonSolver(architecture, grid, my_grid)

    my_model = IncompressibleModel(;
           architecture = architecture,
                   grid = my_grid,
        pressure_solver = pressure_solver,
                   model_kwargs...
    )

    return my_model
end


text:
using MATLAB
using Test

include("engine.jl")
include("matfile.jl")
include("matstr.jl")
include("mxarray.jl")


text:
using Test
include(joinpath("..", "testhelpers.jl"))

@testset "MPIStateArrays reductions" begin
    tests = [
        (1, "basics.jl")
        (1, "broadcasting.jl")
        (1, "reductions.jl")
        (3, "reductions.jl")
        (1, "varsindex.jl")
    ]

    runmpi(tests, @__FILE__)
end


text:
# Ivonne Blaurock
# powered by Julia 0.6.1
# ------------------------------------------------------
# HILL TAIL ESTIMATOR
# ------------------------------------------------------
#

function hill(x::Vector{Float64})
    lx=length(x)
    x=sort(abs.(x))
    ly=Int(round(0.1*lx)) # fixed threshold of 10%
    z=Vector{Float64}(ly)
    for k = 1 : ly
        z[k]=1/((1/k)*sum(log(x[lx-i+1])-log(x[lx-k]) for i in 1:k))
    end
    return z
end

function hill(x::Vector{Float64},y::Float64)
    lx=length(x)
    x=sort(abs.(x))
    ly=Int(round(0.1*lx)) # fixed threshold of 10%
    z=Vector{Float64}(ly)
    lz=Int(round(y*lx))
    for k = 1 : ly
        z[k]=1/((1/k)*sum(log(x[lx-i+1])-log(x[lx-k]) for i in 1:k))
    end
    return z[lz]
end

max_stars_repo_name: luettedeern/HillTailEstimator

max_issues_repo_name: luettedeern/HillTailEstimator

max_forks_repo_name: luettedeern/HillTailEstimator


text:
module Inference

using DataFrames

include("probability.jl")
include("evidence.jl")
include("nlmodel.jl")
include("results.jl")
include("expression.jl")
include("pi.jl")
include("idMap.jl")
include("sampleTissueMap.jl")

function run(configFile, config, verbose)

    if haskey(config, "id-map")
        if verbose
            println("Reading in idMap")
        end
        idMap = IdMap.getIDmap(config["id-map"])
    else
        println("ERROR: Need to specify 'id-map' in config file")
        exit(1)
    end

    tissue = Dict()
    expression = DataFrame()
    if haskey(config, "expression")
        if verbose
            println("Reading in expression")
        end

        if haskey(config, "expression")
            if haskey(config["expression"], "file") == false || haskey(config["expression"], "tissue") == false
                println("ERROR: Expresion section of config needs 'file' and 'tissue' to be specified")
                exit(1)
            else
                expression = Expression.getExpression(config["expression"]["file"])
                tissueParam = config["expression"]["tissue"]
                if haskey(tissueParam, "mapping-file")
                    tissue["map"] = SampleTissueMap.getTissueMap(tissueParam["mapping-file"])
                elseif haskey(tissueParam, "name")
                    tissue["name"] = tissueParam["name"]
                else
                    println("ERROR: Under the expression -> tissue section of the config you need to specify either the key mapping-file or name")
                    exit(1)
                end
            end
        end
    end

    if haskey(config, "upperbound")
        upperbound = config["upperbound"]
    else
        println("ERROR: Need to specify upperbound in config")
        error(1)
    end

    if haskey(config, "lowerbound")
        lowerbound = config["lowerbound"]
    else
        warn("Need to specify lowerbound in config")
        exit(1)
    end

    runID = haskey(config, "id") ? config["id"] : Base.Random.uuid4()
    println("runID: $runID")

    if haskey(config, "results")
        resultsConfig = config["results"]
        if haskey(resultsConfig, "directory")
             dir = resultsConfig["directory"]
             resultsDir = endswith(dir, "/") ? "$dir$runID" : "$dir/$runID"
             pathwayResultsDir = "$resultsDir/pathways"
        else
            error("Need to specify directory in results section of config")
        end
    else
        error("Results section in config required")
    end

    mkpath(resultsDir)
    mkpath(pathwayResultsDir)
    cp(configFile, "$resultsDir/config.yaml", force=true)

    if haskey(config, "evidence")
        evidence = Evidence.getEvidence(config["evidence"], idMap, resultsDir)
    else
        error("Evidence section of config is missing")
    end

    if @isdefined(expression)
        expressionColumns = names(expression)
        for sample in expressionColumns
            if haskey(tissue, "name")
                tissueType = Symbol(tissue["name"])
            elseif haskey(tissue, "map")
                if haskey(tissue["map"], Symbol(sample))
                    tissueType = Symbol(tissue["map"][Symbol(sample)])
                else
                    println("Tissue for sample $sample not found")
                    exit(1)
                end
            end

            if @isdefined(tissueType) && findfirst(isequal(tissueType),expressionColumns) == 0
                println("ERROR tissue type $tissueType not found in expression file")
                exit(1)
            end
        end
    end

    for file in config["pathways"]
        if verbose
            println("Running pathway: $file")
        end

        m = match(r".*/(?<pathway>.*)\.tsv", file)
        pathwayName = m[:pathway]
        pathwayDir = "$pathwayResultsDir/$pathwayName"
        runPathway(file, expression, idMap, evidence, lowerbound, upperbound, tissue, pathwayDir, verbose)
    end
end

function runPathway(file, expression, IDMap, evidence, lowerbound, upperbound, tissue, pathwayDir, verbose)
    pinodes = Pi.readFile(file)
    nodeSampleResults = Dict()
    basePathwayTissueResults = Dict()
    for sample in keys(evidence)
        if haskey(tissue, "name") || haskey(tissue, "map")
            if haskey(tissue, "name")
                tissueType = Symbol(tissue["name"])
            elseif haskey(tissue, "map")
                if haskey(tissue["map"], Symbol(sample))
                    tissueType = Symbol(tissue["map"][Symbol(sample)])
                else
                    println("Tissue for sample $sample not found")
                    exit(1)
                end
            end

            expressionMap = Expression.getTissue(expression, tissueType)
        else
             tissueType = "unspecified"
             expressionMap = Dict()
        end

        if haskey(basePathwayTissueResults, tissueType) == false
            (basePathwayResults, x, x_bar) = NLmodel.runModel(pinodes,
                                                              Dict(),
                                                              lowerbound,
                                                              upperbound,
                                                              expressionMap,
                                                              verbose)
            basePathwayTissueResults[tissueType] = basePathwayResults
        end

        if verbose
            println("Running $sample")
        end

        nodeState = evidence[sample]
        (sampleResults, x, x_bar) = NLmodel.runModel(pinodes,
                                                     nodeState,
                                                     lowerbound,
                                                     upperbound,
                                                     expressionMap,
                                                     verbose)

        for nodeName in keys(sampleResults)
            if length(keys(nodeSampleResults)) == 0 || haskey(nodeSampleResults, nodeName) == false
                nodeSampleResults[nodeName] = Dict()
            end

            calculatedValue = sampleResults[nodeName] - (basePathwayTissueResults[tissueType][nodeName] - 1)
            nodeSampleResults[nodeName][sample] = (calculatedValue < 0.01) ? 0.01 : calculatedValue
        end
    end

    mkpath(pathwayDir)
    resultsPath = "$pathwayDir/results.tsv"

    if verbose
        println("Outputing results: $resultsPath")
    end

    Results.createcsv(nodeSampleResults,
                      keys(evidence),
                      resultsPath)
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
matches(f, x) = f(x) != nothing

function stringindex(x)
  @capture(x, nextind(str_, i_)) && return :($i + 1)
  @capture(x, prevind(str_, i_)) && return :($i - 1)
  return
end

function integer(x)
  x isa Integer && !(x isa Bool) ? rand(typeof(x)) :
  x isa AbstractFloat ? rand(typeof(x)) * x :
  nothing
end

function rmline(x)
  isexpr(x, :block) || return
  is = findall(x -> x isa Expr, x.args)
  isempty(is) && return
  Expr(:block, deleteat!(copy(x.args), rand(is))...)
end

function swapline(x)
  isexpr(x, :block) || return
  is = findall(x -> x isa Expr, x.args)
  length(is) < 2 && return
  i = rand(is)
  j = rand(setdiff(is, i))
  args = copy(x.args)
  args[i], args[j] = args[j], args[i]
  return Expr(:block, args...)
end

function flipcond(x)
  isexpr(x, :if, :elseif) || return
  Expr(x.head, Expr(:call, :!, x.args[1]), x.args[2:end]...)
end

const defaults =
  Any[stringindex, integer, rmline, swapline, flipcond]


text:
# Libz Interfaces
# ===============

mutable struct ZStream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong

    msg::Ptr{UInt8}
    state::Ptr{Cvoid}

    zalloc::Ptr{Cvoid}
    zfree::Ptr{Cvoid}
    opaque::Ptr{Cvoid}

    data_type::Cint

    adler::Culong
    reserved::Culong
end

function ZStream()
    ZStream(
        # input
        C_NULL, 0, 0,
        # output
        C_NULL, 0, 0,
        # message and state
        C_NULL, C_NULL,
        # memory allocation
        C_NULL, C_NULL, C_NULL,
        # data type, adler and reserved
        0, 0, 0)
end

const Z_DEFAULT_COMPRESSION = Cint(-1)

const Z_OK         = Cint(0)
const Z_STREAM_END = Cint(1)
const Z_BUF_ERROR  = Cint(-5)

const Z_NO_FLUSH      = Cint(0)
const Z_SYNC_FLUSH    = Cint(2)
const Z_FINISH        = Cint(4)

# The deflate compression method
const Z_DEFLATED = Cint(8)

const Z_DEFAULT_STRATEGY = Cint(0)

const Z_DEFAULT_MEMLEVEL = Cint(8)
const Z_DEFAULT_WINDOWBITS = Cint(15)

function version()
    return unsafe_string(ccall((:zlibVersion, libz), Ptr{UInt8}, ()))
end

const zlib_version = version()

function deflate_init!(zstream::ZStream, level::Integer, windowbits::Integer)
    return ccall((:deflateInit2_, libz), Cint, (Ref{ZStream}, Cint, Cint, Cint, Cint, Cint, Cstring, Cint), zstream, level, Z_DEFLATED, windowbits, #=default memlevel=#8, #=default strategy=#0, zlib_version, sizeof(ZStream))
end

function deflate_reset!(zstream::ZStream)
    return ccall((:deflateReset, libz), Cint, (Ref{ZStream},), zstream)
end

function deflate_end!(zstream::ZStream)
    return ccall((:deflateEnd, libz), Cint, (Ref{ZStream},), zstream)
end

function deflate!(zstream::ZStream, flush::Integer)
    return ccall((:deflate, libz), Cint, (Ref{ZStream}, Cint), zstream, flush)
end

function inflate_init!(zstream::ZStream, windowbits::Integer)
    return ccall((:inflateInit2_, libz), Cint, (Ref{ZStream}, Cint, Cstring, Cint), zstream, windowbits, zlib_version, sizeof(ZStream))
end

function inflate_reset!(zstream::ZStream)
    return ccall((:inflateReset, libz), Cint, (Ref{ZStream},), zstream)
end

function inflate_end!(zstream::ZStream)
    return ccall((:inflateEnd, libz), Cint, (Ref{ZStream},), zstream)
end

function inflate!(zstream::ZStream, flush::Integer)
    return ccall((:inflate, libz), Cint, (Ref{ZStream}, Cint), zstream, flush)
end

function zerror(zstream::ZStream, code::Integer)
    return throw(ErrorException(zlib_error_message(zstream, code)))
end

function zlib_error_message(zstream::ZStream, code::Integer)
    if zstream.msg == C_NULL
        return "zlib error: <no message> (code: $(code))"
    else
        return "zlib error: $(unsafe_string(zstream.msg)) (code: $(code))"
    end
end


text:
using Base.Test
using DataFrames

let
    #test_group("Confirm that we can read various file types.")

    filenames = ["test/data/blanklines/blanklines.csv",
                 "test/data/compressed/movies.csv.gz",
                 "test/data/newlines/os9.csv",
                 "test/data/newlines/osx.csv",
                 "test/data/newlines/windows.csv",
                 "test/data/newlines/embedded_os9.csv",
                 "test/data/newlines/embedded_osx.csv",
                 "test/data/newlines/embedded_windows.csv",
                 "test/data/padding/space_after_delimiter.csv",
                 "test/data/padding/space_around_delimiter.csv",
                 "test/data/padding/space_before_delimiter.csv",
                 "test/data/quoting/empty.csv",
                 "test/data/quoting/escaping.csv",
                 "test/data/quoting/quotedcommas.csv",
                 "test/data/scaling/10000rows.csv",
                 "test/data/scaling/movies.csv",
                 "test/data/separators/sample_data.csv",
                 "test/data/separators/sample_data.tsv",
                 "test/data/separators/sample_data.wsv",
                 "test/data/typeinference/bool.csv",
                 "test/data/typeinference/standardtypes.csv",
                 "test/data/utf8/corrupt_utf8.csv",
                 "test/data/utf8/short_corrupt_utf8.csv",
                 "test/data/utf8/utf8.csv"]

    for filename in filenames
        try
            df = readtable(filename)
        catch
            error(@sprintf "Failed to read %s\n" filename)
        end
    end

    readtable("test/data/comments/before_after_data.csv", allowcomments = true)
    readtable("test/data/comments/middata.csv", allowcomments = true)
    readtable("test/data/skiplines/skipfront.csv", skipstart = 3)

    # TODO: Implement skipping lines at specified row positions
    # readtable("test/data/skiplines/skipbottom.csv", skiprows = [1, 2, 3])

    # TODO: Implement skipping lines at bottom
    # readtable("test/data/skiplines/skipbottom.csv", skipstartlines = 4)

    #test_group("Confirm that we can read a large file.")

    df = DataFrame()

    nrows, ncols = 100_000, 10

    for j in 1:(ncols - 1)
        df[j] = randn(nrows)
    end
    df[ncols] = "A"

    filename = tempname()

    writetable(filename, df, separator = ',')

    df1 = readtable(filename, separator = ',')

    all(df .== df1)

    rm(filename)

    #
    # Lots of rows
    #

    df = DataFrame()

    nrows, ncols = 1_000_000, 1

    for j in 1:ncols
        df[j] = randn(nrows)
    end

    filename = tempname()

    writetable(filename, df, separator = ',')

    df1 = readtable(filename, separator = ',')

    all(df .== df1)

    rm(filename)

    #test_group("Properties of data frames returned by readtable method.")

    # Readtable with makefactors active should only make factors from columns
    # of strings.
    filename = "test/data/factors/mixedvartypes.csv"
    df = readtable(filename, makefactors = true)

    @assert typeof(df["factorvar"]) == PooledDataArray{UTF8String,Uint32,1}
    @assert typeof(df["floatvar"]) == DataArray{Float64,1}

end


text:
#=
Given a singly linked list and an integer k, remove the kth last element from the list. k is guaranteed to be smaller than the length of the list.

The list is very long, so making more than one pass is prohibitively expensive.

Do this in constant space and in one pass.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
"""
    {{meta}}

Plug in specific meta information for a blog page. The `meta` local page
variable should be given as a list of tuples of pairs like so:
```
@def meta = [("property"=>"og:video", "content"=>"http://example.com/"),
             ("propery"=>"og:title", "content"=>"The Rock")]
```
"""
function hfun_meta()
    m = locvar(:meta)
    io = IOBuffer()
    for tuple in locvar(:meta)
        write(io, "<meta ")
        for (prop, val) in tuple
            write(io, "$prop=\"$val\" ")
        end
        write(io, ">")
    end
    return String(take!(io))
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the nutation and the Equation of the Origins (EO) as in
#   equinox-based IAU-2006 theory.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] ftp://tai.bipm.org/iers/conv2010/chapter5/tab5.2e.txt
#
#   [3] Wallace, P. T., Capitaine, N (2006). Precession-nutation procedures
#       consistent with IAU 2006 resolutions. Astronomy & Astrophysics.
#
#   [4] Capitaine, N., Wallace, P. T (2006). High precision methods for locating
#       the celestial intermediate pole and origin. Astronomy & Astrophysics.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export mean_obliquity_iau2006, nutation_eo_iau2006

"""
    mean_obliquity_iau2006(jd_tt::Number)

Compute the mean obliquity of the ecliptic [rad] using the equinox-based
IAU-2006 theory in the Julian day `jd_tt` [Terrestiral Time].

The algorithm was obtained in **[1]**.

# Reference

- **[1]**: Wallace, P. T., Capitaine, N (2006). Precession-nutation procedures
    consistent with IAU 2006 resolutions. Astronomy & Astrophysics.
"""
function mean_obliquity_iau2006(jd_tt::Number)
    # Compute the Julian Centuries from `jd_tt`.
    t_tt = (jd_tt - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    a2d = 1 / 3600
    d2r = π / 180
    a2r = a2d * d2r

    # Mean obliquity of the ecliptic
    # ==============================

    # Compute the mean obliquity of the ecliptic [s].

    # NOTE: This equation is wrong in [1](p. 216, eq. 3-68)!
    # The one used here was obtained in [3].
    mϵ_2000 = @evalpoly(
        t_tt,
        +84381.406,
        -46.836769,
        -0.0001831,
        +0.00200340,
        -0.000000576,
        -0.0000000434
    )

    # Reduce to the interval [0, 2π].
    mϵ_2000 = mod(mϵ_2000 * a2r, 2π)

    return mϵ_2000
end

"""
    nutation_eo_iau2006(jd_tt::Number)

Compute the nutation parameters and the Equation of Origins (EO) at the Julian
Day `jd_tt` [TT] using the equinox-based 2006 IAU Theory of Nutation. Notice
that one can provide corrections for the nutation in obliquity (`δΔϵ_2000`)
[rad] and in longitude (`δΔψ_2000`) [rad] that are usually obtained from IERS
EOP Data (see [`get_iers_eop`](@ref)).

# Returns

- The mean obliquity of the ecliptic [rad].
- The nutation in obliquity of the ecliptic [rad].
- The nutation in longitude [rad].
- The Equation of Origins (EO) [rad].
"""
function nutation_eo_iau2006(
    jd_tt::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    # Compute the Julian Centuries from `jd_tt`.
    t_tt = (jd_tt - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    a2d = 1 / 3600
    d2r = π / 180
    a2r = a2d * d2r
    r2d = 180 / π
    d2a = 3600
    r2a = r2d * d2a

    # Fundamental arguments
    # =====================

    # Luni-solar part.
    M_s, M_m, u_Mm, D_s, Ω_m = luni_solar_args_iau2006(jd_tt)

    # Planetary part.
    λ_M☿, λ_M♀, λ_Me, λ_M♂, λ_M♃, λ_M♄, λ_M⛢, λ_M♆, p_λ = planetary_args_iau2006(jd_tt)

    # Mean obliquity of the ecliptic
    # ==============================

    mϵ_2000 = mean_obliquity_iau2006(jd_tt)

    # Nutation in the obliquity
    # ==========================================================================

    Δϵ_2000 = _iau2006_sum(
        (
            COEFS_IAU2006_NUT_OBL0,
            COEFS_IAU2006_NUT_OBL1
        ),
        t_tt,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    # Apply the correction given EOP data.
    Δϵ_2000 += δΔϵ_2000 * r2a

    # Nutation in longitude
    # ==========================================================================

    ΔΨ_2000 = _iau2006_sum(
        (
            COEFS_IAU2006_NUT_LON0,
            COEFS_IAU2006_NUT_LON1
        ),
        t_tt,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    # Apply the correction given EOP data.
    ΔΨ_2000 += δΔΨ_2000 * r2a

    # Equation of Origins (EO)
    # ==========================================================================

    ΔEO = _iau2006_sum(
        (
            COEFS_IAU2006_EO_0,
            COEFS_IAU2006_EO_1
        ),
        t_tt,
        M_s,
        M_m,
        u_Mm,
        D_s,
        Ω_m,
        λ_M☿,
        λ_M♀,
        λ_Me,
        λ_M♂,
        λ_M♃,
        λ_M♄,
        λ_M⛢,
        λ_M♆,
        p_λ
    )

    # NOTE: This equation is wrong in [1](p. 216, eq. 3-66)! The coefficients
    # were obtained from [2] and [4](eq. 69).
    EO = @evalpoly(
        t_tt,
        -0.014_506,
        -4612.156_534,
        -1.391_581_7,
        +0.000_000_44,
        +0.000_029_956,
        +3.68e-8
    ) - ΔΨ_2000 * cos(mϵ_2000) - ΔEO

    # Convert to [rad].
    EO *= a2r
    Δϵ_2000 *= a2r
    ΔΨ_2000 *= a2r

    return mϵ_2000, Δϵ_2000, ΔΨ_2000, EO
end


text:
_rescale(x, m, M) = (x .- minimum(x))./(maximum(x)-minimum(x)).*(M-m).+m

"""
    rescale!(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}

Changes the values of the layer given as its first argument, so that it has the
same *range* as the values of the layer given as its second argument.
Modification is done in-place.
"""
function rescale!(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
    return rescale!(layer, extrema(template))
end

"""
    rescale!(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}

Changes the values of the layer given as its first argument, so that it has the
same *range* as the values given as a tuple of values. Modification is done
in-place.
"""
function rescale!(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
    occ = findall(!isnothing, layer.grid)
    layer.grid[occ] .= _rescale(layer.grid[occ], t...)
end

"""
    rescale(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}

Copying version of `rescale!`.
"""
function rescale(layer::TI, template::TJ) where {TI <: SimpleSDMLayer, TJ <: SimpleSDMLayer}
    l = copy(layer)
    return rescale!(l, extrema(template))
end

"""
    rescale(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}

Copying version of `rescale!`.
"""
function rescale(layer::TI, t::Tuple{T,T}) where {TI <: SimpleSDMLayer, T <: Number}
    l = copy(layer)
    return rescale!(l, t)
end

"""
    rescale!(layer::T, p::Vector{Real}) where {T <: SimpleSDMLayer}

Rescale the values of a `layer` so that they match with the quantiles given in
`p`. Internally, this uses the `Statistics.quantile` function.
"""
function rescale!(layer::T, p::Vector{TI}) where {T <: SimpleSDMLayer, TI <: AbstractFloat}
    q = reverse!(quantile(layer, p))
    occupied = findall(!isnothing, layer.grid)
    v = collect(layer)
    for i in 1:length(q)
        layer.grid[occupied[findall(x -> x <= q[i], v)]] .= reverse(p)[i]
    end
    return layer
end

"""
    rescale(layer::T, p::Vector{Real}) where {T <: SimpleSDMLayer}

Copying version of `rescale!`.
"""
function rescale(layer::T, p::Vector{TI}) where {T <: SimpleSDMLayer, TI <: AbstractFloat}
    l = copy(layer)
    return rescale!(l, p)
end


text:
function zoolog(text)
    println("[zooclient] $text")
end

function rand_uniform(rng, lower, upper)
    return rand(rng, Float64) * (upper - lower) + lower
end

function convert_time(second)
    sec = second
    hour = Int64(floor(sec / 3600))
    sec = sec - hour * 3600
    min = Int64(floor(sec / 60))
    sec = Int64(round(sec - min * 60))
    # return "%02d:%02d:%02d" hour min sec
    return "$(hour):$(min):$(sec)"
end

function mydistance(x, y)
    dis = 0
    for i in 1:length(x)
        dis += (x[i] - y[i])^2
    end
    return sqrt(dis)
end


text:
name(ds :: DataSettings) = ds.name;
# Returns the object defining school groups
school_groups(ds :: DataSettings) = ds.sGroups;
# Returns the defined school groups as symbols
CollegeStratBase.ed_symbols(ds :: DataSettings) =
	CollegeStratBase.ed_symbols(school_groups(ds));
wage_regr_settings(ds :: DataSettings) = ds.wageRegressions;
use_parental_dummies(ds :: DataSettings) = use_parental_dummies(ds.wageRegressions);

"""
$(SIGNATURES)

Number of school groups (e.g. HSG, SC, CG). For wage regressions.
"""
CollegeStratBase.n_school(ds :: DataSettings) =
	CollegeStratBase.n_school(school_groups(ds));

"""
	$(SIGNATURES)

Number of college quality groups.
"""
n_colleges(ds :: DataSettings) = ds.nColleges;

"""
	$(SIGNATURES)

Number of two year colleges.
"""
n_2year(ds :: DataSettings) = ds.n2Year;

n_4year(ds :: DataSettings) = n_colleges(ds) - n_2year(ds);

"""
	$(SIGNATURES)

Number of HS GPA or AFQT groups.
"""
n_gpa(ds :: DataSettings) = length(hsgpa_ub(ds));

"""
	$(SIGNATURES)

Number of parental background groups.
"""
n_parental(ds :: DataSettings) = length(parental_ub(ds));

"""
	$(SIGNATURES)

Percentile upper bounds for GPA groups.
"""
hsgpa_ub(ds :: DataSettings) = ds.hsGpaUbV;

"""
	$(SIGNATURES)

Percentile upper bounds for parental groups.
"""
parental_ub(ds :: DataSettings) = ds.parentalUbV;

hsgpa_edges(ds :: DataSettings) = vcat(0.0, hsgpa_ub(ds));
hsgpa_masses(ds :: DataSettings) = diff(hsgpa_edges(ds));
parental_edges(ds :: DataSettings) = vcat(0.0, parental_ub(ds));
parental_masses(ds :: DataSettings) = diff(parental_edges(ds));
mass_gpa_yp(ds :: DataSettings) = hsgpa_masses(ds) * parental_masses(ds)';

# Mean of a Vector by HS GPA
mean_by_gpa(x_pV :: Vector{T1}, ds :: DataSettings) where T1 <: Number =
	sum(x_pV .* hsgpa_masses(ds));
mean_by_parental(x_pV :: Vector{T1}, ds :: DataSettings) where T1 <: Number =
	sum(x_pV .* parental_masses(ds));

function mean_by_qual(x_qV :: AbstractVector{T1},
	ds :: DataSettings) where T1 <: Number

	enroll_qV, _ = frac_enroll_by_qual(ds);
	return sum(x_qV .* enroll_qV);
end

mean_by_gpa_yp(x_gpM :: Matrix{T1}, ds :: DataSettings) where T1 <: Number =
    sum(x_gpM .* mass_gpa_yp(ds));

"""
	$(SIGNATURES)

Given college indices, which ones are two year colleges?
"""
is_two_year(ds :: DataSettings, cIdx :: Integer) = (cIdx .<= n_2year(ds))

"""
	$(SIGNATURES)

Set of two year college indices.
Returns a range, so that assignments of the form
	`m[two_year_colleges(ds)] .= 0.0`
work.
"""
two_year_colleges(ds :: DataSettings) = (1 : n_2year(ds));

"""
	$(SIGNATURES)

Set of four year college indices.
"""
four_year_colleges(ds :: DataSettings) = (n_2year(ds) + 1) : n_colleges(ds);

"""
	$(SIGNATURES)

Can one graduate from this college?
"""
can_graduate(ds :: DataSettings, cIdx :: Integer) = !is_two_year(ds, cIdx);

"""
	$(SIGNATURES)

List of colleges from which one cannot graduate.
"""
no_grad_idx(ds :: DataSettings) = two_year_colleges(ds);

"""
	$(SIGNATURES)

List of colleges from which students can graduate.
"""
grad_idx(ds :: DataSettings) =
	findall(cIdx -> can_graduate(ds, cIdx),  1 : n_colleges(ds))

credits_to_courses(ds :: DataSettings, nCredits) =
	nCredits ./ ds.creditsPerCourse;
courses_to_credits(ds :: DataSettings, nCourses) =
	nCourses .* ds.creditsPerCourse;

## --------------  Files

# Return gpa or afqt suffix for files, whichever is used
afqt_string(ds :: DataSettings) = ds.afqtGpa;
afqt_suffix(ds :: DataSettings) = "_" * afqt_string(ds);
parental_suffix(ds :: DataSettings) = "_inc";
quality_suffix(ds :: DataSettings) = "_q";
grad_drop_suffix(ds :: DataSettings) = "_outcome";

function file_suffix(ds :: DataSettings, suffix :: Symbol)
	if suffix ∈ (:gpa, :afqt)
		return afqt_suffix(ds);
	elseif suffix ∈ (:parental, :yp)
		return parental_suffix(ds);
	elseif suffix ∈ (:qual, :quality)
		return quality_suffix(ds);
	elseif suffix ∈ (:outcome, :gradDrop)
		return grad_drop_suffix(ds);
	else
		error("Invalid suffix: $suffix")
	end
end

function file_suffix(ds :: DataSettings, suffixV :: AbstractVector{Symbol})
	return *([file_suffix(ds, suffixV[j])  for j = 1 : length(suffixV)]...)
end

"""
	$(SIGNATURES)

Construct file name with suffix

# Example:
```julia
file_name(ds, "cumLoans", (:qual, :year), ".dat"; percentile = 90)
```
"""
function file_name(ds :: DataSettings, baseName :: String,
	suffix, fExt :: String = ".dat";
	percentile = nothing)

	fn = baseName * file_suffix(ds, suffix);
	if !isnothing(percentile)
		fn = fn * "_$(percentile)PI";
	end
	fn = fn * fExt;
	return fn
end

## --------------  Individual settings

# List of data settings names. Useful for testing.
data_settings_list() = [:default];		# , :uneven, :twoProfiles];

"""
	$(SIGNATURES)

Make named data settings. Each defines at least:

- `dataSubDir`: sub directory where copied data files are stored.

The optional argument `baseDir` points to the directory where the data files live. This defaults to `DataCollegeStrat.data_dir()` so that the data travel with the main repo.
Other than this setting, the code is location independent.
"""
function make_data_settings(dsName :: Symbol; baseDir :: String = "")
	if isempty(baseDir)
		baseDir = DataCollegeStrat.data_dir();
		# baseDir = joinpath(homedir(), "Documents", "projects", "p2019",
		# 	"college_stratification", "CollegeStrat", "data");
	end
	srcDir = "NLSY 97 moments by AFQT";

	if dsName ∈ (:default, :test)
		subDir = "Updated Types";
	elseif dsName ∈ (:uneven, :twoProfiles)
		error("No longer updated: $dsName")
		subDir = "Uneven Types";
	else
		error("Invalid name: $dsName")
	end

	if dsName == :twoProfiles
		wageRegr = wage_regressions_two();
	else
		wageRegr = default_wage_regressions();
	end

	ds = DataSettings(name = dsName, dataSubDir = joinpath(subDir, srcDir),
		baseDir = baseDir,
		wageRegressions = wageRegr)
	@assert validate_ds(ds);
	return ds
end

test_data_settings() = make_data_settings(:test);

function validate_ds(ds :: DataSettings)
	isValid = true;
	wr = wage_regr_settings(ds);
	sg = school_groups(ds);
	if !isequal(ed_symbols(wr), ed_symbols(sg))
		@warn """
			School groups and wage regression settings do not match
			$sg
			$wr
			"""
		isValid = false;
	end
	return isValid;
end

## --------  Directories

# Data files are stored here
function data_dir(ds :: DataSettings)
    return joinpath(base_dir(ds), data_sub_dir(ds))
end

# This lives in DataCollegeStrat package.
base_dir(ds :: DataSettings) = ds.baseDir;
	# joinpath(homedir(), "Documents", "projects", "p2019", "college_stratification");

data_sub_dir(ds :: DataSettings) = ds.dataSubDir;

pkg_dir() = normpath(joinpath(@__DIR__, ".."));

# Diagnostic output reports go here
out_dir(ds :: DataSettings) = joinpath(pkg_dir(), "out");

## ----------  Show

function Base.show(io :: IO, ds :: DataSettings)
	print(io,  "DataSettings ", name(ds), " with ", n_colleges(ds), " colleges.")
end

function settings_table(ds :: DataSettings)
    # borrowLimitV = borrow_limits(ds, modelUnits = false);
	# borrowLimitV = round.(borrowLimitV, digits = 0);
	gpaUbV = round.(hsgpa_ub(ds), digits = 2);
	pUbV = round.(parental_ub(ds), digits = 2);
    return [
		"DataSettings"  " ";
		"GPA classes"  "$gpaUbV";
		"Parental classes"  "$pUbV"
        # "Borrowing limits"  "$borrowLimitV"
    ]
end

## -----------  Other

function n_groups(ds :: DataSettings, group :: Symbol)
	if group ∈ (:gpa, :afqt)
		n = n_gpa(ds);
	elseif group == :parental
		n = n_parental(ds);
	elseif group == :quality
		n = n_colleges(ds);
	elseif group == :school
		n = n_school(ds);
	else
		error("Unknown group $group");
	end
	return n
end

# --------------


text:
using Genie.Router

route("/") do
  serve_static_file("welcome.html")
end


text:
"""
Calculates the SVD deomposition of a matrix in a type generic manner.
Can for example be used with a `Matrix{BigFloat}`.
Returns a `SVD` factorization object.
"""
genericsvd(x) = genericsvd!(copy(x))
export genericsvd
"""
Same as `genericsvd` but saves space by overwriting the input matrix.
"""
function genericsvd!(x::AbstractArray{T}) where {T<:Union{BigFloat, BigInt}}
  GenericSVD.svd!(x; alg=nothing)
end
function genericsvd!(x::AbstractArray)
  GenericSVD.svd!(x)
end
export genericsvd!