text:
"""Nice, but slower than the ForwardDiff approach"""
function _get_init_derivatives_mtk(prob, order)
    # Output of size order+1
    u0 = prob.u0
    d = length(u0)
    q = order
    out = fill(zero(u0[1]), d*(q+1))

    sys = modelingtoolkitize(prob)
    t = sys.iv()
    @derivatives D'~t
    u = [s(t) for s in sys.states]
    p = [_p() for _p in sys.ps]

    rhs = ModelingToolkit.rhss(sys.eqs)

    substitutions = [u .=> u0; t .=> prob.tspan[1]]
    if length(p) > 0
        substitutions = [substitutions; p .=> prob.p]
    end

    out[1:d] .= u0

    val = substitute.(rhs, (substitutions,))

    out[d+1:2d] .= [v.value for v in val]
    u = D.(u)
    substitutions = [D.(u) .=> val; substitutions]

    for i in 2:(order)
        rhs = ModelingToolkit.jacobian(rhs, [t])[:]
        val = expand_derivatives.(substitute.(rhs, (substitutions,)))
        out[(i*d)+1:(i+1)*d] .= [v.value for v in val]
        u = D.(u)
        substitutions = [u .=> val; substitutions]
    end
    return out
end

function get_initial_states_forwarddiff(prob::ODEProblem, order::Int)
    f = isinplace(prob) ? iip_to_oop(prob.f) : prob.f

    u0 = prob.u0
    t0 = prob.tspan[1]
    p = prob.p

    d = length(u0)
    q = order

    uElType = eltype(u0)
    m0 = zeros(uElType, d*(q+1))
    P0 = zeros(uElType, d*(q+1), d*(q+1))

    m0[1:d] .= u0
    m0[d+1:2d] .= f(u0, p, t0)

    f_derivatives = Function[f]
    for o in 2:q
        _curr_f_deriv = f_derivatives[end]
        dfdu(u, p, t) = ForwardDiff.jacobian(u -> _curr_f_deriv(u, p, t), u)
        dfdt(u, p, t) = ForwardDiff.derivative(t -> _curr_f_deriv(u, p, t), t)
        df(u, p, t) = dfdu(u, p, t) * f(u, p, t) + dfdt(u, p, t)
        push!(f_derivatives, df)
        m0[o*d+1:(o+1)*d] = df(u0, p, t0)
    end

    return m0, P0
end

function iip_to_oop(f!)
    function f(u, p, t)
        du = copy(u)
        f!(du, u, p, t)
        return du
    end
    return f
end

function remake_prob_with_jac(prob)
    try
        p = prob.p isa DiffEqBase.NullParameters ? [] : collect(prob.p)
        prob = remake(prob, p=p)
        sys = modelingtoolkitize(prob)
        jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
        IIP = isinplace(prob.f)
        f = ODEFunction{IIP}(
            prob.f.f,
            jac=jac;
            analytic=prob.f.analytic,
        )
        return remake(prob, f=f)
    catch
        error("Could not generate a jacobian for the problem")
    end
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

"""
    AbstractEvalMetric

The base class for all evaluation metrics. The sub-types should implement the following
interfaces:

* [`update!`](@ref)
* [`reset!`](@ref)
* [`get`](@ref)
"""
abstract type AbstractEvalMetric end

"""
    hasNDArraySupport(metric) -> Val{true/false}

Trait for `_update_single_output` should return `Val{true}() if metric can handle `NDArray`
directly and `Val{false}()` if requires `Array`. Metric that work with NDArrays can be
async, while native Julia arrays require that we copy the output of the network, which is
a blocking operation.
"""
hasNDArraySupport(::AbstractEvalMetric) = Val{true}()

"""
    update!(metric, labels, preds)

Update and accumulate metrics.

# Arguments:
* `metric::AbstractEvalMetric`: the metric object.
* `labels::Vector{NDArray}`: the labels from the data provider.
* `preds::Vector{NDArray}`: the outputs (predictions) of the network.
"""
update!(metric::T, labels::VecOfNDArray, preds::VecOfNDArray) where T<:AbstractEvalMetric =
  _update!(metric, labels, preds, hasNDArraySupport(metric))

