text:
################################################################################
#
#  Map/NfRel.jl : Types for maps with domains of type NfRel
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

mutable struct NfRelToFqMor{T} <: Map{NfRel{T}, FqFiniteField, HeckeMap, NfRelToFqMor}
  header::MapHeader{NfRel{T}, FqFiniteField}

  function NfRelToFqMor{T}() where {T}
    z = new{T}()
    z.header = MapHeader{NfRel{T}, FqFiniteField}()
    return z
  end
end

function _automorphisms(L::NfRel{T}) where T
  if degree(L) == 1
    return morphism_type(L)[id_hom(L)]
  end
  f = L.pol
  Lt, t = PolynomialRing(L, "t", cached = false)
  f1 = change_base_ring(L, f, parent = Lt)
  divpol = Lt([ -gen(L), L(1) ])
  f1 = divexact(f1, divpol)
  lr = roots(f1)
  auts = Vector{morphism_type(L, L)}(undef, length(lr) + 1)
  for i = 1:length(lr)
    auts[i] = hom(L, L, lr[i], check = false)
  end
  auts[end] = id_hom(L)
  return auts
end

function automorphisms(L::T; copy::Bool = true) where {T <: NfRel}
  try
    Aut = _get_automorphisms_nf_rel(L)::Vector{morphism_type(T, T)}
    if copy
      v = Vector{morphism_type(T, T)}(undef, length(Aut))
      for i = 1:length(v)
        v[i] = Aut[i]
      end
      return v::Vector{morphism_type(T, T)}
    else
      return Aut::Vector{morphism_type(T, T)}
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  auts = _automorphisms(L)
  _set_automorphisms_nf_rel(L, auts)
  if copy
    v = Vector{morphism_type(T, T)}(undef, length(auts))
    for i = 1:length(v)
      v[i] = auts[i]
    end
    return v::Vector{morphism_type(T, T)}
  else
    return auts::Vector{morphism_type(T, T)}
  end
end

# Embedding of a relative number field into an algebra over the base field.
# S is the type of the field, T the type of the algebra and Mat the dense matrix
# type of the base ring of either
mutable struct NfRelToAbsAlgAssMor{S, T, Mat} <: Map{S, T, HeckeMap, NfRelToAbsAlgAssMor}
  header::MapHeader{S, T}
  mat::Mat
  t::Mat

  function NfRelToAbsAlgAssMor{S, T, Mat}(K::S, A::T, M::Mat) where { S <: NfRel, T <: AbsAlgAss, Mat <: MatElem }
    @assert base_ring(A) == base_field(K)
    z = new{S, T, Mat}()
    z.mat = M
    z.t = zero_matrix(base_field(K), 1, degree(K))

    function _image(x::NfRelElem)
      for i = 1:degree(K)
        z.t[1, i] = coeff(x, i - 1)
      end
      s = z.t*z.mat
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{S, T}(K, A, _image)
    return z
  end
end

function NfRelToAbsAlgAssMor(K::S, A::T, M::Mat) where { S <: NfRel, T <: AbsAlgAss, Mat <: MatElem }
  return NfRelToAbsAlgAssMor{S, T, Mat}(K, A, M)
end

function haspreimage(m::NfRelToAbsAlgAssMor, a::AbsAlgAssElem)
  A = parent(a)
  t = matrix(base_ring(A), 1, dim(A), coeffs(a))
  b, p = can_solve_with_solution(m.mat, t, side = :left)
  if b
    return true, domain(m)([ p[1, i] for i = 1:nrows(m.mat) ])
  else
    return false, zero(domain(m))
  end
end


text:
using PrimitiveOneHot
using Test

const tests = [
    "primitive",
    "array",
    "op",
]

function should_test_cuda()
    e = get(ENV, "JL_PKG_TEST_CUDA", false)
    e isa Bool && return e
    if e isa String
        x = tryparse(Bool, e)
        return isnothing(x) ? false : x
    else
        return false
    end
end

if should_test_cuda()
    push!(tests, "gpu")
end

@testset "PrimitiveOneHot.jl" begin
    for t in tests
        fp = joinpath(dirname(@__FILE__), "$t.jl")
        @info "Test $(uppercase(t))"
        include(fp)
    end
end

ext: jl


text:
using SparseArrays, Test
using SeDuMi

