ENH: add geomspace function #7268

endolith · 2016-02-17T04:25:53Z

Like logspace, but specifying start and stop directly, without having to
specify the base.

Fixes #7255

njsmith · 2016-02-17T04:39:01Z

Please add link to mailing list discussion for future reference.

endolith · 2016-02-18T00:57:14Z

Added tests, please review

madphysicist · 2016-02-18T01:35:16Z

numpy/core/function_base.py

+    --------
+    logspace : Similar to geomspace, but with endpoints specified using log
+               and base.
+    linspace : Similar to geomspace, but with the samples uniformly distributed


I think "Similar to geomspace but with arithmetic instead of geometric sequence."
would be better.

ok but this was copied from the other functions, should those change too?

I would say let someone else weigh in, preferably a core dev, then make all of them consistent one way or the other.

changed only in geomspace

madphysicist · 2016-02-18T02:24:52Z

I have a few comments, but overall +1 for this. This is much more intuitive than logspace.

endolith · 2016-02-20T04:09:32Z

Mailing list proposals:

Different name?
- I think geomspace is great. (SpacePy has the same function already)
np.linspace(start, stop, steps, spacing="log")
- I really don't like this. Just change the function name to logspace if you want log space.
Specify number of intervals rather than number of steps?
- maybe, but would be confusing coming from existing functions
option for geomspace(start=1, finish=2, ratio=1.03)
- I agree that this is useful, or geomspace(start=1, num=7, ratio=1.03). However...
"call a function with this signature geomrange() instead"
- Agreed. Separate functions is better than a single Swiss army knife function that tries to do everything with mutually-incompatible parameters. This PR will just be about one function to keep it simple.

endolith · 2016-02-20T04:26:29Z

I added support for negative start and stop, but this breaks some complex start and stop, and I'm not sure whether those should be supported anyway. An implementation that doesn't call logspace could support both I guess.

madphysicist · 2016-02-20T14:51:19Z

+1 to two functions. I pretty much agree with the whole explanatory post.

charris · 2016-02-20T17:12:29Z