function _update!(metric::T, labels::VecOfNDArray, preds::VecOfNDArray,
                  ::Val{true}) where T<:AbstractEvalMetric
  if length(labels) != length(preds)
    @warn(
      "The number of labels ($(length(labels))) does not correspond to the " *
      "number of outputs ($(length(preds))). The calculated metric might not be accuracte.",
      maxlog = 1)
  end
  for (label, pred) in zip(labels, preds)
    _update_single_output(metric, label, pred)
  end
end

function _update!(metric::T, labels::VecOfNDArray, preds::VecOfNDArray,
                  ::Val{false}) where {T<:AbstractEvalMetric}
  if length(labels) != length(preds)
    @warn(
      "The number of labels ($(length(labels))) does not correspond to the " *
      "number of outputs ($(length(preds))). The calculated metric might not be accuracte.",
      maxlog = 1)
  end
  for (label, pred) in zip(labels, preds)
     @nd_as_jl ro=(label, pred) begin
       # This is a dynamic dispatch since the conversion from NDArray to
       # Array is not type-stable.
      _update_single_output(metric, label, pred)
    end
  end
end

"""
    reset!(metric)

Reset the accumulation counter.
"""
reset!(metric::AbstractEvalMetric) = throw(MethodError(reset!, (typeof(metric),)))

import Base: get
"""
    get(metric)

Get the accumulated metrics.

Returns `Vector{Tuple{Base.Symbol, Real}}`, a list of name-value pairs.
For example, `[(:accuracy, 0.9)]`.
"""
get(metric::AbstractEvalMetric) = throw(MethodError(get, (typeof(metric),)))

"""
    NullMetric()

A metric that calculates nothing. Can be used to ignore an output during training.
"""
mutable struct NullMetric <: mx.AbstractEvalMetric
end

update!(metric::NullMetric, labels::VecOfNDArray, preds::VecOfNDArray) = nothing

reset!(metric::NullMetric) = nothing

get(metric::NullMetric) = Tuple{Symbol, Float64}[]

"""
    MultiMetric(metrics::Vector{AbstractEvalMetric})

Combine multiple metrics in one and get a result for all of them.

# Usage
To calculate both mean-squared error [`Accuracy`](@ref) and log-loss [`ACE`](@ref):
```julia
  mx.fit(..., eval_metric = mx.MultiMetric([mx.Accuracy(), mx.ACE()]))
```
"""
mutable struct MultiMetric <: AbstractEvalMetric
  metrics :: Vector{mx.AbstractEvalMetric}
end

function update!(metric :: MultiMetric, labels :: Vector{<:NDArray}, preds :: Vector{<:NDArray})
  for m in metric.metrics
    update!(m, labels, preds)
  end
  nothing
end

function reset!(metric :: MultiMetric)
  map(reset!, metric.metrics)
  nothing
end

get(metric::MultiMetric) = mapreduce(get, append!, metric.metrics)

"""
    SeqMetric(metrics::Vector{AbstractEvalMetric})

Apply a different metric to each output. This is especially useful for `mx.Group`.

# Usage
Calculate accuracy [`Accuracy`](@ref) for the first output
and log-loss [`ACE`](@ref) for the second output:
```julia
  mx.fit(..., eval_metric = mx.SeqMetric([mx.Accuracy(), mx.ACE()]))
```
"""
mutable struct SeqMetric <: AbstractEvalMetric
  metrics :: Vector{AbstractEvalMetric}
end

function update!(metric::SeqMetric, labels::VecOfNDArray, preds::VecOfNDArray)
  @assert length(metric.metrics) == length(labels)
  @assert length(metric.metrics) == length(preds)
  for (m, l, p) in zip(metric.metrics, labels, preds)
    update!(m, [l], [p])
  end
  nothing
