SimpleTraits.jl: Simple Traits for Julia

SimpleTraits

NEWS

This package provides a macro-based implementation of traits, using Tim Holy's trait trick. The main idea behind traits is to group types outside the type-hierarchy and to make dispatch work with that grouping. The difference to Union-types is that types can be added to a trait after the creation of the trait, whereas Union types are fixed after creation. The cool thing about Tim's trick is that there is no performance impact compared to using ordinary dispatch. For a bit of background and a quick introduction to traits watch my 10min JuliaCon 2015 talk.

One good example of the use of traits is the abstract array interface in Julia-Base. An abstract array either belongs to the Base.IndexLinear or Base.IndexCartesian trait, depending on how its internal indexing works. The advantage to use a trait there is that one is free to create a type hierarchy independent of this particular "trait" of the array(s).

Tim Holy endorses SimpleTraits, a bit: "I'd say that compared to manually writing out the trait-dispatch, the "win" is not enormous, but it is a little nicer." I suspect that — if you don't write Holy-traits before breakfast — your "win" should be greater ;-)

Manual

Traits are defined with @traitdef:

using SimpleTraits
@traitdef IsNice{X}
@traitdef BelongTogether{X,Y} # traits can have several parameters

All traits have one or more (type-)parameters to specify the type to which the trait is applied. For instance IsNice{Int} signifies that Int is a member of IsNice (although whether that is true needs to be checked with the istrait function). Most traits will be one-parameter traits, however, several parameters are useful when there is a "contract" between several types.

As a style convention, I suggest to use trait names which start with a verb, as above two traits. This makes distinguishing between traits and types easier as type names are usually nouns.

Add types to a trait-group with @traitimpl:

@traitimpl IsNice{Int}
@traitimpl BelongTogether{Int,String}

If there is a function which tests whether a trait is fulfilled then it can be used like so:

@traitimpl IsNice{X} <- isnice(X)
isnice(X) = false # set default

i.e. any type X for which isnice(X)==true belongs to IsNice. Notes:

  • overhead-less (static) dispatch is only possible if isnice is pure: "[A pure method] promises that the result will always be the same constant regardless of when the method is called [for the same input arguments]." (ref).
  • Last note that in above example the @traitimpl IsNice{Int} "wins" over the @traitimpl IsNice{X} <- isnice(X), thus this can be used to define exceptions to a rule.

It can be checked whether a type belongs to a trait with istrait:

using Test
@test istrait(IsNice{Int})
@test !istrait(BelongTogether{Int,Int}) # only BelongTogether{Int,String} was added above

Functions which dispatch on traits are constructed like:

@traitfn f(x::X) where {X; IsNice{X}} = "Very nice!"
@traitfn f(x::X) where {X; !IsNice{X}} = "Not so nice!"

This means that a type X which is part of the trait IsNice will dispatch to the method returning "Very nice!", otherwise to the one returning "Not so nice!":

@test f(5)=="Very nice!"
@test f(5.)=="Not so nice!"

Note that calling a trait-function is just like calling any other function. Thus there is no extra mental gymnastics required for a "user" of a trait-based package.

Similarly for BelongTogether which has two parameters:

@traitfn f(x::X,y::Y) where {X,Y; BelongTogether{X,Y}} = "$x and $y forever!"
@test f(5, "b")=="5 and b forever!"
@test_throws MethodError f(5, 5)

@traitfn f(x::X,y::Y) where {X,Y; !BelongTogether{X,Y}} = "$x and $y cannot stand each other!"
@test f(5, 5)=="5 and 5 cannot stand each other!"

Traitor.jl-like syntax

At JuliaCon 2016 folks suggested an alternate, more compact syntax for trait-functions. However, it only works for single parameter traits. SimpleTraits now supports this. Above function f can be written as:

@traitfn ft(x::::IsNice) = "Very nice!"
@traitfn ft(x::::(!IsNice)) = "Not so nice!"

Note that the parenthesis are needed with negated traits, otherwise a parser error is thrown.

Vararg, default argument and keyword argument trait functions

Vararg, default argument and keyword argument trait functions work. However, with keyword arguments the trait function and negated trait function need both have the same keywords (however different values are allowed). Example:

@traitfn kwfn(x::::Tr1, y...; kw=1) = x+y[1]+kw
@traitfn kwfn(x::::(!Tr1), y...; kw=2) = x+y[1]+kw

For default arguments the rule is slightly different: with default arguments the trait function and negated trait function need both have the same default-argument with the same values.

@traitfn deff(x::::Tr1, y=1) = x+y
@traitfn deff(x::::(!Tr1), y=1) = x+y

Method overwritten warnings

Warnings are issued when methods are overwritten. Due to Tim's trick the @traitfn needs to create two functions the first time it is used for a particular method (see next section for an explanation). But when defining the opposite trait, then better only one method is created or else the warning appears. Some heuristics to avoid the warnings are in-place to check whether a method is defined yet or not but they fail at times (see issue #7). Long story short: define the two methods of a trait and its negation using the same argument names and no warning should be issued. Although note that the warnings are harmless.

Details of method dispatch

