text:
typealias BlockableIO Union{IOStream,LibuvStream,IOBuffer,BlockIO}

type Block{T}
    source::T
    block::Array            # chunk definition
    affinity::Array         # chunk affinity
    filter::Function        # returns a processed block. can be chained. default just passes the chunk definition from "block"
    prepare::Function
end

const no_affinity = []
const loopbackip = IPv4(127,0,0,1)

##
# Pipelined filters
# TODO: relook at IO integration after IO enhancements in base
as_it_is(x) = x

as_io(x) = open(x)
as_io(f::File) = open(f.path)
function as_io(x::Tuple)
    if isa(x[1], LibuvStream)
        aio = x[1]
        maxsize = x[2]
        start_reading(aio)
        Base.wait_readnb(aio, maxsize)
        avlb = min(nb_available(aio.buffer), maxsize)
        return PipeBuffer(read(aio, Array(UInt8,avlb)))
    elseif isa(x[1], BlockableIO)
        io = x[1]
        # using BlockIO to avoid copy.
        # but since BlockIO would skip till the first record delimiter,
        # we must present to it the last read byte (which must be the last read delimiter)
        # therefore we pass the last read position as the start position for BlockIO
        pos = position(io)
        endpos = pos+x[2]
        (0 == pos) && (pos = 1)
        return BlockIO(io, pos:endpos)
    else
        return BlockIO(open(x[1]), x[2])
    end
end
function as_bufferedio(x::BlockableIO)
    buff = read(x, Array(UInt8, nb_available(x)))
    close(x)
    IOBuffer(buff)
end
as_lines(x::BlockableIO) = readlines(x)
as_bytearray(x::BlockableIO) = read(x, Array(UInt8, nb_available(x)))

as_recordio(x::BlockIO, dlm::Char='\n') = BlockIO(x, dlm)
function as_recordio(x::Tuple)
    dlm = (length(x) == 3) ? x[3] : '\n'
    if isa(x[1], LibuvStream)
        aio = x[1]
        minsize = x[2]
        start_reading(aio)
        Base.wait_readnb(aio, minsize)
        tot_avlb = nb_available(aio.buffer)
        avlb = min(tot_avlb, minsize)
        pb = PipeBuffer(read(aio, Array(UInt8,avlb)))
        write(pb, readuntil(aio, dlm))
        return pb
    else
        return BlockIO(as_io(x), dlm)
    end
end

as_wordio(x::BlockIO, dlm::Char=' ') = BlockIO(x, dlm)

|>(b::Block, f::Function) = filter(f, b)
function filter(f::Function, b::Block)
    let oldf = b.filter
        b.filter = (x)->f(oldf(x))
        return b
    end
end

function prepare(b::Block, f::Function)
    let oldf = b.prepare
        b.prepare = (x)->f(oldf(x))
        return b
    end
end

macro prepare(expr)
    local nexpr = expr
    while(isa(nexpr, Expr))
        (nexpr.args[1] == :|>) && (nexpr.args[1] = :prepare)
        nexpr = nexpr.args[2]
    end
    esc(expr)
end

##
# Iterators for blocks and affinities
# There are separate iterators for data and its affinity, because affinity is sometimes common to all chunks.
# But both can be iterated upon together by zip(blocks(b), affinities(b)).
abstract BlockIterator{T}

type BlockAffinityIterator{T} <: BlockIterator{T}
    b::Block{T}
end

type BlockDataIterator{T} <: BlockIterator{T}
    b::Block{T}
end

start{T<:BlockIterator}(bi::T) = 1
done{T<:BlockIterator}(bi::T, status) = (0 == status)

function next{T<:BlockableIO}(bi::BlockAffinityIterator{T}, status)
    blk = bi.b
    data = blk.affinity
    status = eof(blk.source) ? 0 : (status+1)
    (data, status)
end
function next{T<:Any}(bi::BlockAffinityIterator{T}, status)
    blk = bi.b
    data = isempty(blk.affinity) ? blk.affinity : blk.affinity[status]
    status = (status >= length(blk.block)) ? 0 : (status+1)
    (data, status)
end

function next{T<:BlockableIO}(bi::BlockDataIterator{T}, status)
    blk = bi.b
    data = blk.prepare(blk.block[1])
    status = eof(blk.source) ? 0 : (status+1)
    (data, status)
end
function next{T<:Any}(bi::BlockDataIterator{T}, status)
    blk = bi.b
    data = blk.prepare(blk.block[status])
    status = (status >= length(blk.block)) ? 0 : (status+1)
    (data, status)
end

blocks(b::Block) = BlockDataIterator(b)
affinities(b::Block) = BlockAffinityIterator(b)

# Returns blocks that have affinity for local process id.
# If blocks have no affinity, returns a subset after equally dividing the blocks among processors.
# Note: if blocks have affinity with multiple processors, this does not ensure an equal distribution.
localpart(blk::Block) = localpart(blk, myid())
function localpart(blk::Block, pid)
    blks = collect(blocks(blk))
    (nworkers() == 1) && return blks

    local_blocks = Any[]
    affs = collect(affinities(blk))
    if isempty(affs) || isempty(affs[1])
        # pick up only some parts corresponding to pid
        n = pid-1
        while n <= length(blks)
            push!(local_blocks, blks[n])
            n += nworkers()
        end
    else
        for n in 1:length(affs)
            (pid in affs[n]) && push!(local_blocks, blks[n])
        end
    end
    local_blocks
end

##
# Block implementations for common types

function Block(a::Array, by::Int=0, nsplits::Int=0)
    asz = size(a)
    (by == 0) && (by = length(asz))
    (nsplits == 0) && (nsplits = nworkers())
    sz = [1:x for x in asz]
    splits = map((x)->begin sz[by]=x; tuple(sz...); end, Base.splitrange(length(sz[by]), nsplits))
    #blks = [slice(a,x) for x in splits]
    blks = [a[x...] for x in splits]
    Block(a, blks, no_affinity, as_it_is, as_it_is)
end

function Block(A::Array, dims::Array=[])
    isempty(dims) && error("no dimensions specified")
    dimsA = [size(A)...]
    ndimsA = ndims(A)
    alldims = [1:ndimsA;]

    blks = Any[A]
    if dims != alldims
        otherdims = setdiff(alldims, dims)
        idx = fill!(cell(ndimsA), 1)

        Asliceshape = tuple(dimsA[dims]...)
        itershape   = tuple(dimsA[otherdims]...)
        for d in dims
            idx[d] = 1:size(A,d)
        end
        numblocks = *(itershape...)
        blks = cell(numblocks)
        blkid = 1
        if VERSION >= v"0.4.0-dev+3184"
            for idxs in CartesianRange(itershape)
                idx[otherdims] = [idxs[i] for i in 1:length(idxs)]
                blks[blkid] = reshape(A[idx...], Asliceshape)
                blkid += 1
            end
        else
            cartesianmap(itershape) do idxs...
                idx[otherdims] = [idxs...]
                blks[blkid] = reshape(A[idx...], Asliceshape)
                blkid += 1
            end
        end
    end
    Block(A, blks, no_affinity, as_it_is, as_it_is)
end

function Block(f::File, nsplits::Int=0)
    sz = filesize(f.path)
    (sz == 0) && error("missing or empty file: $(f.path)")
    (nsplits == 0) && (nsplits = nworkers())
    splits = Base.splitrange(@compat(Int(sz)), nsplits)
    data = [(f.path, x) for x in splits]
    Block(f, data, no_affinity, as_it_is, as_it_is)
end

function Block(f::File, recurse::Bool=true, nfiles_per_split::Int=0)
    files = AbstractString[]
    all_files(f.path, recurse, files)

    data = Any[]
    (nfiles_per_split == 0) && (nfiles_per_split = iceil(length(files)/nworkers()))

    while length(files) > 0
        np = min(nfiles_per_split, length(files))
        push!(data, files[1:np])
        files = files[(np+1):end]
    end

    Block([f], data, no_affinity, as_it_is, as_it_is)
end

Block(f::File) = isdir(f.path) ? Block(f, true, 0) : Block(f, 0)

Block{T<:BlockableIO}(aio::T, maxsize::Int) = Block(aio, [(aio,maxsize)], no_affinity, as_it_is, as_it_is)
Block{T<:BlockableIO}(aio::T, approxsize::Int, dlm::Char) = Block(aio, [(aio,approxsize,dlm)], no_affinity, as_it_is, as_it_is)

##
# utility methods
function all_files(path::AbstractString, recurse::Bool, list::Array)
    files = readdir(path)

    for file in files
        fullpath = joinpath(path, file)
        if isdir(fullpath) && recurse
            all_files(fullpath, true, list)
        else
            push!(list, fullpath)
        end
    end
    nothing
end


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