end

function reset!(metric::SeqMetric)
  map(reset!, metric.metrics)
  nothing
end

get(metric::SeqMetric) = mapreduce(get, append!, metric.metrics)

"""
    Accuracy

Multiclass classification accuracy.

Calculates the mean accuracy per sample for softmax in one dimension.
For a multi-dimensional softmax the mean accuracy over all dimensions is calculated.
"""
mutable struct Accuracy <: AbstractEvalMetric
  acc_sum  :: Float64
  n_sample :: Int

  Accuracy() = new(0.0, 0)
end

hasNDArraySupport(::Accuracy) = Val{false}()

function _update_single_output(metric::Accuracy, label::Array, pred::Array)
  # Samples are stored in the last dimension
  @assert size(label, ndims(label)) == size(pred, ndims(pred))

  if ndims(pred) == 4 # Multidimensional case
    # Reshape label to be of the same shape as pred.
    # Except for the third dimension where the predictions are stored.
    labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))

    for sample in 1:size(labels, 4)
      for j in 1:size(labels, 2)
        for i in 1:size(labels, 1)
          label = labels[i, j, 1, sample]
          klasses = view(pred, i, j, :, sample)
          klass = argmax(klasses) - 1 # Classes start at 0...k-1

          metric.acc_sum += klass == label
          metric.n_sample += 1
        end
      end
    end
  elseif ndims(pred) == 2 # 1-dimensional case
    for sample in 1:size(label, 1)
      klass = argmax(view(pred, :, sample)) - 1
      metric.acc_sum += klass == label[sample]
      metric.n_sample += 1
    end
  else
    error("Can't handle prediction with dimensions $(ndims(pred)).")
  end
end

get(metric::Accuracy) = [(:accuracy, metric.acc_sum / metric.n_sample)]

function reset!(metric :: Accuracy)
  metric.acc_sum  = 0.0
  metric.n_sample = 0
end

"""
    MSE

Mean Squared Error.

Calculates the mean squared error regression loss.
Requires that label and prediction have the same shape.
"""
mutable struct MSE{N} <: AbstractEvalMetric
  mse_sum  :: Vector{NDArray{MX_float,N}}
  n_sample :: Int

  MSE{N}() where {N} = new(Vector{NDArray{MX_float,N}}(), 0)
end

MSE() = MSE{1}()  # backward compat?

hasNDArraySupport(::MSE) = Val{true}()

function _update_single_output(metric::MSE, label::NDArray{T,N},
                               pred::NDArray{T,N}) where {T,N}
  @assert size(label) == size(pred)
  metric.n_sample += length(label)
  mse_sum = mx.sum((label .- pred).^2)
  push!(metric.mse_sum, mse_sum)
  nothing
end

function get(metric::MSE)
  # Delay copy until last possible moment
  mse_sum = mapreduce(nda->copy(nda)[1], +, 0.0, metric.mse_sum)
  [(:MSE, mse_sum / metric.n_sample)]
end

function reset!(metric::MSE{N}) where N
  metric.mse_sum = Vector{NDArray{Float32,N}}()
  metric.n_sample = 0
end

@doc doc"""
    NMSE

Normalized Mean Squared Error

```math
\sum_i (\frac{label_i - pred_i}{label_i})^2
```

Note that there are various ways to do the *normalization*.
It depends on your own context. Please judge the problem setting you have
first. If the current implementation do not suitable for you,
feel free to file it on GitHub.

Let me show you a use case of this kind of normalization:

Bob is training a network for option pricing. The option pricing problem is
a regression problem (pirce predicting). There are lots of option contracts
on same target stock but different strike price.
For example, there is a stock `S`; it's market price is 1000.
And, there are two call option contracts with different strike price.
Assume Bob obtains the outcome as following table:

```
+--------+----------------+----------------+--------------+
|        | Strike Price   | Market Price   | Pred Price   |
+--------+----------------+----------------+--------------+
| Op 1   | 1500           |  100           | 80           |
+--------+----------------+----------------+--------------+
| Op 2   | 500            |  10            | 8            |
+--------+----------------+----------------+--------------+
```

Now, obviously, Bob will calculate the normalized MSE as:

```math
    (\frac{100 - 80}{100})^2
    \text{ vs }
    (\frac{10 - 8}{10}) ^2
```

Both of the pred prices got the same degree of error.

For more discussion about normalized MSE, please see
[#211](https://github.com/dmlc/MXNet.jl/pull/211) also.

"""
mutable struct NMSE <: AbstractEvalMetric
  nmse_sum  :: Float64
  n_sample :: Int

  NMSE() = new(0.0, 0)