How about expspace as a name? AFAICT, what we have is basically exp(a*linspace(...)). As a generalization, if f is monotone, f(a*linspace(...), so there could be a general underlying function with special cases, the only thing really needed is the ability to invert f to get the endpoints.

endolith · 2016-03-07T13:45:12Z

Ok I modified it to handle negative and complex arguments. With complex log, it actually works out without any special-casing, but I special-cased negative arguments so the output doesn't have negligible imaginary parts. I could also special-case purely-imaginary arguments so they don't have negligible real parts.

Currently:

good:

geomspace(1, 1000, 4)
 array([    1.,    10.,   100.,  1000.])

geomspace(-1, -1000, 4)
 array([   -1.,   -10.,  -100., -1000.])

geomspace(1, 16, 5)
 array([  1.,   2.,   4.,   8.,  16.])

geomspace(1, 1000, 4, dtype=int)
 array([   1,   10,  100, 1000])

geomspace(-1, -1000, 4, dtype=int)
 array([   -1,   -10,  -100, -1000])

reasonable:

geomspace(1, -1, 5)  # dtype derived from real arguments
 array([ nan,  nan,  nan,  nan,  nan])

geomspace(1+0j, -1, 5)  # dtype derived from complex argument
array([  1.00000000e+00 +0.00000000e+00j,
         7.07106781e-01 +7.07106781e-01j,
         6.12323400e-17 +1.00000000e+00j,
        -7.07106781e-01 +7.07106781e-01j,  -1.00000000e+00 +1.22464680e-16j])

geomspace(1, -1+0j, 5)  # dtype derived from complex argument
array([  1.00000000e+00 +0.00000000e+00j,
         7.07106781e-01 +7.07106781e-01j,
         6.12323400e-17 +1.00000000e+00j,
        -7.07106781e-01 +7.07106781e-01j,  -1.00000000e+00 +1.22464680e-16j])

geomspace(1, -1, 5, dtype=complex)  # dtype specified
array([  1.00000000e+00 +0.00000000e+00j,
         7.07106781e-01 +7.07106781e-01j,
         6.12323400e-17 +1.00000000e+00j,
        -7.07106781e-01 +7.07106781e-01j,  -1.00000000e+00 +1.22464680e-16j])

geomspace(1j, 1000j, 4)
Out[250]: 
array([  6.12323400e-17 +1.00000000e+00j,
         6.12323400e-16 +1.00000000e+01j,
         6.12323400e-15 +1.00000000e+02j,   6.12323400e-14 +1.00000000e+03j])

geomspace(-1j, -1000j, 4)
Out[251]: 
array([ -6.12323400e-17 -1.00000000e+00j,
        -6.12323400e-16 -1.00000000e+01j,
        -6.12323400e-15 -1.00000000e+02j,  -6.12323400e-14 -1.00000000e+03j])

bad:

geomspace(1, 16, 5, dtype=int)
 array([ 1,  2,  4,  7, 16])
                     ^

probably should check if dtype is int and use around() first?

madphysicist · 2016-03-07T13:55:36Z

With regards to your "bad" example, you are not obligated to return an array of the same type as the input, and I would say that you definitely should not. Many numpy functions (e.g. percentile comes to mind) will return an array of a reasonable type able to hold the result. In your case, returning a float64 real array for integer inputs is probably not an unreasonable solution. Returning an integer array is unreasonable in almost any case for this function for the reason shown above.

endolith · 2016-03-07T15:01:17Z

the bad example explicitly states dtype=int, though. if dtype is not specified, it's either float or complex output, depending on types of start and stop.

madphysicist · 2016-03-07T15:51:43Z

Looks like you will have to do your own rounding then.

seberg · 2016-03-07T16:56:30Z

You could also just kill the dtype argument ;). This is unsafe casting, ufuncs will nowadays not allow you to use dtype=int, i.e. for log it will tell you that there is no loop for this type and it will tell you it is not possible for float + float = int, unless you also specify casting=unsafe (though it will use float internally then).

madphysicist · 2016-03-07T17:07:05Z

Killing dtype is probably better. Here is another bad case: Inputs are complex, output dtype is real. around won't help with that much.

njsmith · 2016-03-07T17:19:39Z

Output dtype argument is very useful for dtype=float32. Making dtype=int
raise an error is totally reasonable, though.
On Mar 7, 2016 9:07 AM, "Mad Physicist" [email protected] wrote:

Killing dtype is probably better. Here is another bad case: Inputs are
complex, output dtype is real. around won't help with that much.

—
Reply to this email directly or view it on GitHub
#7268 (comment).

endolith · 2016-03-07T17:45:19Z

I don't understand why dtype=int should raise an error. If someone wants a geometric integer progression, why can't they have it? Why shouldn't geomspace(2, 54, 4, dtype=int) output array(2, 6, 18, 54])?

I just need to add a line like:

if np.issubdtype(dtype, np.integer):
    result = _nx.around(result)

madphysicist · 2016-03-07T20:32:55Z

Because someone can also ask for geomspace(2, 7, 4, dtype=int) or geomspace(2, 7+3j, 4, dtype=int) and guessing how they want that is unpythonic. Let the user call astype or around for themselves as they see fit.

njsmith · 2016-03-07T20:53:31Z

If someone wants an integer progression of course they can have it -- the question is whether dtype=int should magically do it, or whether they should just call around themselves if that's what they want :-). I don't think there are any other places in numpy where we give ints a special-case casting mechanism (which is basically what this rounding proposal is), so it's a bit surprising and I doubt people will realize the feature is there...

njsmith · 2016-03-07T20:54:50Z

(To be clear I'm not dead-set against the dtype=int thing, I'm just not 100% convinced and leaning towards erring on the side of conservatism.)

endolith · 2016-03-08T01:53:14Z

I modified it so that purely-imaginary sequences are handled better:

geomspace(1j, 1000j, 4)
Out[376]: array([ 0.   +1.j,  0.  +10.j,  0. +100.j,  0.+1000.j])

geomspace(-1j, -1000j, 4)
Out[377]: array([ 0.   -1.j,  0.  -10.j,  0. -100.j,  0.-1000.j])

I didn't add around().

endolith · 2016-03-08T01:59:18Z

numpy/core/function_base.py

+    result = out_sign * logspace(log_start, log_stop, num=num,
+                                 endpoint=endpoint, base=10.0, dtype=dtype)
+
+    return result.astype(dtype)


is this needed? could probably just return on the above line, since dtype will be handled by logspace and out_sign?

endolith · 2016-03-08T03:55:56Z

Added a bunch more tests, changed the type promotion according to http://stackoverflow.com/q/2875024/125507 so the PhysicalQuantities test passes.

homu · 2016-03-13T15:31:40Z

☔ The latest upstream changes (presumably #7328) made this pull request unmergeable. Please resolve the merge conflicts.

endolith · 2016-03-14T16:43:42Z

I still need to make an example and squash some commits

endolith · 2016-03-19T03:35:34Z

Ok I changed the example and squashed some commits, please check it

charris · 2016-06-17T17:27:48Z

It would be good if the folks involved in this would take it through to completion.

EDIT: otherwise I'm going to drop the 1.12.0 milestone.

endolith · 2016-06-18T00:08:44Z

Anything else I need to do? I added some more examples that will hopefully pass.

charris · 2016-06-18T14:55:49Z

Might squash the commits. Two commits, implementation and tests, or just one would do. I'm not a big fan of the name but no one else seems bothered by it.

Like logspace, but specifying start and stop directly, without having to specify the base. Purely imaginary or purely negative real sequences are converted to positive real, computed, and converted back, so there are no negligible real or imaginary parts. Instead of array([ 6.12323400e-17 +1.00000000e+00j, 6.12323400e-16 +1.00000000e+01j, 6.12323400e-15 +1.00000000e+02j, 6.12323400e-14 +1.00000000e+03j]) it outputs array([ 0. +1.j, 0. +10.j, 0. +100.j, 0.+1000.j]) If dtype is complex it does the math in complex so it can leave the real line and follow a spiral. TST: Added tests for geomspace and more tests for logspace, making PhysicalQuantities tests work for all types of functions PEP8: __all__ after imports, line wrapping

endolith · 2016-06-21T03:14:43Z

@charris Ok I squashed the commits.

@madphysicist

Because someone can also ask for geomspace(2, 7, 4, dtype=int) or geomspace(2, 7+3j, 4, dtype=int) and guessing how they want that is unpythonic.

Yes I agree now. It would be nice if it output exact integer values when it was mathematically supposed to, though. But that could be added in another PR. I documented it in the Examples for now.

charris · 2016-06-21T13:34:45Z

Need to add an entry in doc/source/reference/routines.array-creation.rst where there is already linspace and logspace, and in the 1.12.0 release notes. Otherwise, LGTM, although as @madphysicist says, the more tests the better.

Comment when ready so that github will send out a notification.

charris · 2016-06-21T17:55:11Z

Needs resquash. You can move the update commit and squash it into the initial commit.

endolith · 2016-06-21T18:05:44Z

ok I did that but I can't squash the commits right now

endolith · 2016-06-21T23:19:23Z

ok resquashed

charris · 2016-06-22T00:55:33Z

Thanks @endolith .

endolith · 2017-06-24T04:07:36Z

Also should this be tweaked so that the start and stop points are exactly equal to the function arguments? https://stackoverflow.com/q/44136851/125507

ashwinvis · 2017-11-08T09:52:51Z

@endolith Also this issue #9987

madphysicist reviewed Feb 18, 2016
View reviewed changes

charris added 01 - Enhancement component: numpy.lib labels Feb 19, 2016

endolith reviewed Mar 8, 2016
View reviewed changes

endolith force-pushed the geomspace branch 2 times, most recently from 43bb2b3 to cfb4803 Compare March 19, 2016 13:08

charris added this to the 1.12.0 release milestone Apr 16, 2016

endolith force-pushed the geomspace branch from cfb4803 to f551625 Compare June 17, 2016 23:50

endolith force-pushed the geomspace branch from f551625 to c967361 Compare June 21, 2016 03:13

DOC: Add geomspace to function list, release notes

e53570d

endolith force-pushed the geomspace branch from d377e55 to e53570d Compare June 21, 2016 23:19

charris merged commit 70f50b8 into numpy:master Jun 22, 2016

homu mentioned this pull request Jun 22, 2016

BUG: Let linspace accept input that has an array_interface. #6659

Merged

This was referenced Jun 22, 2016

ENH,TST: Add new warning suppression/filtering context #7763

Closed

Update 5057, ENH: Add a block function for creating stacked block arrays. #7768

Closed

TST: Suppressed warnings #7099

Merged

endolith deleted the geomspace branch June 22, 2016 01:06

endolith mentioned this pull request Jun 22, 2016

DOC: link geomspace from logspace #7771

Merged

endolith mentioned this pull request Nov 8, 2017

Wrong output while using geomspace with dtype=int keyword argument #9987

Open

endolith mentioned this pull request Jan 28, 2024

BUG: branch choices for geomspace on complex arguments #25644

Open

Uh oh!

ENH: add geomspace function #7268

ENH: add geomspace function #7268

Uh oh!

Conversation

endolith commented Feb 17, 2016

Uh oh!

njsmith commented Feb 17, 2016

Uh oh!

endolith commented Feb 18, 2016

Uh oh!

madphysicist Feb 18, 2016

Choose a reason for hiding this comment

Uh oh!

endolith Feb 18, 2016

Choose a reason for hiding this comment

Uh oh!

madphysicist Feb 18, 2016

Choose a reason for hiding this comment

Uh oh!

endolith Feb 20, 2016

Choose a reason for hiding this comment

Uh oh!

madphysicist commented Feb 18, 2016

Uh oh!

endolith commented Feb 20, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

endolith commented Feb 20, 2016

Uh oh!

madphysicist commented Feb 20, 2016

Uh oh!

charris commented Feb 20, 2016

Uh oh!

endolith commented Mar 7, 2016

Uh oh!

madphysicist commented Mar 7, 2016

Uh oh!

endolith commented Mar 7, 2016

Uh oh!

madphysicist commented Mar 7, 2016

Uh oh!

seberg commented Mar 7, 2016

Uh oh!

madphysicist commented Mar 7, 2016

Uh oh!

njsmith commented Mar 7, 2016

Uh oh!

endolith commented Mar 7, 2016

Uh oh!

madphysicist commented Mar 7, 2016

Uh oh!

njsmith commented Mar 7, 2016

Uh oh!

njsmith commented Mar 7, 2016

Uh oh!

endolith commented Mar 8, 2016

Uh oh!

endolith Mar 8, 2016

Choose a reason for hiding this comment

Uh oh!

endolith commented Mar 8, 2016

Uh oh!

homu commented Mar 13, 2016

Uh oh!

endolith commented Mar 14, 2016

Uh oh!

endolith commented Mar 19, 2016

Uh oh!

charris commented Jun 17, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

endolith commented Jun 18, 2016

Uh oh!

charris commented Jun 18, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

endolith commented Jun 21, 2016

Uh oh!

charris commented Jun 21, 2016

endolith commented Feb 20, 2016 •

edited

Loading

charris commented Jun 17, 2016 •

edited

Loading

charris commented Jun 18, 2016 •

edited

Loading