# Example page 4 of SeDuMi_Guide_105R5.pdf
@testset "Linear Programming example" begin
    c = [1.0, -1.0, 0.0, 0.0]
    A = [10.0 -7.0 -1.0 0.0
          1.0  0.5  0.0 1.0]
    b = [5.0, 3.0]
    sol, dual, status = sedumi(A, b, c, fid=0)
    @test sol ≈ sparse([1, 2], [1, 1], [47/24, 25//12], 4, 1)
    @test dual ≈ [1/8, -1/4]
    @test iszero(status["pinf"])
    @test iszero(status["dinf"])
    @test iszero(status["numerr"])
    @test status["feasratio"] == 1
    @test status["iter"] == 4
end


text:
"""
    efficientnetb0(; usedepthwise = false, nl = relu)

"""
function efficientnetb0(; usedepthwise = false, nl = relu)
    # stage 1
    stem = conv(3, 3, 32; nl = nl, stride = 2)

    body = Chain(

        # stage 2
        mbconv(3, 32, 32 * 1, 16, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 3
        mbconv(3, 16, 16 * 6, 24, nl = nl, stride = 2, usedepthwise = usedepthwise),
        mbconv(3, 24, 16 * 6, 24, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 4
        mbconv(5, 24, 24 * 6, 40, nl = nl, stride = 2, usedepthwise = usedepthwise),
        mbconv(5, 40, 40 * 6, 40, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 5
        mbconv(3, 40, 40 * 6, 80, nl = nl, stride = 2, usedepthwise = usedepthwise),
        mbconv(3, 80, 80 * 6, 80, nl = nl, stride = 1, usedepthwise = usedepthwise),
        mbconv(3, 80, 80 * 6, 80, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 6
        mbconv(5, 80, 80 * 6, 112, nl = nl, stride = 1, usedepthwise = usedepthwise),
        mbconv(5, 112, 112 * 6, 112, nl = nl, stride = 1, usedepthwise = usedepthwise),
        mbconv(5, 112, 112 * 6, 112, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 7
        mbconv(5, 112, 112 * 6, 192, nl = nl, stride = 2, usedepthwise = usedepthwise),
        mbconv(5, 192, 192 * 6, 192, nl = nl, stride = 1, usedepthwise = usedepthwise),
        mbconv(5, 192, 192 * 6, 192, nl = nl, stride = 1, usedepthwise = usedepthwise),
        mbconv(5, 192, 192 * 6, 192, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 8
        mbconv(3, 192, 192 * 6, 320, nl = nl, stride = 1, usedepthwise = usedepthwise),

        # stage 9
        conv(1, 320, 1280, nl = nl),
    )

    return Chain(stem, body)
end


text:
Pkg.clone("https://github.com/afniedermayer/InferenceUtilities.jl") # until registered
Pkg.clone(pwd())
Pkg.build("ByteParsers")
Pkg.test("ByteParsers"; coverage=true)


text:
#The key to the non-parametric method describedhere
#is that the sum of squared distances between pointsand
#their centroid is equal to (and can be calculateddirectly from)
#the sum of squared interpoint distancesdivided by the number of points
# https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1442-9993.2001.01070.pp.x

SS(D,N) = sum(D .^2)/N

function F(D,W,N,n,a)
    SST = SS(D,N)
    SSW = SS(W,n)
    SSA = SST - SSW
    return (SSA/(a-1))/(SSW/(N-a))
end

function F(SST ::Real,W,N,n,a)

    SSW = SS(W,n)
    SSA = SST - SSW
    return (SSA/(a-1))/(SSW/(N-a))
end

function getinds(group,N)
    W =Vector{CartesianIndex}(undef,0)
    notW = Vector{CartesianIndex}(undef,0)
    for i in 1:(N-1)
        for j in i+1:N
            if i !==1j
                if group[i] == group[j]
                    push!(W,(CartesianIndex(i,j)))
                else push!(notW,(CartesianIndex(i,j)))
                end
            end
        end
    end
    return W, notW
end

function permutest(D,group, n_perm = 1000)
    N =  length(group)
    a = length(unique(group))
    n = N/a
    Winds ,notWinds =getinds(group,N)

    @inbounds W = D[Winds]
    @inbounds notW= D[notWinds]

    Dvec = vcat(W,notW)

    SST = SS(Dvec,N)

    f = F(SST,W,N,n,a)
    Fs = Vector{Float64}(undef,n_perm)
    inds = 1:length(W)

    @inbounds for i in 1:n_perm
        shuffle!(Dvec)
        Fs[i] = F(SST,view(Dvec,inds),N,n,a)
    end
    P = sum(Fs .>= f)/n_perm

    g = unique(group)
    if a>2
        ppairs = NamedArray(zeros(a,a),( g,g ), ("group","group"))
        tpairs = NamedArray(zeros(a,a),( g,g ), ("group","group"))
        Threads.@threads for i in 1:a-1
            for j in i+1:a
                boolmask = (group .== g[i]) .| (group .== g[j])
                d = D[boolmask,boolmask]
                fstat,pstat =permutest(d,group[boolmask], n_perm)
                ppairs[j,i] = pstat
                tpairs[j,i] = sqrt(fstat)
            end
        end
    return (F = f,P = P, fpairs = tpairs,ppairs = ppairs)
    else
        return (F = f,P = P)
    end

end

X = rand(10,5)
D = pairwise(Euclidean(),X,X,dims = 1)
group = rand(1:2,10)
permutest(D,group)

X = vcat(rand(5,5),rand(5,5) .+3)
D = pairwise(Euclidean(),X,X,dims = 1)
group = [1,1,2,2,2,3,3,3,3,3]
permutest(D,group)

X = vcat(rand(800,500))
D = pairwise(Euclidean(),X,X,dims = 1)
group = vcat(fill(collect(1:10),80)...)
permutest(D,group)

max_stars_repo_name: EvoArt/Hydra

max_issues_repo_name: EvoArt/Hydra

max_forks_repo_name: EvoArt/Hydra


text:
module LogDensityTestSuite

export
    # generic
    samples,
    # primitives
    StandardMultivariateNormal,
    # transformations
    linear, shift, elongate,
    # mixtures
    mix, directional_weight,
    # diagnostics
    quantile_boundaries, bin_counts, two_sided_pvalues, print_ascii_plot

using ArgCheck: @argcheck
using DocStringExtensions: FIELDS, FUNCTIONNAME, SIGNATURES, TYPEDEF
using LinearAlgebra: AbstractTriangular, checksquare, diag, Diagonal, I, dot, logabsdet, lu,
    norm
using LogDensityProblems: LogDensityOrder
import LogDensityProblems: capabilities, dimension, logdensity, logdensity_and_gradient
using MCMCDiagnostics: ess_factor_estimate
using Parameters: @unpack
using Printf: @sprintf
using Sobol: next!, SobolSeq
using Statistics: quantile
using StatsFuns: norminvcdf, normcdf, logaddexp, logistic

include("generic.jl")
include("primitives.jl")
include("transformations.jl")
include("mixtures.jl")
include("diagnostics.jl")

end # module


text:
using ValidatedNumerics
#
# HornerInt initially stores a list of coefficients of a polynomial, a evaluation point,
# and optionally a corresponding list of absolute errors of the coefficients.
# On demand it modifies the values in order to obtain the values of the selected derivatives
# at the given point and the corresponding error value.
#

if ! isdefined(:HornerInt)
type HornerInt{S<:AbstractFloat, T<:AbstractFloat}
    x::Interval{S}						# evaluation spot
    coeff::AbstractVector{Interval{T}}	# polynomial coefficients
    valid::Int						# evaluation status
end
end

# Create Horner table for evaluation at point x.
# Provide polynomial coefficients and optionally absolute errors for x and the coefficients.
function HornerInt{S<:Number,T<:Number,U<:AbstractFloat,V<:AbstractFloat}(xdown::S, xup::T, coeffdown::AbstractVector{U}, coeffup::AbstractVector{V})

    if coeffup == []
        coeffup = copy(coeffdown)
    end
    xdown, xup = min(xdown, xup), max(xdown, xup)
    coeffdown, coeffup = min(coeffdown, coeffup), max(coeffdown, coeffup)
    if length(coeffdown) != length(coeffup) throw(ArgumentError("coeffup length wrong")) end
    R = promote_type(S,T,U,V)
	HornerInt(R(xdown), R(xup), map(R,coeffdown), map(R,coeffup), -1)
end

# Evaluate Horner table up to and including m deviations.
# return the m'th deviation of the polynomial and the error estimation
# If m >= degree, all deviations ot the polynomial have been calculated at x.
# Further calls just return the cached results.
function call{S}(h::HornerInt{S}, m::Int)
    if m < 0 throw(ArgumentError("index $m is not >= 0")) end
	n = length(h.coeffdown) - 1
    if m > h.valid
        if h.xdown != 0 || h.xup != 0
            r = rounding(S)
            setronding(S, RoundDown)
            b, bu = h.coeffdown, h.coeffup
            x, xu = h.xdown, h.xup
            for k = max(h.valid+1,0):min(m, n-1)
                for j = 1:n-k
                    b[j+1], bu[j+1] = intmin((b[j+1], bu[j+1), (b[j], bu[j]), (x, xu))
                end
            end
            setRounding(S, r)
        end
        h.valid = m
    end
    m <= n ? (h.coeffdown[n+1-m], h.coeffup[n+1-m]) : (zero(S), zero(S))
end

#
# evaluate polynomial f and all derivatives at given approximate root-values z
#
function horner_analysis{S<:Number,T<:Number}(z::AbstractVector{S}, f::AbstractVector{T})
    thresh = 1000.0
    m = length(z)
    n = length(z)-1
    if m > length(f)-1
        throw(ArgumentError("polynomial of degree $n with $m zeros"))
    end
	h = Array{Horner}(m)
    for i = 1:m
        h[i] = Horner(z[i], f, eps(abs(z[i])) * 100, map(eps, f) * 100)
    end
	for i = 1:m
        println("root $i: $(z[i])")
        mult = 0
        search = true
		for k = 0:n
            hv, he = h[i](k)
            deviation = abs(hv / (he + eps(zero(real(S)))))
            @printf("%d %d %6.2g %8.8g %8.8g\n", i, k, deviation, abs(hv), he)
            if search && deviation < thresh
                mult = k + 1
            else
                search = false
            end
        end
        println("root $i: $(z[i]) estimated multiplicity: $mult")
    end
end


text:
module graphops

using PyCall
include("utils.jl")

export
    # graph_transformer
    optimize,
    simplify,
    rewrite_einsum_expr,
    # graph_optimizer
    fuse_einsums,
    # graph_inv_optimizer
    optimize_inverse,
    prune_inv_node,
    # graph_generator
    generate_optimal_tree,
    split_einsum,
    generate_optimal_tree_w_constraint,
    # graph_als_optimizer
    generate_sequential_optimal_tree

const gops = PyNULL()

function __init__()
    copy!(gops, pyimport_conda("autohoot.graph_ops", "autohoot"))
end

@func_loader gops.graph_transformer optimize simplify rewrite_einsum_expr
@func_loader gops.graph_optimizer fuse_einsums
@func_loader gops.graph_inv_optimizer optimize_inverse prune_inv_node
@func_loader gops.graph_generator generate_optimal_tree split_einsum generate_optimal_tree_w_constraint
@func_loader gops.graph_als_optimizer generate_sequential_optimal_tree

end


text:
using IRTools

## Get set of freeVars from an IR block
const VarSet = Set{IRTools.Variable}
function union!(s :: Set{T}, ts) where {T}
    for t in ts
        push!(s, t)
    end
end

function freevars(env :: VarSet, x :: Any)
    println("[FREEVARS UNK", typeof(x), ":", x, "]")
end

function freevars(env :: VarSet, x :: Number)
    VarSet()
end

function freevars(env :: VarSet, x :: IRTools.Variable)
    if x in env
        VarSet()
    else
        VarSet([x])
    end
end

function freevars(env :: VarSet, x :: GlobalRef)
    VarSet()
end

function freevars(env :: VarSet, x :: Expr)
    if x.head == :call
        out = VarSet()
        for a in x.args
            union!(out, freevars(env, a))
        end
        out
    else
        println("FREEVARS[EXPR<", x.head, ">:", x, "]")
    end
end

function freevars(env :: VarSet, b :: IRTools.Block)
    args = VarSet(IRTools.arguments(b))
    env = union(env, args)
    out = VarSet()
    for (x,st) in b
        if x != undef
            push!(env, x)
        end
        union!(out, freevars(env, st.expr))
    end
    out
end

freevars(x) = freevars(VarSet(), x)

###################
### Test freevars

using Test

function test_freevars_1(x,y)
    # Block 1: Args (%1,%2,%3), FV: ()
    a = sin(x)        # %4
    b = a+y           # %5

    if x > 3          # %6
        # Block 2: Args (), FV %5
        c = cos(b)    # %7
    else
        # Block 3: Args (), FV %3, %4
        ty = tan(y)   # %8
        c = a + ty    # %9
    end

    # Block 4: Args %10, FV %4
    d = c + a         # %11
    e = atan(d)       # %12
end

ir = @code_ir test_freevars_1(4,5)

set(l) = VarSet(map(IRTools.var, l))
fv(i) = freevars(IRTools.block(ir,i))
@test fv(1) == set([])
@test fv(2) == set([5])
@test fv(3) == set([3,4])
@test fv(4) == set([4])

max_stars_repo_name: tomjaguarpaw/knossos-ksc

max_issues_repo_name: tomjaguarpaw/knossos-ksc

max_forks_repo_name: tomjaguarpaw/knossos-ksc


text:
module SamplingWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Class = @jimport java.lang.Class
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldOrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitFixedStepHandler
const FieldOrekitFixedStepHandlerMultiplexer = @jimport org.orekit.propagation.sampling.FieldOrekitFixedStepHandlerMultiplexer
const FieldOrekitStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitStepHandler
const FieldOrekitStepHandlerMultiplexer = @jimport org.orekit.propagation.sampling.FieldOrekitStepHandlerMultiplexer
const FieldOrekitStepInterpolator = @jimport org.orekit.propagation.sampling.FieldOrekitStepInterpolator
const FieldOrekitStepNormalizer = @jimport org.orekit.propagation.sampling.FieldOrekitStepNormalizer
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const MultiSatStepHandler = @jimport org.orekit.propagation.sampling.MultiSatStepHandler
const Object = @jimport java.lang.Object
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitFixedStepHandlerMultiplexer = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandlerMultiplexer
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const OrekitStepHandlerMultiplexer = @jimport org.orekit.propagation.sampling.OrekitStepHandlerMultiplexer
const OrekitStepInterpolator = @jimport org.orekit.propagation.sampling.OrekitStepInterpolator
const OrekitStepNormalizer = @jimport org.orekit.propagation.sampling.OrekitStepNormalizer
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState

export FieldOrekitFixedStepHandler
export FieldOrekitFixedStepHandlerMultiplexer
export FieldOrekitStepHandler
export FieldOrekitStepHandlerMultiplexer
export FieldOrekitStepInterpolator
export FieldOrekitStepNormalizer
export MultiSatStepHandler
export OrekitFixedStepHandler
export OrekitFixedStepHandlerMultiplexer
export OrekitStepHandler
export OrekitStepHandlerMultiplexer
export OrekitStepInterpolator
export OrekitStepNormalizer
export add
export equals
export get_class
export get_current_state
export get_interpolated_state
export get_previous_state
export handle_step
export hash_code
export init
export is_current_state_interpolated
export is_forward
export is_previous_state_interpolated
export notify
export notify_all
export requires_dense_output
export restrict_step
export to_string
export wait

include("field_orekit_fixed_step_handler.jl")
include("field_orekit_fixed_step_handler_multiplexer.jl")
include("field_orekit_step_handler.jl")
include("field_orekit_step_handler_multiplexer.jl")
include("field_orekit_step_interpolator.jl")
include("field_orekit_step_normalizer.jl")
include("multi_sat_step_handler.jl")
include("orekit_fixed_step_handler.jl")
include("orekit_fixed_step_handler_multiplexer.jl")
include("orekit_step_handler.jl")
include("orekit_step_handler_multiplexer.jl")
include("orekit_step_interpolator.jl")
include("orekit_step_normalizer.jl")

end


text:
function InsertionSort!(arr::Vector{T})where T
    for i in 1:length(arr)-1
        temp=arr[i+1]
        j=i
        while j>0&&arr[j]>temp
            arr[j+1]=arr[j]
            j-=1
        end
        arr[j+1]=temp
    end
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
using AbstractPlotting

function main()
    s1 = slider(LinRange(0.01, 1, 100), raw = true, camera = campixel!, start = 0.3)
    s2 = slider(LinRange(-2pi, 2pi, 100), raw = true, camera = campixel!)
    data = lift(s2[end][:value]) do v
         map(LinRange(0, 2pi, 100)) do x
            4f0 .* Point2f0(sin(x) + (sin(x * v) .* 0.1), cos(x) + (cos(x * v) .* 0.1))
        end
    end
    p = scatter(data)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
__precompile__()

module IndexedArrays

include("delegate.jl")

import Base: copy, in, getindex, findfirst, length, size, isempty, start, done, next, empty!, push!, pop!

abstract AbstractIndexedArray{T} <: AbstractVector{T}

type IndexedArrayError <: Exception # FIXME: or should we just use ArgumentError here?
    msg::AbstractString
end

immutable IndexedArray{T} <: AbstractIndexedArray{T}
    items::Array{T}
    lookup::Dict{T,Int}

    IndexedArray() = new(T[], Dict{T,Int}())
    function IndexedArray{T}(items::Array{T})
        ia = new(items, Dict{T,Int}())
        sizehint!(ia.lookup, length(ia.items))
        for (i, item) in enumerate(ia.items)
            if item in ia
                throw(IndexedArrayError("cannot construct IndexedArray with duplicate items"))
            end
            ia.lookup[item] = i
        end
        @assert length(ia.items) == length(ia.lookup)
        return ia
    end
end

IndexedArray{T}(items::Array{T}) = IndexedArray{T}(items)

@delegate IndexedArray.items [ length, size, isempty, start, done, next ]

getindex(ia::IndexedArray, i::Int) = getindex(ia.items, i)
getindex(ia::IndexedArray, r::UnitRange{Int}) = getindex(ia.items, r)

function empty!(ia::IndexedArray)
    # NOTE: does not provide any exception safety guarantee
    empty!(ia.items)
    empty!(ia.lookup)
    return ia
end

in{T}(item::T, ia::IndexedArray{T}) = haskey(ia.lookup, item)

findfirst{T}(ia::IndexedArray{T}, item::T) = ia.lookup[item] # throws KeyError if not found

function findfirst!{T}(ia::IndexedArray{T}, item::T)
    rv = get!(ia.lookup, item) do
        # NOTE: does not provide any exception safety guarantee
        push!(ia.items, item)
        return length(ia.items)
    end
    @assert length(ia.items) == length(ia.lookup)
    return rv
end

function push!{T}(ia::IndexedArray{T}, item::T)
    if item in ia
        throw(IndexedArrayError("cannot add duplicate item to IndexedArray"))
    end
    # NOTE: does not provide any exception safety guarantee
    push!(ia.items, item)
    ia.lookup[item] = length(ia)
    @assert length(ia.items) == length(ia.lookup)
    return ia
end

function pop!(ia::IndexedArray)
    if isempty(ia.items)
        throw(ArgumentError("array must be non-empty"))
    end
    # NOTE: does not provide any exception safety guarantee
    delete!(ia.lookup, ia.items[end])
    rv = pop!(ia.items)
    @assert length(ia.items) == length(ia.lookup)
    return rv
end

copy{T}(ia::IndexedArray{T}) = IndexedArray{T}(copy(ia.items))

"`swap!(ia::IndexedArray, to::Int, from::Int)` interchange/swap the values on the indices `to` and `from` in the `IndexedArray`"
function swap!(ia::IndexedArray, to::Int, from::Int)
    if to == from
        checkbounds(ia,to)
        return ia
    end
    previous_id  = ia[to]
    future_id    = ia[from]

    ia.items[to]   = future_id
    ia.items[from] = previous_id

    ia.lookup[previous_id] = from
    ia.lookup[future_id]   = to

    return ia
end

export AbstractIndexedArray, IndexedArray, IndexedArrayError, findfirst!, swap!

end # module


text:
# using Optim
# import DifferentialEquations

# ODE system for affine dynamics
function check_sufficient_ode_function(dyn1::AffDyn, dyn2::AffDyn)
    A1 = dyn1.A
    b1 = dyn1.b
    A2 = dyn2.A
    b2 = dyn2.b
    n = length(b1)
    Z = zeros(n, n)
    C = [A1 Z Z; Z A2 Z; A1 -A2 Z]
    d = [b1; b2; b1-b2]
    h(u, p, t) = C * u + d
    return h, n
end

# ODE system for mixed affine and linear dynamics
function check_sufficient_ode_function(dyn1::AffDyn, B::AbstractMatrix)
    n = size(B, 1)
    dyn2 = AffDyn(B, zeros(n))
    return check_sufficient_ode_function(dyn1, dyn2)
end
function check_sufficient_ode_function(A::AbstractMatrix, dyn2::AffDyn)
    n = size(A, 1)
    dyn1 = AffDyn(A, zeros(n))
    return check_sufficient_ode_function(dyn1, dyn2)
end

# ODE system for linear dynamics
function check_sufficient_ode_function(A::AbstractMatrix, B::AbstractMatrix)
    n = size(A, 1)
    Z = zeros(n, n)
    C = [A Z Z; Z B Z; A -B Z]
    h(u, p, t) = C * u
    return h, n
end

function check_sufficient(dyn1, dyn2, x0::AbstractVector, y0::AbstractVector,
                          time; return_solution::Bool=false, log::Bool=false)
    if log
        println("running sufficient check")
    end

    # set up problem
    h, n = check_sufficient_ode_function(dyn1, dyn2)

    z0 = [x0; y0; (x0-y0)]

    problem = DifferentialEquations.ODEProblem(h, z0, time)

    ### solve problem

    sol = DifferentialEquations.solve(problem)

    ### optimize distance

    opt_all = -Inf
    for i in (2*n+1):(3*n)
        F⁺ = (t) -> sol(t, idxs=i)
        F⁻ = (t) -> -sol(t, idxs=i)

        opt⁺ = optimize(F⁺, time[1], time[2])
        opt⁻ = optimize(F⁻, time[1], time[2])
        opt = max(abs(opt⁺.minimum), abs(opt⁻.minimum))

        opt_all = max(opt_all, opt)
    end

    ### plot solution

    if return_solution
        return opt_all, sol
    else
        return opt_all
    end
end

function check_sufficient_points(dyn1, dyn2, x0::AbstractVector, points, time,
                                 ε::Number; log::Bool=false)
    if log
        println("running sufficient check for $(length(points)) points")
    end

    points_sat = Vector{Vector{Float64}}()
    points_unsat = Vector{Vector{Float64}}()
    for p in points
        max_diff = check_sufficient(dyn1, dyn2, x0, p, time)
        if max_diff <= ε
            push!(points_sat, p)
        else
            push!(points_unsat, p)
        end
    end
    return points_sat, points_unsat
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
module jsos

include("sos.jl")
include("utils.jl")

end

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
using PSSFSS
using Documenter

#=
using Literate
cd("literate") do
  include("literate/compile.jl")
end
=#

makedocs(;
    clean=false,
    modules=[PSSFSS],
    authors="Peter Simon <psimon0420@gmail.com> and contributors",
    repo="https://github.com/simonp0420/PSSFSS.jl/blob/{commit}{path}#L{line}",
    sitename="PSSFSS.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://simonp0420.github.io/PSSFSS.jl/stable",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Examples" => "examples.md",
        "Function Reference" => "reference.md",
        "Index" => "function_index.md"
    ],
)

deploydocs(;
    repo="github.com/simonp0420/PSSFSS.jl",
    devbranch = "main"
)


text:
# http://web.mit.edu/puzzle/www/2011/puzzles/civilization/meta/wall_street.html
@testset "wallstreet meta" begin
    f = best_feature(["autumn", "badminton", "trafficpylon", "american", "ingrid", "mercury", "corncake", "gooier", "triskelion", "wandering"], 1)
    @test f.description == "has 1 letters repeated at least 2 times each"
end


text:
using Documenter, BeurlingLasso

makedocs(
    modules = [BeurlingLasso],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Romain Petit",
    sitename = "BeurlingLasso",
    pages = [
        "Introduction" => "index.md",
        "Examples" => "examples.md"
    ]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/rpetit/BeurlingLasso.jl.git",
)


text:
using Caesar, DelimitedFiles, Profile
using Gadfly, Cairo, Fontconfig

fFloor = 250;
fCeil = 1750;
nFFT_czt = 1024;
fSampling = 37500.0
nPhones = 10;
soundSpeed = 1481;
azimuthDivs = 180;
azimuths = range(0,360,length=azimuthDivs)*pi/180;

#Single Test Array of length 10
# posFile = joinpath(ENV["HOME"],"data","sas","test_array_positions.csv");
# dataFile = joinpath(ENV["HOME"],"data","sas","test_array_waveforms.csv");
# arrayElemPos = readdlm(posFile,',',Float64,'\n')
# arrayElemPos = arrayElemPos[:,1:2]

#Load Experimental Data - 10 frames
winstart = 1600;
rawWaveData = zeros(8000,nPhones);
arrayElemPos = zeros(nPhones,2);
for ele in winstart:winstart+nPhones-1
    dataFile = joinpath(ENV["HOME"],"data", "sas", "08_20_sample","waveform$(ele).csv");
    posFile = joinpath(ENV["HOME"],"data", "sas", "08_20_sample","nav$(ele).csv");
    tempRead = readdlm(dataFile,',',Float64,'\n') #first element only
    rawWaveData[:,ele-winstart+1] = adjoint(tempRead[1,:]);
    tempRead = readdlm(posFile,',',Float64,'\n');
    arrayElemPos[ele-winstart+1,:] = tempRead;
end

chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");

FFTfreqs = collect(LinRange(fFloor,fCeil,nFFT_czt))

cfg = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones,azimuths,soundSpeed,FFTfreqs)
myCBF = zeros(Complex{Float64}, getCBFFilter2Dsize(cfg));
lm = zeros(2,1);
dataHolder = zeros(Complex{Float64},nPhones,nFFT_czt)
delaysHolder = zeros(Complex{Float64},nFFT_czt);

@time dataHolder[1,:] = FFTfreqs

@time constructCBFFilter2D!(cfg, arrayElemPos, myCBF, lm, dataHolder)

# MF and CZT
w = exp(-2im*pi*(fCeil-fFloor)/(nFFT_czt*fSampling))
a = exp(2im*pi*fFloor/fSampling)

chirpIn = readdlm(chirpFile,',',Float64,'\n')

#Matched Filter on Data In
nFFT_full = nextpow(2,size(rawWaveData,1))  # MF
mfData = zeros(Complex{Float64}, nFFT_full, nPhones)
mf = prepMF(chirpIn,nFFT_full,nPhones) # MF
mf(rawWaveData,mfData) # MF

cztData = zeros(Complex{Float64}, nFFT_czt,nPhones)
filterCZT = prepCZTFilter(nFFT_full,nPhones,w,nFFT_czt,a)
filterCZT(mfData,cztData)

# CBF step
dataOut = zeros(length(azimuths));
temp1 = zeros(Complex{Float64},nFFT_czt);
temp2 = zeros(Complex{Float64},size(cztData));
@time CBF2D_DelaySum!(cfg, cztData, dataOut,temp1,temp2,myCBF)

pl = Gadfly.plot(
 Gadfly.layer(x=azimuths, y=real(dataOut), Geom.line),
 Gadfly.Coord.cartesian(xmin=0, xmax=2*pi),
 Guide.ylabel(nothing),Guide.xlabel("Azimuths (rad)")
 ) ; pl |> PDF("/tmp/test_CBF.pdf")


text:
function triSolve!(x,tri,b)
# Thomas algorithm, b and x used as intermediate storage
# b is destroyed after the elemination
n = size(tri,2)
x[1] = tri[2,1]
@inbounds @fastmath for i=1:n-1
  x[i+1] = tri[2,i+1]-tri[3,i]/x[i]*tri[1,i+1]
  b[i+1] = b[i+1]-tri[3,i]/x[i]*b[i]
end
x[n] = b[n]/x[n]
@inbounds @fastmath for i=n-1:-1:1
  x[i]=(b[i] - tri[1,i+1]*x[i+1])/x[i]
end
end

function mulUL!(tri,biU,biL)
n = size(biL,2)
tri[2,1] = tri[2,1] - biU[2,1]*biL[1,1] - biU[1,2]*biL[2,1]
tri[3,1] = tri[3,1] - biU[2,2]*biL[2,1]
@inbounds @fastmath for i=2:n-1
  tri[1,i] = tri[1,i] - biU[1,i]*biL[1,i]
  tri[2,i] = tri[2,i] - biU[2,i]*biL[1,i] - biU[1,i+1]*biL[2,i]
  tri[3,i] = tri[3,i] - biU[2,i+1]*biL[2,i]
end
i = n
tri[1,i] = tri[1,i] - biU[1,i]*biL[1,i]
tri[2,i] = tri[2,i] - biU[2,i]*biL[1,i]
end

function mulbiUv!(u,biU,v)
n = size(biU,2)
@inbounds @fastmath for i=1:n-1
  u[i] = u[i] + biU[2,i]*v[i] + biU[1,i+1]*v[i+1]
end
u[n] = u[n] + biU[2,n]*v[n]
end

function mulbiLv!(u,biL,v)
n = size(biL,2)
u[1] = u[1] + biL[1,1]*v[1]
@inbounds @fastmath for i=2:n
  u[i] = u[i] + biL[2,i-1]*v[i-1] + biL[1,i]*v[i]
end
end

function SchurSolve!(k,v,J,fac,Param)
n1=size(v,1);
n2=size(v,2);
n=n1*n2;
tri=J.tri
JWRho=J.JWRho
JRhoW=J.JRhoW
JWTh=J.JWTh
JThW=J.JThW
JWW=J.JWW
rRho=J.rRho
rTh=J.rTh
rw=J.rw
sw=J.sw

@views rRho=reshape(permute(v[:,:,1],[2,1]),n,1);
@views rTh=reshape(permute(v[:,:,5],[2,1]),n,1);
@views rw=reshape(permute(v[:,:,4],[2,1]),n,1);
invfac=1/fac;
invfac2=invfac/fac;
if Param.Damping
# sw=(spdiags(repmat(invfac2,n,1),0,n,n)-invfac*JWW-JWRho*JRhoW-JWTh*JThW)\
#   (invfac*rw+JWRho*rRho+JWTh*rTh);
else
  @views tri[1,:] .= 0
  @views tri[2,:] .= invfac2
  @views tri[3,:] .= 0
  mulUL!(tri,JWRho,JRhoW)
  mulUL!(tri,JWTh,JThW)
  rw .= invfac .* rw
  mulbiUv!(rw,JWRho,rRho)
  mulbiUv!(rw,JWTh,rTh)
  triSolve!(sw,tri,rw)
end
mulbiLv!(rRho,JRhoW,sw)
mulbiLv!(rTh,JThW,sw)

@views k[:,:,1] .= fac .* permute(reshape(rRho,n2,n1),[2,1]);
@views k[:,:,2:3] .= fac*v[:,:,2:3];
@views k[:,:,4] .= permute(reshape(sw,n2,n1),[2,1]);
@views k[:,:,5] .= fac .* permute(reshape(rTh,n2,n1),[2,1]);
end


text:
using FiniteDifferences: rand_tangent, difference

function test_difference(ε::Real, x, dx)
    y = x + ε * dx
    dx_diff = difference(ε, y, x)
    # TODO: `@test isapprox(dx, dx_diff)` once `isapprox` is defined appropriately
    # see https://github.com/JuliaDiff/ChainRulesCore.jl/issues/184
    @test typeof(dx) == typeof(dx_diff)
end

@testset "difference" begin

    @testset "Primal: $(typeof(x))" for (ε, x) in [

        # Test things that don't have tangents.
        (randn(), :a),
        (randn(), 'a'),
        (randn(), "a"),
        (randn(), 0),

        # Test Numbers.
        (randn(Float32), randn(Float32)),
        (randn(Float64), randn(Float64)),
        (randn(Float32), randn(ComplexF32)),
        (randn(Float64), randn(ComplexF64)),

        # Test StridedArrays.
        (randn(), randn(5)),
        (randn(Float32), randn(ComplexF32, 5, 2)),
        (randn(), [randn(1) for _ in 1:3]),
        (randn(), [randn(5, 4), "a"]),

        # Tuples.
        (randn(), (randn(5, 2), )),
        (randn(), (randn(), 4)),
        (randn(), (4, 3, 2)),

        # NamedTuples.
        (randn(), (a=randn(5, 2),)),
        (randn(), (a=randn(), b=4)),
        (randn(), (a=4, b=3, c=2)),

        # Arbitrary structs.
        (randn(), sin),
        (randn(), cos),
        (randn(), Foo(5.0, 4, randn(5, 2))),
        (randn(), Foo(randn(), 1, Foo(randn(), 1, 1))),

        # LinearAlgebra types (also just structs).
        (randn(), UpperTriangular(randn(2, 2))),
        (randn(), Diagonal(randn(4))),
        (randn(), SVector{2, Float64}(1.0, 2.0)),
        (randn(), SMatrix{2, 2, ComplexF64}(1.0, 2.0, 3.0, 4.0)),
        (randn(), Symmetric(randn(2, 2))),
        (randn(), Hermitian(randn(ComplexF64, 1, 1))),
        (randn(), Adjoint(randn(ComplexF64, 3, 3))),
        (randn(), Transpose(randn(3))),
    ]
        test_difference(ε, x, rand_tangent(x))
    end

    # Ensure struct fallback errors for non-struct types.
    @test_throws ArgumentError invoke(difference, Tuple{Float64, T, T} where {T}, 5.0, 5.0, 5.0)
end


text:
using ModelingToolkit, Test
@parameters t a
@variables x(t) y(t)
D = Differential(t)
eqs = [D(x) ~ a*x - x*y,
       D(y) ~ -3y + x*y]
@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
    """
    real[] diffeqf(real t,real[] internal_var___u,real[] internal_var___p,real[] x_r,int[] x_i) {
      real internal_var___du[2];
      internal_var___du[1] = internal_var___p[1] * internal_var___u[1] + -1 * internal_var___u[1] * internal_var___u[2];
      internal_var___du[2] = internal_var___u[1] * internal_var___u[2] + -3 * internal_var___u[2];
      return internal_var___du;
    }
    """

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
  """
  void diffeqf(double* internal_var___du, double* internal_var___u, double* internal_var___p, double t) {
    internal_var___du[0] = internal_var___p[0] * internal_var___u[0] + -1 * internal_var___u[0] * internal_var___u[1];
    internal_var___du[1] = internal_var___u[0] * internal_var___u[1] + -3 * internal_var___u[1];
  }
  """

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
  """
  diffeqf = @(t,internal_var___u) [internal_var___p(1) * internal_var___u(1) + -1 * internal_var___u(1) * internal_var___u(2); internal_var___u(1) * internal_var___u(2) + -3 * internal_var___u(2)];"""

sys = ODESystem(eqs,t,[x,y],[a])

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
    ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                   lhsname=:internal_var___du,
                                   rhsnames=[:internal_var___u,:internal_var___p,:t])

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
      ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.StanTarget())

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
      ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.MATLABTarget())

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t])

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.StanTarget())

@test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.MATLABTarget())