end

hasNDArraySupport(::NMSE) = Val{false}()

function _update_single_output(metric::NMSE, label::Array, pred::Array)
  n_sample = size(pred)[end]
  metric.n_sample += n_sample

  for i = 1:n_sample
    if label[i] == 0.0f0  # in case of batch padding
        continue
    end

    metric.nmse_sum += ((label[i] - pred[i]) / label[i])^2
  end
end

get(metric::NMSE) = [(:NMSE, metric.nmse_sum / metric.n_sample)]

function reset!(metric::NMSE)
  metric.nmse_sum = 0.0
  metric.n_sample = 0
end

"""
    ACE

Calculates the averaged cross-entropy (logloss) for classification.

# Arguments:
* `eps::Float64`: Prevents returning `Inf` if `p = 0`.
"""
mutable struct ACE <: AbstractEvalMetric
  ace_sum  :: Float64
  n_sample :: Int
  eps :: Float64

  ACE(eps=1.0e-8) = new(0.0, 0, eps)
end

get(metric::ACE) = [(:ACE, - metric.ace_sum / metric.n_sample)]

function reset!(metric::ACE)
  metric.ace_sum = 0.0
  metric.n_sample = 0
end

hasNDArraySupport(::ACE) = Val{false}()

function _update_single_output(metric :: ACE, label :: Array{T}, pred :: Array{T}) where T
  eps = convert(T, metric.eps)
  # Samples are stored in the last dimension
  @assert size(label, ndims(label)) == size(pred, ndims(pred))
  if size(label) == size(pred) # simply calculate the cross entropy of the probabilities
    for (q, p) in zip(pred, label)
      # p == true probability
      # q == "unnatural" probability
      metric.ace_sum += p * log(q + eps)
      metric.n_sample += 1
    end
  elseif ndims(pred) == 4
    labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))
    for sample in 1:size(labels, 4)
      for j in 1:size(labels, 2)
        for i in 1:size(labels, 1)
          # Cross-entropy reduces to -(ln(p_1)*0 + ln(p_2)*1) for classification
          # Since we can only target labels right now this is the only thing we can do.
          target = Int(labels[i, j, 1, sample]) + 1 # klasses are 0...k-1 => julia indexing
          p_k = pred[i, j, target, sample]
          metric.ace_sum += log(p_k + eps)
          metric.n_sample += 1
        end
      end
    end
  elseif ndims(pred) == 2 # 1-dimensional case
    for sample in 1:size(label, 1)
      target = Int(label[sample]) + 1    # 0-based indexing => 1-based indexing
      p_k = pred[target, sample]
      metric.ace_sum += log(p_k +eps)
      metric.n_sample += 1
    end
  else
    error("Can't handle prediction with dimensions $(ndims(pred)).")
  end
end

"""
    MultiACE

Calculates the averaged cross-entropy per class and overall (see [`ACE`](@ref)).
This can be used to quantify the influence of different classes on the overall loss.
"""
mutable struct MultiACE <: AbstractEvalMetric
  aces  :: Vector{Float64}
  counts :: Vector{Int}
  eps :: Float64

  MultiACE(nclasses, eps=1.0e-8) = new(Base.zeros(nclasses), Base.zeros(Int, nclasses), eps)