Defining a trait function adds: one new method (or overwrites one) to the generic function, which contains the logic; and one helper method to do the dispatch (Tim's trick), if it has not been defined before.

When calling a generic function which has some trait-methods, dispatch will first work on the types as normal. If the selected method is a trait-method then trait dispatch will kick in too. Example:

@traitdef Tr{X}

fn(x::Integer) = 1 # a normal method
@traitfn fn(x::X) where {X<:AbstractFloat;  Tr{X}} = 2
@traitfn fn(x::X) where {X<:AbstractFloat; !Tr{X}} = 3

@traitimpl Tr{Float32}
@traitimpl Tr{Int} # this does not impact dispatch of `fn`

fn(5) # -> 1; dispatch only happens on the type
fn(Float32(5)) # -> 2; dispatch through traits
fn(Float64(5)) # -> 3; dispatch through traits

Further note that for a particular trait-method dispatch only works on one trait. Continuing above example, this does not work as one may expect:

@traitdef Tr2{X}
@traitfn fn(x::X) where {X<:AbstractFloat; Tr2{X}} = 4

@traitimpl Tr2{Float16}
fn(Float16(5)) # -> 4; dispatch through traits
fn(Float32(5)) # -> MethodError; method defined in previous example
               #    was overwritten above

This last definition of fn just overwrites the definition @traitfn f(x::X) where {X; Tr{X}} = 2 from above.

If you need to dispatch on several traits in a single trait-method, then you're out of luck. But please voice your grievance over in pull request #2.

Performance

There is no performance impact compared to normal functions thanks to Julia's clever design. Continuing the example from above and looking at the native code

julia> @code_native fn(5)
        .text
Filename: REPL[3]
        pushq   %rbp
        movq    %rsp, %rbp
Source line: 1
        movl    $1, %eax
        popq    %rbp
        retq
        nopl    (%rax,%rax)

julia> @code_native fn(Float16(5))
        .text
Filename: SimpleTraits.jl
        pushq   %rbp
        movq    %rsp, %rbp
Source line: 185
        movl    $4, %eax
        popq    %rbp
        retq
        nopl    (%rax,%rax)

shows that the normal method and the trait-method compile down to the same machine instructions.

However, if the trait-grouping function is not constant or a generated function then dispatch may be dynamic. This can be checked with @check_fast_traitdispatch, which checks whether the number of lines of LLVM code is the same for a trait function than a normal one:

checkfn(x) = rand()>0.5 ? true : false # a bit crazy!
@traitdef TestTr{X}
@traitimpl TestTr{X} <- checkfn(X)
# this tests a trait-function with TestTr{Int}:
@check_fast_traitdispatch TestTr
# this tests a trait-function with TestTr{String} and will
# also prints number of LLCM-IR lines of trait vs normal function:
@check_fast_traitdispatch TestTr String true

# Now this is fast:
@traitimpl TestTr{String}
@check_fast_traitdispatch TestTr String true

Advanced features

The macros of the previous section are the official API of the package and should be reasonably stable. What follows in this section is "under the hood" and may well be updated (but still signalled with minor version changes).

Instead of using @traitimpl to add types to traits, it can be programmed. Running @traitimpl IsNice{Int} essentially expands to

SimpleTraits.trait(::Type{IsNice{X1}}) where {X1 <: Int} = IsNice{X1}

I.e. trait is the identity function for a fulfilled trait and returns Not{TraitInQuestion{...}} otherwise (this is the fall-back for <:Any). So instead of using @traitimpl this can be coded directly. Note that anything but a constant function will probably not be inlined away by the JIT and will lead to slower dynamic dispatch (see @check_fast_traitdispatch for a helper to check).

Example leading to static dispatch:

@traitdef IsBits{X}
SimpleTraits.trait(::Type{IsBits{X1}}) where {X1} = isbits(X1) ? IsBits{X1} : Not{IsBits{X1}}
istrait(IsBits{Int}) # true
istrait(IsBits{Array{Int,1}}) # false
struct A
    a::Int
end
istrait(IsBits{A}) # true

Dynamic dispatch can be avoided using a generated function or pure functions (sometimes they need to be annotated with Base.@pure):

@traitdef IsBits{X}
@generated function SimpleTraits.trait(::Type{IsBits{X1}}) where X1
    isbits(X1) ? :(IsBits{X1}) : :(Not{IsBits{X1}})
end

What is allowed in generated functions is heavily restricted, see Julia manual. In particular, no methods which are defined after the generated function are allowed to be called inside the generated function, otherwise this issue is encountered. Generally, try non-generated functions first and only in a pinch generated functions. See also issue 40.

Note that these programmed-traits can be combined with @traitimpl IsBits{XYZ}, i.e. program the general case and add exceptions with @traitimpl IsBits{XYZ}.

Trait-inheritance can also be hand-coded with above trick. For instance, the trait given by (in pseudo syntax) BeautyAndBeast{X,Y} <: IsNice{X}, !IsNice{Y}, BelongTogether{X,Y}:

@traitdef BeautyAndBeast{X,Y}
function SimpleTraits.trait(::Type{BeautyAndBeast{X,Y}}) where {X,Y}
    if istrait(IsNice{X}) && !istrait(IsNice{Y}) && BelongTogether{X,Y}
        BeautyAndBeast{X,Y}
    else
        Not{BeautyAndBeast{X,Y}}
    end
end

Note that this will lead to slower, dynamic dispatch, as the function is not pure (it depends on the global state of which types belong to the traits IsNice and BelongTogether).

Note also that trait functions can be generated functions:

@traitfn @generated fg(x::X) where {X; IsNice{X}} = (println(x); :x)

Innards

The function macroexpand shows the syntax transformations a macro does. Here the edited output of running it for the macros of this package:

julia> macroexpand(:(@traitdef Tr{X}))

struct Tr{X} <: SimpleTraits.Trait
end

julia> macroexpand(:(@traitimpl Tr{Int}))

# this function does the grouping of types into traits:
SimpleTraits.trait(::Type{Tr{X1}}) where X1 <: Int = Tr{X1}
SimpleTraits.istrait(::Type{Tr{X1}}) where X1 <: Int = true # for convenience, really

julia> macroexpand(:(@traitfn g(x::X) where {X; Tr{X}}= x+1))

@inline g(x::X) where {X} = g(trait(Tr{X}), x) # this is Tim's trick, using above grouping-function
g(::Type{Tr{X}},x::X) where {X} = x + 1 # this is the logic

julia> macroexpand(:(@traitfn g(x::X) where {X; !Tr{X}}= x+1000))

# the trait dispatch helper function needn't be defined twice,
# only the logic:
g(::Type{ Not{Tr{X}} }, x::X) where {X} = x + 1000

For a detailed explanation of how Tim's trick works, see Traits.jl: Dispatch on traits. The difference here is I make the methods containing the logic part of the same generic function (there it's in _f).

Base Traits

I started putting some Julia-Base traits together which can be loaded with using SimpleTraits.BaseTraits, see the source for all definitions.

Example, dispatch on whether an argument is immutable or not:

@traitfn f(x::X) where {X; IsImmutable{X}} = X(x.fld+1) # make a new instance
@traitfn f(x::X) where {X; !IsImmutable{X}} = (x.fld += 1; x) # update in-place

# use
mutable struct A; fld end
struct B; fld end
a=A(1)
f(a) # in-place
@assert a.fld == A(2).fld

b=B(1) # out of place
b2 = f(b)
@assert b==B(1)
@assert b2==B(2)

Background