@doc raw"""EventSource contains information for an event.

    IoK8sApiCoreV1EventSource(;
        component=nothing,
        host=nothing,
    )

    - component::String : Component from which the event is generated.
    - host::String : Node name on which the event is generated.
"""
mutable struct IoK8sApiCoreV1EventSource <: SwaggerModel
    component::Any # spec type: Union{ Nothing, String } # spec name: component
    host::Any # spec type: Union{ Nothing, String } # spec name: host

    function IoK8sApiCoreV1EventSource(;component=nothing, host=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EventSource, Symbol("component"), component)
        setfield!(o, Symbol("component"), component)
        validate_property(IoK8sApiCoreV1EventSource, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        o
    end
end # type IoK8sApiCoreV1EventSource

const _property_map_IoK8sApiCoreV1EventSource = Dict{Symbol,Symbol}(Symbol("component")=>Symbol("component"), Symbol("host")=>Symbol("host"))
const _property_types_IoK8sApiCoreV1EventSource = Dict{Symbol,String}(Symbol("component")=>"String", Symbol("host")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1EventSource }) = collect(keys(_property_map_IoK8sApiCoreV1EventSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1EventSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EventSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EventSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EventSource[property_name]

function check_required(o::IoK8sApiCoreV1EventSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1EventSource }, name::Symbol, val)
end

ext: jl


text:
using ..Blink
import Blink: js, id
import JSExpr: JSString, jsstring
import Base: position, size, close

export Window, flashframe, shell, progress, title,
  centre, floating, loadurl, opentools, closetools, tools,
  loadhtml, loadfile, css, front

mutable struct Window
  id::Int
  shell::Shell
  content
  inittask::Union{Nothing, Task}
end

"""
    Window()
    Window(electron_options::Dict; async=true)

Create and open a new Window through Electron.

If `async==false`, this function blocks until the Window is fully initialized
and ready for you to communicate with it via javascript or the Blink API.

The `electron_options` dict is used to initialize the Electron window. See here
for the full set of Electron options:
https://electronjs.org/docs/api/browser-window#new-browserwindowoptions
"""
function Window end

shell(win::Window) = win.shell
id(win::Window) = win.id

const window_defaults = @d(:url => "about:blank",
                           :title => "Julia",
                           "node-integration" => false,
                           "use-content-size" => true,
                           :icon => resolve_blink_asset("deps", "julia.png"))

function raw_window(a::Electron, opts)
  sleep(0.1)
  println(a)
  @js a createWindow($(merge(window_defaults, opts)))
end

function Window(a::Shell, opts::AbstractDict = Dict(); async=false)
  # TODO: Custom urls don't support async b/c don't load Blink.js. (Same as https://github.com/JunoLab/Blink.jl/issues/150)
  return haskey(opts, :url) ?
    Window(raw_window(a, opts), a, nothing, nothing) :
    Window(a, Page(), opts, async=async)
end

function Window(a::Shell, content::Page, opts::AbstractDict = Dict(); async=false)
  id, callback_cond = Blink.callback!()
  url = Blink.localurl(content) * "?callback=$id"
  # Create the window.
  opts = merge(opts, Dict(:url => url))
  println("_")
  rw = raw_window(a, opts)
  w = Window(rw, a, content, nothing)
  # Note: we have to use a task here because of the use of Condition throughout
  # the codebase (it might be better to use Channel or Future which are not
  # edge-triggered). We also need to initialize this after the Window
  # constructor because we have to pass the window into the damn function.
  w.inittask = @async try
    initwindow!(w, callback_cond)
  catch exc
    @error(
      "An error occurred while trying to initialize a Blink window!",
      exception=exc,
    )
  end

  if !async
    wait(w)
  end
  return w
end

function initwindow!(w::Window, callback_cond::Condition)
  initresult = wait(callback_cond)
  if isa(initresult, AbstractDict) && get(initresult, "type", "") == "error"
      throw(JSError(
        get(initresult, "name", "unknown"),
        get(initresult, "message", "blank"),
      ))
  end
  initwebio!(w)
end

Window(args...; kwargs...) = Window(shell(), args...; kwargs...)

dot(a::Electron, win::Integer, code; callback = true) =
  js(a, :(withwin($(win), $(jsstring(code)...))),
     callback = callback)

dot(w::Window, code; callback = true) =
  ifelse(callback, dot(shell(w), id(w), code, callback = callback), w)

dot_(args...) = dot(args..., callback = false)

macro dot(win, code)
  :(dot($(esc(win)), $(esc(Expr(:quote, Expr(:., :this, QuoteNode(code)))))))
end

macro dot_(win, code)
  :(dot_($(esc(win)), $(esc(Expr(:quote, Expr(:., :this, QuoteNode(code)))))))
end

# Base.* methods

function Base.wait(w::Window)
  if w.inittask === nothing
    error("Cannot wait() a \"raw\" window (was this window created with a url arg?).")
  end
  return wait(w.inittask)
end

# Window management APIs

"""
    active(win::Window)::Bool
    active(connection)::Bool

Indicates whether the specified `Window` (or `Page`, `shell`, or other internal component)
is currently "active," meaning it has an open connection to its Electron component.

```julia-repl
julia> w = Window();

julia> active(w)
true

julia> close(w)

julia> active(w)
false
```
"""
function active end

active(s::Electron, win::Integer) =
  @js s windows.hasOwnProperty($win)

active(win::Window) = active(shell(win), id(win))

"""
    flashframe(win::Window, on=true)

Start or stop "flashing" the window to get the user's attention.

In Windows, flashes the window frame. In MacOS, bounces the app in the Dock.
https://github.com/electron/electron/blob/master/docs/api/browser-window.md#winflashframeflag
"""
flashframe(win::Window, on = true) =
  @dot_ win flashFrame($on)

"""
    progress(win::Window, p=-1)

Sets progress value in progress bar. Valid range is [0, 1.0]. Remove progress
bar when progress < 0; Change to indeterminate mode when progress > 1.

https://github.com/electron/electron/blob/master/docs/api/browser-window.md#winsetprogressbarprogress-options
"""
progress(win::Window, p = -1) =
  @dot_ win setProgressBar($p)

"""
    title(win::Window, title)

Set `win`'s title to `title.`
"""
title(win::Window, title) =
  @dot_ win setTitle($title)

"""
    title(win::Window)

Get the window's title.
"""
title(win::Window) =
  @dot win getTitle()

centre(win::Window) =
  @dot_ win center()

position(win::Window, x, y) =
  @dot_ win setPosition($x, $y)

position(win::Window) =
  @dot win getPosition()

size(win::Window, w::Integer, h::Integer) =
  invoke(size, Tuple{Window, Any, Any}, win, w, h)

size(win::Window, w, h) =
  @dot_ win setSize($w, $h)

size(win::Window) =
  @dot win getSize()

floating(win::Window, flag) =
  @dot_ win setAlwaysOnTop($flag)

floating(win::Window) =
  @dot win isAlwaysOnTop()

loadurl(win::Window, url) =
  @dot win loadURL($url)

loadfile(win::Window, f) =
  loadurl(win, "file://$f")

"""
    opentools(win::Window)

Open the Chrome Developer Tools on `win`.

See also: [`closetools`](@ref), [`tools`](@ref)
"""
opentools(win::Window) =
  @dot win openDevTools()

"""
    closetools(win::Window)

Close the Chrome Developer Tools on `win` if open.

See also: [`opentools`](@ref), [`tools`](@ref)
"""
closetools(win::Window) =
  @dot win closeDevTools()

"""
    tools(win::Window)

Toggle the Chrome Developer Tools on `win`.

See also: [`opentools`](@ref), [`closetools`](@ref)
"""
tools(win::Window) =
  @dot win toggleDevTools()

front(win::Window) =
  @dot win showInactive()

close(win::Window) =
  @dot win close()

# Window content APIs

active(::Nothing) = false

handlers(w::Window) = handlers(w.content)

msg(win::Window, m) = msg(win.content, m)

js(win::Window, s::JSString; callback = true) =
  active(win.content) ? js(win.content, s, callback = callback) :
    dot(win, :(this.webContents.executeJavaScript($(s.s))), callback = callback)

const initcss = """
  <style>html,body{margin:0;padding:0;border:0;text-align:center;}</style>
  """

function loadhtml(win::Window, html::AbstractString)
  tmp = string(tempname(), ".html")
  open(tmp, "w") do io
    println(io, initcss)
    println(io, html)
  end
  loadfile(win, tmp)
  @async (sleep(1); rm(tmp))
  return
end


text:
#!/usr/bin/env julia

import GR
const gr3 = GR.gr3

const g = 9.8        # gravitational constant

function rk4(x, h, y, f)
    k1 = h * f(x, y)
    k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
    k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
    k4 = h * f(x + h, y + k3)
    return x + h, y + (k1 + 2 * (k2 + k3) + k4) / 6.0
end

function double_pendulum(theta, length, mass)
    GR.clearws()
    GR.setviewport(0, 1, 0, 1)

    direction = zeros(3, 2)
    position = zeros(3, 3)
    for i=1:2
        direction[:,i] = [ sin(theta[i]) * length[i] * 2,
                          -cos(theta[i]) * length[i] * 2, 0]
        position[:,i+1] = position[:,i] + direction[:,i]
    end

    gr3.clear()
    # draw pivot point
    gr3.drawcylindermesh(1, [0,0.2,0], [0,1,0], [0.4,0.4,0.4], [0.4], [0.05])
    gr3.drawcylindermesh(1, [0,0.2,0], [0,-1,0], [0.4,0.4,0.4], [0.05], [0.2])
    gr3.drawspheremesh(1, [0,0,0], [0.4,0.4,0.4], [0.05])
    # draw rods
    gr3.drawcylindermesh(2, position, direction,
                         [0.6 for i=1:6], [0.05, 0.05], length * 2)
    # draw bobs
    gr3.drawspheremesh(2, position[:,2:3], [1 for i=1:6], mass * 0.2)

    gr3.drawimage(0, 1, 0, 1, 500, 500, gr3.DRAWABLE_GKS)
    GR.updatews()
end

function main()

    function derivs(t, state)
        # The following derivation is from:
        # http://scienceworld.wolfram.com/physics/DoublePendulum.html
        t1, w1, t2, w2 = state
        a = (m1 + m2) * l1
        b = m2 * l2 * cos(t1 - t2)
        c = m2 * l1 * cos(t1 - t2)
        d = m2 * l2
        e = -m2 * l2 * w2^2 * sin(t1 - t2) - g * (m1 + m2) * sin(t1)
        f =  m2 * l1 * w1^2 * sin(t1 - t2) - m2 * g * sin(t2)
        return ([w1, (e*d-b*f) / (a*d-c*b), w2, (a*f-c*e) / (a*d-c*b)])
    end

    l1 = 1.2       # length of rods
    l2 = 1.0
    m1 = 1.0       # weights of bobs
    m2 = 1.5
    t1 = 100.0     # inintial angles
    t2 = -20.0

    w1 = 0.0
    w2 = 0.0
    t = 1
    dt = 0.04
    state = [t1, w1, t2, w2] * pi / 180

    GR.setprojectiontype(2)
    gr3.setcameraprojectionparameters(45, 1, 100)
    gr3.cameralookat(6, -2, 4, 0, -2, 0, 0, 1, 0)
    gr3.setbackgroundcolor(1, 1, 1, 1)
    gr3.setlightdirection(-1, -1, -10)

    now = start = time_ns()

    while t < 30
        t, state = rk4(t, dt, state, derivs)
        t1, w1, t2, w2 = state
        double_pendulum([t1, t2], [l1, l2], [m1, m2])

        now = (time_ns() - start) / 1000000000
        if t > now
            sleep(t - now)
        end
    end
end

main()


text:
include("../server/simulation/stimulus/stimulus.jl")

const SYNAPSE = 3

function test_regular_stream()
    println("Testing regular pattern stream")

    stream = RegularPatternStream()

    # Load pattern
    file = "../../src/data/stim_2.data"

    open(file, "r") do f
        patterns = readlines(f)

        duration = length(patterns[1])
        syn_lanes = length(patterns)
        stimulus = zeros(UInt8, syn_lanes, duration)

        row = 1
        for pat_samples in patterns
            col = 1
            # Format:
            # ..|..|.....| etc
            # println(pat_samples)
            for c in pat_samples
                if c == '|'
                    stimulus[row, col] = UInt8(1)
                end
                col += 1
            end
            # println(samples[row, :])
            row += 1
        end

        config_stream!(stream, stimulus, 50)

        set_synapse!(stream, SYNAPSE)
        for t in 1:50
            print(output(stream))
            step!(stream)
        end
        println("")
    end

    println("Passed")
end

test_regular_stream()


text:
@testset "Parse" begin

    @testset "Parse ID" begin
        @test AAindex.parse_id(test_a1) == "ANDN920101"
        @test AAindex.parse_id(test_a2) == "ALTS910101"
    end

    @testset "Parse Index" begin
        @test AAindex._parse_index(test_index) == test_index_result
    end

    @testset "Parse entry" begin
        parsed = AAindex.parse(test_a1)

        @test parsed.metadata.id == AAindex.parse_id(test_a1)
        @test parsed.data == test_index_result
    end

end


text:
"""
Calculate elastic regression from function data f, for response y

    elastic_regression(f, y, timet; B=None, lambda=0, df=20, max_itr=20,
                       smooth=false)
    :param f: array (M,N) of N functions
    :param y: vector (N) of responses
    :param timet: vector (N) describing time samples
    :param B: matrix describing basis functions (M,N) (default=None generates a
              B-spline basis
    :param lambda: regularization parameter
    :param df: degree of freedom of basis
    :param max_itr: maximum number of iterations
    :param smooth: smooth data

    Returns Dict describing regression
    :return alpha: intercept
    :return beta: regression function
    :return fn: aligned functions
    :return qn: aligned srsfs
    :return gamma: warping functions
    :return q: original srsfs
    :retrun B: basis functions
    :return type: type of regressions
    :return b: coefficients
    :return SSE: sum of squared error
"""
function elastic_regression(f::Array, y::Vector, timet::Vector; B=Union{}, lambda=0, df=20,
                            max_itr=20, smooth=false)
    M, N = size(f);

    if M > 500
        parallel = true;
    elseif N > 100
        parallel = true;
    else
        parallel = false;
    end

    binsize = mean(diff(timet));

    # Create B-spline basis if none provided
    if B == Union{}
        B = bs(timet, df, 4);
    end
    Nb = size(B,2);

    # second derivative for regularization
    Bdiff = zeros(M, Nb);
    for ii in 1:Nb
        Bdiff[:, ii] = gradient(gradient(B[:,ii], binsize), binsize);
    end

    q = f_to_srsf(f, timet, smooth);

    gamma = repmat(linspace(0,1,M), 1, N);
    fn = zeros(M, N);
    qn = zeros(M, N);

    itr = 1;
    SSE = zeros(max_itr);
    alpha = 0;
    b = zeros(Nb+1);
    while itr <= max_itr
        @printf("Iteration: %d\n", itr)
        # align data
        fn = zeros(M, N);
        qn = zeros(M, N);
        for ii in 1:N
            xout = (timet[end] - timet[1]) * gamma[:, ii] + timet[1];
            fn[:,ii] = approx(timet, f[:, ii], xout);
            qn[:,ii] = warp_q_gamma(timet, q[:,ii], gamma[:,ii]);
        end

        # OLS using basis
        Phi = ones(N, Nb+1);
        for ii in 1:N
            for jj in 2:Nb+1
                Phi[ii,jj] = trapz(timet, qn[:,ii].*B[:,jj-1]);
            end
        end

        R = zeros(Nb+1,Nb+1);
        for ii in 2:Nb+1
            for jj in 2:Nb+1
                R[ii,jj] = trapz(timet, Bdiff[:,ii-1].*Bdiff[:,jj-1]);
            end
        end

        xx = Phi.' * Phi;
        inv_xx = inv(xx + lambda * R);
        xy = Phi.' * y;
        b = inv_xx * xy;

        alpha = b[1];
        beta = B * b[2:Nb+1];

        # compute the SSE
        int_X = zeros(N);
        for ii in 1:N
            int_X[ii] = trapz(timet, qn[:,ii].*beta);
        end

        SSE[itr] = sum((y-alpha-int_X).^2);

        # find gamma
        gamma_new = zeros(M, N);
        if parallel
            gamma_new = @parallel (hcat) for i=1:N
                regression_warp(beta, timet, q[:, i], y[i], alpha);
            end
        else

            for i = 1:N
                gamma_new[:,i] = regression_warp(beta, timet, q[:, i], y[i],
                                                 alpha);
            end
        end

        if norm(gamma-gamma_new) < 1e-5
            break
        else
            gamma = copy(gamma_new);
        end

        itr += 1;
    end

    # Last step with centering of gam
    gamI = sqrt_mean_inverse(gamma);
    gamI_dev = gradient(gamI, 1/(M-1));
    timet0 = (timet[end] -timet[1]) .* gamI  + timet[1];
    beta = approx(timet, beta, timet0) .* sqrt(gamI_dev);

    for k in 1:N
        qn[:, k] = approx(timet, qn[:, k], timet0) .* sqrt(gamI_dev);
        fn[:, k] = approx(timet, fn[:, k], timet0);
        gamma[:, k] = approx(timet, gamma[:, k], timet0);
    end

    out = Dict("alpha" => alpha, "beta" => beta, "fn" => fn, "qn" => qn,
               "gamma" => gamma, "q" => q, "B" => B, "type" => "linear",
               "b" => b[2:end], "SSE" => SSE[1:itr-1]);
    return out
end

"""
Calculate elastic logistic regression from function data f, for response y

    elastic_logistic(f, y, timet; B=None, df=20, max_itr=20, smooth=false)
    :param f: array (M,N) of N functions
    :param y: vector (N) of responses
    :param timet: vector (N) describing time samples
    :param B: matrix describing basis functions (M,N) (default=None generates a
              B-spline basis
    :param df: degree of freedom of basis
    :param max_itr: maximum number of iterations
    :param smooth: smooth data

    Returns Dict describing regression
    :return alpha: intercept
    :return beta: regression function
    :return fn: aligned functions
    :return qn: aligned srsfs
    :return gamma: warping functions
    :return q: original srsfs
    :retrun B: basis functions
    :return type: type of regressions
    :return b: coefficients
    :return LL: logistic loss
"""
function elastic_logistic(f, y, timet; B=Union{}, df=20, max_itr=20,
                          smooth=false)

    M, N = size(f);

    if M > 500
        parallel = true;
    elseif N > 100
        parallel = true;
    else
        parallel = false;
    end

    binsize = mean(diff(timet));

    # Create B-spline basis if none provided
    if B == Union{}
        B = bs(timet, df, 4);
    end
    Nb = size(B,2);

    q = f_to_srsf(f, timet, smooth);

    gamma = repmat(linspace(0,1,M), 1, N);
    fn = zeros(M, N);
    qn = zeros(M, N);

    itr = 1;
    LL = zeros(max_itr);
    b = zeros(Nb+1);
    alpha = 0.0;
    beta = zeros(length(timet));
    while itr <= max_itr
        @printf("Iteration: %d\n", itr)
        # align data
        fn = zeros(M, N);
        qn = zeros(M, N);
        for ii in 1:N
            xout = (timet[end] - timet[1]) * gamma[:, ii] + timet[1];
            fn[:,ii] = approx(timet, f[:, ii], xout);
            qn[:,ii] = warp_q_gamma(timet, q[:,ii], gamma[:,ii]);
        end

        Phi = ones(N, Nb+1);
        for ii in 1:N
            for jj in 2:Nb+1
                Phi[ii,jj] = trapz(timet, qn[:,ii].*B[:,jj-1]);
            end
        end

        # find alpha and beta using l_bfgs
        b0 = zeros(Nb+1);
        opt = Opt(:LD_LBFGS, Nb+1);
        maxeval!(opt, 250);
        xtol_abs!(opt, 1e-12);
        ftol_abs!(opt, 1e-30);
        min_objective!(opt, (x,grad) -> logit_optm(x,grad,Phi,y));
        optf, b, ret = optimize(opt, b0);

        alpha = b[1];
        beta = B * b[2:Nb+1];

        # compute the logistic loss
        LL[itr] = logit_loss(b, Phi, y);

        # find gamma
        gamma_new = zeros(M, N);
        if parallel
            gamma_new = @parallel (hcat) for i=1:N
                logistic_warp(beta, timet, q[:, i], y[i]);
            end
        else

            for i = 1:N
                gamma_new[:,i] = logistic_warp(beta, timet, q[:, i], y[i]);
            end
        end

        if norm(gamma-gamma_new) < 1e-5
            break
        else
            gamma = copy(gamma_new);
        end

        itr += 1;
    end

    out = Dict("alpha" => alpha, "beta" => beta, "fn" => fn, "qn" => qn,
               "gamma" => gamma, "q" => q, "B" => B, "type" => "logistic",
               "b" => b[2:end], "Loss" => LL[1:itr-1]);
    return out
end

"""
Calculate elastic m-logistic regression from function data f, for response y

    elastic_mlogistic(f, y, timet; B=None, df=20, max_itr=20, smooth=false)
    :param f: array (M,N) of N functions
    :param y: vector (N) of responses
    :param timet: vector (N) describing time samples
    :param B: matrix describing basis functions (M,N) (default=None generates a
              B-spline basis
    :param df: degree of freedom of basis
    :param max_itr: maximum number of iterations
    :param smooth: smooth data

    Returns Dict describing regression
    :return alpha: intercept
    :return beta: regression function
    :return fn: aligned functions
    :return qn: aligned srsfs
    :return gamma: warping functions
    :return q: original srsfs
    :retrun B: basis functions
    :return type: type of regressions
    :return b: coefficients
    :return n_classes: number of classes
    :return LL: logistic loss
"""
function elastic_mlogistic(f, y, timet; B=Union{}, df=20, max_itr=20,
                           delta=.01, smooth=false)

    M, N = size(f);
    # Code labels
    m = Int(maximum(y));
    Y = zeros(Int32, N, m);
    for ii in 1:N
        Y[ii, y[ii]] = 1;
    end

    if M > 500
        parallel = true;
    elseif N > 100
        parallel = true;
    else
        parallel = false;
    end

    binsize = mean(diff(timet));

    # Create B-spline basis if none provided
    if B == Union{}
        B = bs(timet, df, 4);
    end
    Nb = size(B,2);

    q = f_to_srsf(f, timet, smooth);

    gamma = repmat(linspace(0,1,M), 1, N);
    fn = zeros(M, N);
    qn = zeros(M, N);

    itr = 1;
    LL = zeros(max_itr);
    b = zeros(m*(Nb+1));
    alpha = zeros(m);
    beta = zeros(M, m);
    while itr <= max_itr
        @printf("Iteration: %d\n", itr)
        # align data
        fn = zeros(M, N);
        qn = zeros(M, N);
        for ii in 1:N
            xout = (timet[end] - timet[1]) * gamma[:, ii] + timet[1];
            fn[:,ii] = approx(timet, f[:, ii], xout);
            qn[:,ii] = warp_q_gamma(timet, q[:,ii], gamma[:,ii]);
        end

        Phi = ones(N, Nb+1);
        for ii in 1:N
            for jj in 2:Nb+1
                Phi[ii,jj] = trapz(timet, qn[:,ii].*B[:,jj-1]);
            end
        end

        # find alpha and beta using l_bfgs
        b0 = zeros(m * (Nb+1));
        opt = Opt(:LD_LBFGS, m * (Nb+1));
        maxeval!(opt, 250);
        xtol_abs!(opt, 1e-12);
        ftol_abs!(opt, 1e-30);
        min_objective!(opt, (x,grad) -> mlogit_optm(x,grad,Phi,Y));
        optf, b, ret = optimize(opt, b0);

        B0 = reshape(b, Nb+1, m);
        alpha = B0[1,:];
        beta = zeros(M, m);
        for i in 1:m
            beta[:, i] = B * B0[2:Nb+1, i];
        end

        # compute the logistic loss
        LL[itr] = mlogit_loss(b, Phi, Y);

        # find gamma
        gamma_new = zeros(M, N);
        if parallel
            gamma_new = @parallel (hcat) for i=1:N
                mlogit_warp_grad(alpha, beta, timet, q[:, i], Y[i, :],
                                 delta=delta);
            end
        else

            for i = 1:N
                gamma_new[:,i] = mlogit_warp_grad(alpha, beta, timet, q[:, i],
                                                  Y[i, :], delt=delta);
            end
        end

        if norm(gamma-gamma_new) < 1e-5
            break
        else
            gamma = copy(gamma_new);
        end

        itr += 1;
    end

    out = Dict("alpha" => alpha, "beta" => beta, "fn" => fn, "qn" => qn,
               "gamma" => gamma, "q" => q, "B" => B, "type" => "mlogistic",
               "b" => b[2:end], "n_classes" => m, "Loss" => LL[1:itr-1]);
    return out
end

"""
Prediction from elastic regression model

    elastic_prediction(f, timet, model; y=None, smooth=false)
    :param f: functions to predict
    :param timet: vector describing time samples
    :param model: calculated model (regression, musicologist, mlogistic)
    :param y: true response (default = None)
    :param smooth: smooth data (default = false)

    Returns
    :return y_pred: predicted value
    :return y_labels: labels of predicted value
    :return Perf: Performance metric if truth is supplied
"""
function elastic_prediction(f, timet, model::Dict; y=Union{}, smooth=false)
    q = f_to_srsf(f, timet, smooth);
    n = size(q, 2);

    if model["type"] == "linear" || model["type"] == "logistic"
        y_pred = zeros(n);
    elseif model["type"] == "mlogistic"
        m = model["n_classes"];
        y_pred = zeros(n, m);
    end

    for ii in 1:n
        diff = model["q"] - repmat(q[:, ii], 1, size(model["q"], 2));
        dist = sum(abs(diff).^2, 1).^(1/2);
        q_tmp = warp_q_gamma(timet, q[:, ii],
                             model["gamma"][:, indmin(dist)]);
        if model["type"] == "linear"
            y_pred[ii] = model["alpha"] + trapz(timet, q_tmp.*model["beta"]);
        elseif model["type"] == "logistic"
            y_pred[ii] = model["alpha"] + trapz(timet, q_tmp.*model["beta"]);
        elseif model["type"] == "mlogistic"
            for jj in 1:m
                y_pred[ii,jj] = model["alpha"][jj] + trapz(timet, q_tmp.*model["beta"][:, jj]);
            end
        end
    end

    if y==Union{}
        if model["type"] == "linear"
            Perf = Union{};
            y_labels =  Array(Integer, 0);
        elseif model["type"] == "logistic"
            y_pred = phi(y_pred);
            y_labels = ones(n);
            y_labels[y_pred .< 0.5] = -1;
            Perf = Union{};
        elseif model["type"] == "mlogisistic"
            y_pred = phi(reshape(y_pred,n*m,1));
            y_pred = reshpae(y_pred,n,m);
            max_val, y_labels = findmax(y_pred, 2);
            Perf = Union{};
        end
    else
        if model["type"] == "linear"
            Perf = sum((y-y_pred).^2);
            y_labels = Array(Integer, 0);
        elseif model["type"] == "logistic"
            y_pred = phi(y_pred);
            y_labels = ones(n);
            y_labels[y_pred .< 0.5] = -1;
            TP = sum(y[y_labels .== 1] .== 1);
            FP = sum(y[y_labels .== -1] .== 1);
            TN = sum(y[y_labels .== -1] .== -1);
            FN = sum(y[y_labels .== 1] .== -1);
            Perf = (TP+TN)/(TP+FP+FN+TN);
        elseif model["type"] == "mlogistic"
            y_pred = phi(reshape(y_pred,n*m,1));
            y_pred = reshape(y_pred,n,m);
            max_val, y_labels = findmax(y_pred, 2);
            temp = ind2sub((n,m),vec(y_labels));
            y_labels = temp[2];
            Perf = zeros(m);
            cls_set = collect(1:m+1);
            for ii in 1:m
                cls_sub = setdiff(cls_set, ii);
                TP = sum(y[y_labels .== ii] .== ii);
                FP = sum(y[indexin(y_labels, cls_sub).>0] .== ii)
                TN = sum(y[indexin(y_labels,cls_sub).>0] .==
                         y_labels[indexin(y_labels,cls_sub).>0]);
                FN = sum(indexin(y[y_labels.==ii], cls_sub).>0);
                Perf[ii] = (TP+TN)/(TP+FP+FN+TN);
            end
            Perf = sum(y .==y_labels)/length(y_labels);
        end
    end

    return y_pred, y_labels, Perf
end


text:
using DD
using Test

include("test_bdd.jl")


text:
module TestProductivity
  using Test
  using BioEnergeticFoodWebs

  A = zeros(Int64, (4, 4))
  n = rand(4)

  # Using system-wide regulation, producers with no consumption reach K / n
  parameters = model_parameters(A, productivity=:system)
  s = simulate(parameters, n, start=0, stop=15)
  @test s[:B][end,4] ≈ parameters[:K]/4 atol=0.001

  # Using species-wide regulation, producers with no consumption reach K
  parameters = model_parameters(A, productivity=:species)
  s = simulate(parameters, n, start=0, stop=15)
  @test s[:B][end,4] ≈ parameters[:K] atol=0.001

  # NOTE The following tests start with n = [1 1 1 1]

  # Using competitive regulation with α = 1 is neutral
  parameters = model_parameters(A, productivity=:competitive, α=1.0)
  s = simulate(parameters, ones(4), start=0, stop=15)
  @test s[:B][end,4] ≈ parameters[:K]/4 atol=0.001

  # Using competitive regulation with α > 1 is exclusive
  parameters = model_parameters(A, productivity=:competitive, α=1.05)
  s = simulate(parameters, ones(4), start=0, stop=15)
  @test s[:B][end,4] < parameters[:K]/4

  # Using competitive regulation with α < 1 is overyielding
  parameters = model_parameters(A, productivity=:competitive, α=0.95)
  s = simulate(parameters, ones(4), start=0, stop=15)
  @test s[:B][end,4] > parameters[:K]/4

  # Using competitive regulation with α = 0 is species-level regulation
  parameters = model_parameters(A, productivity=:competitive, α=0.0)
  s = simulate(parameters, ones(4), start=0, stop=15)
  @test s[:B][end,4] ≈ parameters[:K] atol=0.001

end

module TestNutrientIntake
  using Test
  using BioEnergeticFoodWebs

  A = [0 1 ; 0 0]
  b0 = [0.5, 0.5]
  c0 = [2.0, 2.0]
  k1 = [0.2]
  k2 = [0.4]

  p = model_parameters(A, productivity = :nutrients, K1 = k1, K2 = k2)

  @test BioEnergeticFoodWebs.growthrate(p, b0, 2, c = c0)[1] ≈ 2/2.4 atol=0.001
  netgrowth, G = BioEnergeticFoodWebs.get_growth(p, b0, c = c0)
  @test G[2] ≈ 0.41667 atol=0.001
  @test netgrowth[2] ≈ 0.3477 atol=0.001
  dN1dt, dN2dt = BioEnergeticFoodWebs.nutrientuptake(p, b0, c0, G)
  @test dN1dt ≈ 0.0833 atol=0.001
  @test dN2dt ≈ 0.2917 atol=0.001

  # When nutrient supply or turnover is 0, then producers growth is 0 and nutrient growth is 1

  p = model_parameters(A, productivity = :nutrients, K1 = k1, K2 = k2, supply = [0.0])
  s = simulate(p, b0, concentration = c0, stop = 1000)
  @test s[:B][end,1] ≈ .0 atol=1e-6
  @test s[:B][end,2] ≈ .0 atol=1e-6

  p = model_parameters(A, productivity = :nutrients, K1 = k1, K2 = k2, D = 0.0)
  s = simulate(p, b0, concentration = c0, stop = 1000)
  @test s[:B][end,1] ≈ .0 atol=1e-6
  @test s[:B][end,2] ≈ .0 atol=1e-6

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
# FourierLayerState
immutable FourierLayerState <: AbstractScatteredLayerState
    blobs::Vector{Mocha.Blob}
    layer::FourierLayer
end

function FourierLayerState(
        backend::Mocha.CPUBackend,
        layer::FourierLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    blobs = Vector{Mocha.Blob}(length(inputs))
    for idblob in eachindex(inputs)
        innodes = inputs[idblob].nodes
        # TODO: tighten typing
        outnodes = DataStructures.SortedDict{Path,RealFourierNode,
            Base.Order.ForwardOrdering}()
        for path in keys(innodes)
            outnodes[path] = RealFourierNode(
                innodes[path],
                findin(innodes[path].ranges, layer.pathkeys),
                layer.flags,
                layer.timelimit)
        end
        blobs[idblob] = ScatteredBlob(outnodes)
    end
    return FourierLayerState(blobs, layer)
end

function Base.findin{N}(
        ranges::NTuple{N,Pair{PathKey,StepRange{Int,Int}}},
        fourierkeys::Vector{PathKey})
    return find([any(r.first.==fourierkeys) for r in ranges])
end

function Mocha.setup(
        backend::Mocha.Backend,
        layer::FourierLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    return FourierLayerState(backend, layer, inputs, diffs)
end


text:
using JuMP
using GLPK

# Initialize problem constants
const expected_return = [0.1, 0.15, 0.16, 0.08]
const capitalization =
  [ 0.50 0.30 0.25 0.60 ;
    0.30 0.10 0.40 0.20 ;
    0.20 0.60 0.35 0.20
  ]
# const min_capitalization = [0.35, 0.3, 0.15]
const min_capitalization = [0.30, 0.3, 0.15]
const available_capital = 80e3

# Some checks
N_cat, N_funds = size(capitalization)
if length(expected_return) != N_funds
  @show length(expected_return) != N_funds
  error("Dimension error.")
end
if length(min_capitalization) != N_cat
  @show length(min_capitalization) != N_cat
  error("Dimension error.")
end

# Initialize model
model = Model(GLPK.Optimizer)

# Initialize variables
@variable(model, x[1:N_funds] >= 0)

# Set capitalization constraints
c1 = @constraint(
  model,
  [i = 1:N_cat],
  sum(x[j] * capitalization[i, j] for j in 1:N_funds) >=
    min_capitalization[i] * available_capital
)

# Set maximum capital constraint
@constraint(model, sum(x) == available_capital) # could be <=, the Max will take care of using all capital

# Set objective
@objective(model, Max, sum(expected_return[i] * x[i] for i in 1:N_funds))

print(model)

#--- Solve
optimize!(model)
@show solution_summary(model)
@show value.(x)

#--- Do a sensitivity analysis of shadow values
report = lp_sensitivity_report(model);
@show report[x[1]];
@show report[c1[1]]
@show shadow_price(c1[1])

max_stars_repo_name: permanentrance/omif-solutions

max_issues_repo_name: permanentrance/omif-solutions

max_forks_repo_name: permanentrance/omif-solutions


text:
global NodeCounter = 0
global EdgeList = []
global ExprList = Expr[]
global HeaderList = Symbol[]
global ConstList = []

immutable NodeFinder
  ind::Int64
end

step(x::T) where T = x>0 ? one(T) : zero(T)

# Defines conversion for NodeFinder Type
for flt in (Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128,
            UInt8,UInt16,UInt32,UInt64,UInt128)
     @eval function convert(::Type{NodeFinder}, x::$flt)
              global NodeCounter += 1
              push!(ConstList,[NodeCounter,x])
              return NodeFinder(NodeCounter)
           end
     @eval promote_rule(::Type{NodeFinder}, ::Type{$flt}) = NodeFinder
end

for flt in (Int8,Int16,Int32,Int64,Int128,
            UInt8,UInt16,UInt32,UInt64,UInt128)
    # Utility function for updating
    @eval function Update_Tracker_Unary(hd::Symbol,n::$flt)
              push!(EdgeList,[[n,NodeCounter+1]])
              push!(HeaderList,hd)
              global NodeCounter += 1
              return NodeFinder(NodeCounter)
          end

    for flt1 in (Int8,Int16,Int32,Int64,Int128,
                UInt8,UInt16,UInt32,UInt64,UInt128)

        @eval function Update_Tracker_Binary(hd::Symbol,n1::$flt,n2::$flt1)
                  push!(EdgeList,[[n1,NodeCounter+1],[n2,NodeCounter+1]])
                  push!(HeaderList,hd)
                  global NodeCounter += 1
                  return NodeFinder(NodeCounter)
              end
    end
end

# Unary NodeFinder Operation Definitions
exp(x::NodeFinder) = Update_Tracker_Unary(:exp,x.ind)
exp2(x::NodeFinder) = Update_Tracker_Unary(:exp2,x.ind)
exp10(x::NodeFinder) = Update_Tracker_Unary(:exp10,x.ind)
log(x::NodeFinder) = Update_Tracker_Unary(:log,x.ind)
log2(x::NodeFinder) = Update_Tracker_Unary(:log2,x.ind)
log10(x::NodeFinder) = Update_Tracker_Unary(:log10,x.ind)
acosh(x::NodeFinder) = Update_Tracker_Unary(:acosh,x.ind)
cosh(x::NodeFinder) = Update_Tracker_Unary(:cosh,x.ind)
sqrt(x::NodeFinder) = Update_Tracker_Unary(:sqrt,x.ind)
asin(x::NodeFinder) = Update_Tracker_Unary(:asin,x.ind)
sinh(x::NodeFinder) = Update_Tracker_Unary(:sinh,x.ind)
atanh(x::NodeFinder) = Update_Tracker_Unary(:atanh,x.ind)
tan(x::NodeFinder) = Update_Tracker_Unary(:tan,x.ind)
atan(x::NodeFinder) = Update_Tracker_Unary(:atan,x.ind)
acos(x::NodeFinder) = Update_Tracker_Unary(:acos,x.ind)
tanh(x::NodeFinder) = Update_Tracker_Unary(:tanh,x.ind)
asinh(x::NodeFinder) = Update_Tracker_Unary(:asinh,x.ind)
sin(x::NodeFinder) = Update_Tracker_Unary(:sin,x.ind)
cos(x::NodeFinder) = Update_Tracker_Unary(:cos,x.ind)
abs(x::NodeFinder) = Update_Tracker_Unary(:abs,x.ind)
step(x::NodeFinder) = Update_Tracker_Unary(:step,x.ind)
sign(x::NodeFinder) = Update_Tracker_Unary(:sign,x.ind)
inv(x::NodeFinder) = Update_Tracker_Unary(:inv,x.ind)
-(x::NodeFinder) = Update_Tracker_Unary(:-,x.ind)
one(x::NodeFinder) = Update_Tracker_Unary(:one,x.ind)
zero(x::NodeFinder) = Update_Tracker_Unary(:zero,x.ind)

# Binary Nodefinder Operator Definitions
/(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:/,x.ind,y.ind)
-(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:-,x.ind,y.ind)
+(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:+,x.ind,y.ind)
*(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:*,x.ind,y.ind)
min(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:min,x.ind,y.ind)
max(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:max,x.ind,y.ind)
^(x::NodeFinder,y::NodeFinder) = Update_Tracker_Binary(:^,x.ind,y.ind)

# Defines binary interactions with Float64 & Int64
for g in (:+,:-,:+,:*,:min,:max,:^)
    for t in (Float16, Float32, Float64,
              Int8,Int16,Int32,Int64,Int128,
              UInt8,UInt16,UInt32,UInt64,UInt128)
        @eval ($g)(x::NodeFinder,y::$t) = ($g)(x::NodeFinder,convert(NodeFinder,y))
        @eval ($g)(x::$t,y::NodeFinder) = ($g)(convert(NodeFinder,x),y::NodeFinder)
    end
end


text:
# Copyright 2017 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

export GaussianKernel

# using Distributions

struct GaussianKernel <: NaiveBayesKernel
    mean_given_y::Matrix{Float64}
    var_given_y::Matrix{Float64}
end

function GaussianKernel(X::AbstractVector{ItemType}, y::AbstractVector{Int}, nlabels::Int) where ItemType
    dim = length(X[1])
    occ = ones(Float64, dim, nlabels)
    C = zeros(Float64, dim, nlabels)
    V = zeros(Float64, dim, nlabels)
    # α = 1 / length(X)

    @inbounds for i in 1:length(X)
        label = y[i]
        for (j, x) in enumerate(X[i])
            C[j, label] += x
            occ[j, label] += 1
        end
    end
    C = C ./ occ

    @inbounds for i in 1:length(X)
        label = y[i]
        for (j, x) in enumerate(X[i])
           V[j, label] += (x - C[j, label])^2
        end
    end

    V = V ./ occ
    GaussianKernel(C, V)
end

function kernel_prob(nbc::NaiveBayesClassifier, kernel::GaussianKernel, x::AbstractVector{Float64})::Vector{Float64}
    n = length(nbc.le.labels)
    scores = zeros(Float64, n)
    @inbounds for i in 1:n
        pxy = 1.0
        py = nbc.probs[i]
        for j in 1:length(x)
            var2 = 2 * kernel.var_given_y[j, i]
            a = 1 / sqrt(pi * var2)
            # a = 1/sqrt(pi * abs(var2))
            pxy *= a * exp(-(x[j] - kernel.mean_given_y[j, i])^2 / var2)
        end
        scores[i] = py * pxy
    end

    scores
end


text:
using Test
using NL2sol
using Glob

print_steps = false
verbose = false
include("util.jl")

for fn in glob("./optests/*.jl")
    include(fn)
end

@testset "NL2sol tests" begin
    @test rosenbrock(scale=3)
    @test helix(scale=3)
    @test singular(scale=3)
    @test woods(scale=3)
    @test zangwill(scale=3)

    @test engvall(scale=3)
    @test branin(scale=3)
    @test beale(scale=2)
    @test cragg_levy(scale=2)
    @test box(scale=2)
    #@test davidon(scale 1, verbose=true)  # false convergence
    @test freudenstein_roth(scale=3)
    @test watson6(scale=1)
    @test watson9(scale=1)
    @test watson12(scale=1)
    #@test watson20(scale=1, verbose=true) # false convergence
    @test chebyquad(scale=2)
    @test brown_dennis(scale=3)
    @test bard(scale=3)
    @test jennrich_sampson(scale=3)
    @test kowalik_osborne(scale=3)
    @test osborne1(scale=3)
    @test osborne2(scale=2)
    @test madsen(scale=3)
    #@test meyer(verbose=true)             # false convergence
    @test brown5()
    @test brown10()
    @test brown30()
    @test brown40()
end


text:
module RocketTapOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: tap()" begin

    println("Testing: operator tap()")

    run_proxyshowcheck("Tap", tap((_) -> begin end))

    sideeffects1 = []
    sideeffects2 = []
    sideeffects3 = []
    sideeffects4 = []

    run_testset([
        (
            source = from(1:5) |> tap_on_subscribe(() -> sideeffects1 = []) |> tap((d) -> push!(sideeffects1, d)),
            values = @ts([ 1:5, c ])
        ),
        (
            source = from(1:5) |> tap_on_subscribe(() -> sideeffects2 = []) |> tap((d) -> push!(sideeffects2, d)) |> skip_next(),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> tap_on_subscribe(() -> sideeffects3 = []) |> tap((d) -> push!(sideeffects3, d)),
            values = @ts(e(1))
        ),
        (
            source = never() |> tap_on_subscribe(() -> sideeffects4 = []) |> tap((d) -> push!(sideeffects4, d)),
            values = @ts()
        )
    ])

    @test sideeffects1 == [ 1, 2, 3, 4, 5 ]
    @test sideeffects2 == [ 1, 2, 3, 4, 5 ]
    @test sideeffects3 == []
    @test sideeffects4 == []

end

end


text:
#####
##### Player types
#####

export Player, bank_roll
export Human, round_bank_roll
export AbstractAI, Bot5050

abstract type AbstractLifeForm end

struct Human <: AbstractLifeForm end

abstract type AbstractAI <: AbstractLifeForm end
struct Bot5050 <: AbstractAI end

ai_to_use() = Bot5050()

mutable struct Player{LF}
    life_form::LF
    seat_number::Int
    cards::Union{Nothing,Tuple{<:Card,<:Card}}
    bank_roll::Float64
    action_history::Vector
    action_required::Bool
    all_in::Bool
    round_bank_roll::Float64 # bank roll at the beginning of the round
    folded::Bool
    pot_investment::Float64 # accumulation of round_contribution
    checked::Bool
    last_to_raise::Bool
    active::Bool
    round_contribution::Float64
end

function Base.show(io::IO, player::Player, include_type = true)
    include_type && println(io, typeof(player))
    println(io, "$(name(player))        = $(player.cards)")
end

function Player(life_form, seat_number, cards = nothing; bank_roll = 200)
    action_history = []
    action_required = true
    all_in = false
    round_bank_roll = Float64(bank_roll)
    folded = false
    pot_investment = Float64(0)
    checked = false
    active = true
    round_contribution = Float64(0)
    last_to_raise = false
    args = (
        life_form,
        seat_number,
        cards,
        Float64(bank_roll),
        action_history,
        action_required,
        all_in,
        round_bank_roll,
        folded,
        pot_investment,
        checked,
        last_to_raise,
        active,
        round_contribution,
    )
    Player(args...)
end

cards(player::Player) = player.cards
bank_roll(player::Player) = player.bank_roll
seat_number(player::Player) = player.seat_number
name(player::Player{LF}) where {LF <: AbstractLifeForm} = "$(nameof(LF))[$(seat_number(player))]"
folded(player::Player) = player.folded
zero_bank_roll(player::Player) = bank_roll(player) ≈ 0
still_playing(player::Player) = active(player) && !folded(player)
not_playing(player::Player) = !still_playing(player)
action_history(player::Player) = player.action_history
checked(player::Player) = player.checked
last_to_raise(player::Player) = player.last_to_raise
all_in(player::Player) = player.all_in
action_required(player::Player) = player.action_required
active(player::Player) = player.active
inactive(player::Player) = !active(player)
round_bank_roll(player::Player) = player.round_bank_roll
pot_investment(player::Player) = player.pot_investment
round_contribution(player::Player) = player.round_contribution
life_form(player::Player) = player.life_form


text:
# Lambda
unparse(lae::LambdaExpr) = "$(lae.vars...) -> $(unparse(lae.body))"

# AppExpr
isinfix(oe::OExpr) = false
isinfix(s::Symbol) = s in [:+, :-, :*, :/]

function unparse(ae::AppExpr)
  if isinfix(ae.f)
    length(ae.args) == 2 || throw(ArgumentError("Infix Operator but $(length(ae.args)) args"))
    "$(ae.args)[1] $(ae.f) $(ae.args)[2]"
  else
    args =
    "$(unparse(ae.f)) $(map(unparse, ae.args)...)"
  end
end

unparse(le::LetExpr) = "$(le.var) : $(unparse(le.val))\n$(unparse(le.body))"


text:
#=
Analysis methods for finding pole positions and touch locations
=#
#Distance to pole
#(or other landmark)
function calc_p_dist(wx::Array{T,1},wy::Array{T,1},p_x::Real,p_y::Real) where T

    p_dist=1000.0
    p_id=1

    for i=1:length(wx)

        mydist = sqrt((wx[i]-p_x)^2 + (wy[i]-p_y)^2)
        if mydist < p_dist
            p_dist = mydist
            p_id=i
        end

    end
    (p_dist,p_id)
end

function calc_p_dist_interp(wx::Array{T,1},wy::Array{T,1},px::Real,py::Real) where T

end

function interpolate_whisker(x::Array{T,1},y::Array{T,1},interp_res = 1.0) where T

    s_t=WhiskerTracking.total_length(x,y)
    s_out = 1.0:interp_res:s_t
    s_in = zeros(Float64,length(x))
    for i=2:length(x)
        s_in[i] = WhiskerTracking.distance_along(x,y,i) + s_in[i-1]
    end

    x_out = zeros(Float64,length(s_out))
    y_out = zeros(Float64,length(s_out))

    x_out[1] = x[1]
    y_out[1] = y[1]
    x_out[end] = x[end]
    y_out[end] = y[end]

    interp_x = interpolate((s_in,),x,Gridded(Linear()))
    interp_y = interpolate((s_in,),y,Gridded(Linear()))

    for i = 2:(length(s_out)-1)
        x_out[i] = interp_x(s_out[i])
        y_out[i] = interp_y(s_out[i])
    end
    (x_out,y_out)
end

function get_p_distances(wx,wy,p::Array{T,2}) where T
    pd=1000.0.*ones(Float64,length(wx))

    for i=1:length(wx)
        pd[i]=WhiskerTracking.calc_p_dist(wx[i],wy[i],p[i,1],p[i,2])[1]
    end

    pd
end

#Finds the number of unique pole positions (pole positions within threshold radius)
function get_p_pos(p,new_thres=50.0)

    pos_labels=zeros(Int64,size(p,1))

    num_labels=0
    pos_x=Array{Float64,1}()
    pos_y=Array{Float64,1}()

    for i=1:size(p,1)
        if !isnan(p[i,1])

            new_label_flag=true

            for j=1:num_labels

                if sqrt((p[i,1]-pos_x[j])^2+(p[i,2]-pos_y[j])^2) < new_thres
                    new_label_flag=false
                    pos_labels[i] = j
                    break
                end
            end

            if new_label_flag
                num_labels += 1
                push!(pos_x,p[i,1])
                push!(pos_y,p[i,2])
                pos_labels[i]=num_labels
            end
        else

            pos_labels[i]=0

        end
    end

    for i=1:num_labels
        if length(find(pos_labels.==i))<1000
            pos_labels[pos_labels.==i] .= 0
        end
    end

    for i=unique(pos_labels)

        max_in_a_row=0
        temp_in_a_row=0
        in_a_row=false

        for j=1:length(pos_labels)-1

            if pos_labels[j]==i

                if in_a_row
                    temp_in_a_row += 1

                else
                    temp_in_a_row = 1
                    in_a_row = true
                end

            elseif in_a_row

                in_a_row = false
                if temp_in_a_row > max_in_a_row
                    max_in_a_row = temp_in_a_row
                end

            end
        end

        if max_in_a_row < 500

            pos_labels[pos_labels.==i] .= 0

        end
    end

    non_zero_labels=unique(pos_labels)
    non_zero_labels=non_zero_labels[non_zero_labels.!=0]

    (pos_labels, pos_x[non_zero_labels], pos_y[non_zero_labels])
end


text:
function argumentsFull_!(program::Program_, platform::String)

    return Nothing
end

function loopFull_(program::Program_, code::Expr, platform::String)

    code = vectorize_(program.agent, code, program,interaction=true)
    loop = :(
        for ic2_ in 1:N
            $code
        end)
    loop = simpleFirstLoop_(platform, loop)
    loop = subs_(loop,:nnic2_,:ic2_)

    return loop
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Random
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "apply gates to density matrix" begin
    N = 5
    ψ = randn(ComplexF64, 2^N)
    ψ /= norm(ψ)
    # @code_warntype(density_from_statevector(ψ))

    @testset "density matrix apply single qubit gates" begin

        for g in [X, Y, Z, HadamardGate(), SGate(), SdagGate(), TGate(), TdagGate(), RxGate(-1.1), RyGate(0.7), RzGate(0.4), RotationGate([-0.3, 0.1, 0.23]), PhaseShiftGate(0.9), Qaintessent.Proj1Gate()]
            cg = circuit_gate(rand(1:N), g)

            ψs = apply(ψ, cg)

            ρ = density_from_statevector(ψ)
            ρsref = density_from_statevector(ψs)
            ρs1 = apply(ρ, cg)
            ρs2 = apply!(ρ, cg)
            # @code_warntype(apply(cg, ρ))
            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v

            ρ = density_from_statevector(ψ)
            ρsref = density_from_matrix(reshape(0.5 * (kron(conj(ψ), ψs) + kron(conj(ψs), ψ)), 2^N, 2^N))
            ρs1 = Qaintessent.apply_mixed_add(ρ, cg)
            ρs2 = Qaintessent.apply_mixed_add!(ρ, cg)
            # @code_warntype(Qaintessent.apply_mixed_add(cg, ρ))
            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v

            ρ = density_from_statevector(ψ)
            ρsref = density_from_matrix(reshape(0.5im * (kron(conj(ψ), ψs) - kron(conj(ψs), ψ)), 2^N, 2^N))
            ρs1 = Qaintessent.apply_mixed_sub(ρ, cg)
            ρs2 = Qaintessent.apply_mixed_sub!(ρ, cg)
            # @code_warntype(Qaintessent.apply_mixed_sub(cg, ρ))
            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v
        end
    end

    @testset "density matrix apply two qubit gates" begin
        for g in [SwapGate(), EntanglementXXGate(-0.7), EntanglementYYGate(0.2*π), EntanglementZZGate(0.4)]
            i, j = randperm(N)[1:2]
            cg = CircuitGate((i, j), g)
            ψs = apply(ψ, cg)

            ρ = density_from_statevector(ψ)
            ρsref = density_from_statevector(ψs)
            ρs1 = apply(ρ, cg)
            ρs2 = apply!(ρ, cg)

            # @code_warntype(apply(cg, ρ))
            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v

            ρ = density_from_statevector(ψ)
            ρsref = density_from_matrix(reshape(0.5 * (kron(conj(ψ), ψs) + kron(conj(ψs), ψ)), 2^N, 2^N))
            ρs1 = Qaintessent.apply_mixed_add(ρ, cg)
            ρs2 = Qaintessent.apply_mixed_add!(ρ, cg)
            # @code_warntype(Qaintessent.apply_mixed_add(cg, ρ))

            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v

            ρ = density_from_statevector(ψ)
            ρsref = density_from_matrix(reshape(0.5im * (kron(conj(ψ), ψs) - kron(conj(ψs), ψ)), 2^N, 2^N))
            ρs1 = Qaintessent.apply_mixed_sub(ρ, cg)
            ρs2 = Qaintessent.apply_mixed_sub!(ρ, cg)
            # @code_warntype(Qaintessent.apply_mixed_sub(cg, ρ))

            @test ρs1.v ≈ ρsref.v
            @test ρs2.v ≈ ρsref.v
        end
    end

    @testset "density matrix apply controlled gate" begin
        # one target qubit
        iwperm = randperm(N)
        cg = circuit_gate(iwperm[1], RotationGate(rand(3) .- 0.5), (iwperm[2], iwperm[3], iwperm[4]))

        ψs = apply(ψ, cg)
        ρsref = density_from_statevector(ψs)
        ρ = density_from_statevector(ψ)
        ρs1 = apply(ρ, cg)
        ρs2 = apply!(ρ, cg)
        # @code_warntype(apply(cg, ρ))
        @test ρs1.v ≈ ρsref.v
        @test ρs2.v ≈ ρsref.v

        # two target qubits
        iwperm = randperm(N)
        cg = circuit_gate((iwperm[1], iwperm[2]), EntanglementYYGate(2π*rand()), (iwperm[3], iwperm[4]))

        ψs = apply(ψ, cg)

        ρsref = density_from_statevector(ψs)
        ρ = density_from_statevector(ψ)
        ρs1 = apply(ρ, cg)
        ρs2 = apply!(ρ, cg)
        # @code_warntype(apply(cg, ρ))
        @test ρs1.v ≈ ρsref.v
        @test ρs2.v ≈ ρsref.v
    end

    @testset "density matrix apply general unitary gate" begin
        # matrix gate (general unitary gate)
        cg = CircuitGate(NTuple{3,Int64}(randperm(N)[1:3]), MatrixGate(Matrix(qr(randn(ComplexF64, 8, 8)).Q)))
        ψs = apply(ψ, cg)
        ρ = density_from_statevector(ψ)
        ρsref = density_from_statevector(ψs)
        ρs1 = apply(ρ, cg)
        ρs2 = apply!(ρ, cg)
        # @code_warntype(apply(cg, ρ))
        @test ρs1.v ≈ ρsref.v
        @test ρs2.v ≈ ρsref.v
    end

    @testset "density matrix apply multiple unitary gates" begin
        iwperm = randperm(N)

        cgs = [circuit_gate((iwperm[1], iwperm[2]), EntanglementYYGate(2π*rand()), (iwperm[3], iwperm[4])),
            circuit_gate((iwperm[1], iwperm[2]), EntanglementZZGate(2π*rand()), (iwperm[3], iwperm[4])),
            circuit_gate((iwperm[1], iwperm[2]), EntanglementXXGate(2π*rand()), (iwperm[3], iwperm[4])),
        ]
        ψs = apply(ψ, cgs)

        ρsref = density_from_statevector(ψs)
        ρ = density_from_statevector(ψ)
        ρs1 = apply(ρ, cgs)
        ρs2 = apply!(ρ, cgs)
        # @code_warntype(apply(cg, ρ))
        @test ρs1.v ≈ ρsref.v
        @test ρs2.v ≈ ρsref.v
    end
end


text:
using Base: Integer
using LinearAlgebra
using OrdinaryDiffEq
using ForwardDiff
using PyPlot
using ControlSystems
using Distributions
using GaussianProcesses
using BayesianOptimization

abstract type AbstractPolicy <: Function end
abstract type AbstractEnvironment end

struct LQRPolicy{T<:Real} <: AbstractPolicy
    A::Matrix{T}
    B::Vector{T}
    K::AbstractArray
    Q::Matrix{T}
    R::Matrix{T}
    max_ctrl::T
end
(u::LQRPolicy)(x) = clamp(-dot(u.K, x), -u.max_ctrl, u.max_ctrl)
(u::LQRPolicy)(p, x) = clamp(-dot(p, x), -u.max_ctrl, u.max_ctrl)

struct LinearPolicy{T<:Real} <: AbstractPolicy
    K::AbstractArray
    Q::Matrix{T}
    R::Matrix{T}
    max_ctrl::T
end
(u::LinearPolicy)(x) = clamp(-dot(u.K, x), -u.max_ctrl, u.max_ctrl)
(u::LinearPolicy)(p, x) = clamp(-dot(p, x), -u.max_ctrl, u.max_ctrl)

struct IWPParams{T}
    gravity::T
    mp::T
    mr::T
    m::T
    Jp::T
    Jr::T
    J::T
    lp::T
    lr::T
    l::T
    r::T
    max_ctrl::T
end

mutable struct IWPEnv{T<:Real} <: AbstractEnvironment
    params::IWPParams{T}
    policy::AbstractPolicy
    NX::Integer
    NU::Integer
    NT::Integer
    dt::T
    tspan::Tuple{T, T}
    trajectory::Array{T, 2}
    ctrl_input::Vector{T}
    loss::T
end
function Base.show(io::IO, params::IWPParams)
    for f in fieldnames(IWPParams)
        print(io, f)
        print(io, " => ")
        print(io, getfield(params, f))
        println()
    end
end

function IWPEnv(T::Type=Float32;
        gravity=T(9.81),
        mp=T(0.706),
        mr=T(0.223),
        m=T(mp+mr),
        Jp=T(0.00249),
        Jr=T(0.0014615),
        lr=T(0.18217),
        lp=T(lr/2),
        l=T(1/m*(mp*lp+mr*lr)),
        J=T(Jp + mp*lp*lp + mr*lr*lr),
        r=T(45/34*13/3),
        max_ctrl=T(0.25),
        # policy=LQRPolicy(
        #     rand(Float32, 4, 4), rand(Float32, 4), rand(Float32,4),
        #     rand(Float32, 4, 4), Matrix{Float32}(undef, 1, 1), max_ctrl
        # ),
        policy=LinearPolicy(
                Vector{Float32}(undef, 4), Matrix{Float32}(undef, 4, 4),
                Matrix{Float32}(undef, 1, 1), max_ctrl),
        dt=T(1/20),
        tspan=(T(0),T(10)),
    )

    params = IWPParams(
        T(gravity),
        T(mp),
        T(mr),
        T(m),
        T(Jp),
        T(Jr),
        T(J),
        T(lp),
        T(lr),
        T(l),
        T(r),
        T(max_ctrl)
    )

    u = IWPEnv(
        params,
        policy,
        4,
        1,
        length(range(tspan..., step=dt)),
        T(dt),
        T.(tspan),
        Array{T,2}(undef, (4,1)),
        Array{T,1}(undef, 1),
        T(0)
    )
    C = [1 0 0 0; 0 1 0 0];
    Q = C'*C
    R = Matrix{Float32}(undef, 1, 1)
    R[1,1] = 100.0
    # A, B, K, _ = lqr_gains(u; Q=Q, R=R[1])
    # u.policy.A[:] = A[:]
    # u.policy.B[:] = B[:]
    # u.policy.K[:] = K[:]
    u.policy.Q[:] = T.(Q)[:]
    u.policy.R[:] = T.(R)[:]
    u.policy.K[:] = [-3.0, -0.1, -0.5, -0.03]
    return u
end
function Base.show(io::IO, env::IWPEnv)
    print(io, "Reaction Wheel Pendulum Environment: ")
    print(io, "DataType=$(typeof(env).parameters[1]), ")
    print(io, "tspan=$(env.tspan), ")
    print(io, "dt=$(env.dt) ")
end

function _known_dynamics(env::IWPEnv, x, u)
    q1, q2, q1dot, q2dot = x
    J = env.params.J
    Jr = env.params.Jr
    m = env.params.m
    g = env.params.gravity
    l = env.params.l
    r = env.params.r
    return [
        x[3],
        x[4],
        1/(J*0.8) * (m*g*l*sin(q1) - r*u),  # Introduced a 20% error into J
        1/Jr * r * u
    ]
end

function _dynamics(env::IWPEnv, x, u)
    q1, q2, q1dot, q2dot = x
    J = env.params.J
    Jr = env.params.Jr
    m = env.params.m
    g = env.params.gravity
    l = env.params.l
    r = env.params.r
    return [
        x[3],
        x[4],
        1/J * (m*g*l*sin(q1) - r*u),
        1/Jr * r * u
    ]
end

function dynamics(env::IWPEnv, x, u)
    q1, q2, q1dot, q2dot = x
    J = env.params.J
    Jr = env.params.Jr
    m = env.params.m
    g = env.params.gravity
    l = env.params.l
    r = env.params.r
    return [
        x[3],
        x[4],
        1/J * (m*g*l*sin(q1) - r*u[1]),
        1/Jr * r * u[1]
    ]
end

function dynamics(env::IWPEnv, dx, x, u)
    q1, q2, q1dot, q2dot = x
    J = env.params.J
    Jr = env.params.Jr
    m = env.params.m
    g = env.params.gravity
    l = env.params.l
    r = env.params.r

    dx[1] = x[3]
    dx[2] = x[4]
    dx[3] = 1/J * (m*g*l*sin(q1) - r*u[1])
    dx[4] = 1/Jr * r * u[1]
    nothing
end

function lqr_gains(env::IWPEnv{T}; xbar=T[0,0,0,0], Q=I(4), R=1.0) where {T<:Real}
    A = ForwardDiff.jacobian(x->_known_dynamics(env, x, 0), xbar)
    B = ForwardDiff.derivative(u->_known_dynamics(env, xbar, u), 0)
    S = care(A,B,Q,R)
    K = R\B'*S
    A, B, K, S
end

function simulate(env::AbstractEnvironment, x0; tf=last(env.tspan), saveat=env.dt)
    OrdinaryDiffEq.solve(
        ODEProblem(
            ODEFunction((x,p,t) -> dynamics(env, x, env.policy(p,x))),
            x0,
            (first(env.tspan), tf),
            env.policy.K
        ),
        Tsit5(),
        # reltol=1e-8, abstol=1e-6,
        saveat=saveat
    )
end

function simulate(env::AbstractEnvironment, x0, k; tf=last(env.tspan), saveat=env.dt)
    OrdinaryDiffEq.solve(
        ODEProblem(
            ODEFunction((x,p,t) -> dynamics(env, x, env.policy(p,x))),
            x0,
            (first(env.tspan), tf),
            k
        ),
        Tsit5(),
        # reltol=1e-8, abstol=1e-6,
        saveat=saveat
    )
end

function simulate!(env::AbstractEnvironment, x0; tf=last(env.tspan), saveat=env.dt)
    env.trajectory = Array(simulate(env, x0; tf=tf, saveat=saveat))
    env.ctrl_input = mapslices(env.policy, env.trajectory, dims=1) |> vec
    env.loss = lqr_loss(env, env.trajectory)
end

function simulate!(env::AbstractEnvironment, x0, k; tf=last(env.tspan), saveat=env.dt)
    env.trajectory = Array(simulate(env, x0, k; tf=tf, saveat=saveat))
    env.ctrl_input = mapslices(env.policy, env.trajectory, dims=1) |> vec
    env.loss = lqr_loss(env, env.trajectory)
end

function lqr_loss(env::IWPEnv{T}, x) where {T<:Real}
    x = env.trajectory
    u = env.ctrl_input
    return sum(
        1/2*dot(x[:,i], env.policy.Q * x[:,i]) + 1/2*env.policy.R[1,1]*u[i]^2
        for i=1:size(x,2)
    )
end

## Bayesian Optimization to find acceptable linear policy gains from this point options

function find_gains_BO(env::AbstractEnvironment, reset::Bool=false; maxiters::Int=100)

    function f(env::AbstractEnvironment, k, x0=[1*π/180, 0.5, 0., 0])
        simulate!(env, x0, k)
        return env.loss
    end
    g(k) = f(env, k)

    # Choose as a model an elastic GP with input dimensions 2.
    # The GP is called elastic, because data can be appended efficiently.
    model = ElasticGPE(4,                            # 2 input dimensions
        mean = MeanConst(0.),
        kernel = SEArd([0., 0., 0., 0.], 5.),
        logNoise = 0.,
        capacity = 3000)              # the initial capacity of the GP is 3000 samples.
    set_priors!(model.mean, [Normal(1, 2)])

    # Optimize the hyperparameters of the GP using maximum a posteriori (MAP) estimates every 50 steps
    modeloptimizer = MAPGPOptimizer(every = 50, noisebounds = [-4, 3],       # bounds of the logNoise
                kernbounds = [[-1, -1, -1, -1, 0], [4, 4, 4, 4, 10]],  # bounds of the 5 parameters GaussianProcesses.get_param_names(model.kernel)
                maxeval = 40)

    # modeloptimizer = NoModelOptimizer()

    opt = BOpt(g,
            model,
            UpperConfidenceBound(),                   # type of acquisition
            modeloptimizer,
            [-10., -10., -10., -10.], [0., 0., 0., 0.],                     # lowerbounds, upperbounds
            # repetitions = 5,                          # evaluate the function for each input 5 times
            maxiterations = 200,                      # evaluate at 100 input positions
            sense = Min,                              # minimize the function
            acquisitionoptions = (method = :LD_LBFGS, # run optimization of acquisition function with NLopts :LD_LBFGS method
                        restarts = 5,       # run the NLopt method from 5 random initial conditions each time.
                        maxtime = 0.1,      # run the NLopt method for at most 0.1 second each time
                        maxeval = 1000),    # run the NLopt methods for at most 1000 iterations (for other options see https://github.com/JuliaOpt/NLopt.jl)
    verbosity = Progress)

    # x = Array{Array{Float64, 1}, 1}()
    # push!(x, K_estimate[:])
    # push!(x, K_true[:])
    # y = -f.(x)
    # append!(model, hcat(x...), y)

    opt = BOpt(g,
    model,
    #    UpperConfidenceBound(),
    ThompsonSamplingSimple(),
    # ExpectedImprovement(),
    modeloptimizer,
    [-10., -10., -10., -10.], [0., 0., 0., 0.],
    maxiterations = maxiters,
    sense = Min,
    initializer_iterations = 0
    )

    result = boptimize!(opt)
    return result, opt
end


text:
# solve Problem Set 6 Problem 1
#using ArrayViews  # non-copying subarray package
using PETSc

function driver()
  xmin = 0
  xmax = 1
  N = 50  # N+1 = # of grid points per processor
  mpi_rank = MPI.Comm_rank(MPI.COMM_WORLD) + 1
  mpi_size = MPI.Comm_size(MPI.COMM_WORLD)

  Npts_global = (mpi_size -1)*(N+1) + N + 2 - 2
  delta_x = (xmax - xmin)/(Npts_global)
  delta_t = delta_x

  tmax = 4.0

  solve(xmin, xmax, tmax, N, delta_t)
end

function solve(xmin, xmax, tmax, N, delta_t)
# solves the advection-diffusion equation using centered differences in space
# Crank-Nicolson in time
# Dirchlet boundary condition on left, Neumann on the right
# Method of Manufactured Solutions used to verify convergence rate

# Parameters:
# xmin = minimum x coordinate
# xmax = maximum x coordinate
# tmax = maximum time value
# N : N+1 =  number of x points
# r = delta_t/delta_x^2
# sigma = delta_t/delta_x

# Helper Functions:
# ICFunc = initial condition function with signature val = ICFunc(x)
# BCL = left boundary condition function with signature val = BCL(t)
# BCR = right boundary condition function
# source = source function, signature val = source(x, t)
# this function assumes no ghost point on the left, ghost point on the right

# calculate some values
mpi_rank = MPI.Comm_rank(MPI.COMM_WORLD) + 1
mpi_size = MPI.Comm_size(MPI.COMM_WORLD)

if mpi_rank != mpi_size
  mat_size = N+1
else  # last segment has ghost point
  mat_size = N+2
end

Npts_global = (mpi_size -1)*(N+1) + N + 2 - 2
delta_x = (xmax - xmin)/Npts_global

r = delta_t/(delta_x^2)
sigma = delta_t/delta_x
nStep = convert(Int, div(tmax, delta_t))

if mpi_rank == 1
println("tmax = ", tmax)
println("delta_x = ", delta_x)
println("delta_t = ", delta_t)
println("r = ", r)
println("sigma = ", sigma)
println("nStep = ", nStep)
end

# get the Petsc objects
A, u_i, u_ghost, vec_scatter, rhs, ksp, uex, ctx, ghost_offset = createPetscData(mat_size, 3)

# get the ranges to iterate over
idx_range = localpart(u_i)
idx_start = idx_range[1]
idx_end = idx_range[end]
interior_part = (idx_start+1):(idx_end-1)

# apply IC
# Not applying BCL at initial condition
for i in localpart(u_i)
  u_i[i] = ICFunc(xmin + (i-1)*delta_x)
end
# assemble so we can access the values in u_i for the Neumann BC
GridapDistributedPETScWrappers.AssemblyBegin(u_i)
GridapDistributedPETScWrappers.AssemblyEnd(u_i)

if mpi_rank == mpi_size
  u_i[idx_end] = 2*delta_x*BCR(0) + u_i[idx_end-2]
end

# scatter the initial condition u_i to u_ghost, so all procs have the values
# needed to evaluate their stencils
scatter!(vec_scatter, u_i, u_ghost)

# populate LHS matrix A
# do left BC
if mpi_rank == 1
#  println("applying left BC")
  A[1, 1] = 1
end

# do right BC
if mpi_rank == mpi_size
  A[idx_end, idx_end] = 1/(2*delta_x)
  A[idx_end, idx_end-2] = -1/(2*delta_x)  # this BC makes A not Tridiagonal
end

# to matrix interior
stencil_l = -r/2 - sigma/4
stencil_c = 1 + r
stencil_r = -r/2 + sigma/4
for i in interior_part
  A[i, i-1] = stencil_l
  A[i, i] = stencil_c
  A[i, i+1] = stencil_r
end

# do boundary nodes of each proc
if mpi_rank > 1  # all processes except leftmost
#  println("p$mpi_rank doing idx_start")
  A[idx_start, idx_start-1] = stencil_l
  A[idx_start, idx_start] = stencil_c
  A[idx_start, idx_start+1] = stencil_r
end
if mpi_rank != mpi_size # all processes except rightmost
  A[idx_end, idx_end-1] = stencil_l
  A[idx_end, idx_end] = stencil_c
  A[idx_end, idx_end+1] = stencil_r
end

# create rhs stencil
stencil_l = r/2 + sigma/4
stencil_c = 1 - r
stencil_r = r/2 - sigma/4

print("\n")
time = @elapsed for tstep=1:nStep  # loop over timesteps
# advance from timestep tstep to tstep + 1

  if mpi_rank == 1
    println("tstep = ", tstep)
  end

  for i in interior_part
    u_k = u_ghost[ i + ghost_offset ]
    u_k_1 = u_ghost[ i-1 + ghost_offset ]
    u_k_p1 = u_ghost[ i+1 + ghost_offset ]
    src_val = source(xmin + (i-1)*delta_x, (tstep - 0.5)*delta_t)
    rhs[i] = stencil_l*u_k_1 + stencil_c*u_k  + stencil_r*u_k_p1 + delta_t*src_val
  end

  if mpi_rank > 1  # all processes except leftmost
    i = idx_start
    u_k = u_ghost[i + ghost_offset]
    u_k_1 = u_ghost[i-1 + ghost_offset]
    u_k_p1 = u_ghost[i+1 + ghost_offset]
    src_val = source(xmin + (i-1)*delta_x, (tstep - 0.5)*delta_t)
    rhs[i] = stencil_l*u_k_1 + stencil_c*u_k  + stencil_r*u_k_p1 + delta_t*src_val
  end

  if mpi_rank < mpi_size  # all processes except rightmost
    i = idx_end
    u_k = u_ghost[i + ghost_offset]
    u_k_1 = u_ghost[i-1 + ghost_offset]
    u_k_p1 = u_ghost[i+1 + ghost_offset]
    src_val = source(xmin + (i-1)*delta_x, (tstep - 0.5)*delta_t)
    rhs[i] = stencil_l*u_k_1 + stencil_c*u_k  + stencil_r*u_k_p1 + delta_t*src_val
  end

  # apply BC terms to rhs
  if mpi_rank == 1
    rhs[1] = BCL(tstep*delta_t)
  end
  if mpi_rank == mpi_size
    rhs[idx_end] = BCR(tstep*delta_t)
  end

  # solve for next time step, u_i gets overwritten with new solution
  A_ldiv_B!(ksp, rhs, u_i)

  # scatter u_i into u_ghost, which is used to calculate the rhs
  # this is necessary because Petsc does not allow fetching off process values
  scatter!(vec_scatter, u_i, u_ghost)

end

println("time = ", time)

# print numerical solution
GridapDistributedPETScWrappers.AssemblyBegin(u_i)
GridapDistributedPETScWrappers.AssemblyEnd(u_i)
if mpi_rank == 1
  println("u_final = \n")
end
petscview(u_i)

# get exact solution
for i in localpart(uex)
  uex[i] = uExact( (i-1)*delta_x, delta_t*nStep)
end

# print exact solution
GridapDistributedPETScWrappers.AssemblyBegin(uex)
GridapDistributedPETScWrappers.AssemblyEnd(uex)
if mpi_rank == 1
  println("u_ex = \n")
end
petscview(uex)

# calculate error
u_err = uex - u_i
max_err = norm(u_err, Inf)
sleep(2)
if mpi_rank == 1
  println("max error = ", max_err)
end

writeU(u_i, N, nStep)
MPI.Barrier(MPI.COMM_WORLD)

return u_i, delta_t*(nStep)

end

function writeU(u_i, N, tsteps)

  mpi_rank = MPI.Comm_rank(MPI.COMM_WORLD) + 1
  larr = LocalArray(u_i)
  fname = string("u_", tsteps, "_", mpi_rank, ".dat")
  writedlm(fname, larr[1:(N+1)])
end

function ICFunc(x)
  return 2*cos(3*x)
end

function BCZero(x)
  return 0.0
end

function BCL(t)
  return 2*cos(t)
end

function BCR(t)
  return -6*sin(3)*cos(t)
end

function source(x, t)
  return -2*cos(3*x)*sin(t) - 6*sin(3x)*cos(t) + 18*cos(3*x)*cos(t)
end

function uExact(x, t)
  return 2*cos(3*x)*cos(t)
end

function createPetscData(mat_size, stencil_size)
# create all the Petsc objects needed

  mpi_rank = MPI.Comm_rank(MPI.COMM_WORLD) + 1
  mpi_size = MPI.Comm_size(MPI.COMM_WORLD)

  # left hand side matrix
  A = GridapDistributedPETScWrappers.Mat(Float64, GridapDistributedPETScWrappers.C.PETSC_DECIDE, GridapDistributedPETScWrappers.C.PETSC_DECIDE, mlocal=mat_size, nlocal=mat_size, nz=stencil_size, onz=1, comm=MPI.COMM_WORLD)
  A.assembling=true

  # solution at current timestep
  u_i = Vec(Float64, GridapDistributedPETScWrappers.C.VECMPI, comm=MPI.COMM_WORLD);
  u_i.assembling=true
  resize!(u_i, mlocal=mat_size)

  # right hand side
  rhs = Vec(Float64, GridapDistributedPETScWrappers.C.VECMPI, comm=MPI.COMM_WORLD);
  rhs.assembling=true
  resize!(rhs, mlocal=mat_size)

  # exact solution
  uex = Vec(Float64, GridapDistributedPETScWrappers.C.VECMPI, comm=MPI.COMM_WORLD);
  uex.assembling=true
  resize!(uex, mlocal=mat_size)

  # create a ghost vector that includes all points needed for each proc
  # to evaluate the rhs, duplicating the shared points.
  # Set up scattering from u_i to this vector

  # get global indices of local part of u_i
  idx_range = localpart(u_i)
  start_idx = idx_range[1] - 1
  end_idx = idx_range[end] - 1
  # figure out the indices in u_i corresponding to u_ghost
  if mpi_rank == 1  # left boundary (Dirchlet BC)
    idx = collect(start_idx:(end_idx+1))
    ghost_startidx = 0
    ghost_endidx = mat_size
    idx_ghost = collect(ghost_startidx:ghost_endidx)
  elseif mpi_rank == mpi_size  # right boundary (Neumann BC)
    idx = collect( (start_idx-1):end_idx)
    ghost_startidx = mat_size + (mpi_size-2)*(mat_size + 1)
    ghost_endidx = ghost_startidx + mat_size
    idx_ghost = collect(ghost_startidx:ghost_endidx)
  else  # interior
    idx = collect((start_idx-1):(end_idx+1))
    ghost_startidx = mat_size + 1 + (mpi_rank -2)*(mat_size + 2)
    ghost_endidx = ghost_startidx + mat_size + 1
    idx_ghost = collect(ghost_startidx:ghost_endidx)
  end
  # bump everything up to 1-based indexing
  idx += 1
  idx_ghost += 1

  ghost_offset = 2*(mpi_rank-1) # offset from u_i indices to ghost indices
  # set up the vector
  nghost_local = length(idx_ghost)
  u_ghost = Vec(Float64, GridapDistributedPETScWrappers.C.VECMPI, comm=MPI.COMM_WORLD);
  resize!(u_ghost, mlocal=nghost_local)
  u_ghost.assembling=true

  # create the scatter
  is_local = GridapDistributedPETScWrappers.IS(Float64, idx, comm=MPI.COMM_WORLD)
  is_ghost = GridapDistributedPETScWrappers.IS(Float64, idx_ghost, comm=MPI.COMM_WORLD)
  vec_scatter = GridapDistributedPETScWrappers.VecScatter(u_i, is_local, u_ghost, is_ghost)

  ctx = (is_local, is_ghost)  # avoid GC
  ksp = GridapDistributedPETScWrappers.KSP(A, ksp_atol=1e-12, ksp_rtol=1e-14, ksp_monitor="")

  return A, u_i, u_ghost, vec_scatter, rhs, ksp, uex, ctx, ghost_offset
end

# run
driver()


text:
using PyPlot
using LinearAlgebra: norm, dot

a = [ 0.5  -1.0   0.25
      1.0   0.5  -1.0 ]
c = (1/3) * sum(a, dims=2)

figure(1)

function level_sets(a, j, ξ)
    prev = [3, 1, 2]
    succ = [2, 3, 1]
    m = length(ξ)
    start = zeros(2, m)
    finish = zeros(2, m)
    for k = 1:m
        η = -0.3
        start[:,k] = ( a[:,prev[j]] + ξ[k] * ( a[:,j] - a[:,prev[j]] )
                      + η * ( a[:,succ[j]] - a[:,prev[j]] ) )
        η = 1.1 - ξ[k]
        finish[:,k] = ( a[:,prev[j]] + ξ[k] * ( a[:,j] - a[:,prev[j]] )
                      + η * ( a[:,succ[j]] - a[:,prev[j]] ) )
    end
    return start, finish
end

function draw_level_sets(a, j)
    offset = [ 0.1   0.0  -0.1
               0.0   0.1  -0.2 ]
    s = latexstring("\\xi_", j, "=")
    vals = ["0", "1/3", "2/3", "1"]
    ha = ["left", "center", "right"]
    va = ["center", "center", "center"]
    ξ = [0, 1/3, 2/3, 1]
    start, finish = level_sets(a, j, ξ)
    for k = 1:4
        plot([ start[1,k], finish[1,k] ], [ start[2,k], finish[2,k] ], "b",
            linewidth=0.5)
        text(start[1,k]+offset[1,j], start[2,k]+offset[2,j],
             latexstring(s, vals[k]), color="b",
             horizontalalignment=ha[j], verticalalignment=va[j])
    end
end

function draw_triangle(a, c, d)
    J = [1, 2, 3, 1]
    plot(a[1,J], a[2,J], "-k")
    for j = 1:3
        s = latexstring("a_", j)
        offset = a[:,j] - c
        offset *= d / norm(offset)
        text(a[1,j]+offset[1], a[2,j]+offset[2], s, color="r",
            horizontalalignment="center", verticalalignment="center")
    end
    plot(a[1,:], a[2,:], "ro", c[1], c[2], "ko")
end

for j = 1:3
   draw_level_sets(a, j)
end
draw_triangle(a, c, 0.2)
text(c[1]+0.06, c[2]-0.02, L"c",
     horizontalalignment="left", verticalalignment="center")
axis("equal")
axis("off")
savefig("barycentric.pdf")

figure(2)

function draw_vectors(a, b, c)
    prev = [3, 1, 2]
    succ = [2, 3, 1]
    quiver([ c[1], c[1], c[1] ], [ c[2], c[2], c[2] ], b[1,:], b[2,:],
           angles="xy", scale_units="xy", scale=1,
           width=0.004, headwidth=5, color="g")
    offset = [ -0.1   0.1   -0.1
                0.05  0.07   0.0 ]
    for j = 1:3
        text(c[1]+b[1,j]+offset[1,j], c[2]+b[2,j]+offset[2,j],
                 latexstring("b_", j), color="g",
                 horizontalalignment="center", verticalalignment="center")
        bhat = b[:,j] / norm(b[:,j])
        λ = dot(c - a[:,succ[j]], bhat)
        foot = c - λ * bhat
        plot([c[1], foot[1]], [c[2], foot[2]], "--k")
        side = a[:,succ[j]] - a[:,prev[j]]
        v = side / norm(side)
        plot(foot[1] .+ 0.1*[ bhat[1], bhat[1]+v[1], v[1]],
             foot[2] .+ 0.1*[ bhat[2], bhat[2]+v[2], v[2]], linewidth=0.5, "-k")
    end
end

A = [ (a[:,1]-a[:,3]) (a[:,2]-a[:,3]) ]
B = inv(A')
b = [ B[:,1] B[:,2] -sum(B, dims=2) ]

draw_vectors(a, b, c)
draw_triangle(a, c, 0.15)
axis("equal")
axis([-2.4, 1.4, -2.4, 2.4])
axis("off")
savefig("b_vectors.pdf")

figure(3)

function draw_quadratic_triangle(a, c, d, sz=12)
    J = [1, 2, 3, 1]
    plot(a[1,J], a[2,J], "-k")
    m = zeros(2, 3)
    m[:,1] = ( a[:,3] + a[:,2] ) / 2
    m[:,2] = ( a[:,1] + a[:,3] ) / 2
    m[:,3] = ( a[:,2] + a[:,1] ) / 2
    for j = 1:3
        s = latexstring("a_", j)
        offset = a[:,j] - c
        offset *= d / norm(offset)
        text(a[1,j]+offset[1], a[2,j]+offset[2], s, color="r",
            horizontalalignment="center", verticalalignment="center",
            fontsize=sz)
        s = latexstring("m_", j)
        offset = m[:,j] - c
        offset *= d / norm(offset)
        text(m[1,j]+offset[1], m[2,j]+offset[2], s, color="r",
            horizontalalignment="center", verticalalignment="center",
            fontsize=sz)
    end
    plot(a[1,:], a[2,:], "ro", m[1,:], m[2,:], "ro")
end

draw_quadratic_triangle(a, c, 0.15)
axis("equal")
axis("off")
savefig("midpoints.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
module NLSat

export NLSolver, Z3Solver, YicesSolver, SMTSolver
export NLStatus, NLModel
export variables!, constraints!, solve

using DelimitedFiles
using Distributed
using MacroTools: walk, postwalk, @capture, @match, replace
using Dates

const NLModel = Dict{Symbol,Number}

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

abstract type AlgebraicNumber end
export AlgebraicNumber

function AlgebraicNumber(x::Expr, n::Int)
    @warn "Algebraic numbers not yet supported, got $((x, n)), returning $(NaN)"
    NaN
end

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

include("../utils.jl")
include("clauseset.jl")
include("smtlib.jl")
include("lisp.jl")

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

@enum NLStatus sat unsat unknown timeout

abstract type NLSolver end

function variables!(s::NLSolver, d::Dict{Symbol,Type}) end
function solve(s::NLSolver; timeout::Int = -1) end
constraints!(s::NLSolver, cs::ClauseSet) = s.cs &= cs

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

# include("mathematica.jl")
include("smtsolver.jl")

end # module


text:
# Problem 1
println("Problem 1:")
reduce(+, 0, filter(x -> x % 3 == 0 || x % 5 == 0, 1:999))

# Problem 2
# It seams that the problem is wrong.
# The Fib Seq starts with 1, 1, 2, 3, ...
println("Problem 2:")
a, b, s = 1, 1, 0
while a < 4000000
    if iseven(a)
        s += a
    end
    a, b = b, a+b
end
println(s)

# Problem 3
println("Problem 3:")
function  next_prime(primes)
    x = maximum(primes) + 1
    while true
        found_prime = true
        for prime in primes
            if x % prime == 0
                found_prime = false
                break
            end
        end

        if found_prime
            break
        end

        x = x + 1
    end
    return x
end

function p3()
    n = 600851475143
    primes = IntSet([2])
    prime = 2

    while n > 1.0
        if n % prime == 0
            n = n / prime
        else
            prime = next_prime(primes)
            union!(primes, [prime])
        end
    end
    println(prime)
end
@time p3()

max_stars_repo_name: holmescn/practice

max_issues_repo_name: holmescn/practice

max_forks_repo_name: holmescn/practice


text:
###############################################################################
#
# Create LeafOpticals
#
###############################################################################
"""
    create_leaf_opticals(
                sWL::Array{FT,1},
                file::String = OPTI_2021
    ) where {FT<:AbstractFloat}

Create an `AbstractLeafOptiPara` struct, given
- `sWL` Standard wave length
- `opfn` File that saves optical parameters
"""
function create_leaf_opticals(
            sWL::Array{FT,1},
            opfn::String = OPTI_2021
) where {FT<:AbstractFloat}
    N = length(sWL) - 1;

    # reading data
    _opti   = matread(opfn)["optipar"];
    _nr     = _opti["nr"  ];
    _Km     = _opti["Kdm" ];
    _Kab    = _opti["Kab" ];
    _Kant   = _opti["Kant"];
    _Kcar   = _opti["Kca" ];
    _Kw     = _opti["Kw"  ];
    _KBrown = _opti["Ks"  ];
    _phi    = _opti["phi" ];
    _KcaV   = _opti["KcaV"];
    _KcaZ   = _opti["KcaZ"];
    _lambda = _opti["wl"  ];

    # create data to parse
    nr     = zeros(FT, N);
    Km     = zeros(FT, N);
    Kab    = zeros(FT, N);
    Kant   = zeros(FT, N);
    Kcar   = zeros(FT, N);
    Kw     = zeros(FT, N);
    KBrown = zeros(FT, N);
    phi    = zeros(FT, N);
    KcaV   = zeros(FT, N);
    KcaZ   = zeros(FT, N);
    lambda = zeros(FT, N);

    # fill in the data arrays
    # println("Reading Optical Parameters from ", sWL[1], " to ", sWL[end])
    @inbounds for i in 1:N
        wo = findall( (_lambda.>=sWL[i]) .& (_lambda.<sWL[i+1]) );
        if length(wo)==0
            @warn "Warning, some wavelengths out of bounds $(string(sWL[i]))";
        end

        nr[i]     = mean(    _nr[wo]);
        Km[i]     = mean(    _Km[wo]);
        Kab[i]    = mean(   _Kab[wo]);
        Kant[i]   = mean(  _Kant[wo]);
        Kcar[i]   = mean(  _Kcar[wo]);
        Kw[i]     = mean(    _Kw[wo]);
        KBrown[i] = mean(_KBrown[wo]);
        phi[i]    = mean(   _phi[wo]);
        KcaV[i]   = mean(  _KcaV[wo]);
        KcaZ[i]   = mean(  _KcaZ[wo]);
        lambda[i] = mean(_lambda[wo]);
    end

    # return the created struct
    return LeafOpticals{FT}(nr,
                            Km,
                            Kab,
                            Kant,
                            Kcar,
                            Kw,
                            KBrown,
                            phi,
                            KcaV,
                            KcaZ,
                            lambda)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Revise
using AdFem
using LinearAlgebra
using PyPlot
using DelimitedFiles

function calc_residual_and_jacobian(θ, u)
    ψ, J = neo_hookean(u, μ, λ, mmesh)
    B = compute_fem_source_term(zeros(get_ngauss(mmesh)), -0.5*ones(get_ngauss(mmesh)), mmesh)
    ψ = ψ - B
    J, ψ = impose_Dirichlet_boundary_conditions(J, ψ, bdnode, zeros(length(bdnode)))
    ψ, J
end

mmesh = Mesh(8, 8 , 1/8)

E, nu = 10.0, 0.3
μ, λ = E/(2*(1 + nu)) * ones(get_ngauss(mmesh)), E*nu/((1 + nu)*(1 - 2*nu)) * ones(get_ngauss(mmesh))
bdnode = bcnode((x,y)->(x<1e-5) || (x>1-1e-5), mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nr = newton_raphson(calc_residual_and_jacobian, zeros(2mmesh.ndof), missing)

ADCME.options.newton_raphson.verbose = true
sess = Session(); init(sess)
nr = run(sess, nr)

close("all")
# visualize_displacement(nr.x, mmesh)
u = nr.x[1:mmesh.nnode]
v = nr.x[mmesh.ndof+1:mmesh.ndof + mmesh.nnode]
visualize_scalar_on_fem_points(sqrt.(u.^2+v.^2), mmesh)
savefig("test_AdFem.png")


text:
"""
    getnext(x[, ::Val{Tag}])

Get the next node in a list after `x`. Defaults to `x.next`.
`Tag` identifies the list instance, for nodes that may exist in more than one list.
For implementations of nodes that exist in only a single list at a time,
`getnext` can be implemented without the `Tag` parameter.
"""
@inline getnext(x, ::Val{:default}) = getnext(x)
@inline getnext(x) = x.next

"""
    getprev(x[, ::Val{Tag}])

Get the previous node in a list before `x`. Defaults to `x.prev`.
`Tag` identifies the list instance, for nodes that may exist in more than one list.
For implementations of nodes that exist in only a single list at a time,
`getprev` can be implemented without the `Tag` parameter.
"""
@inline getprev(x, ::Val{:default}) = getprev(x)
@inline getprev(x) = x.prev

"""
    setnext!(x, next, [, ::Val{Tag}])

Link `x` to `next` after `x`. Defaults to `x.next = next`.
`Tag` identifies the list instance, for nodes that may exist in more than one list.
For implementations of nodes that exist in only a single list at a time,
`getnext` can be implemented without the `Tag` parameter.
"""
@inline setnext!(x::T, next::T, ::Val{:default}) where {T} = setnext!(x, next)
@inline setnext!(x::T, next::T) where {T} = x.next = next

"""
    setprev!(x, prev, [, ::Val{Tag}])

Link `x` to `prev` before `x`. Defaults to `x.prev = prev`.
`Tag` identifies the list instance, for nodes that may exist in more than one list.
For implementations of nodes that exist in only a single list at a time,
`getnext` can be implemented without the `Tag` parameter.
"""
@inline setprev!(x::T, prev::T, ::Val{:default}) where {T} = setprev!(x, prev)
@inline setprev!(x::T, prev::T) where {T} = x.prev = prev


text:
using Pkg

using YAXArrays
using BenchmarkTools

const SUITE = BenchmarkGroup()
SUITE["mapslices"] = include("bench_mapslices.jl")

max_forks_repo_name: meggart/YAXArrays


text:
abstract type MoistureBC end

"""
    Impermeable() :: MoistureBC

No moisture flux.
"""
struct Impermeable <: MoistureBC end
function atmos_moisture_boundary_state!(
    nf,
    bc_moisture::Impermeable,
    atmos,
    args...,
) end
function atmos_moisture_normal_boundary_flux_diffusive!(
    nf,
    bc_moisture::Impermeable,
    atmos,
    args...,
) end

"""
    PrescribedMoistureFlux(fn) :: MoistureBC

Prescribe the net inward moisture flux across the boundary by `fn`, a function
with signature `fn(state, aux, t)`, returning the flux (in kg/m^2).
"""
struct PrescribedMoistureFlux{FN} <: MoistureBC
    fn::FN
end
function atmos_moisture_boundary_state!(
    nf,
    bc_moisture::PrescribedMoistureFlux,
    atmos,
    args...,
) end
function atmos_moisture_normal_boundary_flux_diffusive!(
    nf,
    bc_moisture::PrescribedMoistureFlux,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    nρd_q_tot = -bc_moisture.fn(state⁻, aux⁻, t)
    fluxᵀn.ρ += nρd_q_tot
    fluxᵀn.ρu += nρd_q_tot / state⁻.ρ .* state⁻.ρu
    # assumes EquilMoist
    fluxᵀn.moisture.ρq_tot += nρd_q_tot
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Test

examples = ["simple_hybrid_ME", "simple_hybrid_CS", "modelica_conference_2021", "advanced_hybrid_ME"]

@testset "FMIFlux.jl Examples" begin
    for example in examples
        @testset "$(example).jl" begin
            path = joinpath(dirname(@__FILE__), "..", "example", example * ".jl")
            @test include(path) == nothing
        end
    end
end


text:
function __init__()
    check_deps()

    Libdl.dlopen(libmkl_core, Libdl.RTLD_GLOBAL)
    Libdl.dlopen(libmkl_rt, Libdl.RTLD_GLOBAL) # maybe only needed on mac
    Libdl.dlopen(libmkl_vml_avx, Libdl.RTLD_GLOBAL)

end

__init__()

struct VMLAccuracy
    mode::UInt
end

const VML_LA = VMLAccuracy(0x00000001)
const VML_HA = VMLAccuracy(0x00000002)
const VML_EP = VMLAccuracy(0x00000003)

Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
                                              m == VML_HA ? "VML_HA" : "VML_EP")
vml_get_mode() = ccall((:_vmlGetMode, libmkl_vml_avx), Cuint, ())
vml_set_mode(mode::Integer) = (ccall((:_vmlSetMode, libmkl_vml_avx), Cuint, (UInt,), mode); nothing)

vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)

vml_set_mode((vml_get_mode() & ~0x0000FF00))
function vml_check_error()
    vml_error = ccall((:_vmlClearErrStatus, libmkl_vml_avx), Cint, ())
    if vml_error != 0
        if vml_error == 1
            throw(DomainError(-1, "This function does not support arguments outside its domain"))
        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
            # Singularity, overflow, or underflow
            # I don't think Base throws on these
        elseif vml_error == 1000
            warn("IntelVectorMath does not support $(vml_get_accuracy); lower accuracy used instead")
        else
            error("an unexpected error occurred in IntelVectorMath ($vml_error)")
        end
    end
end

function vml_prefix(t::DataType)
    if t == Float32
        return "_vmls"
    elseif t == Float64
        return "_vmld"
    elseif t == Complex{Float32}
        return "_vmlc"
    elseif t == Complex{Float64}
        return "_vmlz"
    end
    error("unknown type $t")
end

function def_unary_op(tin, tout, jlname, jlname!, mklname;
        vmltype = tin)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(vmltype))$mklname"))
    exports = Symbol[]
    (@isdefined jlname) || push!(exports, jlname)
    (@isdefined jlname!) || push!(exports, jlname!)
    @eval begin
        function ($jlname!)(out::Array{$tout,N}, A::Array{$tin,N}) where {N}
            size(out) == size(A) || throw(DimensionMismatch())
            ccall(($mklfn, libmkl_vml_avx), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            return out
        end
        $(if tin == tout
            quote
                function $(jlname!)(A::Array{$tin})
                    ccall(($mklfn, libmkl_vml_avx), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, A)
                    vml_check_error()
                    return A
                end
            end
        end)
        function ($jlname)(A::Array{$tin})
            out = similar(A, $tout)
            ccall(($mklfn, libmkl_vml_avx), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            return out
        end
        $(isempty(exports) ? nothing : Expr(:export, exports...))
    end
end

function def_binary_op(tin, tout, jlname, jlname!, mklname, broadcast)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
    exports = Symbol[]
    (@isdefined jlname) || push!(exports, jlname)
    (@isdefined jlname!) || push!(exports, jlname!)
    @eval begin
        $(isempty(exports) ? nothing : Expr(:export, exports...))
        function ($jlname!)(out::Array{$tout,N}, A::Array{$tin,N}, B::Array{$tin,N}) where {N}
            size(out) == size(A) == size(B) || $(broadcast ? :(return broadcast!($jlname, out, A, B)) : :(throw(DimensionMismatch())))
            ccall(($mklfn, libmkl_vml_avx), Nothing, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            return out
        end
        function ($jlname)(A::Array{$tout,N}, B::Array{$tin,N}) where {N}
            size(A) == size(B) || $(broadcast ? :(return broadcast($jlname, A, B)) : :(throw(DimensionMismatch())))
            out = similar(A)
            ccall(($mklfn, libmkl_vml_avx), Nothing, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            return out
        end
    end
end