end

function get(metric :: MultiACE)
  aces = [(Symbol("ACE_$(i-0)"), - metric.aces[i] / metric.counts[i]) for i in 1:length(metric.aces)]
  push!(aces, (:ACE, - Base.sum(metric.aces) / Base.sum(metric.counts)))
  return aces
end

function reset!(metric :: MultiACE)
  metric.aces = Base.zero(metric.aces)
  metric.counts = Base.zero(metric.counts)
end

hasNDArraySupport(::MultiACE) = Val{false}()

function _update_single_output(metric :: MultiACE, label :: Array{T}, pred :: Array{T}) where T
  eps = convert(T, metric.eps)
  # Samples are stored in the last dimension
  @assert size(label, ndims(label)) == size(pred, ndims(pred))
  @assert size(metric.aces) == size(metric.counts)
  if size(label) == size(pred) # simply calculate the cross entropy of the probabilities
    for k in 1:length(metric.aces)
      kpred  = view(pred,  ntuple(d->:, ndims(pred)  - 2)..., k, :)
      klabel = view(label, ntuple(d->:, ndims(label) - 2)..., k, :)
      for (q, p) in zip(kpred, klabel)
        # p == true probability
        # q == "unnatural" probability
        metric.aces[k] += p * log(q + eps)
        metric.counts[k] += 1
      end
    end
  elseif ndims(pred) == 4
    labels = reshape(label, size(pred, 1, 2)..., 1, size(pred, 4))
    for sample in 1:size(labels, 4)
      for j in 1:size(labels, 2)
        for i in 1:size(labels, 1)
          # Cross-entropy reduces to -(ln(p_1)*0 + ln(p_2)*1) for classification
          # Since we can only target labels right now this is the only thing we can do.
          target = Int(labels[i, j, 1, sample]) + 1 # klasses are 0...k-1 => julia indexing
          p_k = pred[i, j, target, sample]

          metric.aces[target] += log(p_k + eps)
          metric.counts[target] += 1
        end
      end
    end
  elseif ndims(pred) == 2
    for sample in 1:size(label, 1)
      target = Int(label[sample]) + 1
      p_k = pred[target, sample]
      metric.aces[target] += log(p_k + eps)
      metric.counts[target] += 1
    end
  else
    error("Can't handle prediction with dimensions $(ndims(pred)).")
  end
end

ext: jl

max_stars_repo_name: paulk-asert/incubator-mxnet

max_issues_repo_name: paulk-asert/incubator-mxnet

max_forks_repo_name: ElaineBao/incubator-mxnet


text:
# Included in benchmark_prepare.jl and in browsertests.jl
# Refers logutils
if !isdefined(:SRCPATH)
    const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") :Base.source_dir()
    const LOGGINGPATH = realpath(joinpath(SRCPATH, "../logutils/"))
    SRCPATH ∉ LOAD_PATH && push!(LOAD_PATH, SRCPATH)
    LOGGINGPATH ∉ LOAD_PATH && push!(LOAD_PATH, LOGGINGPATH)
end
using logutils_ws

"A list of potentially available browsers, to be tried in succession if present"
const BROWSERS = ["chrome", "firefox", "iexplore", "safari", "phantomjs"]
"An complicated browser counter."
mutable struct Countbrowser;value::Int;end
(c::Countbrowser)() =COUNTBROWSER.value += 1
"For next value: COUNTBROWSER(). For current value: COUNTBROWSER.value"
const COUNTBROWSER = Countbrowser(0)
const PORT = [8000]
const PAGE = ["bce.html"]
const URL = ["http://127.0.0.1:$(PORT[1])/$(PAGE[1])"]

"Get application path for developer applications"
function fwhich(s)
    fi = ""
    if Sys.is_windows()
        try
            fi = split(readstring(`where.exe $s`), "\r\n")[1]
            if !isfile(fi)
                fi = ""
            end
        catch
            fi =""
        end
    else
        try
            fi = readchomp(`which $s`)
        catch
            fi =""
        end
    end
    fi
