Ga/Ja_safe flags and mittleff clean-up

Manifest.toml

@@ -8,9 +8,9 @@ version = "0.5.0"
 
 [[Adapt]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe"
+git-tree-sha1 = "f1b523983a58802c4695851926203b36e28f09db"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.2.0"
+version = "3.3.0"
 
 [[ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -64,10 +64,10 @@ uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
 version = "0.8.7"
 
 [[Colors]]
-deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"]
-git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929"
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "82f4e6ff9f847eca3e5ebc666ea2cd7b48e8b47e"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.6"
+version = "0.12.7"
 
 [[CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -338,6 +338,12 @@ git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 version = "0.4.5"
 
+[[MittagLeffler]]
+deps = ["QuadGK", "SpecialFunctions"]
+git-tree-sha1 = "56200ea7b88fda91c7eb9422c52bce78c2fac7b5"
+uuid = "9c257583-4f8f-53fd-abd9-c69d5080dd54"
+version = "0.2.0"
+
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
@@ -358,9 +364,9 @@ uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
 [[MutableArithmetics]]
 deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "6b6bb8f550dc38310afd4a0af0786dc3222459e2"
+git-tree-sha1 = "ff3aa3e4dbc837f80c2031de2f90125c8b3793f3"
 uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-version = "0.2.14"
+version = "0.2.15"
 
 [[NLopt]]
 deps = ["MathOptInterface", "MathProgBase", "NLopt_jll"]
@@ -492,9 +498,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[StatsBase]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
-git-tree-sha1 = "a83fa3021ac4c5a918582ec4721bc0cf70b495a9"
+git-tree-sha1 = "4bc58880426274277a066de306ef19ecc22a6863"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.33.4"
+version = "0.33.5"
 
 [[TOML]]
 deps = ["Dates"]

Project.toml

@@ -13,6 +13,7 @@ FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 InverseLaplace = "e24f45a5-4790-5a38-8b2c-33d770a4e70f"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+MittagLeffler = "9c257583-4f8f-53fd-abd9-c69d5080dd54"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
@@ -28,7 +29,8 @@ FunctionWrappers = "1.1"
 ImageFiltering = "0.6"
 InverseLaplace = "0.3"
 JLD2 = "0.1"
+MittagLeffler = "0.2"
 NLopt = "0.6"
 QuadGK = "2.3"
 SpecialFunctions = "0.10"
-julia = "1.4"
+julia = "1.6"

src/RHEOS.jl

@@ -4,13 +4,17 @@ __precompile__(true)
 module RHEOS
 
 # installed from Julia package repository
-using InverseLaplace
 using NLopt
 using JLD2
 using DataStructures
 using FunctionWrappers: FunctionWrapper
 using Dierckx
 
+# useful for the various model functions
+using InverseLaplace
+import MittagLeffler: mittleff as mittlefforiginal
+import SpecialFunctions: gamma
+
 import DSP.conv
 
 # Base and stdlib imports
@@ -77,8 +81,8 @@ export AFM, Tweezers
 
 ######################################################
 # bundled dependencies from rheos-cambridge forked repos
-MittLeffLiteDir = joinpath(@__DIR__, "..", "deps", "MittLeffLite", "MittLeffLite.jl")
-include(MittLeffLiteDir)
+#MittLeffLiteDir = joinpath(@__DIR__, "..", "deps", "MittLeffLite", "MittLeffLite.jl")
+#include(MittLeffLiteDir)
 
 include("base.jl")
 include("definitions.jl")

src/base.jl

@@ -199,7 +199,7 @@ Preprocessing base functionality
 Generate sigma/std deviation for gaussian smoothing kernel.
 
 Acts as a low pass filter. Information of time scale τ will be half power,
-faster will be increasingly cut. Called by smoothgauss function.
+faster will be increasingly cut. Called by smooth function.
 """
 function getsigma(τ::Real, samplerate::Real)
 

src/datagen.jl

@@ -90,7 +90,7 @@ time as its only argument. The original data's time signal is used.
 Normally used with a RheoTimeData generated using the `timeline` function.
 """
 function strainfunction(data::RheoTimeData, f::T) where T<:Function
-    log = data.log == nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:strainfunction, params=(f=f,), keywords=()),
+    log = data.log === nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:strainfunction, params=(f=f,), keywords=()),
                                     (comment="strain function applied to timeline",) ) ]
 
     return RheoTimeData(data.σ, convert(Vector{RheoFloat}, map(f, data.t)), data.t, log)
@@ -106,7 +106,7 @@ time as its only argument. The original data's time signal is used.
 Normally used with a `RheoTimeData` generated using the `timeline` function.
 """
 function stressfunction(data::RheoTimeData, f::T) where T<:Function
-    log = data.log == nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:stressfunction, params=(f=f,), keywords=()),
+    log = data.log === nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:stressfunction, params=(f=f,), keywords=()),
                                     (comment="stress function applied to timeline",) ) ]
     return RheoTimeData(convert(Vector{RheoFloat}, map(f, data.t)), data.ϵ, data.t, log)
 end

src/definitions.jl

@@ -73,7 +73,7 @@ end
 function RheoTimeData(;ϵ::Vector{T1} = RheoFloat[], σ::Vector{T2} = RheoFloat[], t::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}
     typecheck = check_time_data_consistency(t,ϵ,σ)
     RheoTimeData(convert(Vector{RheoFloat},σ), convert(Vector{RheoFloat},ϵ), convert(Vector{RheoFloat},t),
-                    log == nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
+                    log === nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
 end
 
 @enum TimeDataType invalid_time_data=-1 time_only=0 strain_only=1 stress_only=2 strain_and_stress=3
@@ -176,7 +176,7 @@ function +(d1::RheoTimeData, d2::RheoTimeData)
     @assert (type1!=time_only) "Addition error: time only data cannot be added"
     @assert (d1.t == d2.t) "Addition error: timelines inconsistent"
 
-    log =   if (d1.log == nothing) || (d2.log == nothing)
+    log =   if (d1.log === nothing) || (d2.log === nothing)
                 nothing
             else
                 [ RheoLogItem( (type = :source, funct = :+, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
@@ -207,7 +207,7 @@ function -(d1::RheoTimeData, d2::RheoTimeData)
     @assert (type1!=time_only) "Subtraction error: time only data cannot be added"
     @assert (d1.t == d2.t) "Subtraction error: timelines inconsistent"
 
-    log =   if (d1.log == nothing) || (d2.log == nothing)
+    log =   if (d1.log === nothing) || (d2.log === nothing)
                 nothing
             else
                 [ RheoLogItem( (type = :source, funct = :-, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
@@ -234,7 +234,7 @@ function -(d::RheoTimeData)
     @assert (type!=invalid_time_data) "unary - error: parameter invalid"
     @assert (type!=time_only) "unary - error: time only data cannot be manipulated this way"
 
-    log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:-, params=(), keywords=() ), () ) ]
+    log = d.log === nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:-, params=(), keywords=() ), () ) ]
 
     return RheoTimeData(-d.σ, -d.ϵ, d.t, log)
 
@@ -258,7 +258,7 @@ function *(operand::Real, d::RheoTimeData)
     @assert (type!=invalid_time_data) "* error: parameter invalid"
     @assert (type!=time_only) "* error: time only data cannot be manipulated this way"
 
-    log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:*, params=(operand = operand), keywords=() ), () ) ]
+    log = d.log === nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:*, params=(operand = operand), keywords=() ), () ) ]
 
     return( RheoTimeData(operand .* d.σ, operand .* d.ϵ, d.t, log) )
 end
@@ -280,7 +280,7 @@ function |(d1::RheoTimeData, d2::RheoTimeData)
     # expected case, type1==stress_only, type2==strain_only
     @assert d1.t==d2.t "Data to be combined must have same time-series"
 
-    log = ((d1.log == nothing) || (d2.log == nothing)) ? nothing : [RheoLogItem((type = :source, funct = :|, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ())]
+    log = ((d1.log === nothing) || (d2.log === nothing)) ? nothing : [RheoLogItem((type = :source, funct = :|, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ())]
 
     return RheoTimeData(d1.σ, d2.ϵ, d1.t, log)
 end
@@ -315,7 +315,7 @@ end
 function RheoFreqData(;Gp::Vector{T1} = RheoFloat[], Gpp::Vector{T2} = RheoFloat[], ω::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}
     typecheck = check_freq_data_consistency(ω,Gp,Gpp)
     RheoFreqData(convert(Vector{RheoFloat},Gp), convert(Vector{RheoFloat},Gpp), convert(Vector{RheoFloat},ω),
-    log == nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
+    log === nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
 end
 
 @enum FreqDataType invalid_freq_data=-1 freq_only=0 with_modulus=1
@@ -378,11 +378,11 @@ function rheologrun(rli::RheoLogItem, d=nothing)
 
       if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
          type=rli.action.type
-         if type==:source && d==nothing
+         if type==:source && d===nothing
             return(eval(rli.action.funct)(rli.action.params...;rli.action.keywords...))
-         elseif type==:process && d!=nothing
+         elseif type==:process && d!==nothing
             return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
-         elseif type==:analysis && d!=nothing
+         elseif type==:analysis && d!==nothing
             return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
          end
       end
@@ -397,7 +397,7 @@ execute all actions from the log. It applies them to the data `d` provided, or u
 function rheologrun(arli::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} = nothing)
 
   # check first item is a source item
-  if d == nothing
+  if d === nothing
         @assert typeof(arli[1].action) <: NamedTuple && :type in keys(arli[1].action) && arli[1].action.type == :source
                 "Source missing in RheoLog"
         end
@@ -406,7 +406,7 @@ function rheologrun(arli::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} =
   for rli in arli
       if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
           type=rli.action.type
-          if type==:analysis && d!=nothing
+          if type==:analysis && d!==nothing
               rheologrun(rli, d)
           else  d=rheologrun(rli, d)
           end
@@ -425,7 +425,7 @@ end
 shows the record of operations on a rheological data.
 """
 function showlog(d::Union{RheoTimeData,RheoFreqData})
-    if d.log != nothing
+    if d.log !== nothing
         for idx in 1:length(d.log)
             println(idx)
             print("     action = "); println(d.log[idx].action)
@@ -487,13 +487,16 @@ function Base.show(io::IO, m::RheoModelClass)
     return
 end
 
-rheoconv(t::Real) = RheoFloat(t)
-rheoconv(t::Array{T,1}) where T<:Real = convert(Vector{RheoFloat},t)
 
 
+rheoconv(t::Real) = RheoFloat(t)
+rheoconv(t::RheoFloat) = t
+rheoconv(t::Vector{T}) where T<:Real = convert(Vector{RheoFloat},t)
+rheoconv(t::Vector{RheoFloat}) = t
+
 
 
-invLaplace(f::Function, t::Array{RheoFloat}) = InverseLaplace.talbotarr(f, t)
+invLaplace(f::Function, t::Vector{RheoFloat}) = InverseLaplace.talbotarr(f, t)
 invLaplace(f::Function, t::RheoFloat) = InverseLaplace.talbot(f, t)
 
 
@@ -723,57 +726,10 @@ end
 
 
 
-#
-# These are internal functions to access the array versions of the relaxation and creep moduli
-# They are used in processing.jl as well as in the user level moduli functions
-#
-
-
-
-
-function _Ga(m::RheoModelClass)
-    if true #m.expressions.Ga_safe
-        m._Ga
-#    else
-#        (ta,p) -> [m._G(t,p) for t in ta]
-    end
-end
-
-function _Ja(m::RheoModelClass)
-    if true #m.expressions.Ja_safe
-        m._Ja
-#    else
-#        (ta,p) -> [m._J(t,p) for t in ta]
-    end
-end
-
-
-
-function _Ga(m::RheoModel)
-    if true #m.expressions.Ga_safe
-        m._Ga
-#    else
-#        ta -> [m._G(t) for t in ta]
-    end
-end
-
-function _Ja(m::RheoModel)
-    if true #m.expressions.Ja_safe
-        m._Ja
-#    else
-#        ta -> [m._J(t) for t in ta]
-    end
-end
-
-
-
-
-
-
 
 
 #
-# These are functions enabling users to access the different moduli functions for each model
+# These are exported functions enabling users to access the different moduli functions for each model
 #
 
 
@@ -813,10 +769,10 @@ function relaxmod(m::RheoModel, t::Number)
 end
 
 function relaxmod(m::RheoModel, ta::Vector{T}) where T <: Number
-    _Ga(m)(ta)
+    m._Ga(ta)
 end
 
-relaxmod(m) = x -> relaxmod(m,x)
+relaxmod(m::RheoModel) = x -> relaxmod(m,x)
 
 
 
@@ -827,7 +783,7 @@ end
 
 
 function relaxmod(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
-    _Ga(m)(ta, params)
+    m._Ga(ta, params)
 end
 
 
@@ -895,10 +851,10 @@ function creepcomp(m::RheoModel, t::Number)
 end
 
 function creepcomp(m::RheoModel, ta::Vector{T}) where T <: Number
-    _Ja(m)(ta)
+    m._Ja(ta)
 end
 
-creepcomp(m) = x -> creepcomp(m,x)
+creepcomp(m::RheoModel) = x -> creepcomp(m,x)
 
 
 
@@ -909,7 +865,7 @@ end
 
 
 function creepcomp(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
-    _Ja(m)(ta, params)
+    m._Ja(ta, params)
 end
 
 
@@ -980,7 +936,7 @@ function storagemod(m::RheoModel, ωa::Vector{T}) where T <: Number
     m._Gpa(ωa)
 end
 
-storagemod(m) = x -> storagemod(m,x)
+storagemod(m::RheoModel) = x -> storagemod(m,x)
 
 
 
@@ -1052,7 +1008,7 @@ function lossmod(m::RheoModel, ωa::Vector{T}) where T <: Number
     m._Gppa(ωa)
 end
 
-lossmod(m) = x -> lossmod(m,x)
+lossmod(m::RheoModel) = x -> lossmod(m,x)
 
 
 
@@ -1127,7 +1083,7 @@ function dynamicmod(m::RheoModel, ωa::Vector{T}) where T <: Number
     m._Gpa(ωa) + m._Gppa(ωa) * im
 end
 
-dynamicmod(m) = x -> dynamicmod(m,x)
+dynamicmod(m::RheoModel) = x -> dynamicmod(m,x)