Monotonic interpolations #238

mateuszbaran · 2018-09-26T10:59:57Z

I've prepared an implementation of monotonic interpolations (issue #105 ) based on the code posted there. A few methods are implemented, tested and documented. This code partially covers WIP PR #110 .

Is there anything I should improve?

implemented algorithms: * linear * finite difference * cardinal * Fritsch-Carlson * Fritsch-Butland * Steffen

codecov-io · 2018-09-26T11:09:06Z

Codecov Report

Merging #238 into master will decrease coverage by 0.5%.
The diff coverage is 96.05%.

@@            Coverage Diff             @@
##           master     #238      +/-   ##
==========================================
- Coverage   99.26%   98.76%   -0.51%     
==========================================
  Files          18       19       +1     
  Lines         410      487      +77     
==========================================
+ Hits          407      481      +74     
- Misses          3        6       +3

Impacted Files	Coverage Δ
src/Interpolations.jl	`100% <ø> (ø)`	⬆️
src/io.jl	`94.11% <75%> (-2.44%)`	⬇️
src/monotonic/monotonic.jl	`97.22% <97.22%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 67049ef...e6a38cb. Read the comment docs.

tomasaschan

Nice work - I like what I see here! Thanks for a great contribution!

tomasaschan · 2018-09-26T13:45:03Z

src/monotonic/monotonic.jl

+    MonotonicInterpolation
+
+Monotonic interpolation up to third order represented by type, knots and
+coefficients. Type is any concrete subtype of `MonotonicInterpolationType`.


Type is any concrete subtype of MonotonicInterpolationType.

This sentence seems superfluous, given that it's also specified by the signature. (We don't write out the types of the other arguments, while e.g. K is a much less obvious type argument name...)

I've removed that sentence and renamed K to TKnots,

tomasaschan · 2018-09-26T13:48:08Z

src/monotonic/monotonic.jl

+function MonotonicInterpolation(::Type{TWeights}, it::IType, knots::K, A::AbstractArray{Tel,1},
+    m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, Tel, IType<:MonotonicInterpolationType, K<:AbstractVector{<:Number}}
+
+    isconcretetype(IType) || error("The b-spline type must be a leaf type (was $IType)")


The b-spline type [...]

Should it be

The spline type [...]

?

Right, of course, I've fixed it.

tomasaschan · 2018-09-26T13:53:35Z

src/monotonic/monotonic.jl

+    issorted(knots) || error("knot-vector must be sorted in increasing order")
+end
+
+function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},


These methods could do with some more documentation, that explains the mathematical reasoning behind the algorithm (or at least refers to resources that explain it). The Fritsch & Carlson version has good docs, but the other ones have none...

This code was taken almost straight from the implementation provided by niclasmattsson in discussion for issue #105 . I haven't really looked into the details except for the interpretation of coefficients A, m, c and d. I've added some documentation for these vectors in a new commit.

tomasaschan · 2018-09-26T14:08:29Z

src/monotonic/monotonic.jl

+    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
+end
+
+function derivative(itp::MonotonicInterpolation, x::Number)


Why is this function not named gradient (returning an array of length 1) or gradient1 (returning a scalar)?

I've used the convention from ForwardDiff.jl but gradient1 seems like a more fitting choice. Should I write a new method for the gradient function as well (receiving and returning one element arrays)?

The behavior of gradient should, in order to match how other interpolation schemes work, be

function gradient(itp, x...) g = Vector{eltype(itp)}(undef, lenght(x)) # fill g with values g end

I think (correct me if I'm wrong) that the monotonic implementation here only supports 1D interpolations, which means this would effectively be gradient(itp, x) = [derivative(itp, x)]. gradient1 is just for convenience when the user knows it's a 1D interpolation object, but it can also avoid the array allocation: gradient1(itp, x) = derivative(itp, x). In other words, if you rename derivative to gradient1 and then add gradient(itp, x) = [gradient1(itp, x)] (with typed arguments) you'll be good.

As of #240 there's a fallback defined for gradient! that just calls gradient.

But gradient should return an SVector rather than Vector---there's a ton of overhead for Vector but essentially none for SVector.

I've renamed derivative to gradient1 and introduced new functions: gradient, hessian and hessian1. hessian and gradient return SVectors.

And yes, only 1D -> 1D monotonic interpolation is supported by this code.

* additional documentation for interpolation coefficients * references added to `calcTangents` for Fritsch & Butland and Steffen methods * changed the type parameter name `K` to `TKnots` * changed "b-spline" to "spline" in error string

mateuszbaran · 2018-09-26T18:52:04Z

Thanks! I still don't understand well how things work in this library, so I'm glad this code is good enough :).

tomasaschan

Nice!

After looking at the code again, I see a bunch of other places where the naming of type arguments is inconsistent. Sometimes it's IT, sometimes it's XType, and sometimes is TThing. I prefer the TThing scheme (probably because I'm used to it from writing C# for work all day), but my main concern is that they should at least be consistent.

Also, we should probably have methods for calculating Hessians as well as graidents. For 1D interpolation, the Hessian is simply the 2nd order derivative (in a 1x1 matrix), so it should be quite easy to implement hessian and hessian1 (returning a 1x1 matrix and a scalar, respectively) to fill that gap.

timholy

Looks great, thanks!

timholy · 2018-09-27T09:05:10Z

src/monotonic/monotonic.jl

+    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
+end
+
+function derivative(itp::MonotonicInterpolation, x::Number)


But gradient should return an SVector rather than Vector---there's a ton of overhead for Vector but essentially none for SVector.

mateuszbaran · 2018-09-27T10:47:55Z

Thank you, these are great suggestions. I'll push a new version soon.
I have one question: the getindex function, the way it's defined in the library, does not make much sense here (the nodes/knots are not located at subsequent integers). What should I do about it?

…interpolation

# Conflicts: # test/gradient.jl

tomasaschan

This is coming along very nicely! Almost merge-ready, IMO - big thanks!

tomasaschan · 2018-09-28T07:17:52Z

src/monotonic/monotonic.jl

@@ -180,7 +180,12 @@ function (itp::MonotonicInterpolation)(x::Number)
    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
 end

-function derivative(itp::MonotonicInterpolation, x::Number)
+function gradient(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})


Nice! This is almost what we want, but the type of x should probably be VarArg{<:Number,1} - see e.g. how it's done for Gridded.

I've just checked it and it looks like function gradient(itp::MonotonicInterpolation, x::Number) works. I think gradient isn't really a vararg function when it expects exactly two arguments :)

I tried vararg with hessian first and with function hessian(itp::MonotonicInterpolation, x::Vararg{<:Number,1}) I got stack overflow due to this function but with function hessian(itp::MonotonicInterpolation, x::Number) it just works.

gradient is a varargs function, because the notation args::Vararg{T,N} means "supply N separate arguments of type T." But in your case since all MonotonicInterpolations are 1-dimensional, this is equivalent to gradient(itp::MonotonicInterpolation, x::Number).

x::Vararg{Number,1} should be the same as x::Number, but Vararg{<:Number,1} may not be.

n-dimensional positions are encoded as varargs; like the rest of julia, indexing with vectors means something like this:

julia> A = reshape(1:20, 4, 5) 4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64: 1 5 9 13 17 2 6 10 14 18 3 7 11 15 19 4 8 12 16 20 julia> A[2:3, 3:4] 2×2 Array{Int64,2}: 10 14 11 15

This is not what you mean here, so definitely get rid of the AbstractArray.

I'll change x::AbstractArray{<Number,1} to x::Number then. In this case x::Vararg{Number,1} causes a stack overflow while x::Number does not:

Test threw exception Expression: ≈((Interpolations.hessian(itp, [t]))[1], hessval, atol=1.0e-12) StackOverflowError: Stacktrace: [1] hessian(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,FiniteDifferenceMonotonicInterpolation,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:359 (repeats 80000 times)

My mistake, I checked again and using x::Number leads to the same stack overflow error. I don't know what I should do.

Yes, Interpolations.gradient(itp, 0.3) works, but Interpolations.gradient(itp, [0.3]) and Interpolations.hessian(itp, [0.3]) do not:

julia> using Interpolations [ Info: Recompiling stale cache file /home/mateusz/.julia/compiled/v1.0/Interpolations/VpKVx.ji for Interpolations [a98d9a8b-a2ab-59e6-89dd-64a1c18fca59] julia> itp = interpolate(0.0:0.1:1.0, 0.0:0.1:1.0, LinearMonotonicInterpolation()); julia> itp(0.0) 0.0 julia> itp(0.2) 0.2 julia> Interpolations.gradient(itp, 0.3) 1-element StaticArrays.SArray{Tuple{1},Float64,1,1}: 1.0 julia> Interpolations.gradient(itp, [0.3]) ERROR: Given number 2 is outside of interpolated range [0.0, 1.0] Stacktrace: [1] (::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})(::Int64) at /home/mateusz/.julia/dev/Interpolations/src/monotonic/monotonic.jl:174 [2] getindex at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:383 [inlined] [3] isassigned(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Int64) at ./abstractarray.jl:351 [4] show_delim_array(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Char, ::String, ::Char, ::Bool, ::Int64, ::Int64) at ./show.jl:659 [5] show_delim_array at ./show.jl:649 [inlined] [6] show_vector(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Char, ::Char) at ./arrayshow.jl:442 [7] show_vector at ./arrayshow.jl:432 [inlined] [8] show at ./arrayshow.jl:418 [inlined] [9] print(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}) at ./strings/io.jl:31 [10] #print_to_string#326(::Nothing, ::Function, ::String, ::Vararg{Any,N} where N) at ./strings/io.jl:124 [11] print_to_string at ./strings/io.jl:112 [inlined] [12] string at ./strings/io.jl:143 [inlined] [13] gradient(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:349 [14] top-level scope at none:0 julia> Interpolations.hessian(itp, [0.3]) ERROR: StackOverflowError: Stacktrace: [1] hessian(::Interpolations.MonotonicInterpolation{Float64,Float64,Float64,LinearMonotonicInterpolation,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}) at /home/mateusz/.julia/dev/Interpolations/src/Interpolations.jl:359 (repeats80000 times) julia>

These error messages aren't helpful.

I just pushed some fixes to your branch.

Thanks! Should I do something about that problem with gradient/hessian with array as an argument?

If you want to support Interpolations.gradient(itp, [0.3]) you need a special method that returns an array of gradients, e.g., a Vector{SVector{1,Float64}}. But none of the other methods here support such a syntax:

julia> itp = interpolate(1:10, BSpline(Linear())) 10-element interpolate(::Array{Float64,1}, BSpline(Linear())) with element type Float64: 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 julia> Interpolations.gradient(itp, 2.3) 1-element StaticArrays.SArray{Tuple{1},Float64,1,1}: 1.0 julia> Interpolations.gradient(itp, [2.3]) ERROR: gradient of [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] not supported for position ([2.3],) Stacktrace: [1] error(::String) at ./error.jl:33 [2] gradient(::Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},BSpline{Linear},Tuple{Base.OneTo{Int64}}}, ::Array{Float64,1}) at /home/tim/.julia/dev/Interpolations/src/Interpolations.jl:349 [3] top-level scope at none:0

so I don't think you have to, either.

OK, so I'll just remove the failing tests for gradient/hessian with array as an argument.

tomasaschan · 2018-09-28T07:18:23Z

src/monotonic/monotonic.jl

@@ -180,7 +180,12 @@ function (itp::MonotonicInterpolation)(x::Number)
    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
 end

-function derivative(itp::MonotonicInterpolation, x::Number)
+function gradient(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})
+    length(x) == 1 || error("Given vector x ($x) should have exactly one element")


With x::VarArg{<:Number,1}, I don't think this check is needed (but I'm not sure).

This check isn't needed with function gradient(itp::MonotonicInterpolation, x::Number) and it works.

Yeah; for monotonic interpolations it doesn't matter very much.

My main concern here is that I want the interface to match gradient for other interpolation types exactly, with the only exception being the type of itp. The idea is that you should be able to have code that initializes an interpolation object of some type, and then (perhaps a lot later in the code) use any part of the common API without caring which type you used. So if you e.g. start out with regular B-splines, and then realize you want monotonic splines, you switch out the initialization part and nothing else.

It might be that x::Number accomplishes that well enough, and if so it's totally fine by me. I just want to make sure that the signature of this method (as well as the others, e.g. itp(<whatever>) and hessian(itp, <whatever>)) is designed with the common API in mind.

Yes, I understand, but Vararg{<:Number,1} didn't work (I wrote about it here: #238 (comment) ). gradient(itp, [1.0]) still works (UnexpandedIndexTypes union contains AbstractVector). Is there another use case I should keep in mind?

Hum. I'll ping in @timholy here...

tomasaschan · 2018-09-28T07:18:59Z

src/monotonic/monotonic.jl

@@ -189,13 +194,27 @@ function derivative(itp::MonotonicInterpolation, x::Number)
    return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*xdiff
 end

+function hessian(itp::MonotonicInterpolation, x::AbstractArray{<:Number, 1})
+    length(x) == 1 || error("Given vector x ($x) should have exactly one element")


Same thing here: x should be VarArg, and the length check is then probably redundant.

No AbstractArray here.

mateuszbaran · 2018-09-28T08:08:18Z

Thank you too for a great library!

Before merging, I'd still like to know what should I do about getindex because right now printing an object of type MonotonicInterpolation may results in an error due to evaluation beyond interpolated range. Gridded just prints #undef in such cases but I don't see how it's done with evaluation outside interpolated range still throwing an error.

tomasaschan · 2018-09-28T08:51:29Z

@timholy Two outstanding questions here:

How should getindex handle out-of-bounds evaluation?
What's the proper way to declare gradient's signature here?

I really don't know what the best way to think here is, so would appreciate your input :)

timholy · 2018-09-28T12:18:50Z

I fixed the getindex problem by throwing a BoundsError. It's not great for indicating the natural bounds of the object, unfortunately; I'm not quite sure what to do about this, short of creating a special exception type.

tomasaschan · 2018-09-28T12:25:24Z

I'd argue that now that we've decided that the "correct" semantics of interpolation is itp(x), rather than itp[x], I think we should carefully evaluate how much value itp[x] gives at all, and whether we should just decide not to inherit AbstractArray at all. As a consequence, I don't think we should spend too much effort on it here.

timholy · 2018-09-28T12:28:02Z

My vote for itp[x] is to interpret x with respect to the parent array indexes, see the last paragraph of #227. But I don't think we change the meaning of [] without a nice long deprecation period.

tomasaschan · 2018-09-28T12:33:56Z

Ah, true. I filed #242 to continue this discussion - it's a little OT here...

They test gradient/hessian with array as an argument.

tomasaschan

I'd say this LGTM now, but I'd like to see @timholy approve it as well before merging.

timholy · 2018-09-28T13:44:32Z

test/gradient.jl

@@ -194,7 +194,6 @@ using Test, Interpolations, DualNumbers, LinearAlgebra
                for t in grid
                    gradval = epsilon(itp(dual(t, 1.0)))
                    @test Interpolations.gradient1(itp, t) ≈ gradval atol = 1.e-12
-                    @test Interpolations.gradient(itp, [t])[1] ≈ gradval atol = 1.e-12


Rather than remove these can you change [t] to t?

timholy · 2018-09-28T13:45:36Z

Once those gradient/hessian tests are re-inserted without the [] then I approve merging this. Thanks in advance, @mateuszbaran!

mateuszbaran · 2018-09-28T14:03:47Z

Tests are re-inserted. Thank you for help 👍 .

timholy · 2018-09-28T14:07:41Z

Just to make sure you understand @mateusbaran, in ForwardDiff a multidimensional position is encoded by a vector but here it is encoded by varargs. This package follows indexing rules where vectors along each dimension means a rectangular sub-region of an array, i.e. the Cartesian product of all positions along each supplied axis.

Whether we stick with that will depend on whether we continue to make these subtypes of AbstractArray.

mateuszbaran · 2018-09-28T14:39:28Z

I understand now. I didn't pay attention to it earlier because I've used this library only for one-dimensional interpolation.

timholy · 2018-09-29T08:37:38Z

Merged as #243. Thanks for the great contribution @mateuszbaran!

tomasaschan · 2018-10-01T08:07:33Z

Great work here, @mateuszbaran. Big thanks!

mateuszbaran added 4 commits September 25, 2018 20:03

basic implementation of monotonic interpolation methods

c3d3d0b

implemented algorithms: * linear * finite difference * cardinal * Fritsch-Carlson * Fritsch-Butland * Steffen

argument checks for monotonic interpolation

5854513

derivatives for monotonic interpolation

28ad194

Documentation for monotonic interpolation

1bfd8d2

faster knot searching in monotonic interpolation

9d65102

tomasaschan approved these changes Sep 26, 2018

View reviewed changes

minor changes in monotonic interpolation

9c687b4

* additional documentation for interpolation coefficients * references added to `calcTangents` for Fritsch & Butland and Steffen methods * changed the type parameter name `K` to `TKnots` * changed "b-spline" to "spline" in error string

tomasaschan suggested changes Sep 27, 2018

View reviewed changes

more coherent type parameter naming in monotonic interpolation

c3763e3

timholy reviewed Sep 27, 2018

View reviewed changes

mateuszbaran added 2 commits September 27, 2018 19:11

new gradient, gradient1, hessian and hessian1 function for monotonic …

89dec3b

…interpolation

Merge branch 'master' into monotonic-interpolations

a19eff7

# Conflicts: # test/gradient.jl

tomasaschan reviewed Sep 28, 2018

View reviewed changes

Small fixes for Monotonic indexing rules

514185a

tomasaschan mentioned this pull request Sep 28, 2018

Drop AbstractArray inheritance? #242

Open

Removal of two incorrect tests for monotonic interpolation.

a2144cf

They test gradient/hessian with array as an argument.

tomasaschan approved these changes Sep 28, 2018

View reviewed changes

timholy reviewed Sep 28, 2018

View reviewed changes

Re-insertion of fixed gradient/hessian tests for monotonic interpolation

e6a38cb

timholy mentioned this pull request Sep 28, 2018

basic implementation of monotonic interpolation methods #243

Merged

timholy closed this Sep 29, 2018

Monotonic interpolations #238

Monotonic interpolations #238

Conversation

mateuszbaran commented Sep 26, 2018

codecov-io commented Sep 26, 2018 • edited Loading

Codecov Report

tomasaschan left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mateuszbaran Sep 27, 2018 • edited Loading

Choose a reason for hiding this comment

mateuszbaran commented Sep 26, 2018

tomasaschan left a comment

Choose a reason for hiding this comment

timholy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mateuszbaran commented Sep 27, 2018 • edited Loading

tomasaschan left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mateuszbaran Sep 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy Sep 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mateuszbaran commented Sep 28, 2018

tomasaschan commented Sep 28, 2018 • edited Loading

timholy commented Sep 28, 2018

tomasaschan commented Sep 28, 2018

timholy commented Sep 28, 2018

tomasaschan commented Sep 28, 2018

tomasaschan left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy commented Sep 28, 2018

mateuszbaran commented Sep 28, 2018

timholy commented Sep 28, 2018

mateuszbaran commented Sep 28, 2018

timholy commented Sep 29, 2018

tomasaschan commented Oct 1, 2018

codecov-io commented Sep 26, 2018 •

edited

Loading

mateuszbaran Sep 27, 2018 •

edited

Loading

mateuszbaran commented Sep 27, 2018 •

edited

Loading

mateuszbaran Sep 28, 2018 •

edited

Loading

timholy Sep 28, 2018 •

edited

Loading

tomasaschan commented Sep 28, 2018 •

edited

Loading