end
function browser_path_unix_apple(shortname)
    trypath = ""
    if shortname == "chrome"
        if Base.Sys.is_apple()
            return "Google Chrome"
        else
            return "google-chrome"
        end
    end
    if shortname == "firefox"
        return "firefox"
    end
    if shortname == "safari"
        if Base.Sys.is_apple()
            return "safari"
        else
            return ""
        end
    end
    if shortname == "phantomjs"
        return fwhich(shortname)
    end
    return ""
end
function browser_path_windows(shortname)
    # windows accepts English paths anywhere.
    # Forward slash is acceptable too.
    # Searches C and D drives....
    trypath = ""
    homdr = ENV["HOMEDRIVE"]
    path32 = homdr * "/Program Files (x86)/"
    path64 = homdr * "/Program Files/"
    if shortname == "chrome"
        trypath = path64 * "Chrome/Application/chrome.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Google/Chrome/Application/chrome.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "firefox"
        trypath = path64 * "Mozilla Firefox/firefox.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Mozilla Firefox/firefox.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "safari"
        trypath = path64 * "Safari/Safari.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Safari/Safari.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "iexplore"
        trypath = path64 * "Internet Explorer/iexplore.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Internet Explorer/iexplore.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "phantomjs"
        return fwhich(shortname)
    end
    return ""
end
"Constructs launch command"
function launch_command(shortbrowsername)
    if Sys.is_windows()
        pt = browser_path_windows(shortbrowsername)
    else
        pt = browser_path_unix_apple(shortbrowsername)
    end
    pt == "" && return ``
    if shortbrowsername == "iexplore"
        prsw = "-private"
    else
        prsw = "--incognito"
    end
    if shortbrowsername == "phantomjs"
        if isdefined(:SRCPATH)
            script = joinpath(SRCPATH, "phantom.js")
        else
            script = "phantom.js"
        end
        return Cmd(`$pt $script $URL`)
    else
        if Sys.is_windows()
            return Cmd( [ pt, URL, prsw])
        else
            if Sys.is_apple()
                return Cmd(`open --fresh -n $URL -a $pt --args $prsw`)
            elseif Sys.is_linux() || Sys.is_bsd()
                return Cmd(`xdg-open $(URL) $pt`)
            end
        end
    end
end

function open_testpage(shortbrowsername)
    id = "open_testpage"
    dmc = launch_command(shortbrowsername)
    if dmc == ``
        clog(id, "Could not find " * shortbrowsername)
        return false
    else
        try
            if shortbrowsername == "phantomjs"
                # Run enables text output of phantom messages in the REPL. In Windows
                # standalone REPL, run will freeze the main thread if not run async.
                @async run(dmc)
            else
                spawn(dmc)
            end
        catch
            clog(id, :red, "Failed to spawn " * shortbrowsername)
            return false
        end
    end
    return true
end

"Try to open one browser from BROWSERS.
In some cases we expect an immediate indication
of failure, for example when the corresponding browser
is not found on the system. In other cases, we will
just wait in vain. In those cases,
call this function again after a reasonable timeout.
The function remembers which browsers were tried before.
"
function open_a_browser()
    id = "open_next_browser"
    if COUNTBROWSER.value > length(BROWSERS)
        return false
    end
    success = false
    b = ""
    while COUNTBROWSER.value < length(BROWSERS) && !success
        b = BROWSERS[COUNTBROWSER()]
        clog(id, "Trying to launch browser no. ", COUNTBROWSER.value, ", ", :yellow, b)
        success = open_testpage(b)
    end
    success && clog(id, "seems to work:", :yellow, b, :normal, " on ", URL)
    success, b
end


text:
using JuliaChem

using Documenter, JuliaChem.JCInput, JuliaChem.JCBasis, JuliaChem.JCRHF

