text:
function CC.abstract_call_gf_by_type(interp::XInterpreter, @nospecialize(f),
                                     fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any}, @nospecialize(atype),
                                     sv::InferenceState, max_methods::Int = InferenceParams(interp).MAX_METHODS)
    ret = @invoke abstract_call_gf_by_type(interp::AbstractInterpreter, @nospecialize(f),
                                           fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any}, @nospecialize(atype),
                                           sv::InferenceState, max_methods::Int)

    return ret
end

function CC.bail_out_toplevel_call(interp::XInterpreter, @nospecialize(sig), sv)
    ret = @invoke bail_out_toplevel_call(interp::AbstractInterpreter, @nospecialize(sig), sv)

    return ret
end

function CC.bail_out_call(interp::XInterpreter, @nospecialize(t), sv)
    ret = @invoke bail_out_call(interp::AbstractInterpreter, @nospecialize(t), sv)

    return ret
end

function CC.add_call_backedges!(interp::XInterpreter,
                                @nospecialize(rettype),
                                edges::Vector{MethodInstance},
                                fullmatch::Vector{Bool}, mts::Vector{Core.MethodTable}, @nospecialize(atype),
                                sv::InferenceState)
    ret = @invoke add_call_backedges!(interp::AbstractInterpreter,
                                      @nospecialize(rettype),
                                      edges::Vector{MethodInstance},
                                      fullmatch::Vector{Bool}, mts::Vector{Core.MethodTable}, @nospecialize(atype),
                                      sv::InferenceState)

    return ret
end

function CC.abstract_call_method_with_const_args(interp::XInterpreter, result::MethodCallResult,
                                                 @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                                 sv::InferenceState, va_override::Bool)
    ret = @invoke abstract_call_method_with_const_args(interp::AbstractInterpreter, result::MethodCallResult,
                                                       @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                                       sv::InferenceState, va_override::Bool)

    return ret
end

function CC.maybe_get_const_prop_profitable(interp::XInterpreter, result::MethodCallResult,
                                         @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                         sv::InferenceState)
    ret = @invoke maybe_get_const_prop_profitable(interp::AbstractInterpreter, result::MethodCallResult,
                                                  @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                                  sv::InferenceState)

    return ret
end

function CC.const_prop_entry_heuristic(interp::XInterpreter, result::MethodCallResult, sv::InferenceState)
    ret = @invoke const_prop_entry_heuristic(interp::AbstractInterpreter, result::MethodCallResult, sv::InferenceState)

    return ret
end

function CC.const_prop_argument_heuristic(interp::XInterpreter, argtypes::Vector{Any})
    ret = @invoke const_prop_argument_heuristic(interp::AbstractInterpreter, argtypes::Vector{Any})

    return ret
end

function CC.const_prop_rettype_heuristic(interp::XInterpreter, @nospecialize(rettype))
    ret = @invoke const_prop_rettype_heuristic(interp::AbstractInterpreter, @nospecialize(rettype))

    return ret
end

function CC.force_const_prop(interp::XInterpreter, @nospecialize(f), method::Method)
    ret = @invoke force_const_prop(interp::AbstractInterpreter, @nospecialize(f), method::Method)

    return ret
end

function CC.const_prop_function_heuristic(interp::XInterpreter, @nospecialize(f), argtypes::Vector{Any}, nargs::Int, allconst::Bool)
    ret = @invoke const_prop_function_heuristic(interp::AbstractInterpreter, @nospecialize(f), argtypes::Vector{Any}, nargs::Int, allconst::Bool)

    return ret
end

function CC.const_prop_methodinstance_heuristic(interp::XInterpreter, method::Method, mi::MethodInstance)
    ret = @invoke const_prop_methodinstance_heuristic(interp::AbstractInterpreter, method::Method, mi::MethodInstance)

    return ret
end

function CC.abstract_call_method(interp::XInterpreter, method::Method, @nospecialize(sig), sparams::SimpleVector, hardlimit::Bool, sv::InferenceState)
    ret = @invoke abstract_call_method(interp::AbstractInterpreter, method::Method, sig, sparams::SimpleVector, hardlimit::Bool, sv::InferenceState)

    return ret
end

function CC.precise_container_type(interp::XInterpreter, @nospecialize(itft), @nospecialize(typ), sv::InferenceState)
    ret = @invoke precise_container_type(interp::AbstractInterpreter, itft, typ, sv::InferenceState)

    return ret
end

function CC.abstract_iteration(interp::XInterpreter, @nospecialize(itft), @nospecialize(itertype), sv::InferenceState)
    ret = @invoke abstract_iteration(interp::AbstractInterpreter, itft, itertype, sv::InferenceState)

    return ret
end

function CC.abstract_apply(interp::XInterpreter, @nospecialize(itft), @nospecialize(aft),    aargtypes::Vector{Any}, sv::InferenceState,
                        max_methods::Int = InferenceParams(interp).MAX_METHODS)
    ret = @invoke abstract_apply(interp::AbstractInterpreter, itft, aft, aargtypes::Vector{Any}, sv::InferenceState,
                                 max_methods::Int)

    return ret
end

function CC.bail_out_apply(interp::XInterpreter, @nospecialize(t), sv)
    ret = @invoke bail_out_apply(interp::AbstractInterpreter, @nospecialize(t), sv)

    return ret
end

function CC.abstract_call_builtin(interp::XInterpreter, f::Builtin, fargs::Union{Nothing,Vector{Any}},
                               argtypes::Vector{Any}, sv::InferenceState, max_methods::Int)
    ret = @invoke abstract_call_builtin(interp::AbstractInterpreter, f::Builtin, fargs::Union{Nothing,Vector{Any}},
                                        argtypes::Vector{Any}, sv::InferenceState, max_methods::Int)

    return ret
end

function CC.abstract_invoke(interp::XInterpreter, argtypes::Vector{Any}, sv::InferenceState)
    ret = @invoke abstract_invoke(interp::AbstractInterpreter, argtypes::Vector{Any}, sv::InferenceState)

    return ret
end

function CC.abstract_call_known(interp::XInterpreter, @nospecialize(f),
                                fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any},
                                sv::InferenceState,
                                max_methods::Int = InferenceParams(interp).MAX_METHODS)
    ret = @invoke abstract_call_known(interp::AbstractInterpreter, f,
                                      fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any},
                                      sv::InferenceState,
                                      max_methods::Int)
end

function CC.abstract_call(interp::XInterpreter, fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any},
                       sv::InferenceState, max_methods::Int = InferenceParams(interp).MAX_METHODS)
    ret = @invoke abstract_call(interp::AbstractInterpreter, fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any},
                                sv::InferenceState, max_methods::Int)

    return ret
end