text:
using NearestNeighborModels
using StableRNGs
using Test
using MLJBase
using OffsetArrays

import NearestNeighborModels: check_onebased_indexing, combine_weights,
    err_if_given_invalid_K, Fill, get_weights, KNNResult, NeighDistsMatrix,
    NeighIndsMatrix, RankOneMatrix, _replace!, scale, sort_idxs, _sum

include("testutils.jl")

@testset "NearestNeighborModels.jl" begin
    include("models.jl")
    include("utils.jl")
    include("kernels.jl")
end


text:
import Base: getindex

##### getindex #####
struct GetindexPosthookCtx end

@wrapped function getindex_posthook(f, x::AbstractArray)
    if hasmetadata(x, GetindexPosthookCtx)
        g = getmetadata(x, GetindexPosthookCtx)
        setmetadata(x,
                    GetindexPosthookCtx,
                    (res, args...) -> f(g(res, args...), args...))
    else
        setmetadata(x, GetindexPosthookCtx, f)
    end
end
function Base.getindex(x::SymArray, idx::CartesianIndex)
    return x[Tuple(idx)...]
end

function Base.getindex(x::SymArray, idx...)
    if all(i->symtype(i) <: Integer, idx)
        res = Term{eltype(symtype(x))}(getindex, [x, idx...])
    else
        input_idx = []
        output_idx = []
        ranges = Dict{Sym, AbstractRange}()
        subscripts = makesubscripts(length(idx))
        for (j, i) in enumerate(idx)
            if symtype(i) <: Integer
                push!(input_idx, i)
            elseif i isa Colon
                push!(output_idx, subscripts[j])
                push!(input_idx, subscripts[j])
            elseif i isa AbstractVector
                isym = subscripts[j]
                push!(output_idx, isym)
                push!(input_idx, isym)
                ranges[isym] = i
            else
                error("Don't know how to index by $i")
            end
        end

        term = Term{Any}(getindex, [x, idx...])
        T = eltype(symtype(x))
        N = ndims(x) - count(i->symtype(i) <: Integer, idx)
        res = ArrayOp(atype(symtype(x)){T,N},
                      (output_idx...,),
                      x[input_idx...],
                      +,
                      term,
                      ranges)
    end

    if hasmetadata(x, GetindexPosthookCtx)
        f = getmetadata(x, GetindexPosthookCtx, identity)
        f(res, x, idx...)
    else
        res
    end