makedocs(
    authors="David Poole",
    sitename="JuliChem Documentation",
    modules = [JCInput,
               JCBasis,
               JCRHF],
    #format = LaTeX()
)


text:
export
    dissect,
    Diagonal,
    partialdervativeof,
    Jacobian

@doc raw"""

    dissect(p::MPolyElem)

Iterator for the coefficients and exponent vectors of the given polynomial.

"""
function dissect(p::MPolyElem)
    return zip(coeffs(p), exponent_vectors(p))
end

function Nemo.nullspace_right_rational(N::AbstractArray{T}) where {T<:Integer}
    Nnemo = Nemo.matrix(Nemo.FlintZZ, N)
    r, U = Nemo.nullspace_right_rational(Nnemo)
    return convert_to_array(U[:,1:r],T(1))
end

"""

    matrix(R, V::AbstractVector{T}) where {T}

Given a Julia vector `V` of entries, construct the corresponding AbstractAlgebra.jl one-column matrix over the given ring `R`, assuming all the entries can be coerced into `R`.
"""
function Nemo.matrix(R::AbstractAlgebra.Ring, V::AbstractVector)
    return Nemo.matrix(R, reshape(V, (:, 1)))
end

@doc raw"""

    Diagonal(R, V::AbstractVector{T}) where {T}

Given a Julia vector `V` of entries, construct the corresponding AbstractAlgebra.jl diagonal matrix over the given ring `R`, assuming all the entries can be coerced into `R`.
"""
function Diagonal(R::AbstractAlgebra.Ring, V::AbstractVector{T}) where {T}
    n = size(V, 1)
    D = Nemo.zero_matrix(R, n, n)
    for i in 1:n D[i,i] = V[i] end
    return D
end

function partialdervativeof(i::Integer)
    return p -> Nemo.derivative(p, i)
end

"""

    Jacobian(R, M)
    Jacobian(R, M, vars=1:length(gens(R)))

Returns the Jacobian matrix of a one-column matrix of polynomials `M` with respect to the generators of `R` indexed by `vars`. When `vars` is omitted all the generators of `R` are used.

# Examples
```jldoctest; setup = :(using CRNT, Nemo)
julia> using Nemo

julia> R, vars = PolynomialRing(ZZ, vcat(["k\$i" for i in 1:5], ["x\$i" for i in 1:4]));

julia> M = [vars[1]*vars[2]*vars[6]-vars[8]; vars[3]*vars[9]+2*vars[7]]
[k1*k2*x1-x3]
[ k3*x4+2*x2]

julia> Jacobian(R, M, 6:9)
[k1*k2  0  -1   0]
[    0  2   0  k3]

julia> Jacobian(R,M)
[k2*x1  k1*x1   0  0  0  k1*k2  0  -1   0]
[    0      0  x4  0  0      0  2   0  k3]
```
"""
function Jacobian(R::MPolyRing, M::AbstractVector, vars=1:length(gens(R)))
    rows = size(M, 1)
    cols = length(vars)
    J = fill(zero(R), rows, cols)
    for j in 1:cols
        J[:,j] = partialdervativeof(vars[j]).(M)
    end
    return J
end


text:
# test rotatation tools

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
Nwin = 10
cc_len = 40.96 # length of noise window
cc_step = 40.96 # step between windows
azi = 320 / 180 * π
baz = 140 / 180 * π
maxlag = 20. # seconds of coda to save
starttime = d2u(DateTime(Date(now()))) # starttime for SeisChanel / SeisData