function CC.abstract_eval_cfunction(interp::XInterpreter, e::Expr, vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_cfunction(interp::AbstractInterpreter, e::Expr, vtypes::VarTable, sv::InferenceState)

    return ret
end

function CC.abstract_eval_value_expr(interp::XInterpreter, e::Expr, vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_value_expr(interp::AbstractInterpreter, e::Expr, vtypes::VarTable, sv::InferenceState)

    return ret
end

function CC.abstract_eval_special_value(interp::XInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_special_value(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    return ret
end

function CC.abstract_eval_value(interp::XInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_value(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    return ret
end

function CC.abstract_eval_statement(interp::XInterpreter, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_statement(interp::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    return ret
end

function CC.typeinf_local(interp::XInterpreter, frame::InferenceState)
    ret = @invoke typeinf_local(interp::AbstractInterpreter, frame::InferenceState)

    return ret
end

function CC.typeinf_nocycle(interp::XInterpreter, frame::InferenceState)
    ret = @invoke typeinf_nocycle(interp::AbstractInterpreter, frame::InferenceState)

    return ret
end


text:
# Dual number with multiple epsilon components.
# Fast inlined storage (will switch to tuples when they're fast)
# This is useful because we can avoid recomputing expensive functions and share computations across the epsilon components.
# Interface will improve, so not exported yet.

immutable Dual4{T<:Real} <: Number
    re::T
    du1::T
    du2::T
    du3::T
    du4::T
end
Dual4(x::Real, y1::Real, y2::Real, y3::Real, y4::Real) = Dual4(promote(x,y1,y2,y3,y4)...)
Dual4(x::Real) = Dual4(x, zero(x), zero(x), zero(x), zero(x))

real(z::Dual4) = z.re
epsilon1(z::Dual4) = z.du1
epsilon2(z::Dual4) = z.du2
epsilon3(z::Dual4) = z.du3
epsilon4(z::Dual4) = z.du4

eps(z::Dual4) = eps(real(z))
eps{T}(::Type{Dual4{T}}) = eps(T)
one(z::Dual4) = Dual4(one(real(z)))
one{T}(::Type{Dual4{T}}) = Dual4(one(T))
@deprecate inf{T}(::Type{Dual4{T}}) convert(Dual4{T}, Inf)
@deprecate nan{T}(::Type{Dual4{T}}) convert(Dual4{T}, NaN)
isnan(z::Dual4) = isnan(real(z))

convert{T<:Real}(::Type{Dual4{T}}, x::Real) =
  Dual4{T}(convert(T, x), zero(T), zero(T), zero(T), zero(T))
convert{T<:Real}(::Type{Dual4{T}}, z::Dual4{T}) = z
convert{T<:Real}(::Type{Dual4{T}}, z::Dual4) =
  Dual4{T}(real(z), epsilon1(z), epsilon2(z), epsilon3(z), epsilon4(z))

convert{T<:Real}(::Type{T}, z::Dual4) =
  (epsilon1(z)==0 && epsilon2(z) == 0 && epsilon3(z) == 0 && epsilon4(z) == 0 ? convert(T, real(z)) : throw(InexactError()))

promote_rule{T<:Real, S<:Real}(::Type{Dual4{T}}, ::Type{Dual4{S}}) =
    Dual4{promote_type(T, S)}
# these promotion rules shouldn't be used for scalar operations -- they're slow
promote_rule{T<:Real}(::Type{Dual4{T}}, ::Type{T}) = Dual4{T}
promote_rule{T<:Real, S<:Real}(::Type{Dual4{T}}, ::Type{S}) =
  Dual4{promote_type(T, S)}

isdual(x::Dual4) = true

real_valued{T<:Real}(z::Dual4{T}) = epsilon1(z) == 0 && epsilon2(z) == 0 && epsilon3(z) == 0 && epsilon4(z) == 0

isfinite(z::Dual4) = isfinite(real(z))

## Generic functions of dual numbers ##

convert(::Type{Dual4}, z::Dual) = z
convert(::Type{Dual4}, x::Real) = Dual(x)

==(z::Dual4, w::Dual4) = real(z) == real(w)
==(z::Dual4, x::Real) = real(z) == x
==(x::Real, z::Dual4) = real(z) == x

isequal(z::Dual4, w::Dual4) =
  isequal(real(z),real(w)) && isequal(epsilon1(z), epsilon1(w)) &&
  isequal(epsilon2(z), epsilon2(w)) && isequal(epsilon3(z),epsilon3(w)) &&
  isequal(epsilon4(z), epsilon4(w))
isequal(z::Dual4, x::Real) = real_valued(z) && isequal(real(z), x)
isequal(x::Real, z::Dual4) = real_valued(z) && isequal(real(z), x)

isless(z::Dual4,w::Dual4) = real(z) < real(w)
isless(z::Real,w::Dual4) = z < real(w)
isless(z::Dual4,w::Real) = real(z) < w

#hash(z::Dual) =
#  (x = hash(real(z)); real_valued(z) ? x : bitmix(x,hash(epsilon(z))))

# we don't support Dual{Complex}, so conj is a noop
conj(z::Dual4) = z
abs(z::Dual4)  = (real(z) >= 0) ? z : -z
abs2(z::Dual4) = z*z

# algebraic definitions
#conjdual(z::Dual) = Dual(real(z),-epsilon(z))
#absdual(z::Dual) = abs(real(z))
#abs2dual(z::Dual) = abs2(real(z))

+(z::Dual4, w::Dual4) = Dual4(real(z)+real(w), epsilon1(z)+epsilon1(w), epsilon2(z) + epsilon2(w), epsilon3(z) + epsilon3(w), epsilon4(z) + epsilon4(w))
+(z::Real, w::Dual4) = Dual4(z+real(w), epsilon1(w), epsilon2(w), epsilon3(w), epsilon4(w))
+(z::Dual4, w::Real) = w+z

-(z::Dual4) = Dual4(-real(z), -epsilon1(z), -epsilon2(z), -epsilon3(z), -epsilon4(z))
-(z::Dual4, w::Dual4) = Dual4(real(z)-real(w), epsilon1(z)-epsilon1(w), epsilon2(z)-epsilon2(w), epsilon3(z)-epsilon3(w), epsilon4(z)-epsilon4(w))
-(z::Real, w::Dual4) = Dual4(z-real(w), -epsilon1(w), -epsilon2(w), -epsilon3(w), -epsilon4(w))
-(z::Dual4, w::Real) = Dual4(real(z)-w, epsilon1(z), epsilon2(z), epsilon3(z), epsilon4(z))

# avoid ambiguous definition with Bool*Number
*(x::Bool, z::Dual4) = ifelse(x, z, ifelse(signbit(real(z))==0, zero(z), -zero(z)))
*(x::Dual4, z::Bool) = z*x

*(z::Dual4, w::Dual4) = Dual4(real(z)*real(w), epsilon1(z)*real(w)+real(z)*epsilon1(w),
epsilon2(z)*real(w)+real(z)*epsilon2(w), epsilon3(z)*real(w)+real(z)*epsilon3(w),
epsilon4(z)*real(w)+real(z)*epsilon4(w))
*(x::Real, z::Dual4) = Dual4(x*real(z), x*epsilon1(z), x*epsilon2(z), x*epsilon3(z), x*epsilon4(z))
*(z::Dual4, x::Real) = x*z

function /(z::Real, w::Dual4)
    df = -z/real(w)^2
    Dual4(z/real(w), epsilon1(w)*df, epsilon2(w)*df, epsilon3(w)*df, epsilon4(w)*df)
end
/(z::Dual4, x::Real) = Dual4(real(z)/x, epsilon1(z)/x, epsilon2(z)/x, epsilon3(z)/x, epsilon4(z)/x)
function /(z::Dual4, w::Dual4)
    denom = 1/(real(w)*real(w))
    Dual4(real(z)/real(w), (epsilon1(z)*real(w)-real(z)*epsilon1(w))*denom,
        (epsilon2(z)*real(w)-real(z)*epsilon2(w))*denom,
        (epsilon3(z)*real(w)-real(z)*epsilon3(w))*denom,
        (epsilon4(z)*real(w)-real(z)*epsilon4(w))*denom)
end

for f in [:^, :(NaNMath.pow)]
    @eval function ($f)(z::Dual4, w::Dual4)
        re = $f(real(z),real(w))
        powval = real(w)*(($f)(real(z),real(w)-1))
        logval = ($f)(real(z),real(w))*log(real(z))
        du1 = epsilon1(z)*powval+epsilon1(w)*logval
        du2 = epsilon2(z)*powval+epsilon2(w)*logval
        du3 = epsilon3(z)*powval+epsilon3(w)*logval
        du4 = epsilon4(z)*powval+epsilon4(w)*logval
        Dual4(re, du1, du2, du3, du4)
    end
end

# these two definitions are needed to fix ambiguity warnings
function ^(z::Dual4, n::Integer)
    zn1 = n*real(z)^(n-1)
    Dual4(real(z)^n, epsilon1(z)*zn1, epsilon2(z)*zn1, epsilon3(z)*zn1, epsilon4(z)*zn1)
end
^(z::Dual4, n::Rational) = invoke(^, Tuple{Dual4,Real}, z,n)

function ^(z::Dual4, n::Real)
    zn1 = n*real(z)^(n-1)
    Dual4(real(z)^n, epsilon1(z)*zn1, epsilon2(z)*zn1, epsilon3(z)*zn1, epsilon4(z)*zn1)
end

function NaNMath.pow(z::Dual4, n::Real)
    powval = n*NaNMath.pow(real(z),n-1)
    Dual4(NaNMath.pow(real(z),n), epsilon1(z)*powval, epsilon2(z)*powval, epsilon3(z)*powval, epsilon4(z)*powval)
end
function NaNMath.pow(z::Real, w::Dual4)
    logval = NaNMath.pow(z,real(w))*log(z)
    Dual4(NaNMath.pow(z,real(w)), epsilon1(w)*logval, epsilon2(w)*logval, epsilon3(w)*logval, epsilon4(w)*logval)
end

for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
    funsym == :exp && continue
    funsym == :abs2 && continue
    @eval function $(funsym)(z::Dual4)
        x = real(z)
        df = $exp
        Dual4($(funsym)(x),epsilon1(z)*df,epsilon2(z)*df,epsilon3(z)*df,epsilon4(z)*df)
    end
    # extend corresponding NaNMath methods
    if funsym in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(funsym))
        @eval function $(funsym)(z::Dual4)
            x = real(z)
            df = $(to_nanmath(exp))
            Dual4($(funsym)(x),epsilon1(z)*df,epsilon2(z)*df,epsilon3(z)*df,epsilon4(z)*df)
        end
    end
end

# only need to compute exp once
function exp(z::Dual4)
    df = exp(real(z))
    return Dual4(df, epsilon1(z)*df, epsilon2(z)*df, epsilon3(z)*df, epsilon4(z)*df)
end

ext: jl


text:
##########################################################################
# General utilities
stringify_method(f::Function) = string(f.env.name)
@compat stringify_method(f::AbstractString) = f

function clean_simulation_params_for_jld(simulation_params)
    cleaned_simulation_params = copy(simulation_params)

    # Precoding method function pointers
    if haskey(simulation_params, "precoding_methods")
        cleaned_precoding_methods = Array(ASCIIString, 0)
        for method in simulation_params["precoding_methods"]
            push!(cleaned_precoding_methods, stringify_method(method))
        end
        cleaned_simulation_params["precoding_methods"] = cleaned_precoding_methods
    end
    if haskey(simulation_params, "precoding_methods_nosweep")
        cleaned_precoding_methods_nosweep = Array(ASCIIString, 0)
        for method in simulation_params["precoding_methods_nosweep"]
            push!(cleaned_precoding_methods_nosweep, stringify_method(method))
        end
        cleaned_simulation_params["precoding_methods_nosweep"] = cleaned_precoding_methods_nosweep
    end

    # Assignment method function pointers
    if haskey(simulation_params, "assignment_methods")
        cleaned_assignment_methods = Array(ASCIIString, 0)
        for method in simulation_params["assignment_methods"]
            push!(cleaned_assignment_methods, stringify_method(method))
        end
        cleaned_simulation_params["assignment_methods"] = cleaned_assignment_methods
    end

    # Independent variable function pointer
    if haskey(simulation_params, "independent_variable")
        cleaned_independent_variable = (string(simulation_params["independent_variable"][1]), simulation_params["independent_variable"][2])
        cleaned_simulation_params["independent_variable"] = cleaned_independent_variable
    end

    # Auxiliary independent variable function pointers
    if haskey(simulation_params, "aux_independent_variables")
        cleaned_aux_independent_variables = Array((@compat Tuple{ASCIIString, Any}), 0)
        for aux_idp in simulation_params["aux_independent_variables"]
            push!(cleaned_aux_independent_variables, (string(aux_idp[1]), aux_idp[2]))
        end
        cleaned_simulation_params["aux_independent_variables"] = cleaned_aux_independent_variables
    end

    return cleaned_simulation_params
end

##########################################################################
# Parameter dict defaultization
macro defaultize_param!(params, key, default_val)
    # Note that we are intentionally violating macro hygiene here...
    return esc(:(haskey($params, $key) || setindex!($params, $default_val, $key)))
end


text:
using Plots
using Reproduce
using JLD2
using FileIO
using Statistics
using ProgressMeter

include("plot_utils.jl")

previous_GKSwstype = get(ENV, "GKSwstype", "")
ENV["GKSwstype"] = "100"

save_file = "plotting/chunlok/generated_plots/state_vistation_heatmap_reproduce.gif"

num_steps = 40000
log_interval = 100

p = plot()

ic = ItemCollection("./OneDTMaze_NEW_DYNAMICS_NEW_EVAL/")
best_ic = ic

data = load_results(best_ic, :oned_tmaze_state_visitation; return_type="array")

state_visitation_data = cat(data..., dims=4)

num_log = size(state_visitation_data)[3]
num_frames = 50
frames_between = floor(Int32, num_log / num_frames)

start_frame = 1
prog = Progress(num_frames, 1)
anim = @animate for i ∈ start_frame:num_frames
    next!(prog)
    init_frame = i * frames_between
    last_frame = min((i+1) * frames_between, num_log)
    frame_info = state_visitation_data[:, :, init_frame:last_frame, 12:12]
    # println(size(frame_info))

    # println(size(sum(frame_info, dims=(3, 4))))

    frame_info = dropdims(sum(frame_info, dims=(3,4)), dims=(3, 4))

    # println(sum(frame_info))
    frame_info = frame_info ./ sum(frame_info)

    # println(size(frame_info))

    frame_info = permutedims(frame_info, [2, 1])
    cgradient = cgrad(:rainbow, scale=:exp)
    heatmap!(frame_info, title="Step: $(i * frames_between * log_interval)", c=cgradient, clim=(0, 1), zscale=:log10)
end

gif(anim, save_file, fps = 2)
ProgressMeter.finish!(prog)

# gif(anim, save_file, fps=15)
# savefig(save_file)

ENV["GKSwstype"] = previous_GKSwstype

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using SafeTestsets

@safetestset "Atomic Solutions                                                                " begin
    include("atomic_solutions_tests.jl")
end
@safetestset "Deterministic Solutions                                                         " begin
    include("deterministic_solutions_tests.jl")
end
@safetestset "IO                                                                              " begin
    include("io_tests.jl")
end


text:
println("===============================================")
println("Testing base.jl")
println("===============================================")

function _trapz(tol)
    x = [i for i in -5.0:0.01:5.0]
    y = 3*x.^2

    isapprox(RHEOS.trapz(y, x), 250.0, atol=tol)
end
@test _trapz(tol)

function _derivCD(tol)
    x = Vector(0.0:0.001:1.5)
    y = x.^2
    dy = 2*x
    dy_numeric = RHEOS.derivCD(y, x)

    isapprox(dy_numeric, dy, atol=tol)
end
@test _derivCD(tol)

function _derivBD(tol)
    x = collect(0.0:0.001:1.5)
    y = x.^2
    dy = 2*x
    dy_numeric = RHEOS.derivBD(y, x)
    # note that the backwards difference method has a weaker
    # numerical test than derivCD as it is less accurate
    all(i -> isapprox(dy[i], dy_numeric[i], atol=tol), eachindex(dy_numeric))
end
@test _derivBD(tol)

function _doublederivCD_quadratic(tol)
    x = Vector(0.0:0.001:1.5)
    y = x.^2
    dy = [2.0 for xi in x]
    dy_numeric = RHEOS.doublederivCD(y, x)

    isapprox(dy_numeric, dy, atol=tol)
end
@test _doublederivCD_quadratic(tol)

function _doublederivCD_cubic(tol)
    x = Vector(0.0:0.001:1.5)
    y = 50*x.^2 .- 2*x.^3 .+ x .- 100.0
    dy = [100.0 - 12*xi for xi in x]
    dy_numeric = RHEOS.doublederivCD(y, x)

    N = length(x)
    test1 = all(i -> isapprox(dy[i], dy_numeric[i], atol=tol), 2:(N-1))
    test2 = isapprox(dy[1], dy_numeric[1], atol=1.1*tol)
    test3 = isapprox(dy[N], dy_numeric[N], atol=1.1*tol)

    test1 && test2 && test3
end
@test _doublederivCD_cubic(tol)

function _doublederivCD_cubic_oscillatory(tol)
    x = Vector(0.0:0.001:20.0)
    y = 50*x.^2 .- 2*x.^3 .+ x .- 100.0 .+ 50.0*sin.(x)
    dy = [100.0 - 12*xi - 50*sin(xi) for xi in x]
    dy_numeric = RHEOS.doublederivCD(y, x)

    N = length(x)
    test1 = all(i -> isapprox(dy[i], dy_numeric[i], atol=tol), 2:(N-1))
    test2 = isapprox(dy[1], dy_numeric[1], atol=10*tol)
    test3 = isapprox(dy[N], dy_numeric[N], atol=10*tol)

    test1 && test2 && test3
end
@test _doublederivCD_cubic_oscillatory(tol)

@test RHEOS.constantcheck([1.0, 2.0, 3.0])==true
@test RHEOS.constantcheck([1.0, 2.0, 4.0])==false

@test RHEOS.closestindex([1.0, 2.0, 3.0], 1.7)==2

@test RHEOS.closestindices([1.0, 2.0, 3.0], [1.7, 3.4])==[2, 3]

@test RHEOS.getsigma(10.0, 1000.0)≈1873.9062512927758

# function _downsample_1region()
#     x0 = collect(0.0:0.01:1.0)
#     y0 = x0.^2
#
#     x1 = collect(0.0:0.1:1.0)
#     y1 = x1.^2
#
#     xout, yout = RHEOS.fixed_resample(x0, y0, [1, length(x0)], [-10])
#
#     xout==x1 && yout==y1
# end
# @test _downsample_1region()
#
# function _upsample_1region()
#     # must be linear as interpolation for upsampling is linear
#     x0 = collect(0.0:0.1:10.0)
#     y0 = 2*x0
#
#     x1 = collect(0.0:0.01:10.0)
#     y1 = 2*x1
#
#     xout, yout = RHEOS.fixed_resample(x0, y0, [1, length(x0)], [10])
#
#     # note that we expect x (generated from range) to be exactly
#     # the same as target higher-sampled data. y will only be the
#     # same within floating point error, hence the approximate
#     # comparison
#     xout==x1 && yout≈y1
# end
# @test _upsample_1region()
#
# function _upanddown_multipleregions()
#     # data to be resampled
#     x0 = collect(-10.0:0.5:10.0)
#     y0 = 2*x0
#
#     # data for comparison composed
#     # of multiple sections
#     x1a = collect(-10.0:0.1:-5.1)
#     x1b = collect(-5.0:1.0:0.0)
#     x1c = collect(0.05:0.05:4.95)
#     x1d = collect(5.0:2.0:10.0)
#     x1e = [x0[end]]
#
#     x1 = vcat(x1a, x1b, x1c, x1d, x1e)
#     y1 = 2*x1
#
#     indices = [1, RHEOS.closestindex(x0, -5.0), RHEOS.closestindex(x0, 0.0), RHEOS.closestindex(x0, 5.0), length(x0)]
#
#     xout, yout = RHEOS.fixed_resample(x0, y0, indices, [5, -2, 10, -4])
#
#     xout==x1 && yout≈y1
# end
# @test _upanddown_multipleregions()

@test RHEOS.singularitytest(x->1/x)
@test RHEOS.singularitytest(x->1/(x-5.0), t1 = 5.0)
@test RHEOS.singularitytest(x->NaN)

function _boltzintegral_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, ramp_loading_derivative)

    all(i -> isapprox(exact_response[i], ramp_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzintegral_nonsing_ramp(tol)

function _boltzintegral_nonsing_ramp_2sections(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, ramp_loading_derivative)

    all(i -> isapprox(exact_response[i], ramp_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzintegral_nonsing_ramp_2sections(tol)

function _boltzintegral_nonsing_linear_3sections(tol)
    t1 = Vector{RheoFloat}(0.0:0.01:9.99)
    t2 = Vector{RheoFloat}(10.0:0.1:16.9)
    t3 = Vector{RheoFloat}(17.0:0.001:20.0)
    t = [t1;t2;t3]
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, ramp_loading_derivative)

    all(i -> isapprox(exact_response[i], ramp_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzintegral_nonsing_linear_3sections(tol)

function _boltzintegral_step(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = exp.(-t)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)
    step_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, step_loading_deriv)

    all(i -> isapprox(exact_response[i], step_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzintegral_step(tol)

function _boltzintegral_linearcombo(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)
    step_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, step_loading_deriv)

    ramp_loading = t
    ramp_loading_deriv = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, ramp_loading_deriv)

    combined_loading = t .+ ones(length(t))
    combined_loading_deriv = RHEOS.derivBD(combined_loading, t)
    combined_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, combined_loading_deriv)

    all(i -> isapprox(combined_response[i], (step_response[i] + ramp_response[i]), atol=tol), eachindex(combined_response))
end
@test _boltzintegral_linearcombo(tol)

function _boltzintegral_nonsing_parabolic(tol)
    # response of Maxwell model to
    # a parabola: 2500 - (t-50)^2
    t = Vector{RheoFloat}(0.0:0.001:20.0)
    exact_response = 102 .- 102*exp.(-t) .- 2t

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzintegral_nonsing(x->exp.(-x), t, loading_derivative)

    # note that tol is 5x higher here (and sample rate is higher)
    # as trapezoidal method is relatively innacurate.
    all(i -> isapprox(exact_response[i], integration_response[i], atol=5*tol), eachindex(exact_response))
end
@test _boltzintegral_nonsing_parabolic(tol)

function _boltzintegral_sing_linear(tol)
    # response of a power-law model
    # to a linear loading: t
    t = Vector{RheoFloat}(0.0:0.01:20.0)
    β = 0.5
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzintegral_sing(x->x.^(-β), t, loading_derivative)

    all(i -> isapprox(exact_response[i], integration_response[i], atol=4*tol), eachindex(exact_response))
end
@test _boltzintegral_sing_linear(tol)

function _boltzintegral_sing_linear_2sections(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    β = 0.5
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzintegral_sing(x->x.^(-β), t, loading_derivative)

    all(i -> isapprox(exact_response[i], integration_response[i], atol=4*tol), eachindex(exact_response))
end
@test _boltzintegral_sing_linear_2sections(tol)

# function _boltzintegral_sing_linear_3sections(tol)
#     # note this test identified that going from higher sample period to lower sample period
#       could cause innacuracies due to singularity approximation
#     t1 = Vector{RheoFloat}(0.0:0.1:4.9)
#     t2 = Vector{RheoFloat}(5.0:0.1:9.9)
#     t3 = Vector{RheoFloat}(10.0:0.01:15.0)
#     t = [t1; t2; t3]
#     β = 0.5
#     exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

#     loading = t
#     loading_derivative = RHEOS.derivBD(loading, t)

#     integration_response = RHEOS.boltzintegral_sing(x->x.^(-β), t, loading_derivative)

#     all(i -> isapprox(exact_response[i], integration_response[i], atol=4*tol), eachindex(exact_response))
# end
# @test _boltzintegral_sing_linear_3sections(tol)

function _boltzintegral_sing_step(tol)
    # response of a power-law model
    # to a step loading
    t = Vector{RheoFloat}(0.0:0.01:20.0)
    β = 0.5
    exact_response = t.^(-0.5)

    loading = ones(length(t))
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzintegral_sing(x->x.^(-β), t, loading_derivative)

    # note that first element is skipped due to singularity
    all(i -> isapprox(exact_response[i], integration_response[i], atol=tol), 2:length(t))
end
@test _boltzintegral_sing_step(tol)

function _boltzintegral_sing_parabolic(tol)
    # response of power-law model to
    # a parabola: 2500 - (t-50)^2
    t = Vector{RheoFloat}(0.0:0.001:20.0)
    β = 0.5
    exact_response = (100/(1-β))*t.^(1-β) .- (2/((1-β)*(2-β)))*t.^(2-β)

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzintegral_sing(x->x.^(-β), t, loading_derivative)
    # note that first element is skipped due to singularity
    # and the higher tolerance and skipping of many elements
    # this is one of the hardest cases for the trapezoidal
    # method of hereditary integration to handle when then there
    # is a singularity.
    all(i -> isapprox(exact_response[i], integration_response[i], atol=1.0), 250:length(t))
end
@test _boltzintegral_sing_parabolic(tol)

function _obj_var_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)

    cost = RHEOS.obj_var_nonsing(nothing, nothing, (x, params)->exp.(-x), t, ramp_loading_derivative, exact_response)

    cost < length(t)*tol^2
end
@test _obj_var_nonsing_ramp(tol)

function _obj_var_nonsing_step(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = exp.(-t)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)

    cost = RHEOS.obj_var_nonsing(nothing, nothing, (x, params)->exp.(-x), t, step_loading_deriv, exact_response)

    cost < length(t)*tol^2
end
@test _obj_var_nonsing_step(tol)

function _obj_var_nonsing_parabolic()
    # response of Maxwell model to
    # a parabola: 2500 - (t-50)^2
    t = Vector{RheoFloat}(0.0:0.001:20.0)
    exact_response = 102 .- 102*exp.(-t) .- 2t

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    cost = RHEOS.obj_var_nonsing(nothing, nothing, (x, params)->exp.(-x), t, loading_derivative, exact_response)
    # note that cost is very high for parabolic as
    # hereditary integral approximation is not
    # good for this case.
    cost < 15*length(t)*tol
end
@test _obj_var_nonsing_parabolic()

function _obj_var_sing_linear(tol)
    # response of a power-law model
    # to a linear loading: t
    t = Vector{RheoFloat}(0.0:0.01:20.0)
    β = 0.5
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    cost = RHEOS.obj_var_sing(nothing, nothing, (x, params)->x.^(-β), t, loading_derivative, exact_response)

    cost < 3*length(t)*tol
end
@test _obj_var_sing_linear(tol)

function _boltzconvolve_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, ramp_loading_derivative)

    all(i -> isapprox(exact_response[i], ramp_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzconvolve_nonsing_ramp(tol)

function _boltzconvolve_step(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = exp.(-t)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)
    step_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, step_loading_deriv)

    all(i -> isapprox(exact_response[i], step_response[i], atol=tol), eachindex(exact_response))
end
@test _boltzconvolve_step(tol)

function _boltzconvolve_linearcombo(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)
    step_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, step_loading_deriv)

    ramp_loading = t
    ramp_loading_deriv = RHEOS.derivBD(ramp_loading, t)
    ramp_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, ramp_loading_deriv)

    combined_loading = t .+ ones(length(t))
    combined_loading_deriv = RHEOS.derivBD(combined_loading, t)
    combined_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, combined_loading_deriv)

    all(i -> isapprox(combined_response[i], (step_response[i] + ramp_response[i]), atol=tol), eachindex(combined_response))
end
@test _boltzconvolve_linearcombo(tol)

function _boltzconvolve_nonsing_parabolic(tol)
    # response of Maxwell model to
    # a parabola: 2500 - (t-50)^2
    dt = 0.001
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 102 .- 102*exp.(-t) .- 2t

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    integration_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, loading_derivative)

    # note that tol is 5x higher here (and sample rate is higher)
    # as trapezoidal method is relatively innacurate.
    all(i -> isapprox(exact_response[i], integration_response[i], atol=5*tol), eachindex(exact_response))
end
@test _boltzconvolve_nonsing_parabolic(tol)

function _boltzconvolve_sing_linear(tol)
    # response of a power-law model
    # to a linear loading:
    dt = 0.001
    t = Vector{RheoFloat}(0.0:dt:20.0)
    β = 0.5
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    t[1] = dt/10.0
    integration_response = RHEOS.boltzconvolve(x->x.^(-β), t, dt, loading_derivative)

    all(i -> isapprox(exact_response[i], integration_response[i], atol=5*tol), eachindex(exact_response))
end
@test _boltzconvolve_sing_linear(tol)

function _boltzconvolve_sing_step(tol)
    # response of a power-law model
    # to a step loading
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    β = 0.5
    exact_response = t.^(-0.5)

    loading = ones(length(t))
    loading_derivative = RHEOS.derivBD(loading, t)

    t[1] = dt/10.0
    integration_response = RHEOS.boltzconvolve(x->x.^(-β), t, dt, loading_derivative)

    # note that first element is skipped due to singularity
    all(i -> isapprox(exact_response[i], integration_response[i], atol=tol), 2:length(t))
end
@test _boltzconvolve_sing_step(tol)

function _boltzconvolve_sing_parabolic(tol)
    # response of power-law model to
    # a parabola: 2500 - (t-50)^2
    dt = 0.001
    t = Vector{RheoFloat}(0.0:dt:20.0)
    β = 0.5
    exact_response = (100/(1-β))*t.^(1-β) .- (2/((1-β)*(2-β)))*t.^(2-β)

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    t[1] = dt/10.0
    integration_response = RHEOS.boltzconvolve(x->x.^(-β), t, dt, loading_derivative)
    # note that first element is skipped due to singularity
    # and the higher tolerance and skipping of many elements
    # this is one of the hardest cases for the trapezoidal
    # method of hereditary integration to handle when then there
    # is a singularity.
    all(i -> isapprox(exact_response[i], integration_response[i], atol=7.0), 250:length(t))
end
@test _boltzconvolve_sing_parabolic(tol)

function _obj_const_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)

    cost = RHEOS.obj_const_nonsing(nothing, nothing, (x, params)->exp.(-x), t, dt, ramp_loading_derivative, exact_response)

    cost < length(t)*tol^2
end
@test _obj_const_nonsing_ramp(tol)

function _obj_const_nonsing_step(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = exp.(-t)
    step_loading = ones(length(t))
    step_loading_deriv = RHEOS.derivBD(step_loading,t)

    cost = RHEOS.obj_const_nonsing(nothing, nothing, (x, params)->exp.(-x), t, dt, step_loading_deriv, exact_response)

    cost < length(t)*tol^2
end
@test _obj_const_nonsing_step(tol)

function _obj_const_nonsing_parabolic()
    # response of Maxwell model to
    # a parabola: 2500 - (t-50)^2
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 102 .- 102*exp.(-t) .- 2t

    loading = 2500.0 .- (t .- 50).^2
    loading_derivative = RHEOS.derivBD(loading, t)

    cost = RHEOS.obj_const_nonsing(nothing, nothing, (x, params)->exp.(-x), t, dt, loading_derivative, exact_response)
    # note that cost is very high for parabolic as
    # hereditary integral approximation is not
    # good for this case.
    cost < 15*length(t)*tol
end
@test _obj_const_nonsing_parabolic()

function _obj_const_sing_linear(tol)
    # response of a power-law model
    # to a linear loading: t
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    β = 0.5
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    t[1] = dt/10.0
    cost = RHEOS.obj_const_sing(nothing, nothing, (x, params)->x.^(-β), t, dt, loading_derivative, exact_response)

    cost < 3*length(t)*tol
end
@test _obj_const_sing_linear(tol)

function _leastsquares_init_const_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    #ramp_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, ramp_loading_derivative)
    modulus = (t, params)->(params[1]*exp.(-t/params[2]))

    init_params = [1.0, 1.0]
    results = RHEOS.leastsquares_init(init_params, [0.90, 0.90], [1.1, 1.1], modulus, t, dt, ramp_loading_derivative, exact_response; constant_sampling = true)

    found_params = results[2]

    isapprox(found_params, init_params, atol=tol)
end
@test _leastsquares_init_const_nonsing_ramp(tol)

function _leastsquares_init_var_nonsing_ramp(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t
    ramp_loading_derivative = RHEOS.derivBD(ramp_loading, t)
    #ramp_response = RHEOS.boltzconvolve(x->exp.(-x), t, dt, ramp_loading_derivative)
    modulus = (t, params)->(params[1]*exp.(-t/params[2]))

    init_params = [1.0, 1.0]
    results = RHEOS.leastsquares_init(init_params, [0.90, 0.90], [1.1, 1.1], modulus, t, dt, ramp_loading_derivative, exact_response; constant_sampling = false)

    found_params = results[2]

    isapprox(found_params, init_params, atol=tol)
end
@test _leastsquares_init_var_nonsing_ramp(tol)

function _leastsquares_init_const_sing_ramp(tol)
    # response of a power-law model
    # to a linear loading:
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    β = 0.5
    exact_response = t.^(1.0 - β) / (1.0 - β)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    modulus = (t, params)->t.^(-params[1])
    init_params = [0.5]
    results = RHEOS.leastsquares_init(init_params, [0.3], [0.7], modulus, t, dt, loading_derivative, exact_response; constant_sampling = true, singularity=true)
    found_params = results[2]

    isapprox(found_params, init_params, atol=tol)
end
@test _leastsquares_init_const_sing_ramp(tol)

function _leastsquares_init_var_sing_ramp(tol)
    # response of a power-law model
    # to a linear loading:
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    β = 0.5
    exact_response = t.^(1.0 - β) / (1.0 - β)

    loading = t
    loading_derivative = RHEOS.derivBD(loading, t)

    modulus = (t, params)->t.^(-params[1])
    init_params = [0.5]
    results = RHEOS.leastsquares_init(init_params, [0.3], [0.7], modulus, t, dt, loading_derivative, exact_response; constant_sampling = false, singularity=true)
    found_params = results[2]

    isapprox(found_params, init_params, atol=tol)
end
@test _leastsquares_init_var_sing_ramp(tol)

function _obj_step_nonsing()
    t = Vector{RheoFloat}(0.0:0.001:20.0)
    params = [2.0, 0.5]
    modulus = (x, p)->p[1]*exp.(-x/p[2])
    loading = 1.0
    expected_response = params[1]*exp.(-t/params[2])
    cost = RHEOS.obj_step_nonsing(params, nothing, modulus, t, loading, expected_response)

    cost < tol
end
@test _obj_step_nonsing()

function _obj_step_weighted()
    t = Vector{RheoFloat}(0.0:0.5:20.0)
    params = [2.0, 0.5]
    modulus = (x, p)->p[1]*x.^(-p[2])
    loading = 1.0
    expected_response = loading*params[1]*t.^(-params[2])
    cost = RHEOS.obj_step_weighted(params, nothing, modulus, t, loading, expected_response[2:end], collect(Integer, 2:length(t)))

    cost < tol
end
@test _obj_step_weighted()

function _leastsquares_stepinit_nonsing_ramp(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    params = [2.0, 0.5]
    loading = 1.0
    expected_response = loading*params[1]*exp.(-t/params[2])

    modulus = (x, p)->p[1]*exp.(-x/p[2])
    init_params = [1.0, 1.0]
    results = RHEOS.leastsquares_stepinit(init_params, [0.70, 0.1], [3.0, 1.5], modulus, t, loading, expected_response)

    found_params = results[2]

    isapprox(params, found_params, atol=tol)
end
@test _leastsquares_stepinit_nonsing_ramp(tol)

function _leastsquares_stepinit_sing_ramp(tol)
    dt = 0.1
    t = Vector{RheoFloat}(0.0:dt:20.0)
    params = [2.0, 0.5]
    loading = 1.0
    expected_response = loading*params[1]*t.^(-params[2])

    modulus = (x, p) -> p[1]*x.^(-p[2])
    init_params = [1.0, 1.0]
    results = RHEOS.leastsquares_stepinit(init_params, [0.70, 0.1], [3.0, 1.5], modulus, t, loading, expected_response; singularity=true)

    found_params = results[2]

    isapprox(params, found_params, atol=tol)
end
@test _leastsquares_stepinit_sing_ramp(tol)


text:
#=
  Bayesian Data Analysis of Titanic data set

  This is also a test of the ARFFFiles package for reading ARFF Files
  (the format used in the data mining system Weka, which happens to
  be one of my favorite data mining system).

  The file used is from my Weka page: http://www.hakank.org/weka/titanic.arff

  This is the simple version of the Titanic dataset with 2201 rows and the following
  attributes, their values (and the transformed values):
  - class: 1st, 2nd, 3rd, crew (=> 1, 2, 3, 4)
  - age  : adult, child        (=> 1, 2)
  - sex  : male, female        (=> 1, 2)
  - survived: yes, no          (=> 1, 0)

Output

 Row │ variable  mean      min    median   max    nmissing  eltype
     │ Symbol    Float64   Int64  Float64  Int64  Int64     DataType
─────┼───────────────────────────────────────────────────────────────
   1 │ class     2.97728       1      3.0      4         0  Int64
   2 │ age       1.95048       1      2.0      2         0  Int64
   3 │ sex       1.21354       1      1.0      2         0  Int64
   4 │ survived  0.323035      0      0.0      1         0  Int64

Mean value for each categories: [0, 1]
 (0, [3.175838926174497 1.9651006711409396 1.0845637583892618])
 (1, [2.561181434599156 1.919831223628692 1.4838255977496484])

Distribution of the categories
0 : 1490 (67.69650159018627%)
1 : 711 (32.30349840981372%)

mean(x): [7.231329977749316e-16 -1.0330471396784736e-16 1.9369633868971382e-17]
std(x): [0.9999999999999998 0.9999999999999993 0.9999999999999996]

Model:
Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           a    0.3229    0.0088     0.0003    0.0004   1294.2108    1.0011      338.9761
   params[1]   -0.0530    0.0093     0.0003    0.0002   1056.5414    1.0010      276.7264
   params[2]   -0.0208    0.0089     0.0003    0.0002   1789.4310    1.0019      468.6828
   params[3]    0.1908    0.0095     0.0003    0.0003   1073.7710    1.0016      281.2391
       sigma    0.4133    0.0062     0.0002    0.0002   1006.8319    0.9996      263.7066

Predictions:
Num correct: 1708.0
Num incorrect: 493.0
Accuracy: 0.7760109041344844
Confusion matrix:
2×2 Matrix{Float64}:
 1364.0  126.0
  367.0  344.0

=#
using Turing, StatsPlots, Distributions
using DataFrames, StatsBase
using CSV
using ARFFFiles
include("jl_utils.jl")

#
# Read dataset and preprocessing
#
df = DataFrame(ARFFFiles.load(download("http://www.hakank.org/weka/titanic.arff")))
display(describe(df))
println()

# Convert categories -> numeric values
transform!(df, [:class]    => ByRow(i -> Dict("1st"=>1,"2nd"=>2,"3rd"=>3,"crew"=>4)[i]) => :class)
transform!(df, [:age]      => ByRow(i -> Dict("child"=>1,"adult"=>2)[i]) => :age)
transform!(df, [:sex]      => ByRow(i -> Dict("male"=>1,"female"=>2)[i]) => :sex)
transform!(df, [:survived] => ByRow(i -> Dict("no"=>0,"yes"=>1)[i]) => :survived)

display(describe(df))
println()

# Convert to matrix and extract x and y.
num_rows,num_cols = size(df)
println("Num rows: $num_rows Num columns: $num_cols")
data = Matrix(df)
x = data[:,1:end-1]
y = data[:,end]

categories = unique(y) |> sort
println("\nMean value for each categories: ", categories)
display(categories .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

println("Distribution of the categories:")
for (k,v) in make_hash(y)
  println("$k : $v ($(100*v/num_rows)%)")
end

# Standardize the dataset (better and faster)
means = mean(x,dims=1)
stds = std(x,dims=1)
x = (x .- means) ./ stds

println("\nmean(x): ", mean(x,dims=1))
println("std(x): ", std(x,dims=1))
println()

# This is defined in jl_utils.jl
model = linear_regression_model(x,y)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

chns_orig = chns
display(chns)

# display(plot(chns))
println("\nPredictions:")

mpred = linear_regression_model(x,missing)
pred = predict(mpred,chns)
# display(pred)

cats = Dict(0=>1,1=>2)
confusion, bads = confusion_matrix(y,pred,cats,false)
println("\nConfusion matrix:")
display(confusion)
# println("\nBad classfied:")
# display(bads)
println("\nnum bad: $(length(bads))")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

function configure(version; experimental::Bool=false)
    name = "LLVMLibUnwind"

    hash = Dict(
        v"11.0.0" => "8455011c33b14abfe57b2fd9803fb610316b16d4c9818bec552287e2ba68922f",
        v"11.0.1" => "6db3b173d872911c9ce1f2779ea4463b3b7e582b4e5fda9d3a005c1ed5ec517f",
    )

    # Collection of sources required to complete build
    sources = [
        ArchiveSource("https://github.com/llvm/llvm-project/releases/download/llvmorg-$(version)/libunwind-$(version).src.tar.xz", hash[version]),
        DirectorySource("./bundled"; follow_symlinks=true),
    ]

    # Bash recipe for building across all platforms
    script = raw"""
cd $WORKSPACE/srcdir/libunwind*

CMAKE_FLAGS=()
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=$prefix)
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=MinSizeRel)  # Note: Disables debug messages
CMAKE_FLAGS+=(-DLIBUNWIND_INCLUDE_DOCS=OFF)
CMAKE_FLAGS+=(-DLIBUNWIND_ENABLE_PEDANTIC=OFF)

if [[ ${target} == x86_64-w64-mingw32 ]]; then
    # Support for threading requires Windows Vista.
    export CXXFLAGS="-D_WIN32_WINNT=0x0600"
fi

# Apply all our patches
if [ -d $WORKSPACE/srcdir/patches ]; then
    for f in $WORKSPACE/srcdir/patches/*.patch; do
        echo "Applying patch ${f}"
        atomic_patch -p2 ${f}
    done
fi

mkdir build && cd build
cmake "${CMAKE_FLAGS[@]}" ..
make -j${nprocs}
make install

# Install header files. Required to access patched in functions
cp -aR ../include ${prefix}/

# Move over the DLL. TODO: There may be a CMAKE flag for this.
if [[ ${target} == *mingw32* ]]; then
    mkdir -p "${libdir}"
    mv -v lib/libunwind.dll "${libdir}"
fi
"""

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = supported_platforms(; experimental=experimental)

    # The products that we will ensure are always built
    products = [
        LibraryProduct("libunwind", :libunwind),
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = Dependency[
    ]

    return name, version, sources, script, platforms, products, dependencies
end

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
import DataFrames
import DataArrays: isna

function test_dataframes()
    df = connect(Postgres, "localhost", "postgres", "", "julia_test") do conn
        stmt = prepare(conn, "SELECT 4::integer as foo, 4.0::DOUBLE PRECISION as bar, " *
            "NULL::integer as foobar;")
        result = execute(stmt)
        return fetchdf(result)
    end

    @test isa(df, DataFrames.DataFrame)
    @test df[:foo] == Int32[4]
    @test df[:bar] == Float64[4.0]
    @test isna(df[:foobar][1])

    df2 = connect(Postgres, "localhost", "postgres", "", "julia_test") do conn
        run(conn, """CREATE TEMPORARY TABLE dftable (foo integer, bar double precision,
            foobar integer);""")
        testdberror(conn, PostgreSQL.CONNECTION_OK)
        # stmt = prepare(conn, "INSERT INTO dftable (foo, bar, foobar) VALUES (\$1, \$2, \$3);")
        run(conn, "INSERT INTO dftable (foo, bar, foobar) VALUES (4, 4.0, NULL);")
        testdberror(conn, PostgreSQL.CONNECTION_OK)
        # executemany(stmt, Vector[{4, 4.0, 4}])
        # testdberror(stmt, PostgreSQL.PGRES_COMMAND_OK)
        stmt = prepare(conn, "SELECT * FROM dftable;")
        result = execute(stmt)
        return fetchdf(result)
    end

    @test names(df) == names(df2)
    for col in names(df)
        @test all(isna(df[col]) == isna(df2[col]))
        @test all(i -> (isna(df[col][i]) && isna(df2[col][i])) || (df[col][i] == df2[col][i]),
            1:length(df[col]))
    end
end

test_dataframes()

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using HDF5

Base.@kwdef mutable struct ObsInfo <: Info
    obsid::String = "unknown:unknown:YYmmddTHHMMSSZ"
    freq_array::AbstractVector{Float64} = Vector{Float64}(undef, 0)
    phase_center_ra::Float64 = NaN
    phase_center_dec::Float64 = NaN
    instrument_name::String = "Unknown"
end # mutable struct ObsInfo

function to_hdf5(group::HDF5.Group, info::ObsInfo)
    invoke(to_hdf5, Tuple{HDF5.Group, Info}, group, info;
           obsid=(layout=HDF5.H5D_COMPACT,),
           phase_center_ra=(layout=HDF5.H5D_COMPACT,),
           phase_center_dec=(layout=HDF5.H5D_COMPACT,),
           instrument_name=(layout=HDF5.H5D_COMPACT,)
          )
end


text:
using PyPlot
using HDF5, JLD
using Statistics

using SubmodularMaximization

close("all")

name = "connectivity_example"

#num_agents = 50
num_agents = 10

num_sensors = 10
nominal_area = 2.0

num_hops = 1
num_partitions = 8

sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 2 * sensor_radius

# works out to the maximum radius for all interactions
communication_radius = 3*station_radius

agent_specification = CircleAgentSpecification(sensor_radius, station_radius,
                                               num_sensors)

f(x) = mean_area_coverage(x, 100)

problem = generate_connected_problem(agent_specification, num_agents,
                                     communication_radius=communication_radius,
                                     objective=f)

# Leave an adjacency matrix around for fun
adjacency = make_adjacency_matrix(problem, communication_radius)

function evaluate_solver(make_solver, name)
  println("$name solver running")
  solver = make_solver(problem)

  @time solution = solve_problem(solver, problem, threaded=true)

  figure()
  xlim([0, 1])
  ylim([0, 1])
  colors = generate_colors(problem.partition_matroid)
  visualize_agents(problem.partition_matroid, colors)
  visualize_solution(problem, solution, colors)

  plot_adjacency(problem, communication_radius)

  gca().set_aspect("equal")
  #savefig("$(fig_path)/$(to_file(name)).png", pad_inches=0.00, bbox_inches="tight")

  coverage = solution.value

  title("$name Solver ($coverage)")

  println("Coverage: ", coverage)
  println("Span: ", communication_span(solver))
  println("Messages: ", communication_messages(solver))
  println("Volume: ", communication_volume(solver))
end

function make_hop_solver(p)
  num_agents = length(p.partition_matroid)

  partition_solver = generate_by_global_partition_size(num_agents,
                                                       num_partitions)

  MultiHopSolver(p, solver=partition_solver,
                 communication_range=communication_radius,
                 num_hops=num_hops)
end

function make_sequential_solver(p)
  SequentialCommunicationSolver(p, communication_radius)
end

# Converging auction solver
# (should converge in the same time as the sequential solver or faster)
function make_auction_solver(p)
  AuctionSolver(p, communication_radius)
end

# Converging auction solver
# (should converge in the same time as the sequential solver or faster)
function make_local_auction_solver(p)
  LocalAuctionSolver(p, communication_radius)
end

for (make_solver, name) in [
                            (make_hop_solver, "RSP Hops"),
                            (make_sequential_solver, "Sequential"),
                            (make_auction_solver, "Auction"),
                            (make_local_auction_solver, "Local Auction")
                           ]
  evaluate_solver(make_solver, name)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
###
### Subject
###
"""
    subject(::BIDSDataType)

Returns the subject label. Required for *ALL* BIDS formats.
"""
subject(p::BIDSDataType) = subject(getproperty(p, :subject))

"""
    substr(::BIDSDataType) -> "sub-<label>"

Returns the subject label as a string for use in a BIDS file name.
"""
substr(p::BIDSDataType) = substr(subject(p))
substr(p::AbstractString) = "sub-$(p)"
substr(p::Any) = substr(String(p))

###
### Session
###
"""
    session(::BIDSDataType)
"""
session(p::BIDSDataType) = session(getproperty(p, :session))

"""
    sesstr(::BIDSDataType) -> "ses-<label>"
"""
sesstr(p::BIDSDataType) = sesstr(session(p))
sesstr(p::Nothing) = ""
sesstr(p::AbstractString) = "ses-$(p)"
sesstr(p::Any) = sesstr(String(p))

###
### Task
###
"""
    tasklabel(::BIDSDataType)
"""
tasklabel(p::BIDSDataType) = getproperty(p, :tasklabel)

"""
    taskstr(::BIDSDataType) -> "task-<label>"
"""
taskstr(p::BIDSDataType) = taskstr(tasklabel(p))
taskstr(p::Nothing) = ""
taskstr(p::AbstractString) = "task-$(p)"
taskstr(p::Any) = acqstr(String(p))

###
### Acquisition
###
"""
    acquisitiontype(::BIDSDataType)
"""
acquisitiontype(p::BIDSDataType) = getproperty(p, :acquisitiontype)

"""
    acqstr(::BIDSDataType) -> "acq-<label>"
"""
acqstr(p::BIDSDataType) = acqstr(acquisitiontype(p))
acqstr(p::AbstractString) = "acq-$(p)"
acqstr(p::Any) = acqstr(String(p))

###
### Contrast enhancing agent
###
"""
    contrastagent(::BIDSDataType)
"""
contrastagent(p::BIDSDataType) = getproperty(p, :contrastagent)

"""
cestr
"""
cestr(p::BIDSDataType) = cestr(contrastagent(p))
cestr(p::Nothing) = ""
cestr(p::AbstractString) = "ce-$p"
cestr(p::Any) = cestr(String(p))

###
### Reconstruction
###
"""
    reconstruction(::BIDSDataType)
"""
reconstruction(p::BIDSDataType) = getproperty(p, :reconstruction)

"""
    recstr(::BIDSDataType) -> \"rec-<label>\"
"""
recstr(p::BIDSDataType) = cestr(reconstruction(p))
recstr(p::Nothing) = ""
recstr(p::AbstractString) = "rec-$p"
recstr(p::Any) = cestr(String(p))

###
### Phase direction
###
"""
    phasedirection
"""
phasedirection(p::BIDSDataType) = getproperty(p, :phasedirection)

"pdstr(x) -> \"dir-<label>\""
pdstr(p::BIDSDataType) = pdstr(phasedirection(p))
pdstr(p::Nothing) = ""
pdstr(p::AbstractString) = "dir-$p"
pdstr(p::Any) = pdstr(String(p))

###
### Run index
###
"runindex(x)"
runindex(x::BIDSDataType) = getproperty(p, :runindex)

"ristr(x) -> \"run-<index>\""
ristr(x::BIDSDataType) = ristr(runindex(x))
ristr(x::Nothing) = ""
ristr(x::AbstractString) = "run-$x"
ristr(x::Any) = ristr(String(x))

###
### Corresponding modality
###
# TODO I do NOT like this name
"""
    corrmodality(x)

Returns the modality associated with the file. Note that this is _not_ the
file's modality. For example, a brain mask for fMRI data would not have the
same "modality" as "corrmodality".
"""
corrmodality(x::BIDSDataType) = getproperty(x, :corrmodality)

"cmstr(x) -> \"mod-<label>\""
cmstr(x::BIDSDataType) = cmstr(corrmodality(x))
cmstr(x::Nothing) = ""
cmstr(x::AbstractString) = "mod-$x"
cmstr(x::Any) = cmstr(String(x))

###
### Echo index
###
"""
    echoindex(x)
"""
echoindex(x::BIDSDataType) = getproperty(x, :echoindex)

"echostr(x) -> \"echo-<label>\""
echostr(x::BIDSDataType) = echostr(echoindex(x))
echostr(x::Nothing) = ""
echostr(x::AbstractString) = "echo-$x"
echostr(x::Any) = echostr(String(x))

###
### Recording label
###
"""
    recordinglabel(x)
"""
recordinglabel(x::BIDSDataType) = getproperty(x, :recordinglabel)

"echostr(x) -> \"echo-<label>\""
recordingstr(x::BIDSDataType) = recordingstr(recordinglabel(x))
recordingstr(x::Nothing) = ""
recordingstr(x::AbstractString) = "recording-$x"
recordingstr(x::Any) = recordingstr(String(x))

###
### Processed (on device)
###
"""
"""
isprocessed(x::BIDSDataType)

###
### Space
###


text:
module ClassicalLaminateTheory

# package code goes here
using Compat

include("clt.jl")

if !isdefined(Main, :JULIA_SVG_BROWSER)
  JULIA_SVG_BROWSER = ""
  try
    JULIA_SVG_BROWSER = ENV["JULIA_SVG_BROWSER"]
  catch e
    println("Environment variable JULIA_SVG_BROWSER not found.")
    JULIA_SVG_BROWSER = ""
  end
end

export
  # From clt.jl
  qmat,
  qbarmat,
  createLaminate!,

  # From ClassicalLaminateTheory.jl
  JULIA_SVG_BROWSER

end # module


text:
T = Float64

## Pendulum
xf = Problems.pendulum.xf
x0 = Problems.pendulum.x0
N = Problems.pendulum.N
n = Problems.pendulum.model.n; m = Problems.pendulum.model.m

# options
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,constraint_tolerance=1.0e-5,
    cost_tolerance=1.0e-5,cost_tolerance_intermediate=1.0e-5,opts_uncon=opts_ilqr,iterations=30,
    penalty_scaling=10.0)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,opts_al=opts_al,R_inf=1.0,resolve_feasible_problem=false)
opts_altro_resolve = ALTROSolverOptions{T}(verbose=verbose,opts_al=opts_al,R_inf=1.0,resolve_feasible_problem=true)

X0 = line_trajectory(x0,xf,N)

# unconstrained infeasible solve
prob = update_problem(copy(Problems.pendulum))
copyto!(prob.X,X0)
solve!(prob,opts_altro)
@test norm(prob.X[end] - xf) < 1.0e-3

prob_resolve = update_problem(copy(Problems.pendulum))
copyto!(prob_resolve.X,X0)
solve!(prob_resolve,opts_altro_resolve)
@test norm(prob_resolve.X[end] - xf) < 1.0e-3

@test norm(prob.X[end] - prob_resolve.X[end]) < 1.0e-5

# constrained infeasible solve
u_bnd = 3.
bnd = BoundConstraint(n,m,u_min=-u_bnd,u_max=u_bnd,trim=true)
goal = goal_constraint(xf)
constraints = Constraints([bnd],N)

prob = update_problem(copy(Problems.pendulum),constraints=copy(constraints))
prob.constraints[N] += goal
copyto!(prob.X,X0)
solve!(prob,opts_altro)

@test max_violation(prob) < opts_al.constraint_tolerance

prob_resolve = update_problem(copy(Problems.pendulum),constraints=copy(constraints))
prob_resolve.constraints[N] += goal

copyto!(prob_resolve.X,X0)
solve!(prob_resolve,opts_altro_resolve)
@test max_violation(prob_resolve) < opts_al.constraint_tolerance

@test norm(prob.X[end] - prob_resolve.X[end]) < 1.0e-5

# ## Quadrotor in Maze
# verbose=false
#
# opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,iterations=300,live_plotting=:off)
#
# opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,opts_uncon=opts_ilqr,
#     iterations=40,cost_tolerance=1.0e-5,cost_tolerance_intermediate=1.0e-4,constraint_tolerance=1.0e-3,penalty_scaling=10.,penalty_initial=1.)
#
# opts_altro = ALTROSolverOptions{T}(verbose=verbose,resolve_feasible_problem=false,opts_al=opts_al,R_inf=0.001);
# opts_altro_resolve = ALTROSolverOptions{T}(verbose=verbose,resolve_feasible_problem=false,opts_al=opts_al,R_inf=0.001);
#
# prob = copy(Problems.quadrotor_maze)
# solve!(prob,opts_altro)
# @test max_violation(prob) < opts_al.constraint_tolerance
#
# prob_resolve = copy(Problems.quadrotor_maze)
#
# solve!(prob_resolve,opts_altro_resolve)
# @test max_violation(prob_resolve) < opts_al.constraint_tolerance
# @test norm(prob.X[end] - prob_resolve.X[end]) < 1.0e-3

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
"""
    MatDeforModule

Module to support general operations for deformation material models.
"""
module MatDeforModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.DeforModelRedModule: DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D

"""
    MatDefor

Abstract type that represents deformable materials.
"""
abstract type MatDefor end

"""
    strain2x2tto3v!(v::FFltVec, t::FFltMat)

Convert a matrix of 2x2 strain components  into a 3-component vector.
"""
function strain2x2tto3v!(v::FFltVec, t::FFltMat)
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = t[1,2] + t[2,1];
    return v
end

"""
    strain3vto2x2t!(t::FFltMat, v::FFltVec)

Convert a strain 3-vector to a  matrix of 2x2 strain components (symmetric tensor).
"""
function strain3vto2x2t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3]/2.;
    t[2,1] = v[3]/2.;
    return t
end

"""
    strain3x3tto6v!(v::FFltVec, t::FFltMat)

Convert a matrix of 3x3 strain components to a 6-component vector.
"""
function strain3x3tto6v!(v::FFltVec, t::FFltMat)
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = t[3,3];
    v[4] = t[1,2] + t[2,1];
    v[5] = t[1,3] + t[3,1];
    v[6] = t[3,2] + t[2,3];
    return v
end

"""
    strain6vto3x3t!(t::FFltMat, v::FFltVec)

Convert a strain 6-vector to a  matrix of 3x3 strain components (symmetric tensor)..
"""
function strain6vto3x3t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[3,3] = v[3];
    t[1,2] = v[4]/2.;
    t[2,1] = v[4]/2.;
    t[1,3] = v[5]/2.;
    t[3,1] = v[5]/2.;
    t[3,2] = v[6]/2.;
    t[2,3] = v[6]/2.;
    return t
end

"""
    strain3x3tto9v!(v::FFltVec, t::FFltMat)

Convert a matrix of 3x3 strain components to a 9-component vector.

The  strain components are in the order
     ex, ey, ez, gxy, gyx, gyz, gzy, gxz, gzx
"""
function strain3x3tto9v!(v::FFltVec, t::FFltMat)
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = t[3,3];
    v[4] = t[1,2];
    v[5] = t[2,1];
    v[6] = t[2,3];
    v[7] = t[3,2];
    v[8] = t[1,3];
    v[9] = t[3,1];
    return v
end

"""
    strain9vto3x3t!(v::FFltVec, t::FFltMat)

Convert a matrix of 3x3 strain components to a 9-component vector.

The  strain components are in the order
     ex, ey, ez, gxy, gyx, gyz, gzy, gxz, gzx
"""
function strain9vto3x3t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[3,3] = v[3];
    t[1,2] = v[4];
    t[2,1] = v[5];
    t[2,3] = v[6];
    t[3,2] = v[7];
    t[1,3] = v[8];
    t[3,1] = v[9];
    return t
end

"""
    stress2x2to3v!(v::FFltVec, t::FFltMat)

Convert a symmetric matrix of 2x2 stress components to a 3-component vector.
"""
function stress2x2to3v!(v::FFltVec, t::FFltMat)
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = 1.0/2.0*(t[1,2] + t[2,1]);
    return v
end

"""
    stress3vto2x2t!(t::FFltMat, v::FFltVec)

Convert a 3-vector to a  matrix of 2x2 stress components (symmetric tensor).
"""
function stress3vto2x2t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3];
    t[2,1] = v[3];
    return t
end

"""
    stress3vto3x3t!(t::FFltMat, v::FFltVec)

Convert a 3-vector to a matrix of 3x3 stress components (symmetric tensor).
"""
function stress3vto3x3t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3];
    t[2,1] = v[3];
    return t
end

"""
    stress4vto3x3t!(t::FFltMat, v::FFltVec)

Convert a 4-vector to a  matrix of 3x3 stress components (tensor).

Convert a 4-vector to a *symmetric*
matrix of 3x3 stress components (tensor).  This is
conversion routine that would be useful for plane strain or
axially symmetric conditions.
The stress vector components need to be ordered as:
    sigmax, sigmay, tauxy, sigmaz,
which is the ordering used for the plane-strain model reduction.
Therefore, for axially symmetric analysis the components need to be
reordered, as from the constitutive equation they come out
as sigmax, sigmay, sigmaz, tauxy.
"""
function stress4vto3x3t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3];
    t[2,1] = v[3];
    t[3,3] = v[4];
    return t
end

"""
    stress6vto3x3t!(t::FFltMat, v::FFltVec)

Convert a 6-vector to a  matrix of 3x3 stress components (symmetric tensor).
"""
function stress6vto3x3t!(t::FFltMat, v::FFltVec)
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[3,3] = v[3];
    t[1,2] = v[4];
    t[2,1] = v[4];
    t[1,3] = v[5];
    t[3,1] = v[5];
    t[3,2] = v[6];
    t[2,3] = v[6];
    return t
end

"""
    stress3x3tto6v!(v::FFltVec, t::FFltMat)

Convert a matrix of 3x3 stress components to a 6-component vector.
"""
function stress3x3tto6v!(v::FFltVec, t::FFltMat)
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = t[3,3];
    v[4] = 1/2.0*(t[1,2] + t[2,1]);
    v[5] = 1/2.0*(t[1,3] + t[3,1]);
    v[6] = 1/2.0*(t[3,2] + t[2,3]);
    return v
end

"""
    strain9vto6v!(t::FFltVec, v::FFltVec)

Convert a strain 9-vector to a  strain 6-vector components (tensor).
"""
function strain9vto6v!(t::FFltVec, v::FFltVec)
    t[1] = v[1];
    t[2] = v[2];
    t[3] = v[3];
    t[4] = v[4]+v[5];
    t[5] = v[8]+v[9];
    t[6] = v[6]+v[7];
    return t
end

"""
    strain6vto9v!(t::FFltVec, v::FFltVec)

Convert a strain 6-vector to a  strain 9-vector components (tensor).

The  strain components are in the order
     ex, ey, ez, gxy/2, gxy/2, gyz/2, gyz/2, gxz/2, gxz/2
"""
function strain6vto9v!(t::FFltVec, v::FFltVec)
    t[1] = v[1];
    t[2] = v[2];
    t[3] = v[3];
    t[4] = v[4]/2.;
    t[5] = v[4]/2.;
    t[6] = v[6]/2.;
    t[7] = v[6]/2.;
    t[8] = v[5]/2.;
    t[9] = v[5]/2.;
    return t
end

"""
    stress9vto6v!(t::FFltVec, v::FFltVec)

Convert a stress 9-vector (tensor) to a  stress 6-vector components.

The  stress components are in the order
     sigx, sigy, sigz, tauxy, tauxy, tauyz, tauyz, tauxz, tauxz
"""
function stress9vto6v!(t::FFltVec, v::FFltVec)
    t[1] = v[1];
    t[2] = v[2];
    t[3] = v[3];
    t[4] = v[4];
    t[5] = v[8];
    t[6] = v[6];
    return t
end

"""
    stress6vto9v!(t::FFltVec, v::FFltVec)

Convert a stress 6-vector to a  stress 9-vector components (tensor).

The  stress components are in the order
     sigx, sigy, sigz, tauxy, tauxy, tauyz, tauyz, tauxz, tauxz
"""
function stress6vto9v!(t::FFltVec, v::FFltVec)
    t[1] = v[1];
    t[2] = v[2];
    t[3] = v[3];
    t[4] = v[4];
    t[5] = v[4];
    t[6] = v[6];
    t[7] = v[6];
    t[8] = v[5];
    t[9] = v[5];
    return t
end

################################################################################
# 3-D model

"""
    rotstressvec(::Type{DeforModelRed3D},  outstress::FFltVec,
      instress::FFltVec, Rm::FFltMat)

Rotate the stress vector by the supplied rotation matrix.

Calculate the rotation of the stress vector to the
'bar' coordinate system given by the columns of the rotation matrix `Rm`.

`Rm` = columns are components of 'bar' basis vectors on the 'plain'
     basis vectors
"""
function  rotstressvec(::Type{DeforModelRed3D},  outstress::FFltVec,  instress::FFltVec, Rm::AbstractArray{FFlt, 2})
      # # Derivation of the transformation matrix [T]
    # #This is from Barbero''s  book Finite element analysis of composite
    # #materials  using Abaqus.  Note that his matrix "a"  is the transpose of
    # #the FinEALE matrix "Rm".
    # #  We also use the FinEALE numbering of the strains.
    #
    # syms T alpha R
    # syms a a11 a12 a13 a21 a22 a23 a31 a32 a33 real
    # a = [a11,a12,a13;
    #     a21,a22,a23;
    #     a31,a32,a33];
    # a = a';#his matrix "a"  is the transpose of the FinEALE matrix "Rm".
    # # # it can be done in terms of l,m,n's as well
    # # syms a l1 m1 n1 l2 m2 n2 l3 m3 n3
    # # a = [l1,m1,n1;l2,m2,n2;l3,m3,n3]
    # #  We also use the FinEALE numbering of the strains.
    # Numbering =[1,4,5;
    #     4,2,6;
    #     5,6,3];
    # T(1:6,1:6) = 0;
    # for i=1:1:3
    #     for j=1:1:3
    #         #if i==j; alpha = j; else alpha = 9-i-j; end
    #         alpha = Numbering(i,j);
    #         for p=1:1:3
    #             for q=1:1:3
    #                 #   if p==q beta = p; else beta = 9-p-q; end
    #                 beta = Numbering(p,q);
    #                 T(alpha,beta) = 0;
    #                 if alpha<=3 & beta<= 3; T(alpha,beta)=a(i,p)*a(i,p); end
    #                 if alpha> 3 & beta<= 3; T(alpha,beta)=a(i,p)*a(j,p); end
    #                 if alpha<=3 & beta>3; T(alpha,beta)=a(i,q)*a(i,p)+a(i,p)*a(i,q);end
    #                 if alpha>3 & beta>3; T(alpha,beta)=a(i,p)*a(j,q)+a(i,q)*a(j,p);end
    #             end
    #         end
    #     end
    # end
    # T
    # R = eye(6,6); R(4,4)=2; R(5,5)=2; R(6,6)=2; # Reuter matrix
    # Tbar = R*T*R^(-1)
    a11=Rm[1,1]; a12=Rm[1,2]; a13=Rm[1,3];
    a21=Rm[2,1]; a22=Rm[2,2]; a23=Rm[2,3];
    a31=Rm[3,1]; a32=Rm[3,2]; a33=Rm[3,3];

    outstress[1] =  (a11^2)*instress[1] + (a21^2)*instress[2] +
                    (a31^2)*instress[3] + (2*a11*a21)*instress[4] +
                    (2*a11*a31)*instress[5] + (2*a21*a31)*instress[6]
    outstress[2] =  (a12^2)*instress[1] +  (a22^2)*instress[2] +
                    (a32^2)*instress[3] + (2*a12*a22)*instress[4] +
                    (2*a12*a32)*instress[5] + (2*a22*a32)*instress[6]
    outstress[3] =  (a13^2)*instress[1] + (a23^2)*instress[2] +
                    (a33^2)*instress[3] + (2*a13*a23)*instress[4] +
                    (2*a13*a33)*instress[5] + (2*a23*a33)*instress[6]
    outstress[4] =  a11*a12*instress[1] +  a21*a22*instress[2] +
                    a31*a32*instress[3] + (a11*a22 + a12*a21)*instress[4] +
                    (a11*a32 + a12*a31)*instress[5] + (a21*a32 + a22*a31)*instress[6]
    outstress[5] =  (a11*a13)*instress[1] + (a21*a23)*instress[2] +
                    (a31*a33)*instress[3] + (a11*a23 + a13*a21)*instress[4] +
                    (a11*a33 + a13*a31)*instress[5] + (a21*a33 + a23*a31)*instress[6]
    outstress[6] =  (a12*a13)*instress[1] + (a22*a23)*instress[2] +
                    (a32*a33)*instress[3] + (a12*a23 + a13*a22)*instress[4] +
                    (a12*a33 + a13*a32)*instress[5] + (a22*a33 + a23*a32)*instress[6]
    return outstress
end

################################################################################
# 2-D plane strain model

"""
    rotstressvec(::Type{DeforModelRed2DStrain},  outstress::FFltVec,
      instress::FFltVec,  Rm::FFltMat)

Rotate the stress vector by the supplied rotation matrix.

Calculate the rotation of the stress vector to the
'bar' coordinate system given by the columns of the rotation matrix `Rm`.

`Rm` = columns are components of 'bar' basis vectors on the 'plain'
     basis vectors
"""
function  rotstressvec(::Type{DeforModelRed2DStrain},  outstress::FFltVec,  instress::FFltVec,  Rm::AbstractArray{FFlt, 2})
    a11=Rm[1,1]; a12=Rm[1,2]; a13=0.0;
    a21=Rm[2,1]; a22=Rm[2,2]; a23=0.0;
    a31=0.0; a32=0.0; a33=1.0;
    # Note the special arrangement  of the components for plane strain
    outstress[1] =  (a11^2)*instress[1] + (a21^2)*instress[2] + (a31^2)*instress[4] + (2*a11*a21)*instress[3]
    outstress[2] =  (a12^2)*instress[1] +  (a22^2)*instress[2] + (a32^2)*instress[4] + (2*a12*a22)*instress[3]
    outstress[4] =  (a13^2)*instress[1] + (a23^2)*instress[2] + (a33^2)*instress[4] + (2*a13*a23)*instress[3]
    outstress[3] =  a11*a12*instress[1] +  a21*a22*instress[2] + a31*a32*instress[4] + (a11*a22 + a12*a21)*instress[3]
    return outstress
end

################################################################################
# 2-D plane stress model

"""
    rotstressvec(::Type{DeforModelRed2DStress},  outstress::FFltVec,
      instress::FFltVec,  Rm::FFltMat)

Rotate the stress vector by the supplied rotation matrix.

Calculate the rotation of the stress vector to the
'bar' coordinate system given by the columns of the rotation matrix `Rm`.

`Rm` = columns are components of 'bar' basis vectors on the 'plain'
     basis vectors
"""
function  rotstressvec(::Type{DeforModelRed2DStress},  outstress::FFltVec,  instress::FFltVec,  Rm::AbstractArray{FFlt, 2})
    a11=Rm[1,1]; a12=Rm[1,2];
    a21=Rm[2,1]; a22=Rm[2,2];
    outstress[1] =  (a11^2)*instress[1] + (a21^2)*instress[2] + (2*a11*a21)*instress[3]
    outstress[2] =  (a12^2)*instress[1] +  (a22^2)*instress[2] + (2*a12*a22)*instress[3]
    outstress[3] =  (a11*a12)*instress[1] + (a21*a22)*instress[2] + (a11*a22+a12*a21)*instress[3]
    return outstress
end

################################################################################
# 2-D axially symmetric stress model

"""
    rotstressvec(::Type{DeforModelRed2DAxisymm},  outstress::FFltVec,
      instress::FFltVec,  Rm::FFltMat)

Rotate the stress vector by the supplied rotation matrix.

Calculate the rotation of the stress vector to the
'bar' coordinate system given by the columns of the rotation matrix `Rm`.

`Rm` = columns are components of 'bar' basis vectors on the 'plain'
     basis vectors
"""
function  rotstressvec(::Type{DeforModelRed2DAxisymm},  outstress::FFltVec,  instress::FFltVec,  Rm::AbstractArray{FFlt, 2})
    a11=Rm[1,1]; a12=Rm[1,2];
    a21=Rm[2,1]; a22=Rm[2,2];
    outstress[1] =  (a11^2)*instress[1] + (a21^2)*instress[2] + (0.0)*instress[3] + (2*a11*a21)*instress[4]
    outstress[2] =  (a12^2)*instress[1] + (a22^2)*instress[2] + (0.0)*instress[3] + (2*a12*a22)*instress[4]
    outstress[3] =  (0.0)*instress[1] + (0.0)*instress[2] + (1.0)*instress[3] + (0.0)*instress[4]
    outstress[4] =  (a11*a12)*instress[1] + (a21*a22)*instress[2] + (0.0)*instress[3] + (a11*a22 + a12*a21)*instress[4]
    return outstress
end

################################################################################
# 1-D stress model

"""
    rotstressvec(::Type{DeforModelRed2DAxisymm},  outstress::FFltVec,
      instress::FFltVec,  Rm::FFltMat)

Rotate the stress vector by the supplied rotation matrix.

Calculate the rotation of the stress vector to the
'bar' coordinate system given by the columns of the rotation matrix `Rm`.

`Rm` = columns are components of 'bar' basis vectors on the 'plain'
     basis vectors
"""
function  rotstressvec(::Type{DeforModelRed1D},  outstress::FFltVec,  instress::FFltVec,  Rm::AbstractArray{FFlt, 2})
    copyto!(outstress, instress)
    return outstress
end

end

# function stressvectorrotation{MR<:DeforModelRed2DStress}(::Type{MR},
#   Rm::FFltMat)
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # switch self.reduction
#     #     case {'axisymm','strain'}
#     #         T =[[   a11^2,   a21^2, 0,         2*a11*a21]
#     #             [   a12^2,   a22^2, 0,         2*a12*a22]
#     #             [       0,       0, 1,                 0]
#     #             [ a11*a12, a21*a22, 0, a11*a22 + a12*a21]];
#     #             case 'stress'
#     T =[[   a11^2   a21^2       2*a11*a21]
#         [   a12^2   a22^2       2*a12*a22]
#         [ a11*a12 a21*a22 a11*a22+a12*a21]];
#     return T
# end
#
# function strainvectorrotation{MR<:DeforModelRed2DStress}(::Type{MR},
#   Rm::FFltMat)
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # switch self.reduction
#     #     case {'axisymm','strain'}
#     #         Tbar =[ [     a11^2,     a21^2, 0,           a11*a21]
#     #                 [     a12^2,     a22^2, 0,           a12*a22]
#     #                 [         0,         0, 1,                 0]
#     #                 [ 2*a11*a12, 2*a21*a22, 0, a11*a22 + a12*a21]];
#     #   case 'stress'
#     Tbar =[ [     a11^2,     a21^2,           a11*a21]
#            [     a12^2,     a22^2,           a12*a22]
#            [ 2*a11*a12, 2*a21*a22, a11*a22 + a12*a21]];
#     return Tbar
# end
#
#
# ################################################################################
# # 2-D plane axially symmetric model
#
# function stressvectorrotation(::Type{MR},
#   Rm::FFltMat) where {MR<:DeforModelRed2DAxisymm}
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     T =[[   a11^2    a21^2  0          2*a11*a21]
#         [   a12^2    a22^2  0          2*a12*a22]
#         [       0        0  1                  0]
#         [ a11*a12  a21*a22  0  a11*a22 + a12*a21]];
#     return T
# end
#
# function strainvectorrotation(::Type{MR},
#   Rm::FFltMat) where {MR<:DeforModelRed2DAxisymm}
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     Tbar =[ [     a11^2      a21^2  0            a11*a21]
#            [     a12^2      a22^2  0            a12*a22]
#            [         0          0  1                  0]
#            [ 2*a11*a12  2*a21*a22  0  a11*a22 + a12*a21]];
#     return Tbar
# end
#
# ################################################################################
# # Generic versions of rotations of stiffness and compliance matrices
#
# function rotatestiffness!(::Type{MR}, D::FFltMat,
#                                                            Rm::FFltMat) where {MR<:DeforModelRed}
#     # Rotate constitutive stiffness matrix of the material.
#     #
#     #         function D=transform_stiffness(self,D,Rm)
#     #
#     # Rotate constitutive stiffness matrix of the material to the
#     # coordinate system given by the columns of the rotation matrix Rm.
#     T =stressvectorrotation(MR,Rm);
#     D = T*D*T';
#     return D
# end

#
# function rotatecompliance!(::Type{MR}, C::FFltMat,
#   Rm::FFltMat) where {MR<:DeforModelRed}
#     # Rotate constitutive compliance matrix of the material.
#     #
#     #   C = rotate_compliance(self,C,Rm)
#     #
#     # Rotate constitutive compliance matrix of the material to the
#     # coordinate system given by the columns of the rotation matrix Rm.
#     Tbar =strainvectorrotation(MR,Rm);
#     C = Tbar*C*Tbar';
#     return C
# end
# function stressvectorrotation(::Type{DeforModelRed2DStrain}, Rm::FFltMat)
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # sigmax, sigmay, sigmaz, tauxy
#     T =[[   a11^2    a21^2  0          2*a11*a21]
#         [   a12^2    a22^2  0          2*a12*a22]
#         [       0        0  1                  0]
#         [ a11*a12  a21*a22  0  a11*a22 + a12*a21]];
#     return T
# end
#
# function strainvectorrotation(::Type{DeforModelRed2DStrain},  Rm::FFltMat)
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     Tbar =[ [     a11^2      a21^2  0            a11*a21]
#            [     a12^2      a22^2  0            a12*a22]
#            [         0          0  1                  0]
#            [ 2*a11*a12  2*a21*a22  0  a11*a22 + a12*a21]];
#     return Tbar
# end
# function strainvectorrotation(::Type{DeforModelRed3D}, Rm::FFltMat)
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2]; a13=Rm[1,3];
#     a21=Rm[2,1]; a22=Rm[2,2]; a23=Rm[2,3];
#     a31=Rm[3,1]; a32=Rm[3,2]; a33=Rm[3,3];
#     Tbar =[
#            [     a11^2     a21^2     a31^2           a11*a21           a11*a31           a21*a31]
#            [     a12^2     a22^2     a32^2           a12*a22           a12*a32           a22*a32]
#            [     a13^2     a23^2     a33^2           a13*a23           a13*a33           a23*a33]
#            [ 2*a11*a12 2*a21*a22 2*a31*a32 a11*a22 + a12*a21 a11*a32 + a12*a31 a21*a32 + a22*a31]
#            [ 2*a11*a13 2*a21*a23 2*a31*a33 a11*a23 + a13*a21 a11*a33 + a13*a31 a21*a33 + a23*a31]
#            [ 2*a12*a13 2*a22*a23 2*a32*a33 a12*a23 + a13*a22 a12*a33 + a13*a32 a22*a33 + a23*a32]];
#     return Tbar
# end


text:
function optimize!(structure::VerticalStructure, optimizer::AbstractStructuredOptimizer, x₀::AbstractVector)
    # Sanity checks
    supports_structure(optimizer, structure) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    check_loadable(optimizer, structure)
    # Load structure
    load_structure!(optimizer, structure, x₀)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function set_master_optimizer!(structure::VerticalStructure, optimizer)
    set_optimizer(structure.first_stage, optimizer)
    return nothing
end

function set_master_optimizer_attribute!(structure::VerticalStructure, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(backend(structure.first_stage), attr, value)
    return nothing
end

function set_subproblem_optimizer!(structure::VerticalStructure, optimizer)
    set_optimizer!(scenarioproblems(structure), optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(structure::VerticalStructure, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(scenarioproblems(structure), attr, value)
    return nothing
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

using Documenter, AbaqusReader

makedocs(modules=[AbaqusReader],
         format = :html,
         checkdocs = :all,
         sitename = "AbaqusReader.jl",
         pages = ["index.md"]
        )


text:
# Julia wrapper for header: portaudio.h
# Automatically generated using Clang.jl

function Pa_GetVersion()
    ccall((:Pa_GetVersion, libportaudio), Cint, ())
end

function Pa_GetVersionText()
    ccall((:Pa_GetVersionText, libportaudio), Cstring, ())
end

function Pa_GetVersionInfo()
    ccall((:Pa_GetVersionInfo, libportaudio), Ptr{PaVersionInfo}, ())
end

function Pa_GetErrorText(errorCode)
    ccall((:Pa_GetErrorText, libportaudio), Cstring, (PaError,), errorCode)
end

function Pa_Initialize()
    ccall((:Pa_Initialize, libportaudio), PaError, ())
end

function Pa_Terminate()
    ccall((:Pa_Terminate, libportaudio), PaError, ())
end

function Pa_GetHostApiCount()
    ccall((:Pa_GetHostApiCount, libportaudio), PaHostApiIndex, ())
end

function Pa_GetDefaultHostApi()
    ccall((:Pa_GetDefaultHostApi, libportaudio), PaHostApiIndex, ())
end

function Pa_GetHostApiInfo(hostApi)
    ccall((:Pa_GetHostApiInfo, libportaudio), Ptr{PaHostApiInfo}, (PaHostApiIndex,), hostApi)
end

function Pa_HostApiTypeIdToHostApiIndex(type)
    ccall((:Pa_HostApiTypeIdToHostApiIndex, libportaudio), PaHostApiIndex, (PaHostApiTypeId,), type)
end

function Pa_HostApiDeviceIndexToDeviceIndex(hostApi, hostApiDeviceIndex)
    ccall((:Pa_HostApiDeviceIndexToDeviceIndex, libportaudio), PaDeviceIndex, (PaHostApiIndex, Cint), hostApi, hostApiDeviceIndex)
end

function Pa_GetLastHostErrorInfo()
    ccall((:Pa_GetLastHostErrorInfo, libportaudio), Ptr{PaHostErrorInfo}, ())
end

function Pa_GetDeviceCount()
    ccall((:Pa_GetDeviceCount, libportaudio), PaDeviceIndex, ())
end

function Pa_GetDefaultInputDevice()
    ccall((:Pa_GetDefaultInputDevice, libportaudio), PaDeviceIndex, ())
end

function Pa_GetDefaultOutputDevice()
    ccall((:Pa_GetDefaultOutputDevice, libportaudio), PaDeviceIndex, ())
end

function Pa_GetDeviceInfo(device)
    ccall((:Pa_GetDeviceInfo, libportaudio), Ptr{PaDeviceInfo}, (PaDeviceIndex,), device)
end

function Pa_IsFormatSupported(inputParameters, outputParameters, sampleRate)
    ccall((:Pa_IsFormatSupported, libportaudio), PaError, (Ptr{PaStreamParameters}, Ptr{PaStreamParameters}, Cdouble), inputParameters, outputParameters, sampleRate)
end

function Pa_OpenStream(stream, inputParameters, outputParameters, sampleRate, framesPerBuffer, streamFlags, streamCallback, userData)
    ccall((:Pa_OpenStream, libportaudio), PaError, (Ptr{Ptr{PaStream}}, Ptr{PaStreamParameters}, Ptr{PaStreamParameters}, Cdouble, Culong, PaStreamFlags, Ptr{PaStreamCallback}, Ptr{Cvoid}), stream, inputParameters, outputParameters, sampleRate, framesPerBuffer, streamFlags, streamCallback, userData)
end

function Pa_OpenDefaultStream(stream, numInputChannels, numOutputChannels, sampleFormat, sampleRate, framesPerBuffer, streamCallback, userData)
    ccall((:Pa_OpenDefaultStream, libportaudio), PaError, (Ptr{Ptr{PaStream}}, Cint, Cint, PaSampleFormat, Cdouble, Culong, Ptr{PaStreamCallback}, Ptr{Cvoid}), stream, numInputChannels, numOutputChannels, sampleFormat, sampleRate, framesPerBuffer, streamCallback, userData)
end

function Pa_CloseStream(stream)
    ccall((:Pa_CloseStream, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_SetStreamFinishedCallback(stream, streamFinishedCallback)
    ccall((:Pa_SetStreamFinishedCallback, libportaudio), PaError, (Ptr{PaStream}, Ptr{PaStreamFinishedCallback}), stream, streamFinishedCallback)
end

function Pa_StartStream(stream)
    ccall((:Pa_StartStream, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_StopStream(stream)
    ccall((:Pa_StopStream, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_AbortStream(stream)
    ccall((:Pa_AbortStream, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_IsStreamStopped(stream)
    ccall((:Pa_IsStreamStopped, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_IsStreamActive(stream)
    ccall((:Pa_IsStreamActive, libportaudio), PaError, (Ptr{PaStream},), stream)
end

function Pa_GetStreamInfo(stream)
    ccall((:Pa_GetStreamInfo, libportaudio), Ptr{PaStreamInfo}, (Ptr{PaStream},), stream)
end

function Pa_GetStreamTime(stream)
    ccall((:Pa_GetStreamTime, libportaudio), PaTime, (Ptr{PaStream},), stream)
end

function Pa_GetStreamCpuLoad(stream)
    ccall((:Pa_GetStreamCpuLoad, libportaudio), Cdouble, (Ptr{PaStream},), stream)
end

function Pa_ReadStream(stream, buffer, frames)
    ccall((:Pa_ReadStream, libportaudio), PaError, (Ptr{PaStream}, Ptr{Cvoid}, Culong), stream, buffer, frames)
end

function Pa_WriteStream(stream, buffer, frames)
    ccall((:Pa_WriteStream, libportaudio), PaError, (Ptr{PaStream}, Ptr{Cvoid}, Culong), stream, buffer, frames)
end

function Pa_GetStreamReadAvailable(stream)
    ccall((:Pa_GetStreamReadAvailable, libportaudio), Clong, (Ptr{PaStream},), stream)
end

function Pa_GetStreamWriteAvailable(stream)
    ccall((:Pa_GetStreamWriteAvailable, libportaudio), Clong, (Ptr{PaStream},), stream)
end

function Pa_GetSampleSize(format)
    ccall((:Pa_GetSampleSize, libportaudio), PaError, (PaSampleFormat,), format)
end

function Pa_Sleep(msec)
    ccall((:Pa_Sleep, libportaudio), Cvoid, (Clong,), msec)
end