end

# Wrapped array should wrap the elements too
function Base.getindex(x::Arr, idx...)
    wrap(unwrap(x)[idx...])
end
function Base.getindex(x::Arr, idx::Symbolic{<:Integer}...)
    wrap(unwrap(x)[idx...])
end

function propagate_ndims(::typeof(getindex), x, idx...)
    ndims(x) - count(x->symtype(x) <: Integer, idx)
end

function propagate_shape(::typeof(getindex), x, idx...)
    @oops axes = shape(x)

    idx1 = to_indices(CartesianIndices(axes), axes, idx)
    ([1:length(x) for x in idx1 if !(symtype(x) <: Number)]...,)
end

propagate_eltype(::typeof(getindex), x, idx...) = geteltype(x)

function SymbolicUtils.promote_symtype(::typeof(getindex), X, ii...)
    @assert all(i-> i <:Integer, ii) "user arrterm to create arr term."

    eltype(X)
end

#### Broadcast ####
#

using Base.Broadcast

Base.broadcastable(s::SymArray) = s
struct SymBroadcast <: Broadcast.BroadcastStyle end
Broadcast.BroadcastStyle(::Type{<:SymArray}) = SymBroadcast()
Broadcast.result_style(::SymBroadcast) = SymBroadcast()
Broadcast.BroadcastStyle(::SymBroadcast, ::Broadcast.BroadcastStyle) = SymBroadcast()

