From 6a228b56107105739e324f855108ef45a5691c11 Mon Sep 17 00:00:00 2001 From: Dan Date: Tue, 5 May 2020 15:34:34 -0700 Subject: [PATCH 01/10] removing coefficient_functions.jl and the corresponding test these do not adhere to the documented interface for coeff_func, and also contain a uninitialized memory bug --- .../coefficient_functions.jl | 67 -------------- test/coefficient_functions.jl | 88 ------------------- 2 files changed, 155 deletions(-) delete mode 100644 src/derivative_operators/coefficient_functions.jl delete mode 100644 test/coefficient_functions.jl diff --git a/src/derivative_operators/coefficient_functions.jl b/src/derivative_operators/coefficient_functions.jl deleted file mode 100644 index b2e816e58..000000000 --- a/src/derivative_operators/coefficient_functions.jl +++ /dev/null @@ -1,67 +0,0 @@ -""" -``` -compute_coeffs(coeff_func, current_coeffs) -``` -calculates the coefficients for the stencil of UpwindDifference operators. -""" -function compute_coeffs!(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number} - return current_coeffs .+= coeff_func -end - -function compute_coeffs!(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number} - return current_coeffs[:] += coeff_func -end - -# Coefficient functions when coeff_func is a Function and current_coeffs exists -function compute_coeffs!(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number} - if hasmethod(coeff_func, (Vector{T},)) - current_coeffs[:] = coeff_func(current_coeffs) - else - map!(coeff_func, current_coeffs, current_coeffs) - end - return current_coeffs - # if hasmethod(coeff_func, (Int,)) # assume we want to provide the index of coefficients to coeff_func - # for i = 1:length(current_coeffs) - # current_coeffs[i] += coeff_func(i) - # end - # return current_coeffs - # elseif hasmethod(coeff_func, (Vector{Int},)) # assume we want to provide the index of coefficients to coeff_func - # current_coeffs[:] += coeff_func(collect(1:length(current_coeffs))) - # return current_coeffs - # elseif hasmethod(coeff_func, (UnitRange{Int},)) # assume we want to provide the index of coefficients to coeff_func - # current_coeffs[:] += coeff_func(1:length(current_coeffs)) - # return current_coeffs - # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients - # map!(coeff_func, current_coeffs, current_coeffs) - # return current_coeffs - # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients - # current_coeffs[:] = coeff_func(current_coeffs) - # return current_coeffs - # else - # error("Coefficient functions with the arguments of $coeff_func have not been implemented.") - # end -end - -compute_coeffs(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number} = current_coeffs .+ coeff_func -compute_coeffs(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number} = coeff_func + current_coeffs - -function compute_coeffs(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number} - if hasmethod(coeff_func, (Vector{T},)) - return coeff_func(current_coeffs) - else - return map(coeff_func, current_coeffs) - end - # if hasmethod(coeff_func, (Int,)) - # return current_coeffs + map(coeff_func, collect(1:length(current_coeffs))) - # elseif hasmethod(coeff_func, (Vector{Int},)) - # return current_coeffs + coeff_func(collect(1:length(current_coeffs))) - # elseif hasmethod(coeff_func, (UnitRange{Int},)) - # return current_coeffs + coeff_func(1:length(current_coeffs)) - # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients - # return map(coeff_func, current_coeffs) - # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients - # return coeff_func(current_coeffs) - # else - # error("Coefficient functions with the arguments of coeff_func $coeff_func have not been implemented.") - # end -end diff --git a/test/coefficient_functions.jl b/test/coefficient_functions.jl deleted file mode 100644 index 5e7a5cf23..000000000 --- a/test/coefficient_functions.jl +++ /dev/null @@ -1,88 +0,0 @@ -using Test, DiffEqOperators - -# Set up coefficient functions to test with -vec_fcn = Vector{Function}(undef,0) -f1(x::Float64) = sin(x) -f2(x::Vector{Float64}) = sin.(x) -push!(vec_fcn, f1) -push!(vec_fcn, f2) - -@testset "Test coefficient functions when current_coeffs exists" begin - - vec_fcn_ans = Vector{Vector{Float64}}(undef,0) - for i = 1:2 - push!(vec_fcn_ans, [sin(1.1), sin(2.3), sin(4.5)]) - end - - current_coeffs = [1.1, 2.3, 4.5] - - # Test inplace versions - current_coeffs1 = Vector{Float64}(undef,3) - @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ 2.54 .* ones(3) - current_coeffs1 = Vector{Float64}(undef,3) - @test DiffEqOperators.compute_coeffs!([1.1, 2.3, 4.5], current_coeffs1) ≈ current_coeffs - @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ [3.64, 4.84, 7.04] - @test all(current_coeffs .!= current_coeffs1) - current_coeffs1 = copy(current_coeffs) - @test DiffEqOperators.compute_coeffs!([2.54, 1.03, .18], current_coeffs1) ≈ [3.64, 3.33, 4.68] - @test all(current_coeffs .!= current_coeffs1) - for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans) - current_coeffs1 = copy(current_coeffs) - @test DiffEqOperators.compute_coeffs!(fcn, current_coeffs1) ≈ fcn_ans - @test all(current_coeffs1 .!= current_coeffs) - end - - # Test not in place versions - current_coeffs1 = Vector{Float64}(undef,3) - @test DiffEqOperators.compute_coeffs(2.54, current_coeffs) ≈ [3.64, 4.84, 7.04] - @test DiffEqOperators.compute_coeffs([2.54, 1.03, 0.18], current_coeffs) ≈ [3.64, 3.33, 4.68] - for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans) - @test DiffEqOperators.compute_coeffs(fcn, current_coeffs) ≈ fcn_ans - end -end - -@testset "Check coefficients of DerivativeOperators are properly computed" begin - - # Set up operators to test construction with different coeff_func - vec_op = Vector{DerivativeOperator}(undef,0) - dxv = [1.4, 1.1, 2.3, 3.6, 1.5] - push!(vec_op, UpwindDifference(1, 1, 1., 3, 1.)) - push!(vec_op, UpwindDifference(1, 1, 1., 3, [1., 2.25, 3.5])) - push!(vec_op, UpwindDifference(1, 1, 1., 3, cos)) - push!(vec_op, UpwindDifference(1, 1, dxv, 3, 1.)) - push!(vec_op, UpwindDifference(1, 1, dxv, 3, [1., 2.25, 3.5])) - push!(vec_op, UpwindDifference(1, 1, dxv, 3, cos)) - - - vec_op_ans = Vector{Any}(undef,0) - push!(vec_op_ans, ones(3)) - push!(vec_op_ans, [1., 2.25, 3.5]) - push!(vec_op_ans, ones(3)) - push!(vec_op_ans, ones(3)) - push!(vec_op_ans, [1., 2.25, 3.5]) - push!(vec_op_ans, ones(3)) - - # Check constructors - @test UpwindDifference(1, 1, 1., 3).coefficients + ones(3) - ones(3) ≈ zeros(3) - @test CenteredDifference(2, 2, 1., 3, 0.).coefficients + ones(3) - ones(3) ≈ zeros(3) - for (L, op_ans) in zip(vec_op, vec_op_ans) - @test L.coefficients ≈ op_ans - end - push!(vec_op, CenteredDifference(2, 2, 1., 3)) == nothing - - # Compute answers to coeff_func * operator - func_mul_op_ans1 = Vector{Any}(undef,0) - for i = 1:2 - push!(func_mul_op_ans1, zeros(3)) - end - func_mul_op_ans2 = Vector{Any}(undef,0) - for i = 1:2 - push!(func_mul_op_ans2, sin(1) .* ones(3)) - end - func_mul_op_ans3 = Vector{Any}(undef,0) - for i = 1:2 - push!(func_mul_op_ans3, [0., sin(1.5), 0.]) - end -end - -nothing From 7530b13167af88ee7497f60c961a11afe94433f8 Mon Sep 17 00:00:00 2001 From: Dan Date: Tue, 5 May 2020 15:39:50 -0700 Subject: [PATCH 02/10] whitespace --- src/derivative_operators/derivative_operator.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index c32d91bdd..03962d398 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -161,13 +161,13 @@ function UpwindDifference{N}(derivative_order::Int, low_boundary_x = 0.0:(boundary_stencil_length-1) L_boundary_deriv_spots = 1.0:boundary_stencil_length - 2.0 - _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots] - low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) - high_boundary_x = 0.0:-1.0:-(boundary_stencil_length-1.0) + high_boundary_x = 0.0:-1.0:-(boundary_stencil_length-1.0) R_boundary_deriv_spots = -1.0:-1.0:-(boundary_stencil_length-2.0) - _high_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots] - high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) + _high_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots] + high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) coefficients = Vector{T}(undef,len) if coeff_func != nothing From 965bea4e75a0f950c2b21b78ac74f84a60fd3c78 Mon Sep 17 00:00:00 2001 From: Dan Date: Tue, 5 May 2020 15:42:28 -0700 Subject: [PATCH 03/10] remove references to coefficient_functions.jl and related test --- src/DiffEqOperators.jl | 1 - test/runtests.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 5438a85d1..5712c2052 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -33,7 +33,6 @@ include("derivative_operators/convolutions.jl") include("derivative_operators/concretization.jl") include("derivative_operators/ghost_derivative_operator.jl") include("derivative_operators/derivative_operator_functions.jl") -include("derivative_operators/coefficient_functions.jl") ### Composite Operators include("composite_operators.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 8c6da1a5c..bab478b5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,7 +29,6 @@ if GROUP == "All" || GROUP == "Interface" @time @safetestset "Convolutions" begin include("convolutions.jl") end @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end @time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end - @time @safetestset "Coefficient Functions" begin include("coefficient_functions.jl") end @time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end @time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end @time @safetestset "Basic SDO Examples" begin include("BasicSDOExamples.jl") end From 5c24d99f6c4053fc164c35cc74dcfae598133c93 Mon Sep 17 00:00:00 2001 From: Dan Date: Tue, 5 May 2020 19:27:28 -0700 Subject: [PATCH 04/10] initialize DerivativeOperator.coefficients to ones instead of undef when coeff_func is provided --- .../derivative_operator.jl | 26 +++++++------------ 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 03962d398..c329aa846 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,6 +26,7 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end + struct CenteredDifference{N} end function CenteredDifference{N}(derivative_order::Int, @@ -53,7 +54,7 @@ function CenteredDifference{N}(derivative_order::Int, high_boundary_coefs = convert(SVector{boundary_point_count},reverse(map(reverse, _low_boundary_coefs*(-1)^derivative_order))) - coefficients = coeff_func isa Nothing ? nothing : Vector{T}(undef,len) + coefficients = coeff_func isa Nothing ? nothing : ones(T, len) DerivativeOperator{T,N,false,T,typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), @@ -103,7 +104,7 @@ function CenteredDifference{N}(derivative_order::Int, calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1] high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) - coefficients = coeff_func isa Nothing ? nothing : Vector{T}(undef,len) + coefficients = coeff_func isa Nothing ? nothing : ones(T, len) DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), @@ -128,12 +129,12 @@ constructs a DerivativeOperator that automatically implements upwinding. ### Inputs * `dx::T` or `dx::Vector{T}`: grid spacing -* `coeff_func`: function mapping index in the grid to coefficient at that grid location +* `coeff_func`: Same as for DerivativeOperator. ### Examples julia> drift = [1., 1., -1.] -julia> L1 = UpwindDifference(1, 1, 1., 3, i -> drift[i]) -julia> L2 = UpwindDifference(1, 1, 1., 3, i -> 1.) +julia> L1 = UpwindDifference(1, 1, 1., 3, drift) +julia> L2 = UpwindDifference(1, 1, 1., 3, 1.) julia> Q = Neumann0BC(1, 1.) julia> Array(L1 * Q)[1] 3×3 Array{Float64,2}: @@ -169,13 +170,10 @@ function UpwindDifference{N}(derivative_order::Int, _high_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots] high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) - coefficients = Vector{T}(undef,len) - if coeff_func != nothing - compute_coeffs!(coeff_func, coefficients) - end + coefficients = coeff_func isa Nothing ? nothing : ones(T, len) DerivativeOperator{T,N,true,T,typeof(stencil_coefs), - typeof(low_boundary_coefs),Vector{T}, + typeof(low_boundary_coefs),typeof(coefficients), typeof(coeff_func)}( derivative_order, approximation_order, dx, len, stencil_length, stencil_coefs, @@ -214,14 +212,10 @@ function UpwindDifference{N}(derivative_order::Int, _downwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2:i+1])) for i in len-boundary_point_count+1:len]) high_boundary_coefs = [_upwind_coefs ; _downwind_coefs] - # Compute coefficients - coefficients = Vector{T}(undef,len) - if coeff_func != nothing - compute_coeffs!(coeff_func, coefficients) - end + coefficients = coeff_func isa Nothing ? nothing : ones(T, len) DerivativeOperator{T,N,true,typeof(dx),typeof(stencil_coefs), - typeof(low_boundary_coefs),Vector{T}, + typeof(low_boundary_coefs),typeof(coefficients), typeof(coeff_func)}( derivative_order, approximation_order, dx, len, stencil_length, stencil_coefs, From fcdf7752ace3890255fa7cf325668a7dc7e2db1d Mon Sep 17 00:00:00 2001 From: Dan Date: Tue, 5 May 2020 19:45:24 -0700 Subject: [PATCH 05/10] add function to initialize coefficients based on coeff_func --- .../derivative_operator.jl | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index c329aa846..e59b51b9c 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,6 +26,18 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end +init_coefficients(coeff_func::Nothing, len::Int) = nothing + +init_coefficients(coeff_func::Number, len::Int) = coeff_func * ones(typeof(coeff_func), len) + +function init_coefficients(coeff_func::AbstractVector{T}, len::Int) where T <: Number + coeff_func +end + +function init_coefficients(coeff_func::Function, len::Int) + ones(Float64, len) +end + struct CenteredDifference{N} end @@ -54,7 +66,7 @@ function CenteredDifference{N}(derivative_order::Int, high_boundary_coefs = convert(SVector{boundary_point_count},reverse(map(reverse, _low_boundary_coefs*(-1)^derivative_order))) - coefficients = coeff_func isa Nothing ? nothing : ones(T, len) + coefficients = init_coefficients(coeff_func, len) DerivativeOperator{T,N,false,T,typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), @@ -104,7 +116,7 @@ function CenteredDifference{N}(derivative_order::Int, calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1] high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) - coefficients = coeff_func isa Nothing ? nothing : ones(T, len) + coefficients = init_coefficients(coeff_func, len) DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), @@ -170,7 +182,7 @@ function UpwindDifference{N}(derivative_order::Int, _high_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots] high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs) - coefficients = coeff_func isa Nothing ? nothing : ones(T, len) + coefficients = init_coefficients(coeff_func, len) DerivativeOperator{T,N,true,T,typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), @@ -212,7 +224,7 @@ function UpwindDifference{N}(derivative_order::Int, _downwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2:i+1])) for i in len-boundary_point_count+1:len]) high_boundary_coefs = [_upwind_coefs ; _downwind_coefs] - coefficients = coeff_func isa Nothing ? nothing : ones(T, len) + coefficients = init_coefficients(coeff_func, len) DerivativeOperator{T,N,true,typeof(dx),typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(coefficients), From f80886e537bfeb71b1fb78f312440270036d59ce Mon Sep 17 00:00:00 2001 From: Dan Date: Wed, 6 May 2020 15:04:07 -0700 Subject: [PATCH 06/10] consolidate coefficient handling to coeffcients.jl --- src/DiffEqOperators.jl | 1 + src/derivative_operators/coefficients.jl | 32 +++++++++++++++++ src/derivative_operators/convolutions.jl | 34 +++++++++---------- .../derivative_operator.jl | 13 ------- 4 files changed, 50 insertions(+), 30 deletions(-) create mode 100644 src/derivative_operators/coefficients.jl diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 5712c2052..b1d6aa3a9 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -27,6 +27,7 @@ include("derivative_operators/multi_dim_bc_operators.jl") ### Derivative Operators include("derivative_operators/fornberg.jl") +include("derivative_operators/coefficients.jl") include("derivative_operators/derivative_operator.jl") include("derivative_operators/abstract_operator_functions.jl") include("derivative_operators/convolutions.jl") diff --git a/src/derivative_operators/coefficients.jl b/src/derivative_operators/coefficients.jl new file mode 100644 index 000000000..bb49f77b4 --- /dev/null +++ b/src/derivative_operators/coefficients.jl @@ -0,0 +1,32 @@ +""" +``` +init_coefficients(coeff_func, len::Int) +``` + +Return the initial value of an operator's `coefficients` field based on the type +of `coeff_func`. +""" +init_coefficients(coeff_func::Nothing, len::Int) = nothing + +init_coefficients(coeff_func::Number, len::Int) = coeff_func * ones(typeof(coeff_func), len) + +function init_coefficients(coeff_func::AbstractVector{T}, len::Int) where T <: Number + coeff_func +end + +init_coefficients(coeff_func::Function, len::Int) = ones(Float64, len) + + + +""" +``` +get_coefficient(coefficients, index) +``` +""" +get_coefficient(coefficients::AbstractVector, index::Int) = coeff[i] + +# FIXME: I don't think this case is used anymore +get_coefficient(coefficients::Number, index::Int) = coefficients + +# FIXME: Why use "true" here for the value 1? +get_coefficient(coefficients::Nothing, index::Int) = true diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index fe8f17808..dc44f4569 100644 --- a/src/derivative_operators/convolutions.jl +++ b/src/derivative_operators/convolutions.jl @@ -33,7 +33,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count) xtempi = zero(T) cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) for idx in 1:A.stencil_length xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx] end @@ -43,7 +43,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A for i in [(1+A.boundary_point_count):(A.boundary_point_count+offset); (length(x_temp)-A.boundary_point_count-offset+1):(length(x_temp)-A.boundary_point_count)] xtempi = zero(T) cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) for idx in 1:A.stencil_length xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx] end @@ -58,7 +58,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A: coeff = A.coefficients for i in 1 : A.boundary_point_count cur_stencil = stencil[i] - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) xtempi = cur_coeff*stencil[i][1]*x[1] for idx in 2:A.boundary_stencil_length xtempi += cur_coeff * cur_stencil[idx] * x[idx] @@ -73,7 +73,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A coeff = A.coefficients for i in 1 : A.boundary_point_count cur_stencil = stencil[i] - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) xtempi = cur_coeff*stencil[i][end]*x[end] for idx in (A.boundary_stencil_length-1):-1:1 xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx] @@ -93,7 +93,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A:: for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count) xtempi = zero(T) cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) cur_stencil = cur_coeff >= 0 ? cur_stencil : A.derivative_order % 2 == 0 ? reverse(cur_stencil) : -1*reverse(cur_stencil) for idx in 1:A.stencil_length x_idx = cur_coeff < 0 ? x[i - A.stencil_length + 1 + idx] : x[i + idx] @@ -163,7 +163,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A:: coeff = A.coefficients for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.stencil_coefs[1,i-bpc] @@ -190,7 +190,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::D coeff = A.coefficients for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.low_boundary_coefs[1,i] @@ -218,7 +218,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A:: coeff = A.coefficients for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff < 0 xtempi = zero(T) cur_stencil = A.high_boundary_coefs[2,i-len+bpc] @@ -250,7 +250,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector for i in (2+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)-1 xtempi = zero(T) cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i - A.boundary_point_count) @inbounds for idx in 1:A.stencil_length xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1] end @@ -264,7 +264,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector, coeff = A.coefficients for i in 1 : A.boundary_point_count cur_stencil = stencil[i] - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) xtempi = cur_coeff*cur_stencil[1]*_x.l @inbounds for idx in 2:A.boundary_stencil_length xtempi += cur_coeff * cur_stencil[idx] * _x.u[idx-1] @@ -277,7 +277,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector, i = 1 + A.boundary_point_count xtempi = zero(T) cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i - A.boundary_point_count) xtempi = cur_coeff*cur_stencil[1]*_x.l @inbounds for idx in 2:A.stencil_length xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1] @@ -296,7 +296,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector i = length(x_temp)-A.boundary_point_count xtempi = zero(T) cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i - A.boundary_point_count) xtempi = cur_coeff*cur_stencil[end]*_x.r @inbounds for idx in 1:A.stencil_length-1 xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1] @@ -304,7 +304,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector x_temp[i] = xtempi + !overwrite*x_temp[i] for i in 1 : A.boundary_point_count cur_stencil = stencil[i] - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[bc_start + i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, bc_start + i) xtempi = cur_coeff*cur_stencil[end]*_x.r @inbounds for idx in A.stencil_length:-1:1 xtempi += cur_coeff * cur_stencil[end-idx] * _x.u[end-idx+1] @@ -331,7 +331,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count) xtempi = zero(T) cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil - cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true + cur_coeff = get_coefficient(coeff, i) cur_stencil = cur_coeff >= 0 ? cur_stencil : A.derivative_order % 2 == 0 ? reverse(cur_stencil) : -1*reverse(cur_stencil) for idx in 1:A.stencil_length index = cur_coeff < 0 ? i - A.stencil_length + 1 + idx : i + idx @@ -419,7 +419,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.stencil_coefs[1,i-bpc] @@ -450,7 +450,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, _x_len = length(_x) for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.low_boundary_coefs[1,i] @@ -482,7 +482,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, _x_len = length(_x) for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(i) if cur_coeff < 0 xtempi = zero(T) cur_stencil = A.high_boundary_coefs[2,i-len+bpc] diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index e59b51b9c..1fd3fb6db 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,19 +26,6 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end -init_coefficients(coeff_func::Nothing, len::Int) = nothing - -init_coefficients(coeff_func::Number, len::Int) = coeff_func * ones(typeof(coeff_func), len) - -function init_coefficients(coeff_func::AbstractVector{T}, len::Int) where T <: Number - coeff_func -end - -function init_coefficients(coeff_func::Function, len::Int) - ones(Float64, len) -end - - struct CenteredDifference{N} end function CenteredDifference{N}(derivative_order::Int, From 514a99b520c2bf7d26bf2344864332fa5e2ac282 Mon Sep 17 00:00:00 2001 From: Dan Date: Wed, 6 May 2020 16:02:32 -0700 Subject: [PATCH 07/10] forgot two convolution function in last commit --- src/derivative_operators/convolutions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index dc44f4569..288807526 100644 --- a/src/derivative_operators/convolutions.jl +++ b/src/derivative_operators/convolutions.jl @@ -355,7 +355,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, x = _x.u for i in 1:A.boundary_point_count - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) xtempi = 0.0 if cur_coeff >= 0 && i+A.stencil_length <= length(_x) cur_stencil = eltype(upwind_stencils) <: AbstractVector ? upwind_stencils[i] : upwind_stencils @@ -383,7 +383,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, _x_len = length(_x) for i in 1:A.boundary_point_count - cur_coeff = coeff[x_temp_len-A.boundary_point_count+i] + cur_coeff = get_coefficient(coeff, x_temp_len-A.boundary_point_count+i) xtempi = 0.0 if cur_coeff < 0 && _x_len-A.stencil_length - A.boundary_point_count + i >= 1 cur_stencil = eltype(downwind_stencils) <: AbstractVector ? downwind_stencils[i] : downwind_stencils From c0b85cd0677bd14ec0c1a1337f9149e8fcec007d Mon Sep 17 00:00:00 2001 From: Dan Date: Wed, 6 May 2020 16:15:18 -0700 Subject: [PATCH 08/10] use get_coefficient in concretization routines --- src/derivative_operators/concretization.jl | 38 +++++++++++----------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl index 30adbbcd3..6968ab5ec 100644 --- a/src/derivative_operators/concretization.jl +++ b/src/derivative_operators/concretization.jl @@ -13,7 +13,7 @@ function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) wh coeff = A.coefficients get_coeff = if coeff isa AbstractVector - i -> coeff[i] + i = get_coefficient(coeff, i) elseif coeff isa Number i -> coeff else @@ -314,7 +314,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh stencils = A.stencil_coefs for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 cur_stencil = stencils L[i,i+1:i+stl] = cur_coeff*cur_stencil @@ -325,7 +325,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) cur_stencil = stencils cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil) if cur_coeff >= 0 @@ -336,7 +336,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 cur_stencil = stencils cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil) @@ -358,7 +358,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) coeff = A.coefficients for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i] else @@ -367,7 +367,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc] else @@ -376,7 +376,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc] else @@ -403,7 +403,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int= stencils = A.stencil_coefs for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 cur_stencil = stencils L[i,i+1:i+stl] = cur_coeff*cur_stencil @@ -414,7 +414,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int= end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) cur_stencil = stencils cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil) if cur_coeff >= 0 @@ -425,7 +425,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int= end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 cur_stencil = stencils cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil) @@ -447,7 +447,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In coeff = A.coefficients for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i] else @@ -456,7 +456,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc] else @@ -465,7 +465,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc] else @@ -492,7 +492,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A stencils = A.stencil_coefs for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 cur_stencil = stencils L[i,i+1:i+stl] = cur_coeff*cur_stencil @@ -503,7 +503,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) cur_stencil = stencils cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil) if cur_coeff >= 0 @@ -514,7 +514,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 cur_stencil = stencils cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil) @@ -537,7 +537,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int L = BandedMatrix{T}(Zeros(len, len+2), (stl-2, stl)) for i in 1:bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i] else @@ -546,7 +546,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int end for i in bpc+1:len-bpc - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc] else @@ -555,7 +555,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int end for i in len-bpc+1:len - cur_coeff = coeff[i] + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc] else From db02dba12ef1a112caef90c145720d5c66aeae5b Mon Sep 17 00:00:00 2001 From: Dan Date: Fri, 8 May 2020 16:04:11 -0700 Subject: [PATCH 09/10] small fixes. tests passing --- src/derivative_operators/coefficients.jl | 2 +- src/derivative_operators/convolutions.jl | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/derivative_operators/coefficients.jl b/src/derivative_operators/coefficients.jl index bb49f77b4..4b4b9cfa4 100644 --- a/src/derivative_operators/coefficients.jl +++ b/src/derivative_operators/coefficients.jl @@ -23,7 +23,7 @@ init_coefficients(coeff_func::Function, len::Int) = ones(Float64, len) get_coefficient(coefficients, index) ``` """ -get_coefficient(coefficients::AbstractVector, index::Int) = coeff[i] +get_coefficient(coefficients::AbstractVector, index::Int) = coefficients[index] # FIXME: I don't think this case is used anymore get_coefficient(coefficients::Number, index::Int) = coefficients diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index 288807526..51bd26712 100644 --- a/src/derivative_operators/convolutions.jl +++ b/src/derivative_operators/convolutions.jl @@ -163,7 +163,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A:: coeff = A.coefficients for i in bpc+1:len-bpc - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.stencil_coefs[1,i-bpc] @@ -190,7 +190,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::D coeff = A.coefficients for i in 1:bpc - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.low_boundary_coefs[1,i] @@ -218,7 +218,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A:: coeff = A.coefficients for i in len-bpc+1:len - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 xtempi = zero(T) cur_stencil = A.high_boundary_coefs[2,i-len+bpc] @@ -419,7 +419,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, for i in bpc+1:len-bpc - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.stencil_coefs[1,i-bpc] @@ -450,7 +450,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, _x_len = length(_x) for i in 1:bpc - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff >= 0 xtempi = zero(T) cur_stencil = A.low_boundary_coefs[1,i] @@ -482,7 +482,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, _x_len = length(_x) for i in len-bpc+1:len - cur_coeff = get_coefficient(i) + cur_coeff = get_coefficient(coeff, i) if cur_coeff < 0 xtempi = zero(T) cur_stencil = A.high_boundary_coefs[2,i-len+bpc] From edd7abdaf5bcfda63c69a2902930dc648d654180 Mon Sep 17 00:00:00 2001 From: Dan Date: Mon, 14 Dec 2020 23:48:49 +0100 Subject: [PATCH 10/10] small fix in Base.copyto! in concretization.jl (credit simenhu) --- src/derivative_operators/concretization.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl index 32e91c719..3b0a0380e 100644 --- a/src/derivative_operators/concretization.jl +++ b/src/derivative_operators/concretization.jl @@ -12,12 +12,12 @@ function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) wh bstl = A.boundary_stencil_length coeff = A.coefficients - get_coeff = if coeff isa AbstractVector - i = get_coefficient(coeff, i) + get_coeff(i) = if coeff isa AbstractVector + get_coefficient(coeff, i) elseif coeff isa Number - i -> coeff + coeff else - i -> true + true end for i in 1:bl