This package grew out of an attempt to reduce the complexity of Traits.jl, but at the same time staying compatible (but which it isn't). Compared to Traits.jl, it drops support for:

  • Trait definition in terms of methods and constraints. Instead the user needs to assign types to traits manually. This removes the most complex part of Traits.jl: the checking whether a type satisfies a trait definition.
  • trait functions which dispatch on more than one trait. This allows to remove the need for generated functions, as well as removing the rules for trait-dispatch.

The reason for splitting this away from Traits.jl are:

  • creating a more reliable and easier to maintain package than Traits.jl
  • exploring inclusion in Base (see #13222).

My JuliaCon 2015 talk gives a 10 minute introduction to Traits.jl and SimpleTraits.jl.

The Future

The future of traits in Julia-Base: According to Stefan Karpinski's JuliaCon 2016 talk, Julia 1.0, traits are scheduled to land after Julia 1.0. My crystal ball tells me that all or most of the functionality of this package will be supported in the future trait system (multiparameter-traits may not be). Thus I expect the transition will be mostly a matter of a syntax update and less of a semantic update. Also, an advantage to using this package versus hand-coding Holy-traits will be that all occurrences of trait usage are clearly marked and thus easier to update.

The future of this package: I see it as light-weight package focusing on letting functions use dispatch based on traits. This dispatch is currently fairly limited, see section "Gotcha" above, but may be expanded in the future: either through something like in PR m3/multitraits.

In the unlikely event that I find myself with too much time on my hands, I may try to develop a companion package to allow the specification of a trait in terms of interfaces. The combination of the two packages would then have similar functionality to my experimental package Traits.jl. If anyone fancies a go at writing this companion package, I would be very happy to help and contribute. After the type-system overhaul lands, this should be much less hackish than what's in Traits.jl.

References

  • Traits.jl and its references. In particular here is an in-depth discussion on limitations of Holy-Traits, which this package implements.

To ponder

  • There is a big update sitting in the branch m3/multitraits; but I never quite finished it. It would also address the next point:
  • Could type inheritance be used for sub-traits (Jutho's idea)? In particular could it be used in such a way that it is compatible with the multiple inheritance used in Traits.jl?

Download Details:

Author: Mauro3
Source Code: https://github.com/mauro3/SimpleTraits.jl 
License: MIT license

#julia #simple 

What is GEEK

Buddha Community

SimpleTraits.jl: Simple Traits for Julia

SimpleTraits.jl: Simple Traits for Julia

SimpleTraits

NEWS

This package provides a macro-based implementation of traits, using Tim Holy's trait trick. The main idea behind traits is to group types outside the type-hierarchy and to make dispatch work with that grouping. The difference to Union-types is that types can be added to a trait after the creation of the trait, whereas Union types are fixed after creation. The cool thing about Tim's trick is that there is no performance impact compared to using ordinary dispatch. For a bit of background and a quick introduction to traits watch my 10min JuliaCon 2015 talk.

One good example of the use of traits is the abstract array interface in Julia-Base. An abstract array either belongs to the Base.IndexLinear or Base.IndexCartesian trait, depending on how its internal indexing works. The advantage to use a trait there is that one is free to create a type hierarchy independent of this particular "trait" of the array(s).

Tim Holy endorses SimpleTraits, a bit: "I'd say that compared to manually writing out the trait-dispatch, the "win" is not enormous, but it is a little nicer." I suspect that — if you don't write Holy-traits before breakfast — your "win" should be greater ;-)

Manual

Traits are defined with @traitdef:

using SimpleTraits
@traitdef IsNice{X}
@traitdef BelongTogether{X,Y} # traits can have several parameters

All traits have one or more (type-)parameters to specify the type to which the trait is applied. For instance IsNice{Int} signifies that Int is a member of IsNice (although whether that is true needs to be checked with the istrait function). Most traits will be one-parameter traits, however, several parameters are useful when there is a "contract" between several types.

As a style convention, I suggest to use trait names which start with a verb, as above two traits. This makes distinguishing between traits and types easier as type names are usually nouns.

Add types to a trait-group with @traitimpl:

@traitimpl IsNice{Int}
@traitimpl BelongTogether{Int,String}

If there is a function which tests whether a trait is fulfilled then it can be used like so:

@traitimpl IsNice{X} <- isnice(X)
isnice(X) = false # set default

i.e. any type X for which isnice(X)==true belongs to IsNice. Notes:

  • overhead-less (static) dispatch is only possible if isnice is pure: "[A pure method] promises that the result will always be the same constant regardless of when the method is called [for the same input arguments]." (ref).
  • Last note that in above example the @traitimpl IsNice{Int} "wins" over the @traitimpl IsNice{X} <- isnice(X), thus this can be used to define exceptions to a rule.

It can be checked whether a type belongs to a trait with istrait:

using Test
@test istrait(IsNice{Int})
@test !istrait(BelongTogether{Int,Int}) # only BelongTogether{Int,String} was added above

Functions which dispatch on traits are constructed like:

@traitfn f(x::X) where {X; IsNice{X}} = "Very nice!"
@traitfn f(x::X) where {X; !IsNice{X}} = "Not so nice!"

This means that a type X which is part of the trait IsNice will dispatch to the method returning "Very nice!", otherwise to the one returning "Not so nice!":

@test f(5)=="Very nice!"
@test f(5.)=="Not so nice!"

Note that calling a trait-function is just like calling any other function. Thus there is no extra mental gymnastics required for a "user" of a trait-based package.

Similarly for BelongTogether which has two parameters:

@traitfn f(x::X,y::Y) where {X,Y; BelongTogether{X,Y}} = "$x and $y forever!"
@test f(5, "b")=="5 and b forever!"
@test_throws MethodError f(5, 5)

@traitfn f(x::X,y::Y) where {X,Y; !BelongTogether{X,Y}} = "$x and $y cannot stand each other!"
@test f(5, 5)=="5 and 5 cannot stand each other!"

Traitor.jl-like syntax

At JuliaCon 2016 folks suggested an alternate, more compact syntax for trait-functions. However, it only works for single parameter traits. SimpleTraits now supports this. Above function f can be written as:

@traitfn ft(x::::IsNice) = "Very nice!"
@traitfn ft(x::::(!IsNice)) = "Not so nice!"

Note that the parenthesis are needed with negated traits, otherwise a parser error is thrown.

Vararg, default argument and keyword argument trait functions

Vararg, default argument and keyword argument trait functions work. However, with keyword arguments the trait function and negated trait function need both have the same keywords (however different values are allowed). Example:

@traitfn kwfn(x::::Tr1, y...; kw=1) = x+y[1]+kw
@traitfn kwfn(x::::(!Tr1), y...; kw=2) = x+y[1]+kw

For default arguments the rule is slightly different: with default arguments the trait function and negated trait function need both have the same default-argument with the same values.

@traitfn deff(x::::Tr1, y=1) = x+y
@traitfn deff(x::::(!Tr1), y=1) = x+y

Method overwritten warnings

Warnings are issued when methods are overwritten. Due to Tim's trick the @traitfn needs to create two functions the first time it is used for a particular method (see next section for an explanation). But when defining the opposite trait, then better only one method is created or else the warning appears. Some heuristics to avoid the warnings are in-place to check whether a method is defined yet or not but they fail at times (see issue #7). Long story short: define the two methods of a trait and its negation using the same argument names and no warning should be issued. Although note that the warnings are harmless.

Details of method dispatch

Defining a trait function adds: one new method (or overwrites one) to the generic function, which contains the logic; and one helper method to do the dispatch (Tim's trick), if it has not been defined before.

When calling a generic function which has some trait-methods, dispatch will first work on the types as normal. If the selected method is a trait-method then trait dispatch will kick in too. Example:

@traitdef Tr{X}

fn(x::Integer) = 1 # a normal method
@traitfn fn(x::X) where {X<:AbstractFloat;  Tr{X}} = 2
@traitfn fn(x::X) where {X<:AbstractFloat; !Tr{X}} = 3

@traitimpl Tr{Float32}
@traitimpl Tr{Int} # this does not impact dispatch of `fn`

fn(5) # -> 1; dispatch only happens on the type
fn(Float32(5)) # -> 2; dispatch through traits
fn(Float64(5)) # -> 3; dispatch through traits

Further note that for a particular trait-method dispatch only works on one trait. Continuing above example, this does not work as one may expect:

@traitdef Tr2{X}
@traitfn fn(x::X) where {X<:AbstractFloat; Tr2{X}} = 4

@traitimpl Tr2{Float16}
fn(Float16(5)) # -> 4; dispatch through traits
fn(Float32(5)) # -> MethodError; method defined in previous example
               #    was overwritten above

This last definition of fn just overwrites the definition @traitfn f(x::X) where {X; Tr{X}} = 2 from above.

If you need to dispatch on several traits in a single trait-method, then you're out of luck. But please voice your grievance over in pull request #2.

Performance

There is no performance impact compared to normal functions thanks to Julia's clever design. Continuing the example from above and looking at the native code

julia> @code_native fn(5)
        .text
Filename: REPL[3]
        pushq   %rbp
        movq    %rsp, %rbp
Source line: 1
        movl    $1, %eax
        popq    %rbp
        retq
        nopl    (%rax,%rax)

julia> @code_native fn(Float16(5))
        .text
Filename: SimpleTraits.jl
        pushq   %rbp
        movq    %rsp, %rbp
Source line: 185
        movl    $4, %eax
        popq    %rbp
        retq
        nopl    (%rax,%rax)

shows that the normal method and the trait-method compile down to the same machine instructions.

However, if the trait-grouping function is not constant or a generated function then dispatch may be dynamic. This can be checked with @check_fast_traitdispatch, which checks whether the number of lines of LLVM code is the same for a trait function than a normal one:

checkfn(x) = rand()>0.5 ? true : false # a bit crazy!
@traitdef TestTr{X}
@traitimpl TestTr{X} <- checkfn(X)
# this tests a trait-function with TestTr{Int}:
@check_fast_traitdispatch TestTr
# this tests a trait-function with TestTr{String} and will
# also prints number of LLCM-IR lines of trait vs normal function:
@check_fast_traitdispatch TestTr String true

# Now this is fast:
@traitimpl TestTr{String}
@check_fast_traitdispatch TestTr String true

Advanced features

The macros of the previous section are the official API of the package and should be reasonably stable. What follows in this section is "under the hood" and may well be updated (but still signalled with minor version changes).

Instead of using @traitimpl to add types to traits, it can be programmed. Running @traitimpl IsNice{Int} essentially expands to

SimpleTraits.trait(::Type{IsNice{X1}}) where {X1 <: Int} = IsNice{X1}

I.e. trait is the identity function for a fulfilled trait and returns Not{TraitInQuestion{...}} otherwise (this is the fall-back for <:Any). So instead of using @traitimpl this can be coded directly. Note that anything but a constant function will probably not be inlined away by the JIT and will lead to slower dynamic dispatch (see @check_fast_traitdispatch for a helper to check).

Example leading to static dispatch:

@traitdef IsBits{X}
SimpleTraits.trait(::Type{IsBits{X1}}) where {X1} = isbits(X1) ? IsBits{X1} : Not{IsBits{X1}}
istrait(IsBits{Int}) # true
istrait(IsBits{Array{Int,1}}) # false
struct A
    a::Int
end
istrait(IsBits{A}) # true

Dynamic dispatch can be avoided using a generated function or pure functions (sometimes they need to be annotated with Base.@pure):

@traitdef IsBits{X}
@generated function SimpleTraits.trait(::Type{IsBits{X1}}) where X1
    isbits(X1) ? :(IsBits{X1}) : :(Not{IsBits{X1}})
end

What is allowed in generated functions is heavily restricted, see Julia manual. In particular, no methods which are defined after the generated function are allowed to be called inside the generated function, otherwise this issue is encountered. Generally, try non-generated functions first and only in a pinch generated functions. See also issue 40.

Note that these programmed-traits can be combined with @traitimpl IsBits{XYZ}, i.e. program the general case and add exceptions with @traitimpl IsBits{XYZ}.

Trait-inheritance can also be hand-coded with above trick. For instance, the trait given by (in pseudo syntax) BeautyAndBeast{X,Y} <: IsNice{X}, !IsNice{Y}, BelongTogether{X,Y}:

@traitdef BeautyAndBeast{X,Y}
function SimpleTraits.trait(::Type{BeautyAndBeast{X,Y}}) where {X,Y}
    if istrait(IsNice{X}) && !istrait(IsNice{Y}) && BelongTogether{X,Y}
        BeautyAndBeast{X,Y}
    else
        Not{BeautyAndBeast{X,Y}}
    end
end

Note that this will lead to slower, dynamic dispatch, as the function is not pure (it depends on the global state of which types belong to the traits IsNice and BelongTogether).

Note also that trait functions can be generated functions:

@traitfn @generated fg(x::X) where {X; IsNice{X}} = (println(x); :x)

Innards

The function macroexpand shows the syntax transformations a macro does. Here the edited output of running it for the macros of this package:

julia> macroexpand(:(@traitdef Tr{X}))

struct Tr{X} <: SimpleTraits.Trait
end

julia> macroexpand(:(@traitimpl Tr{Int}))

# this function does the grouping of types into traits:
SimpleTraits.trait(::Type{Tr{X1}}) where X1 <: Int = Tr{X1}
SimpleTraits.istrait(::Type{Tr{X1}}) where X1 <: Int = true # for convenience, really

julia> macroexpand(:(@traitfn g(x::X) where {X; Tr{X}}= x+1))

@inline g(x::X) where {X} = g(trait(Tr{X}), x) # this is Tim's trick, using above grouping-function
g(::Type{Tr{X}},x::X) where {X} = x + 1 # this is the logic

julia> macroexpand(:(@traitfn g(x::X) where {X; !Tr{X}}= x+1000))

# the trait dispatch helper function needn't be defined twice,
# only the logic:
g(::Type{ Not{Tr{X}} }, x::X) where {X} = x + 1000

For a detailed explanation of how Tim's trick works, see Traits.jl: Dispatch on traits. The difference here is I make the methods containing the logic part of the same generic function (there it's in _f).

Base Traits

I started putting some Julia-Base traits together which can be loaded with using SimpleTraits.BaseTraits, see the source for all definitions.

Example, dispatch on whether an argument is immutable or not:

@traitfn f(x::X) where {X; IsImmutable{X}} = X(x.fld+1) # make a new instance
@traitfn f(x::X) where {X; !IsImmutable{X}} = (x.fld += 1; x) # update in-place

# use
mutable struct A; fld end
struct B; fld end
a=A(1)
f(a) # in-place
@assert a.fld == A(2).fld

b=B(1) # out of place
b2 = f(b)
@assert b==B(1)
@assert b2==B(2)

Background

This package grew out of an attempt to reduce the complexity of Traits.jl, but at the same time staying compatible (but which it isn't). Compared to Traits.jl, it drops support for:

  • Trait definition in terms of methods and constraints. Instead the user needs to assign types to traits manually. This removes the most complex part of Traits.jl: the checking whether a type satisfies a trait definition.
  • trait functions which dispatch on more than one trait. This allows to remove the need for generated functions, as well as removing the rules for trait-dispatch.

The reason for splitting this away from Traits.jl are:

  • creating a more reliable and easier to maintain package than Traits.jl
  • exploring inclusion in Base (see #13222).

My JuliaCon 2015 talk gives a 10 minute introduction to Traits.jl and SimpleTraits.jl.

The Future

The future of traits in Julia-Base: According to Stefan Karpinski's JuliaCon 2016 talk, Julia 1.0, traits are scheduled to land after Julia 1.0. My crystal ball tells me that all or most of the functionality of this package will be supported in the future trait system (multiparameter-traits may not be). Thus I expect the transition will be mostly a matter of a syntax update and less of a semantic update. Also, an advantage to using this package versus hand-coding Holy-traits will be that all occurrences of trait usage are clearly marked and thus easier to update.

The future of this package: I see it as light-weight package focusing on letting functions use dispatch based on traits. This dispatch is currently fairly limited, see section "Gotcha" above, but may be expanded in the future: either through something like in PR m3/multitraits.

In the unlikely event that I find myself with too much time on my hands, I may try to develop a companion package to allow the specification of a trait in terms of interfaces. The combination of the two packages would then have similar functionality to my experimental package Traits.jl. If anyone fancies a go at writing this companion package, I would be very happy to help and contribute. After the type-system overhaul lands, this should be much less hackish than what's in Traits.jl.

References

  • Traits.jl and its references. In particular here is an in-depth discussion on limitations of Holy-Traits, which this package implements.

To ponder

  • There is a big update sitting in the branch m3/multitraits; but I never quite finished it. It would also address the next point:
  • Could type inheritance be used for sub-traits (Jutho's idea)? In particular could it be used in such a way that it is compatible with the multiple inheritance used in Traits.jl?

Download Details:

Author: Mauro3
Source Code: https://github.com/mauro3/SimpleTraits.jl 
License: MIT license

#julia #simple 

Traits.jl: Exploration Of Traits in Julia

Traits.jl

This package is currently abandoned, although with the new type system in Julia 0.6 it should become viable to make this package maintainable; also 23117 should help. SimpleTraits.jl package provides traits, although without the interface specification present in this package.

Traits.jl allows to:

define traits/interfaces with @traitdef

implement interfaces with @traitimpl

make functions which dispatch on traits with @traitfn

It's based on what I think traits should be:

contracts on a type or a tuple of types. The contract can contain required methods but also other assertions. (Assertions could be that certain fields are present or that it has some storage structure, etc.)

they needn't be declared explicitly, but can be.

they allow dispatch to work with them

Julia's generic functions are very good to set up contracts as mentioned in (1). But Julia does not support (2) or (3) yet. (2) is fairly easy to implement. However, dispatch on a "contract" is not easily possible, but Tim Holy recently came up with a trick. The cool thing about that trick is that the generated machine-code for a trait-dispatch function should be identical to a duck-typed function, i.e. there is no loss in performance.

Traits.jl adds those kind of traits to Julia, using Tim's trick combined with stagedfunctions and extensive facilities to define traits. See also the Julia-issue 6975 concerning interfaces/traits.

My JuliaCon 2015 talk gives a 10 minute introduction to Traits.jl. Also, Jeff mentioned Traits.jl during Q&A in his JuliaCon talk, suggesting that trait functionality may well be added to Julia-Base (but different to Traits.jl).

Note that there is also the SimpleTraits.jl package which is a lot simpler than Traits.jl, essentially just macro-sugar for Holy-traits. Thus SimpleTraits.jl does not contain the possibility to define traits by specifying required methods. Instead types need to be added to a trait manually; the rest is almost identical to Traits.jl.

Example examples/ex1.jl:

using Traits
# Check Cmp-trait (comparison) which is implemented in Traits.jl/src/commontraits.jl
@assert istrait(Cmp{Int,Float64})        # Int and Float64 can be compared
@assert istrait(Cmp{Int,String})==false  # Int and String cannot be compared

# make a new trait and add a type to it:
@traitdef MyTr{X,Y} begin
    foobar(X,Y) -> Bool # All type-tuples for which there is a method foo
                        # with that signature belong to MyTr
end
type A
    a::Int
end
@assert istrait(MyTr{A,A})==false  # foobar not implement yet
foobar(a::A, b::A) = a.a==b.a      # implement it
@assert istrait(MyTr{A,A})         # voila!
@assert istrait(MyTr{Int,Int})==false

# make a function which dispatches on traits:
@traitfn ft1{X,Y; Cmp{X,Y}}(x::X,y::Y)  = x>y ? 5 : 6
@traitfn ft1{X,Y; MyTr{X,Y}}(x::X,y::Y) = foobar(x,y) ? -99 : -999

ft1(4,5)        # ==6    i.e. dispatches to first definition
ft1(A(5), A(6)) # ==-999 i.e. dispatches to second definition

ft1("asdf", 6)
# -> ERROR: TraitException("No matching trait found for function ft1")

Package status

New features are documented in NEWS as they are added. I keep some notes, musings and plans in dev_notes.md.

This is a fairly experimental package and I will not try to keep backwards compatibility as I move on. Please try it out and give me feedback, issues or pull requests!

Update April 2016: I will not have much time to work on this in the near-ish future and it will be a low priority for me to update Traits.jl once Julia 0.5 is out.

Syntax

The source of below examples is in examples/ex2.jl. Most of the important functions are documented and will respond to ? in the REPL.

Trait definition (for details see traitdef.md):

using Traits
# simple
@traitdef Tr1{X} begin
    fun1(X) -> Number   # this means a method with signature fun1(::X)
                        # returning a Number
end
@traitdef Tr2{X,Y} begin
    fun2(X,Y) -> Number
end
# subtrait
@traitdef Tr3{X,Y} <: Tr1{X}, Tr2{X,Y} begin
    fun3(X,Y,Int)
end
# with additional constraint on the types
@traitdef Tr4{X,Y} begin
    fun4(X,Y)
    @constraints begin
        # both Types need to start with the same letter:
        string(X.name)[1]==string(Y.name)[1]
    end
end

# using associated types
@traitdef Tr5{X,Y} begin
    Z = promote_type(X,Y) # calculates Z from X and Y
    fun5(X,Y) -> Z
end

Note that return-type checking is quite experimental. It can be turned off with check_return_types(false).

Trait implementation and checking with istrait:

# manual definiton, i.e. just define the functions
fun1(x::Int) = 5x
@assert istrait(Tr1{Int})

# using @traitimpl
@traitimpl Tr1{Float64} begin
    fun1(x::Float64) = 7x # the explicit "::Float64" is needed at the moment
end
@assert istrait(Tr1{Float64})

# wrong usage of @traitimpl
try
    @traitimpl Tr1{Float32} begin
        fun1(x::Float64) = 7x # if the explicit type is wrong, it may error
    end
catch e
    println(e)
end

# This gives an error because supertypes have not been defined yet:
try
    eval(:(
    @traitimpl Tr3{Int, Int} begin
        fun3(x::Int, y::Int, t::Int) = x+y+t
    end))
catch e
    println(e)
end

# this works:
@traitimpl Tr2{Int, Int} begin
    fun2(x::Int, y::Int) = x+y
end
@traitimpl Tr3{Int, Int} begin
    fun3(x::Int, y::Int, t::Int) = x+y+t
end
@traitimpl Tr4{Int, Int} begin
    fun4(x::Int, y::Int) = x+y
end

# This gives an error because constraints are not satisfied:
# Int starts with an "I" whereas Float64 with an "F":
try
    eval(:(
    @traitimpl Tr4{Int, Float64} begin
        fun4(x::Int, y::Float64) = x+y
    end))
catch e
    println(e)  # ErrorException("assertion failed: istrait(Tr4{Int,Float64})")
end

Trait functions & dispatch (for details see traitfns.md):

@traitfn tf1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y) = fun1(a) + fun1(b)             # I
@traitfn tf1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y, c::Int) = fun1(a) + fun1(b) + c # II
@traitfn tf1{X, Y; Tr2{X,Y}}(a::X, b::Y) = fun2(a,b)                           # III
# Note that all the type-parameters are in the {} and that all
# arguments need a type parameter (a limitation of the
# macro-parser). This doesn't work:
#
# julia> @traitfn ttt1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y, c) = fun1(a) + fun1(b) + c
# ERROR: type Symbol has no field args
#
# But this works:
#
# julia> @traitfn ttt1{X, Y, Z; Tr1{X}, Tr1{Y}}(a::X, b::Y, c::Z) = fun1(a) + fun1(b) + c
# ttt1 (generic function with 6 methods)


# tf1 now dispatches on traits
@assert tf1(5.,6.)==77. # -> 77 ; dispatches to I because istrait(Tr1{Float64})
                        #         but not istrait(Tr2{Float64,Float64})
@assert tf1(5.,6.,77)==154. # -> 154. ; dispatches to II because of the extra argument

# Errors because of dispatch ambiguity:
try
    tf1(5,6)  # istrait(Tr1{Int}) and istrait(Tr2{Int,Int}) are both true!
catch e
    println(e)
end

# Implementing Tr1 for a type will make it work with tf1:
type MyType
    a::Int
end
try
    tf1(MyType(8), 9) # not implemented yet
catch e
    println(e)
end
@traitimpl Tr1{MyType} begin
    fun1(x::MyType) = x.a+9
end

@assert tf1(MyType(8), 9)==62 # -> 62 ; dispatches to I

Generated code

Continuing the example from last section, let's have a look at the llvm code:

f(x,y) = 7x + 7y
@code_llvm f(5.,6.)
@code_llvm tf1(5.,6.)

both produces

define double @"julia_f;41342"(double, double) {
top:
  %2 = fmul double %0, 7.000000e+00, !dbg !1388
  %3 = fmul double %1, 7.000000e+00, !dbg !1388
  %4 = fadd double %2, %3, !dbg !1388
  ret double %4, !dbg !1388
}

However, for more complicated functions code is not quite the same, see test/perf/perf.jl.

Inner workings

Julia is already good at specifying function-interfaces for types/families of types with its ability of multiple dispatch. However, for a particular type these function-interfaces are implicitly defined by defining the particular functions needed for that "interface". For example, to support iteration a type needs to implement the functions start, next, and done with a specific calling convention. What Julia is lacking is a way to formally define and implement an interface, and, probably more importantly, to dispatch on interfaces. This package implements both of these abilities.

Probably easiest to figure out what Traits.jl is doing with its macros, is to have a look at the files tests/manual-trait*.jl. There I manually coded what the macros do.

Dispatch on traits

In Julia dispatch works on types, to extend this to traits I use @timholy's trick. His trick uses a method to assign its arguments to a trait. That trait-function is then used for dispatch in another function. Example of Tim's trick (examples/ex_tims_traits.jl):

type Trait1 end
type Trait2 end
type Trait3 end
# Now define function f which should dispatch on those traits:
f(x,y) = _f(x,y, traitfn(x,y))
# Logic which dispatches on trait:
_f(x,y,::Type{Trait1}) = x+y
_f(x,y,::Type{Trait2}) = x-y
_f(x,y,::Type{Trait3}) = x*y
# Association of types with traits through method definitions:
# Throw error as default
traitfn{T,S}(x::T,y::S) = error("Function f not implemented for type ($T,$S)")
# Add types-tuples to Trait1, Trait2 or Trait3
traitfn(::Int, ::Int) = Trait1
traitfn(::Int, ::FloatingPoint) = Trait2
traitfn(::FloatingPoint, ::FloatingPoint) = Trait3
# use
@assert f(3,4)==7      # Trait1
@assert f(3,4.)==-1.0  # Trait2
@assert f(3.,4.)==12.0 # Trait3
# add another type-tuple to Trait3
traitfn(::String, ::String) = Trait3
@assert f("Lorem ", "Ipsum")=="Lorem Ipsum"

What does this add compared to what we had before with usual dispatch? When a new type, say A, is created it can made to work with the function f without needing to re-define f for that particular type. Instead all that is needed is to add it to the traitfn, and choosing the exact behavior of f by the type traitfn returns:

traitfn(::A, ::Int) = Trait1()
traitfn(::Int, ::A) = Trait1()

Therefore traitfn is in effect a function that groups type-tuples into different Traits (via method definitions) and returns the appropriate type when called (which is then used inside f for dispatch). However, the limitation of this approach is that traitfn is married to f as can be seen from trying to reuse it for another function g which wants to dispatch in different traits:

g(x,y) = _g(x,y, traitfn(x,y))
# Logic which dispatches on trait:
_g(x,y,::Type{Trait1}) = 2x+2y
_g(x,y,::Type{Trait4}) = 2x-2y  # g doesn't care about Trait2&3 but about 4

# However, say Trait4 should also be implemented by {Int, FloatingPoint} just
# like Trait2:
traitfn(::Int, ::FloatingPoint) = Trait4 # this will overwrite the
                                         # Trait2 definition above!
g(5, 6.) # doesn't work

This limitation can be overcome having a different traitfn for each function which uses trait dispatch. However, it becomes rather tricky to remember to update all different traitfns if a type-tuple is added to a certain trait! This problem is solved in Traits.jl by de-coupling the trait definition from the trait dispatch helper function, both of which was done above by the traitfn.

Whether a trait is defined is checked by the istrait function (completely independent of any functions doing trait-dispatch). For instance istrait(Tr1{Int,Float64}) will check whether Tr1 is implemented by Tuple{Int,Float64}.

For the trait dispatch of a function, say f1, a generated-method is used (which also belongs to the generic function f1, so I needn't worry about scopes). The first time the generated method is called with arguments of a specific type, it figures out which trait or traits-tuple featuring in the method definitions of f1 that type satisfies, and constructs a constant method returning that trait. This trait is then used for dispatch. Time for an example!

For methods definition like so

@traitfn f1{X,Y<:Integer; D1{Y}, D4{X,Y}}(x::X,y::Y) = x + sin(y)
@traitfn f1{S,T<:Integer; D1{S}, D1{T}  }(s::S,t::T) = sin(s) - sin(t)
@traitfn f1{X,Y<:FloatingPoint; D1{X}, D1{Y}  }(x::X,y::Y) = cos(x) - cos(y)

the underlying definitions are:

f1{X,Y<:Integer}(x::X, y::Y)       = f1(f1(_TraitDispatch,x, y), x, y)
f1{S,T<:Integer}(s::S, t::T)       = f1(f1(_TraitDispatch,s, t), s, t)
f1{X,Y<:FloatingPoint}(x::X, y::Y) = f1(f1(_TraitDispatch,x, y), x, y)

# The logic for different traits is:
@inline f1{X,Y<:Integer}(::Type{(D1{Y}, D4{X,Y})}, x::X, y::Y) = x + sin(y)
@inline f1{S,T<:Integer}(::Type{(D1{S}, D1{T})}, s::S, t::T) = sin(s) - sin(t)
@inline f1{X,Y<:FloatingPoint}(::Type{(D1{X}, D1{Y})}, x::X, y::Y) = cos(x) - cos(y)

# Trait dispatch happens in these generated functions
@generated function f1{X1,X2<:Integer}(::Type{_TraitDispatch}, x1::X1, x2::X2)
    # Figure out which trait matches.  Note below list is updated as more
    # trait-dispatched methods are added to f1.
    traittypes = [(D1{X2}, D4{X1,X2}), (D1{X1}, D1{X2})]

    # errors if not a single match is found:
    traittyp = Traits.traitdispatch(traittypes, $(fn.name))

    out = :(())
    for s in traittyp
        push!(out.args, :($s))
    end
    return out
end
# For each type signature there is a trait-dispatch function
@generated function f1{X1,X2<:FloatingPoint}(::Type{_TraitDispatch}, x1::X1, x2::X2)
...
end

Dispatch, happening in the function Traits.traitdispatch is quite simple taking trait-hierarchies into account. Although, note that it is easily possible to have unsolvable ambiguities with trait-dispatch as traits do not have a strict hierarchy like types.

Other trait implementations

See the Julia-issue 6975 for a discussion about interfaces/traits.

My simpler traits implementation: SimpleTraits.jl Similar to Traits.jl but:

  • types are added to traits manually
  • some restrictions on @traitfn

@Rory-Finnegan's Interfaces.jl

  • does dispatch on traits
  • only single parameter traits
  • uses a new type of mutable Union coded in C

Jason Morton's package Typeclass.jl

  • does multiple parameters
  • no dispatch

@pao's https://gist.github.com/pao/2432554

  • simple
  • no dispatch on trait

https://github.com/JuliaLang/julia/pull/7025 (and https://gist.github.com/tknopp/ed53dc22b61062a2b283)

  • @tknopp
  • interfaces are just added to types
  • no dispatch on interfaces

https://gist.github.com/abe-egnor/503661eb4cc0d66b4489

  • @abe-egnor
  • no dispatch

https://github.com/abeschneider/TypeTraits.jl

  • only does fields of types, as far as I can tell

Graphs.jl: http://graphsjl-docs.readthedocs.org/en/latest/interface.html

  • interface verification

@timholy's trick https://github.com/JuliaLang/julia/issues/2345#issuecomment-54537633

  • does limited dispatch: a function returns a true/false type depending on the input types
  • Jeff suggested some additions to it.

Download Details:

Author: Mauro3
Source Code: https://github.com/mauro3/Traits.jl 
License: View license

#julia 

LsqFit.jl: Simple Curve Fitting in Julia

LsqFit.jl

The LsqFit package is a small library that provides basic least-squares fitting in pure Julia under an MIT license. The basic functionality was originally in Optim.jl, before being separated into this library. At this time, LsqFit only utilizes the Levenberg-Marquardt algorithm for non-linear fitting.

Basic Usage

There are top-level methods curve_fit() and margin_error() that are useful for fitting data to non-linear models. See the following example. Let's first define the model function:

using LsqFit

# a two-parameter exponential model
# x: array of independent variables
# p: array of model parameters
# model(x, p) will accept the full data set as the first argument `x`.
# This means that we need to write our model function so it applies
# the model to the full dataset. We use `@.` to apply the calculations
# across all rows.
@. model(x, p) = p[1]*exp(-x*p[2])

The function applies the per observation function p[1]*exp(-x[i]*p[2]) to the full dataset in x, with i denoting an observation row. We simulate some data and chose our "true" parameters.

# some example data
# xdata: independent variables
# ydata: dependent variable
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

Now, we're ready to fit the model.

fit = curve_fit(model, xdata, ydata, p0)
# fit is a composite type (LsqFitResult), with some interesting values:
#    dof(fit): degrees of freedom
#    coef(fit): best fit parameters
#    fit.resid: residuals = vector of residuals
#    fit.jacobian: estimated Jacobian at solution
lb = [1.1, -0.5]
ub = [1.9, Inf]
p0_bounds = [1.2, 1.2] # we have to start inside the bounds
# Optional upper and/or lower bounds on the free parameters can be passed as an argument.
# Bounded and unbouded variables can be mixed by setting `-Inf` if no lower bounds
# is to be enforced for that variable and similarly for `+Inf`
fit_bounds = curve_fit(model, xdata, ydata, p0_bounds, lower=lb, upper=ub)

# We can estimate errors on the fit parameters,
# to get standard error of each parameter:
sigma = stderror(fit)
# to get margin of error and confidence interval of each parameter at 5% significance level:
margin_of_error = margin_error(fit, 0.05)
confidence_inter = confidence_interval(fit, 0.05)

# The finite difference method is used above to approximate the Jacobian.
# Alternatively, a function which calculates it exactly can be supplied instead.
function jacobian_model(x,p)
    J = Array{Float64}(undef, length(x), length(p))
    @. J[:,1] = exp(-x*p[2])     #dmodel/dp[1]
    @. @views J[:,2] = -x*p[1]*J[:,1] #dmodel/dp[2], thanks to @views we don't allocate memory for the J[:,1] slice
    J
end
fit = curve_fit(model, jacobian_model, xdata, ydata, p0)

Multivariate regression

There's nothing inherently different if there are more than one variable entering the problem. We just need to specify the columns appropriately in our model specification.

using LsqFit

x = collect(range(0, stop=200, length=201))
y = collect(range(0, stop=200, length=201))

xy = hcat(x, y)

function twoD_Gaussian(xy, p)
    amplitude, xo, yo, sigma_x, sigma_y, theta, offset = p
    a = (cos(theta)^2)/(2*sigma_x^2) + (sin(theta)^2)/(2*sigma_y^2)
    b = -(sin(2*theta))/(4*sigma_x^2) + (sin(2*theta))/(4*sigma_y^2)
    c = (sin(theta)^2)/(2*sigma_x^2) + (cos(theta)^2)/(2*sigma_y^2)

    # creating linear meshgrid from xy
    x = xy[:, 1]
    y = xy[:, 2]
    g = offset .+ amplitude .* exp.( - (a.*((x .- xo).^2) + 2 .* b .* (x .- xo) .* (y .- yo) + c * ((y .- yo).^2)))
    return g[:]
end

p0 = Float64.([3, 100, 100, 20, 40, 0, 10])
data = twoD_Gaussian(xy, p0)

# Noisy data
data_noisy = data + 0.2 * randn(size(data))

fit = LsqFit.curve_fit(twoD_Gaussian, xy, data_noisy, p0)

Evaluating the Jacobian and using automatic differentiation

The default is to calculate the Jacobian using a central finite differences scheme if no Jacobian function is provided. The default is to use central differences because it can be more accurate than forward finite differences, but at the expense of computational cost. It is possible to switch to forward finite differences, like MINPACK uses for example, by specifying

fit = curve_fit(model, xdata, ydata, p0; autodiff=:finiteforward)

It is also possible to use forward mode automatic differentiation as implemented in ForwardDiff.jl by using the autodiff=:forwarddiff keyword.

fit = curve_fit(model, xdata, ydata, p0; autodiff=:forwarddiff)

Here, you have to be careful not to manually restrict any types in your code to, say, Float64, because ForwardDiff.jl works by passing a special number type through your functions, to automagically calculate the value and gradient with one evaluation.

In-place model and Jacobian

It is possible to either use an in-place model, or an in-place model and an in-place Jacobian. It might be pertinent to use this feature when curve_fit is slow, or consumes a lot of memory.

model_inplace(F, x, p) = (@. F = p[1] * exp(-x * p[2]))

function jacobian_inplace(J::Array{Float64,2},x,p)
        @. J[:,1] = exp(-x*p[2])
        @. @views J[:,2] = -x*p[1]*J[:,1]
    end
fit = curve_fit(model_inplace, jacobian_inplace, xdata, ydata, p0; inplace = true)

Geodesic acceleration

This package implements optional geodesic acceleration, as outlined by this paper. To enable it, one needs to specify the function computing the directional second derivative of the function that is fitted, as the avv! parameter. It is also preferable to set lambda and min_step_qualityto 0:

curve_fit(model, xdata, ydata, p0; avv! = Avv!,lambda=0, min_step_quality = 0)

Avv! must have the following form:

  • p is the array of parameters
  • vis the direction in which the direction is taken
  • dir_deriv is the output vector (the function is necessarily in-place)
function Avv!(dir_deriv,p,v)
        v1 = v[1]
        v2 = v[2]
        for i=1:length(xdata)
            #compute all the elements of the Hessian matrix
            h11 = 0
            h12 = (-xdata[i] * exp(-xdata[i] * p[2]))
            #h21 = h12
            h22 = (xdata[i]^2 * p[1] * exp(-xdata[i] * p[2]))

            # manually compute v'Hv. This whole process might seem cumbersome, but
            # allocating temporary matrices quickly becomes REALLY expensive and might even
            # render the use of geodesic acceleration terribly inefficient
            dir_deriv[i] = h11*v1^2 + 2*h12*v1*v2 + h22*v2^2

        end
end

Typically, if the model to fit outputs [y_1(x),y_2(x),...,y_m(x)], and that the input data is xdata then Avv!should output an array of size m, where each element is v'*H_i(xdata,p)*v, where H_iis the Hessian matrix of the output y_iwith respect to the parameter vector p.

Depending on the size of the dataset, the complexity of the model and the desired tolerance in the fit result, it may be worthwhile to use automatic differentiation (e.g. via Zygote.jl or ForwardDiff.jl) to determine the directional derivative. Although this is potentially less efficient than calculating the directional derivative manually, this additional information will generally lead to more accurate results.

An example of such an implementation is given by:

using LinearAlgebra, Zygote

function Avv!(dir_deriv,p,v)
    for i=1:length(xdata)
        dir_deriv[i] = transpose(v) * Zygote.hessian(z->model(xdata[i],z),p) * v
    end
end

Existing Functionality

fit = curve_fit(model, [jacobian], x, y, [w,] p0; kwargs...):

  • model: function that takes two arguments (x, params)
  • jacobian: (optional) function that returns the Jacobian matrix of model
  • x: the independent variable
  • y: the dependent variable that constrains model
  • w: (optional) weight applied to the residual; can be a vector (of length(x) size or empty) or matrix (inverse covariance matrix)
  • p0: initial guess of the model parameters
  • kwargs: tuning parameters for fitting, passed to levenberg_marquardt, such as maxIter, show_trace or lower and upper bounds
  • fit: composite type of results (LsqFitResult)

This performs a fit using a non-linear iteration to minimize the (weighted) residual between the model and the dependent variable data (y). The weight (w) can be neglected (as per the example) to perform an unweighted fit. An unweighted fit is the numerical equivalent of w=1 for each point (although unweighted error estimates are handled differently from weighted error estimates even when the weights are uniform).


sigma = stderror(fit; atol, rtol):

  • fit: result of curve_fit (a LsqFitResult type)
  • atol: absolute tolerance for negativity check
  • rtol: relative tolerance for negativity check

This returns the error or uncertainty of each parameter fit to the model and already scaled by the associated degrees of freedom. Please note, this is a LOCAL quantity calculated from the Jacobian of the model evaluated at the best fit point and NOT the result of a parameter exploration.

If no weights are provided for the fits, the variance is estimated from the mean squared error of the fits. If weights are provided, the weights are assumed to be the inverse of the variances or of the covariance matrix, and errors are estimated based on these and the Jacobian, assuming a linearization of the model around the minimum squared error point.

margin_of_error = margin_error(fit, alpha=0.05; atol, rtol):

  • fit: result of curve_fit (a LsqFitResult type)
  • alpha: significance level
  • atol: absolute tolerance for negativity check
  • rtol: relative tolerance for negativity check

This returns the product of standard error and critical value of each parameter at alpha significance level.

confidence_interval = confidence_interval(fit, alpha=0.05; atol, rtol):

  • fit: result of curve_fit (a LsqFitResult type)
  • alpha: significance level
  • atol: absolute tolerance for negativity check
  • rtol: relative tolerance for negativity check

This returns confidence interval of each parameter at alpha significance level.


covar = estimate_covar(fit):

  • fit: result of curve_fit (a LsqFitResult type)
  • covar: parameter covariance matrix calculated from the Jacobian of the model at the fit point, using the weights (if specified) as the inverse covariance of observations

This returns the parameter covariance matrix evaluated at the best fit point.

Download Details:

Author: JuliaNLSolvers
Source Code: https://github.com/JuliaNLSolvers/LsqFit.jl 
License: View license

#julia #curve #fitting 

AudioIO.jl: Simple Audio IO in Julia

AudioIO.jl 

AudioIO interfaces to audio streams, including real-time recording, audio processing, and playback through your sound card using PortAudio. It also supports reading and writing audio files in a variety of formats. It is under active development and the low-level API could change, but the basic functionality (reading and writing files, the play function, etc.) should be stable and usable by the general Julia community.

File I/O

File I/O is handled by libsndfile, so we can support a wide variety of file and sample formats. Use the AudioIO.open function to open a file. It has the same API as the built-in Base.open, but returns an AudioFile type. Opening an audio file and reading its contents into an array is as simple as:

f = AudioIO.open("data/never_gonna_give_you_up.wav")
data = read(f)
close(f)

Or to hand closing the file automatically (including in the case of unexpected exceptions), we support the do block syntax:

data = AudioIO.open("data/never_gonna_let_you_down.wav") do f
    read(f)
end

By default the returned array will be in whatever format the original audio file is (Float32, UInt16, etc.). We also support automatic conversion by supplying a type:

data = AudioIO.open("data/never_gonna_run_around.wav") do f
    read(f, Float32)
end

Basic Array Playback

Arrays in various formats can be played through your soundcard. Currently the native format that is delivered to the PortAudio backend is Float32 in the range of [-1, 1]. Arrays in other sizes of float are converted. Arrays in Signed or Unsigned Integer types are scaled so that the full range is mapped to [-1, 1] floating point values.

To play a 1-second burst of noise:

julia> v = rand(44100) * 0.1
julia> play(v)

AudioNodes

In addition to the basic play function you can create more complex networks of AudioNodes in a render chain. In fact, when using the basic play to play an Array, behind the scenes an instance of the ArrayPlayer type is created and added to the master AudioMixer inputs. Audionodes also implement a stop function, which will remove them from the render graph. When an implicit AudioNode is created automatically, such as when using play on an Array, the play function should return the audio node that is playing the Array, so it can be stopped if desired.

To explictly do the same as above:

julia> v = rand(44100) * 0.1
julia> player = ArrayPlayer(v)
julia> play(player)

To generate 2 sin tones:

julia> osc1 = SinOsc(440)
julia> osc2 = SinOsc(660)
julia> play(osc1)
julia> play(osc2)
julia> stop(osc1)
julia> stop(osc2)

All AudioNodes must implement a render function that can be called to retreive the next block of audio.

AudioStreams

AudioStreams represent an external source or destination for audio, such as the sound card. The play function attaches AudioNodes to the default stream unless a stream is given as the 2nd argument.

AudioStream is an abstract type, which currently has a PortAudioStream subtype that writes to the sound card, and a TestAudioStream that is used in the unit tests.

Currently only 1 stream at a time is supported so there's no reason to provide an explicit stream to the play function. The stream has a root mixer field which is an instance of the AudioMixer type, so that multiple AudioNodes can be heard at the same time. Whenever a new frame of audio is needed by the sound card, the stream calls the render method on the root audio mixer, which will in turn call the render methods on any input AudioNodes that are set up as inputs.

Installation

To install the latest release version, simply run

julia> Pkg.add("AudioIO")

If you want to install the lastest master, it's almost as easy:

julia> Pkg.clone("AudioIO")
julia> Pkg.build("AudioIO")

This package is no longer being maintained. I suggest you check out the JuliaAudio family of packages for audio-related work in Julia

Download Details:

Author: ssfrr
Source Code: https://github.com/ssfrr/AudioIO.jl 
License: MIT license

#julia #audio 

Monty  Boehm

Monty Boehm

1659656880

DThree.jl: Simple interface to D3.js From Julia

Simple interface to d3 (http://d3js.org) and other JavaScript libraries for chart making.

This package for Julia provides a simple interface for using d3 syntax within julia.

It isn't very clever, basically it takes a d3 command like:

d3.selectAll("p").style("color", "white")

And turns it into a julia call like:

d3.selectAll("p").style("color", "white")

(Only after)

using DThree
d3 = D3()

Okay, you might guess the style. This just pieces together a string of JavaScript that will get inserted into a web page.

The implementation is pretty stupid, it just makes fields named after the main methods and creates a function when the object is instantiated. The functions return a D3 instance so they can be chained, as above.

If the field isn't present, then the interface can look one of two ways:

d3[:selectAll]("p")._("color", "white") ## two ways [:symbol](args...) or _("meth", args...)

By default, strings are quoted. To stop that, wrap the string in asis (like R's I operator). This is necessary when the argument refers to a JavaScript object.

using DThree
style = """
.chart div {
  font: 10px sans-serif;
  background-color: steelblue;
  text-align: right;
  padding: 3px;
  margin: 1px;
  color: white;
  }
"""
  
w = figure()
DThree.empty_page(w, style=style) # loads D3 libraries

d3 = D3()
d3.select("body").append("div").attr("class", "chart") |> js

data = [4,8,16,23,42]

d3.select(".chart").selectAll("div").
    data(data).
      enter().append("div").
       style("width", asis"""function(d) { return d * 10 + "px"; }""").
         text(asis"""function(d) { return d; }""")  |> js
    

The js call comes from Blink, as does the figure object.

Blink

This package also borrows the figure manipulation tools of Immerse and the HTML windows of Blink to create canvases to manipulate. The basic idea would follow from this example

using DThree, Plots; plotly()
w = figure()
Blink.title(w, "title")
plot(sin, 0, 2pi)
plot!(cos)

Author: jverzani
Source Code: https://github.com/jverzani/DThree.jl 
License: MIT license

#julia #d3