*Matlab behavior:*
>> format long
>> x = 0:pi/100:pi;
>> length(x)
ans =
101
>> x(1) - 0
ans =
0
>> x(end) - pi
ans =
0
>> diff([max(diff(x)), min(diff(x))])
ans =
-4.440892098500626e-16
*Julia behavior:*
julia> x = 0:pi/100:pi;
julia> length(x)
100
julia> x[1]-0
0.0
julia> x[end]-pi
-0.031415926535897754
julia> diff([maximum(diff(x)), minimum(diff(x))])
1-element Array{Float64,1}:
-4.44089e-16
I highlighted the important differences in red. IMO the Matlab behavior is
more intuitive. If you choose an increment that "very nearly" (apparently
within a few ulp, as stated by Stefan) evenly divides the difference
between the first and last entries in the colon notation, Matlab assumes
that you want the resulting vector to begin and end with exactly these
values. This makes sense in most situations, IMO. E.g. when using a range
to specify sample points in a numerical integration scheme, or when using a
range to specify x-coordinates in an interpolation scheme. In both of
these cases you don't want the last point to be adjusted away from what you
specified using the colon notation.
--Pete
--Peter
On Thursday, April 24, 2014 5:28:19 AM UTC-7, andrew cooke wrote:
>
>
> how does matlab differ from julia here? (i though the original problem
> was related to "fundamental" issues of how you represent numbers on a
> computer, not language specific).
>
>
> On Thursday, 24 April 2014 02:24:59 UTC-3, Peter Simon wrote:
>>
>> Thanks, Simon, that construct works nicely to solve the problem I posed.
>>
>> I have to say, though, that I find Matlab's colon range behavior more
>> intuitive and generally useful, even if it isn't as "exact" as Julia's.
>>
>> --Peter
>>
>> On Wednesday, April 23, 2014 7:17:23 PM UTC-7, Simon Kornblith wrote:
>>>
>>> pi*(0:0.01:1) or similar should work.
>>>
>>> On Wednesday, April 23, 2014 7:12:58 PM UTC-4, Peter Simon wrote:
>>>>
>>>> Thanks for the explanation--it makes sense now. This question arose
>>>> for me because of the example presented in
>>>> https://groups.google.com/d/msg/julia-users/CNYaDUYog8w/QH9L_Q9Su9YJ :
>>>>
>>>> x = [0:0.01:pi]
>>>>
>>>> used as the set of x-coordinates for sampling a function to be
>>>> integrated (ideally over the interval (0,pi)). But the range defined in x
>>>> has a last entry of 3.14, which will contribute to the inaccuracy of the
>>>> integral being sought in that example. As an exercise, I was trying to
>>>> implement the interpolation solution described later in that thread by
>>>> Cameron McBride: "BTW, another possibility is to use a spline
>>>> interpolation on the original data and integrate the spline evaluation
>>>> with quadgk()". It seems that one cannot use e.g. linspace(0,pi,200) for
>>>> the x values, because CoordInterpGrid will not accept an array as its
>>>> first
>>>> argument, so you have to use a range object. But the range object has a
>>>> built-in error for the last point because of the present issue. Any
>>>> suggestions?
>>>>
>>>> Thanks,
>>>>
>>>> --Peter
>>>>
>>>> On Wednesday, April 23, 2014 3:24:10 PM UTC-7, Stefan Karpinski wrote:
>>>>>
>>>>> The issue is that float(pi) < 100*(pi/100). The fact that pi is not
>>>>> rational – or rather, that float64(pi) cannot be expressed as the
>>>>> division
>>>>> of two 24-bit integers as a 64-bit float – prevents rational lifting of
>>>>> the
>>>>> range from kicking in. I worried about this kind of issue when I was
>>>>> working on FloatRanges, but I'm not sure what you can really do, aside
>>>>> from
>>>>> hacks where you just decide that things are "close enough" based on some
>>>>> ad
>>>>> hoc notion of close enough (Matlab uses 3 ulps). For example, you can't
>>>>> notice that pi/(pi/100) is an integer – because it isn't:
>>>>>
>>>>> julia> pi/(pi/100)
>>>>> 99.99999999999999
>>>>>
>>>>>
>>>>> One approach is to try to find a real value x such that float64(x/100)
>>>>> == float64(pi)/100 and float64(x) == float64(pi). If any such value
>>>>> exists,
>>>>> it makes sense to do a lifted FloatRange instead of the default naive
>>>>> stepping seen here. In this case there obviously exists such a real
>>>>> number
>>>>> – π itself is one such value. However, that doesn't quite solve the
>>>>> problem
>>>>> since many such values exist and they don't necessarily all produce the
>>>>> same range values – which one should be used? In this case, π is a good
>>>>> guess, but only because we know that's a special and important number.
>>>>> Adding in ad hoc special values isn't really satisfying or acceptable. It
>>>>> would be nice to give the right behavior in cases where there is only one
>>>>> possible range that could have been intended (despite there being many
>>>>> values of x), but I haven't figured out how determine if that is the case
>>>>> or not. The current code handles the relatively straightforward case
>>>>> where
>>>>> the start, step and stop values are all rational.
>>>>>
>>>>>
>>>>> On Wed, Apr 23, 2014 at 5:59 PM, Peter Simon <psimo...@gmail.com>wrote:
>>>>>
>>>>>> The first three results below are what I expected. The fourth result
>>>>>> surprised me:
>>>>>>
>>>>>> julia> (0:pi:pi)[end]
>>>>>> 3.141592653589793
>>>>>>
>>>>>> julia> (0:pi/2:pi)[end]
>>>>>> 3.141592653589793
>>>>>>
>>>>>> julia> (0:pi/3:pi)[end]
>>>>>> 3.141592653589793
>>>>>>
>>>>>> julia> (0:pi/100:pi)[end]
>>>>>> 3.1101767270538954
>>>>>>
>>>>>> Is this behavior correct?
>>>>>>
>>>>>> Version info:
>>>>>> julia> versioninfo()
>>>>>> Julia Version 0.3.0-prerelease+2703
>>>>>> Commit 942ae42* (2014-04-22 18:57 UTC)
>>>>>> Platform Info:
>>>>>> System: Windows (x86_64-w64-mingw32)
>>>>>> CPU: Intel(R) Core(TM) i7 CPU 860 @ 2.80GHz
>>>>>> WORD_SIZE: 64
>>>>>> BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
>>>>>> LAPACK: libopenblas
>>>>>> LIBM: libopenlibm
>>>>>>
>>>>>>
>>>>>> --Peter
>>>>>>
>>>>>>
>>>>>