isonedim(x, i) = shape(x) == Unknown() ? false : isone(size(x, i))

function Broadcast.materialize(bc::Broadcast.Broadcasted{SymBroadcast})
    # Do the thing here
    ndim = mapfoldl(ndims, max, bc.args, init=0)
    subscripts = makesubscripts(ndim)

    expr_args′ = map(bc.args) do x
        if ndims(x) != 0
            subs = map(i-> isonedim(x, i) ?
                       1 : subscripts[i], 1:ndims(x))
            x[subs...]
        else
            x
        end
    end

    expr = term(bc.f, expr_args′...) # Imagine x .=> y -- if you don't have a term
                                     # then you get pairs, and index matcher cannot
                                     # recurse into pairs
    Atype = propagate_atype(broadcast, bc.f, bc.args...)
    ArrayOp(Atype{symtype(expr), ndim},
            (subscripts...,),
            expr,
            +,
            Term{Any}(broadcast, [bc.f, bc.args...]))
end

# On wrapper:
struct SymWrapBroadcast <: Broadcast.BroadcastStyle end

Base.broadcastable(s::Arr) = s

Broadcast.BroadcastStyle(::Type{<:Arr}) = SymWrapBroadcast()

Broadcast.result_style(::SymWrapBroadcast) = SymWrapBroadcast()

Broadcast.BroadcastStyle(::SymWrapBroadcast,
                         ::Broadcast.BroadcastStyle) = SymWrapBroadcast()

function Broadcast.materialize(bc::Broadcast.Broadcasted{SymWrapBroadcast})
    args = map(bc.args) do arg
        if arg isa Broadcast.Broadcasted
            return Broadcast.materialize(arg)
        else
            return arg
        end
    end
    wrap(broadcast(bc.f, map(unwrap, args)...))
end

#################### TRANSPOSE ################
#
@wrapped function Base.adjoint(A::AbstractMatrix)
    @syms i::Int j::Int
    @arrayop A' (i, j) A[j, i]
end

@wrapped function Base.adjoint(b::AbstractVector)
    @syms i::Int
    @arrayop b' (1, i) b[i]
end

import Base: *, \

using LinearAlgebra

isdot(A, b) = isadjointvec(A) && ndims(b) == 1

isadjointvec(A::Adjoint) = ndims(parent(A)) == 1
isadjointvec(A::Transpose) = ndims(parent(A)) == 1

function isadjointvec(A::Term)
    (operation(A) === (adjoint) ||
     operation(A) == (transpose)) && ndims(arguments(A)[1]) == 1
end

isadjointvec(A::ArrayOp) = isadjointvec(A.term)

isadjointvec(A) = false

# TODO: add more such methods
function getindex(A::AbstractArray, i::Symbolic{<:Integer}, ii::Symbolic{<:Integer}...)
    Term{eltype(A)}(getindex, [A, i, ii...])
end

function getindex(A::AbstractArray, i::Int, j::Symbolic{<:Integer})
    Term{eltype(A)}(getindex, [A, i, j])
end

function getindex(A::AbstractArray, j::Symbolic{<:Integer}, i::Int)
    Term{eltype(A)}(getindex, [A, j, i])
end

function getindex(A::Arr, i::Int, j::Symbolic{<:Integer})
    wrap(unwrap(A)[i,j])
end

function getindex(A::Arr, j::Symbolic{<:Integer}, i::Int)
    wrap(unwrap(A)[j,i])
end

function _matmul(A,B)
    @syms i::Int j::Int k::Int
    if isadjointvec(A)
        op = operation(A.term)
        return op(op(B) * first(arguments(A.term)))
    end
    return @arrayop (A*B) (i, j) A[i, k] * B[k, j]
end

@wrapped (*)(A::AbstractMatrix, B::AbstractMatrix) = _matmul(A, B)
@wrapped (*)(A::AbstractVector, B::AbstractMatrix) = _matmul(A, B)

function _matvec(A,b)
    @syms i::Int k::Int
    if isdot(A, b)
        make_shape((i,), A[i, k] * b[k]) # This is a dimension check
        T = SymbolicUtils.promote_symtype(*, eltype(A), eltype(b))
        S = SymbolicUtils.promote_symtype(+, T,T)
        return Term{S}(*, [A, b])
    end
    @arrayop (A*b) (i,) A[i, k] * b[k]
end
@wrapped (*)(A::AbstractMatrix, b::AbstractVector) = _matvec(A, b)
(*)(A::AbstractArray, b::Union{AbstractVector, SymVec}) = _matvec(A, b)

#################### MAP-REDUCE ################
#

@wrapped Base.map(f, x::AbstractArray) = _map(f, x)
@wrapped Base.map(f, x::AbstractArray, xs...) = _map(f, x, xs...)
@wrapped Base.map(f, x, y::AbstractArray, z...) = _map(f, x, y, z...)
@wrapped Base.map(f, x, y, z::AbstractArray, w...) = _map(f, x, y, z, w...)

function _map(f, x, xs...)
    N = ndims(x)
    idx = makesubscripts(N)

    expr = f(map(a->a[idx...], [x, xs...])...)

    Atype = propagate_atype(map, f, x, xs...)
    ArrayOp(Atype{symtype(expr), N},
            (idx...,),
            expr,
            +,
            Term{Any}(map, [f, x, xs...]))
end

@inline _mapreduce(f, g, x, dims, kw) = mapreduce(f, g, x; dims=dims, kw...)

function scalarize_op(::typeof(_mapreduce), t)
    f,g,x,dims,kw = arguments(t)
    # we wrap and unwrap to make things work smoothly.
    # we need the result unwrapped to allow recursive scalarize to work.
    unwrap(_mapreduce(f, g, collect(wrap(x)), dims, kw))
end

@wrapped function Base.mapreduce(f, g, x::AbstractArray; dims=:, kw...)
    idx = makesubscripts(ndims(x))
    out_idx = [dims == (:) || i in dims ? 1 : idx[i] for i = 1:ndims(x)]
    expr = f(x[idx...])
    T = symtype(g(expr, expr))
    if dims === (:)
        return Term{T}(_mapreduce, [f, g, x, dims, (kw...,)])
    end

    Atype = propagate_atype(_mapreduce, f, g, x, dims, (kw...,))
    ArrayOp(Atype{T, ndims(x)},
            (out_idx...,),
            expr,
            g,
            Term{Any}(_mapreduce, [f, g, x, dims, (kw...,)]))
end

for (ff, opts) in [sum => (identity, +, false),
                  prod => (identity, *, true),
                  any => (identity, (|), false),
                  all => (identity, (&), true)]

    f, g, init = opts
    @eval @wrapped function (::$(typeof(ff)))(x::AbstractArray;
                                     dims=:, init=$init)
        mapreduce($f, $g, x, dims=dims, init=init)
    end
    @eval @wrapped function (::$(typeof(ff)))(f::Function, x::AbstractArray;
                                     dims=:, init=$init)
        mapreduce(f, $g, x, dims=dims, init=init)
    end
end


text:
module ConformingFESpacesTests

using Test
using Gridap
using Gridap.CellValuesGallery
using Gridap.ConformingFESpaces: CellEqClass

model = CartesianDiscreteModel(domain=(0.0,1.0,-1.0,2.0), partition=(2,2))

cell_to_nfaces = [[1,2,4,5,7,8,9,10,14],[2,3,5,6,11,12,10,13,15]]
nface_to_dofs = [[i*10,] for i in 1:15]

cell_to_nfaces = CellValueFromArray(cell_to_nfaces)
nface_to_dofs = CellValueFromArray(nface_to_dofs)

D = 2
order = 2
orders = fill(order,D)
polytope = Polytope(fill(HEX_AXIS,D)...)
fe = LagrangianRefFE(Float64,polytope, orders)

cell_to_dofs = CellEqClass(cell_to_nfaces,nface_to_dofs,fe)

r = [
  [10, 20, 40, 50, 70, 80, 90, 100, 140],
  [20, 30, 50, 60, 110, 120, 100, 130, 150]]
test_index_cell_array(cell_to_dofs,r)

model = CartesianDiscreteModel(domain=(0.0,1.0,-1.0,2.0), partition=(2,2))

D = pointdim(model)

grid = Grid(model,D)
trian = Triangulation(grid)
graph = GridGraph(model)
labels = FaceLabels(model)
tags = [1,2,3,4]

order = 1
orders = fill(order,D)
polytope = Polytope(fill(HEX_AXIS,D)...)
fe = LagrangianRefFE(Float64,polytope, orders)

fespace = ConformingFESpace(fe,trian,graph,labels,tags)

@test num_free_dofs(fespace) == 5
@test num_diri_dofs(fespace) == 4

@test diri_tags(fespace) === tags

r = [[-1, 1, 2, 3], [1, -2, 3, 4], [2, 3, -3, 5], [3, 4, 5, -4]]

@test r == collect(fespace.cell_eqclass)

order = 2
orders = fill(order,D)
polytope = Polytope(fill(HEX_AXIS,D)...)
fe = LagrangianRefFE(Float64,polytope, orders)

tags = [1,2,3,4,6,5]
fespace = ConformingFESpace(fe,trian,graph,labels,tags)

@test num_free_dofs(fespace) == 15
@test num_diri_dofs(fespace) == 10

r = [[-1, -2, 1, 2, -7, 4, 5, 6, 12],
  [-2, -3, 2, 3, -8, 7, 6, 8, 13],
  [1, 2, -4, -5, 4, -9, 9, 10, 14],
  [2, 3, -5, -6, 7, -10, 10, 11, 15]]

@test r == collect(fespace.cell_eqclass)

fun(x) = sin(x[1])*cos(x[2])

free_vals, diri_vals = interpolate_values(fespace,fun)

rf = [
  0.0, 0.420735, 0.73846, 0.217117, 0.0, 0.464521,
  0.598194, 0.815312, 0.0, 0.151174, 0.265335,
  0.239713, 0.660448, 0.078012, 0.214936]

rd = [
  0.0, 0.259035, 0.454649, -0.0, -0.199511,
  -0.350175, 0.133673, 0.368291, -0.102956, -0.283662]

@test isapprox(free_vals,rf,rtol=1.0e-5)
@test isapprox(diri_vals,rd,rtol=1.0e-5)

diri_vals = interpolate_diri_values(fespace,fun)
@test isapprox(diri_vals,rd,rtol=1.0e-5)

diri_vals = interpolate_diri_values(fespace,fill(fun,length(tags)))
@test isapprox(diri_vals,rd,rtol=1.0e-5)

uh = FEFunction(fespace,free_vals,diri_vals)

@test free_dofs(uh) === free_vals
@test diri_dofs(uh) === diri_vals
@test FESpace(uh) === fespace

zh = zero(fespace)

@test free_dofs(zh) == zeros(Float64,num_free_dofs(fespace))
@test diri_dofs(zh) == zeros(Float64,num_diri_dofs(fespace))
@test FESpace(zh) === fespace

U = TrialFESpace(fespace,fun)

zh = zero(U)
@test free_dofs(zh) == zeros(Float64,num_free_dofs(U))
@test diri_dofs(zh) === U.diri_dofs
@test FESpace(zh) === U

cellbasis = CellBasis(fespace)

quad = CellQuadrature(trian,degree=2)

a(v,u) = varinner(v,u)

bfun(x) = x[2]

b(v) = varinner(v,CellField(trian,bfun))

mmat = integrate(a(cellbasis,cellbasis),trian,quad)

bvec = integrate(b(cellbasis),trian,quad)

cellids = IdentityCellNumber(Int,length(bvec))
bvec2 = apply_constraints(fespace,bvec,cellids)
dofs = celldofids(fespace)

@test bvec2 === bvec

@test dofs == fespace.cell_eqclass

mmat2 = apply_constraints_rows(fespace,mmat,cellids)
dofs = celldofids(fespace)

@test mmat2 === mmat

@test dofs == fespace.cell_eqclass

mmat3 = apply_constraints_cols(fespace,mmat,cellids)
dofs = celldofids(fespace)

@test mmat3 === mmat

@test dofs == fespace.cell_eqclass

uh = interpolate(fespace,fun)
@test isa(uh,FEFunction)

q = coordinates(quad)
uhq = evaluate(uh,q)

grad_uh = gradient(uh)
grad_uhq = evaluate(grad_uh,q)

v = collect(uhq)
g = collect(grad_uhq)

test_index_cell_field(uh,q,v,g)

test_fe_space(fespace,15,10,mmat,bvec,fun)

fespace = H1ConformingFESpace(Float64,model,order,tags)

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))

order = 1
fespace = H1ConformingFESpace(Float64,model,order,tags)

grid = Grid(model,D)
trian = Triangulation(grid)

fun1(x) = x[1]
uh1 = interpolate(fespace,fun1)

fun2(x) = x[1]*x[2]
uh2 = interpolate(fespace,fun2)

fun3(x) = sin(x[1])*cos(x[2])
uh3 = interpolate(fespace,fun3)

quad = CellQuadrature(trian,degree=2)
q = coordinates(quad)
uhq = evaluate(uh1,q)

# Vector valued

order = 1
T = VectorValue{2,Float64}
tags = [1,2,3,4]
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))
fespace = H1ConformingFESpace(T,model,order,tags)

ufun(x) = VectorValue(x[2],x[1])
uh = interpolate(fespace,ufun)

@test free_dofs(uh) == [0.0, 0.5, 0.5, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 0.5]
@test diri_dofs(uh) == [0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0]

end # module


text:
# Basic test of DERelative

using Test
using DifferentialEquations
using IncrementalInference
using Dates
using Statistics

## plotting functions

# using Plots
# using Cairo, RoMEPlotting
# Gadfly.set_default_plot_size(25cm,20cm)

##

@testset "First order DERelative" begin

##

function firstOrder!(dstate, state, force, t)
  β = -0.2
  dstate[1] = β*state[1] + force(t)
  nothing
end

# testing function parameter version (could also be array of data)
tstForce(t) = 0

## build a representative factor graph with ODE built inside

fg = initfg()
# the starting points and "0 seconds"
addVariable!(fg, :x0, ContinuousScalar, timestamp=DateTime(2000,1,1,0,0,0))
# pin with a simple prior
addFactor!(fg, [:x0], Prior(Normal(1,0.01)))

doautoinit!(fg, :x0)

prev = :x0

for i in 1:3

  nextSym = Symbol("x$i")

  # another point in the trajectory 5 seconds later
  addVariable!(fg, nextSym, ContinuousScalar, timestamp=DateTime(2000,1,1,0,0,5*i))
  oder = IIF.DERelative(fg, [prev; nextSym],
                        ContinuousEuclid{1},
                        firstOrder!,
                        tstForce,
                        dt=0.05,
                        problemType=ODEProblem )
  #
  addFactor!( fg, [prev;nextSym], oder, graphinit=false )
  initManual!(fg, nextSym, zeros(1,100))

  prev = nextSym
end

## basic sample test

meas = sampleFactor(fg, :x0x1f1, 10)
@test size(meas[1],1) == 1
@test size(meas[1],2) == 10

## do all forward solutions

pts, = sampleFactor(fg, :x0f1, 100)
initManual!(fg, :x0, pts)
pts = approxConv(fg, :x0x1f1, :x1)
@test 0.3 < Statistics.mean(pts) < 0.4

## check that the reverse solve also works

initManual!(fg, :x1, pts)
pts = approxConv(fg, :x0x1f1, :x0)

# check the reverse solve to be relatively accurate
@test (pts - (getBelief(fg, :x0) |> getPoints)) |> norm < 1e-4

##

oder_ = DERelative( fg, [:x0; :x3],
                    ContinuousEuclid{1},
                    firstOrder!,
                    tstForce,
                    dt=0.05,
                    problemType=ODEProblem )

oder_.forwardProblem.u0 .= [1.0]
sl = DifferentialEquations.solve(oder_.forwardProblem)

##

# Plots.plot(sl,linewidth=2,xaxis="unixtime [s]",layout=(1,1))

# for lb in [:x0; :x1;:x2;:x3]
#   x = getTimestamp(getVariable(fg, lb)) |> DateTime |> datetime2unix
#   xx = [x;x]
#   yy = [0;1]
#   Plots.plot!(xx, yy, show=true)
# end

##

tfg = initfg()
pts = approxConv(fg, :x0f1, :x3, setPPE=true, tfg=tfg)
# initManual!(tfg, :x3, pts)

##

@test getPPE(tfg, :x0).suggested - sl(getVariable(fg, :x0) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test getPPE(tfg, :x1).suggested - sl(getVariable(fg, :x1) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test getPPE(tfg, :x2).suggested - sl(getVariable(fg, :x2) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test       Statistics.mean(pts) - sl(getVariable(fg, :x3) |> getTimestamp |> DateTime |> datetime2unix)[1] < 0.1

##

# plotKDE(tfg, [:x0;:x1;:x2;:x3])

## Now test a full solve

solveTree!(fg);

##

@test getPPE(fg, :x0).suggested - sl(getVariable(fg, :x0) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test getPPE(fg, :x1).suggested - sl(getVariable(fg, :x1) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test getPPE(fg, :x2).suggested - sl(getVariable(fg, :x2) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1
@test getPPE(fg, :x3).suggested - sl(getVariable(fg, :x3) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.1

##

end

##

@testset "Damped Oscillator DERelative" begin

## setup some example dynamics

# Lets build an damped oscillator to demonstrate the process in state space
# https://en.wikipedia.org/wiki/Harmonic_oscillator
# ddx/ddt = β dx/dt  -  ω x  +  force[t]
# dx/dt   = dx/dt
function dampedOscillator!(dstate, state, force, t)
  ω = 0.7
  β = -0.3
  dstate[2] = β*state[2] - ω*state[1] + force(t)
  dstate[1] = state[2]
  nothing
end

# testing function parameter version (could also be array of data)
tstForce(t) = 0

## build a representative factor graph with ODE built inside

fg = initfg()

# the starting points and "0 seconds"
addVariable!(fg, :x0, ContinuousEuclid{2}, timestamp=DateTime(2000,1,1,0,0,0))
# pin with a simple prior
addFactor!(fg, [:x0], Prior(MvNormal([1;0],0.01*diagm(ones(2)))))

##

prev = :x0
DT = 2

for i in 1:7

  nextSym = Symbol("x$i")

  # another point in the trajectory 5 seconds later
  addVariable!(fg, nextSym, ContinuousEuclid{2}, timestamp=DateTime(2000,1,1,0,0,DT*i))
  oder = DERelative( fg, [prev; nextSym],
                      ContinuousEuclid{2},
                      dampedOscillator!,
                      tstForce,
                      # (state, var)->(state[1] = var[1]),
                      # (var, state)->(var[1] = state[1]),
                      dt=0.05,
                      problemType=ODEProblem )
  #
  addFactor!( fg, [prev;nextSym], oder )

  prev = nextSym
end

## check forward and backward solving

pts = approxConv(fg, :x0f1, :x0)
@test norm(Statistics.mean(pts, dims=2) - [1;0]) < 0.3

initManual!(fg, :x0, pts)
X0_ = deepcopy(pts)

pts = approxConv(fg, :x0x1f1, :x1)
@test norm(Statistics.mean(pts, dims=2) - [0;-0.6]) < 0.4

# now check the reverse direction solving
initManual!(fg, :x1, pts)
pts = approxConv(fg, :x0x1f1, :x0)

@test (X0_ - pts) |> norm < 1e-4

##

tfg = initfg()
for s in ls(fg)
  initManual!(fg, s, zeros(2,100))
end

pts = approxConv(fg, :x0f1, :x7, setPPE=true, tfg=tfg)
# initManual!(tfg, :x7, pts)

##

# plotKDE(tfg, ls(fg) |> sortDFG, dims=[1] )

##

oder_ = DERelative( fg, [:x0; :x7],
                    ContinuousEuclid{2},
                    dampedOscillator!,
                    tstForce,
                    # (state, var)->(state[1] = var[1]),
                    # (var, state)->(var[1] = state[1]),
                    dt=0.05,
                    problemType=ODEProblem )

oder_.forwardProblem.u0 .= [1.0;0.0]
sl = DifferentialEquations.solve(oder_.forwardProblem)

## check the solve values are correct

for sym = ls(tfg)
  @test getPPE(tfg, sym).suggested - sl(getVariable(fg, sym) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.2
end

##

# Plots.plot(sl,linewidth=2,xaxis="unixtime [s]",label=["ω [rad/s]" "θ [rad]"],layout=(2,1))

# for lb in sortDFG(ls(fg))
#   x = getTimestamp(getVariable(tfg, lb)) |> DateTime |> datetime2unix
#   xx = [x;x]
#   yy = [-1;1]
#   Plots.plot!(xx, yy, show=true)
# end

##

@error "Disabling useMsgLikelihood for DERelative test, follow fix on #1010 as rough guide"
getSolverParams(fg).useMsgLikelihoods = false

solveTree!(fg);

##

for sym = ls(fg)
  @test getPPE(fg, sym).suggested - sl(getVariable(fg, sym) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.2
end

##

end

##

@testset "Parameterized Damped Oscillator DERelative" begin

## setup some example dynamics

# Lets build an damped oscillator to demonstrate the process in state space
# https://en.wikipedia.org/wiki/Harmonic_oscillator
# ddx/ddt = β dx/dt  -  ω x  +  force[t]
# dx/dt   = dx/dt
# force_ωβ = (data, ωβ)
function dampedOscillatorParametrized!(dstate, state, force_ωβ, t)
  # 3rd variable in this factor graph test example
  force = force_ωβ[1]
  ω     = force_ωβ[2][1]
  β     = force_ωβ[2][2]
  # classic ODE between first and second fg variables
  dstate[2] = β*state[2] - ω*state[1] + force(t)
  dstate[1] = state[2]
  nothing
end

# testing function parameter version (could also be array of data)
tstForce(t) = 0

## build a representative factor graph with ODE built inside

fg = initfg()

# the starting points and "0 seconds"
addVariable!(fg, :x0, ContinuousEuclid{2}, timestamp=DateTime(2000,1,1,0,0,0))
# pin with a simple prior
addFactor!(fg, [:x0], Prior(MvNormal([1;0],0.01*diagm(ones(2)))))
doautoinit!(fg, :x0)

# and the new parameterized variable
ω = 0.7
β = -0.3

# these are the stochastic parameters
addVariable!(fg, :ωβ, ContinuousEuclid{2}) # timestamp should not matter
# pin with a simple prior
addFactor!(fg, [:ωβ], Prior(MvNormal([ω;β],0.0001*diagm(ones(2)))))
doautoinit!(fg, :ωβ)

##

prev = :x0
DT = 2

for i in 1:7

  nextSym = Symbol("x$i")

  # another point in the trajectory 5 seconds later
  addVariable!(fg, nextSym, ContinuousEuclid{2}, timestamp=DateTime(2000,1,1,0,0,DT*i))
  oder = DERelative( fg, [prev; nextSym; :ωβ],
                      ContinuousEuclid{2},
                      dampedOscillatorParametrized!,
                      tstForce, # this is passed in as `force_ωβ[1]`
                      # (state, var)->(state[1] = var[1]),
                      # (var, state)->(var[1] = state[1]),
                      # dt=0.05,
                      problemType=ODEProblem )
  #
  addFactor!( fg, [prev; nextSym; :ωβ], oder, graphinit=false )

  prev = nextSym
end

## check forward and backward solving

pts = approxConv(fg, :x0f1, :x0)
@test norm(Statistics.mean(pts, dims=2) - [1;0]) < 0.3

initManual!(fg, :x0, pts)
X0_ = deepcopy(pts)

pts = approxConv(fg, :x0x1ωβf1, :x1)
@test norm(Statistics.mean(pts, dims=2) - [0;-0.6]) < 0.4

# now check the reverse direction solving
initManual!(fg, :x1, pts)

# failing here
pts = approxConv(fg, :x0x1ωβf1, :x0)

@test (X0_ - pts) |> norm < 1e-2

##

tfg = initfg()
for s in ls(fg)
  initManual!(fg, s, zeros(2,100))
end

# must initialize the parameters
pts = approxConv(fg, :ωβf1, :ωβ)
initManual!(fg, :ωβ, pts)

# project forward
forcepath = [:x0f1;]
push!(forcepath, :x0)
push!(forcepath, :x0x1ωβf1)
push!(forcepath, :x1)
push!(forcepath, :x1x2ωβf1)
push!(forcepath, :x2)
push!(forcepath, :x2x3ωβf1)
push!(forcepath, :x3)
push!(forcepath, :x3x4ωβf1)
push!(forcepath, :x4)
push!(forcepath, :x4x5ωβf1)
push!(forcepath, :x5)
push!(forcepath, :x5x6ωβf1)
push!(forcepath, :x6)
push!(forcepath, :x6x7ωβf1)
push!(forcepath, :x7)
pts = approxConv(fg, :x0f1, :x7, setPPE=true, tfg=tfg, path=forcepath)

##

# plotKDE(tfg, ls(tfg) |> sortDFG, dims=[1] )

##

# getBelief(fg, :ωβ) |> getPoints

# plotKDE(tfg, :ωβ)

##

oder_ = DERelative( fg, [:x0; :x7; :ωβ],
                    ContinuousEuclid{2},
                    dampedOscillatorParametrized!,
                    tstForce,
                    # (state, var)->(state[1] = var[1]),
                    # (var, state)->(var[1] = state[1]),
                    dt=0.05,
                    problemType=ODEProblem )

oder_.forwardProblem.u0 .= [1.0;0.0]
oder_.data[2] .= [ω;β]
sl = DifferentialEquations.solve(oder_.forwardProblem)

## check the approxConv is working right

for sym in setdiff(ls(tfg), [:ωβ])
  @test getPPE(tfg, sym).suggested - sl(getVariable(fg, sym) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.2
end

##

# Plots.plot(sl,linewidth=2,xaxis="unixtime [s]",label=["ω [rad/s]" "θ [rad]"],layout=(2,1))

# for lb in sortDFG(ls(fg))
#   x = getTimestamp(getVariable(tfg, lb)) |> DateTime |> datetime2unix
#   xx = [x;x]
#   yy = [-1;1]
#   Plots.plot!(xx, yy, show=true)
# end

## test convolution to the parameter (third) variable

# easy test with good starting points
pts = approxConv(fg, :ωβf1, :ωβ)
initManual!(fg, :ωβ, pts)

# make sure the other variables are in the right place
@test Statistics.mean(getBelief(fg, :x0) |> getPoints, dims=2) - [1;0] |> norm < 0.1
@test Statistics.mean(getBelief(fg, :x1) |> getPoints, dims=2) - [0;-0.6] |> norm < 0.2

pts = approxConv(fg, :x0x1ωβf1, :ωβ)

@test Statistics.mean(pts, dims=2) - [0.7;-0.3] |> norm < 0.1

# repeat with more difficult starting point

initManual!(fg, :ωβ, zeros(2,100))

pts = approxConv(fg, :x0x1ωβf1, :ωβ)

@test Statistics.mean(pts, dims=2) - [0.7;-0.3] |> norm < 0.1

@warn "n-ary DERelative test on :ωβ requires issue #1010 to be resolved first before being reintroduced."
# ## do a complete solve (must first resolve #1010)

# solveTree!(fg);

# ## Solve quality might not yet be good enough for this particular test case

# @test getPPE(fg, :ωβ).suggested - [0.7;-0.3] |> norm < 0.2

# for sym in setdiff(ls(tfg), [:ωβ])
#   @test getPPE(fg, sym).suggested - sl(getVariable(fg, sym) |> getTimestamp |> DateTime |> datetime2unix) |> norm < 0.2
# end

##

end

@error "DERelative not tested for `multihypo=` case yet, see issue #1025"

#


text:
function annotate(vcf_path)
    shell('table_annovar.pl anno_tmp.vcf ~/tools/annovar-170715/humandb '
            '-buildver hg19 --remove --otherinfo --outfile annotated '
            '-operation g,f,f,f '
            '-protocol refGene,cosmic70,1000g2014oct_all,esp6500si_alltable_annovar.pl anno_tmp.vcf ~/tools/annovar-170715/humandb '
            '-buildver hg19 --remove --otherinfo --outfile annotated '
            '-operation g,f,f,f '
            '-protocol refGene,cosmic70,1000g2014oct_all,esp6500si_all')
end

#=
function annotate(vcf_path, method, dbs,)
    shell('table_annovar.pl anno_tmp.vcf ~/tools/annovar-170715/humandb '
            '-buildver hg19 --remove --otherinfo --outfile annotated '
            '-operation g,f,f,f '
            '-protocol refGene,cosmic70,1000g2014oct_all,esp6500si_alltable_annovar.pl anno_tmp.vcf ~/tools/annovar-170715/humandb '
            '-buildver hg19 --remove --otherinfo --outfile annotated '
            '-operation g,f,f,f '
            '-protocol refGene,cosmic70,1000g2014oct_all,esp6500si_all')
end
=#


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

using FMI

pathToFMU = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

display("Event Indicators: $(myFMU.modelDescription.numberOfEventIndicators)")

fmiInstantiate!(myFMU; loggingOn=true)

fmiSetupExperiment(myFMU, 0.0)

fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

dt = 0.01
t_start = 0.0
t_stop = 8.0

ts = t_start:dt:t_stop
ss = []

solution = fmiSimulateME(myFMU, dt, t_start, t_stop)
fmiPlot(myFMU, solution)

fmiUnload(myFMU)


text:
# Handle faces.

struct FaceMesh{T}
end


text:
using Test
using DECAES
using DECAES.NNLS
using LinearAlgebra

function verify_solution(m, n)
    A = randn(m, n)
    b = randn(m)
    work = NNLSWorkspace(A, b)
    nnls!(work)

    GC.@preserve work begin
        x  = NNLS.solution(work)
        w  = NNLS.dual(work)
        i₊ = NNLS.components(work)
        m₊ = NNLS.ncomponents(work)
        m₀ = m - m₊
        i₀ = setdiff(1:length(x), i₊)
        x₊, x₀ = x[i₊], x[i₀]
        w₀, w₋ = w[i₊], w[i₀]
        A₊, A₀ = A[:,i₊], A[:,i₀]
        L  = NNLS.choleskyfactor(work, Val(:L))
        U  = NNLS.choleskyfactor(work, Val(:U))

        # Internals
        @test isperm(work.idx)

        # Solution
        @test all(>(0), x₊)
        @test all(==(0), x₀)
        @test x₊ ≈ (A₊'A₊) \ (A₊'b)
        @test x₊ ≈ A₊ \ b

        # Dual (i.e. gradient)
        if m₊ < m
            # Less components than rows; in general solution is not exact and gradient has negative components
            @test all(<(0), w₋)
            @test all(==(0), w₀)
            @test w₋ ≈ A₀'b - A₀'A₊ * x₊
        else
            # At least as many components as rows; solution is exact, gradient is zero
            @test all(==(0), w)
        end

        # Cholesky factors
        @test A₊'A₊ ≈ L'L
        @test A₊'A₊ ≈ U*U'
    end
    work
end

@testset "NNLS" begin
    for m in [1,2,5,8,16,25,32,50], n in [1,2,5,8,16,25,32,50]
        verify_solution(m, n)
    end
end


text:
function bracket_minimum(f, x=0; s=1e-2, k=2.0)
    a, ya = x, f(x)
    b, yb = a + s, f(a + s)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        s = -s
    end
    while true
        c, yc = b + s, f(b + s)
        if yc > yb
            return a < c ? (a,b, c) : (c,b, a)
        end
        a, ya, b, yb = b, yb, c, yc
        s *= k
    end
end

function quadratic_fit_search(f, a, b, c, n)
    ya, yb, yc = f(a), f(b), f(c)
    for i in 1:n-3
        x = 0.5*(ya*(b^2-c^2)+yb*(c^2-a^2)+yc*(a^2-b^2)) /
            (ya*(b-c) +yb*(c-a) +yc*(a-b))
        yx = f(x)
        if x > b
            if yx > yb
                c, yc = x, yx
            else
                a, ya, b, yb = b, yb, x, yx
            end
        elseif x < b
            if yx > yb
                a, ya = x, yx
            else
                c, yc, b, yb = b, yb, x, yx
            end
        end
    end
    return (a, b, c)
end

φ = MathConstants.φ
function golden_section_search(f, a, b, n)
    ρ = φ-1
    d = ρ * b + (1 - ρ)*a
    yd = f(d)
    for i = 1 : n-1
        c = ρ*a + (1 - ρ)*b
        yc = f(c)
        if yc < yd
            b, d, yd = d, c, yc
        else
            a, b = b, c
        end
    end
    return a < b ? (a, b) : (b, a)
end

function line_search(f, x, d)
    objective = α -> f(x + α*d)
    a, b,c = bracket_minimum(objective)

    #numero de iteracoes
    n=(c-a)/(10^-3*log(φ))

    #pelo metodo golden section search
    α,beta = golden_section_search(objective, a, c,n)

    #pelo metodo quadratic fit search
    #α,beta,gamma = quadratic_fit_search(objective, a, b, c,n)
    return (x + α*d,α)
end

#funcao multivariavel
f(x) = sin(x[1]*x[2])+exp(x[2]+x[3])-x[3]

#ponto inicial
x0 = [1,2,3]

#direcao de busca
d = [0,-1,-1]

(x,α) = line_search(f, x0, d)
println("x* = ", x)


text:
module Potato

using JuMP
using XLSX
using DataFrames
using Distances
using UUIDs
using Statistics

include("struct.jl")
include("settings.jl")
include("base.jl")
include("core.jl")
include("search.jl")

#struct
export vtx,veh,col,node,dv,β,bound,S

#settings
export set_optimizer!
export get_optimizer
export reset_optimizer
export set_slack_coeff!
export set_surp_coeff!
export su_C
export sl_C
export silence!
export silent

#base
export extract!
export V,dist,K,T,d
export root

#core
export passes
export callSub,buildSub!
export master,getDuals
export sub,getCols
export colGen
export origin
export Q,s,f
export prioritize_vertex!,i_priority
export prioritize_vehicle!,k_priority
export fractional_column
export find_branch
export separate
export find_separator

#search
export leaf,traverse
export createBranch

end