@testset "rotation" begin

    Ch = SeisIO.RandSeis.randSeisChannel(c=false,s=true)
    ungap!(Ch)
    Ch.x = rand(T,Int(cc_len*Nwin*fs +1)) .- T(0.5)
    Ch.fs = fs
    Ch.t = [1 starttime * 1e6;length(Ch.x) 0]
    Ch.loc = GeoLoc()
    R = RawData(Ch,cc_len,cc_step)
    F1s = [rfft(R) for ii in 1:3]
    F2s = [rfft(R) for ii in 1:3]

    # change FFTData names to ENZ naming convention
    F1s[1].name = "NN.STA1.00.HHE"
    F1s[2].name = "NN.STA1.00.HHN"
    F1s[3].name = "NN.STA1.00.HHZ"
    F2s[1].name = "NN.STA2.00.HHE"
    F2s[2].name = "NN.STA2.00.HHN"
    F2s[3].name = "NN.STA2.00.HHZ"

    # cross-correlate
    Cs = Array{CorrData}(undef,9)
    let  kk = 1
        for ii = 1:3
            for jj = 1:3
                Cs[kk] = correlate(F1s[ii],F2s[jj],maxlag)
                kk += 1
            end
        end
    end

    # rotate
    Crot = rotate(Cs,azi,baz)
    comps = [Crot[ii].comp for ii = 1:9]
    @test Set(comps) == Set(["RR","RT","RZ","TR","TT","TZ","ZR","ZT","ZZ"])
    @test all([Crot[ii].rotated == true for ii = 1:9])
    @test all([eltype(Crot[ii].corr) == T for ii = 1:9])

    # NEEDED: test with real data

    # test in-place
    rotate!(Cs,azi,baz)
    comps = [Cs[ii].comp for ii = 1:9]
    @test Set(comps) == Set(["RR","RT","RZ","TR","TT","TZ","ZR","ZT","ZZ"])
    @test all([Cs[ii].rotated == true for ii = 1:9])
    @test all([Cs[ii].corr == Crot[ii].corr for ii = 1:9])
    @test all([eltype(Cs[ii].corr) == T for ii = 1:9])
end


text:
function to_real_field_bivariate_function(obj::FieldBivariateFunction, arg0::Field)
    return jcall(obj, "toRealFieldBivariateFunction", RealFieldBivariateFunction, (Field,), arg0)
end

function value(obj::FieldBivariateFunction, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end


text:
parallel_execution_cache(::NotParallel) = NotParallel()

workers_mean!(target, data, par_type::NotParallel) =
    mean!(target, data)

workers_mean!(target, par_type::NotParallel) = target

workers_mean(data, par_type::NotParallel) = mean(data)

workers_sum!(target, par_type::NotParallel) = target
workers_sum!(target, data, par_type::NotParallel) = target .= data

num_workers(::NotParallel) = 1
worker_local_seed(seed, ::NotParallel) = seed
worker_block(iter, ::NotParallel) = iter
my_block(stuff, ::NotParallel) = stuff[1]

worker_allgatherv!(data, intervals, ::NotParallel) =
    data


text:
using LinearAlgebra

function duplication_matrix(n)
    mm = ones(n,n)
    tt = trues(n,n)
    uptt = UpperTriangular(tt)
    lott = LowerTriangular(tt)
    dd = Diagonal(tt)
    uptt = uptt .⊻ dd
    ind = 1:(n+1)*n/2
    mm[lott] = ind

    lott = lott .⊻ dd

    mm[uptt] = mm[lott]'
    mm = reshape(mm,n^2,1)
    m_fim = [mm[i] == ind[j] for i = 1:size(mm,1),j = 1:size(ind,1)]
    return m_fim
end

function vech(A)
    if !issymmetric(A)
        error("Not Symmetric")
    end
    n = size(A,1)
    tt = trues(n,n)
    n = Int(n*(n+1)/2)
    lott = LowerTriangular(tt)
    A_vec = A[lott]
    return reshape(A_vec,n,1)
end

function vecc(A)
    n = size(A,1)*size(A,2)
    return reshape(A',n,1)
end

function commutation_matrix(m,n)
    mat = zeros(m*n,m*n)
    seq = 1:(m*n)
    ind = mod.(seq .-1,m)*n + floor.((seq .-1)/m) .+1
    ind = Int.(ind)
    for k = 1:(m*n)
        mat[k,ind[k]] = 1
    end
    return mat
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC