These are chat archives for SHTOOLS/SHTOOLS

5th
Aug 2016
MMesch
@MMesch
Aug 05 2016 07:52
@MarkWieczorek ok. Tell me if you want me to merge it
MMesch
@MMesch
Aug 05 2016 08:09
for reference, here is explained how to upload the stuff to pypi: https://packaging.python.org/distributing/#uploading-your-project-to-pypi
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 08:30
It will be done by noon. Just cleaning things up...
MMesch
@MMesch
Aug 05 2016 09:08
ok. I need to go through it again anyway because I want to add a few things
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 09:10
@MMesch: I think you also need to add examples/notebooks in MANIFEST.in
MMesch
@MMesch
Aug 05 2016 10:04
yes
MMesch
@MMesch
Aug 05 2016 10:15
shtools_logo.svg.png
something like this for the logo ?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 11:13
@MMesch did you modifiy something in src/SHBias.f95?
its just white space. no worries.
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 11:21
ok. just put in a pull request. I'm going to lunch now, so if you are happy, feel free to merge it. Otherwise, I will do it in about 45 minutes. I made some modifications to the README.md file, but I didn't have time to preview it. Feel free to update it.
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 12:33
who made the logo?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 12:50
I have two comments about the setup.py file. First, I just tried it, and it only installs the python2 version, and not the python3 version.
Second, is there a "clean" command that will remove all the junk files (like in build and pyshtools.egg-info)?
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 12:52
What command do you use for setup.py?
Because it installs
the version from which it called
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:01
ok, it works fine with pip3. And is there an uninstall? pip uninstall . doesn't seem to do anything.
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 13:03
Yes, you need pip uninstall pyshtools
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:06
If I understand, that just removes " /usr/local/lib/python3.5/site-packages/pyshtools.egg-link". Is there a way to bring the original folder back to its initial state?
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 13:08
Did you I install with -e flag?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:10
yeah.
The problem is that "make clean" doesn't look for the files created by pip. I could add them to make clean if I have too.
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 13:13
It can be added to setup.py as well. setup.py has clean command, but it removes only temporary file in build/tmp
And it shouldn't, actually, because egg file is still in the system
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 13:19
An linked to local dir if you using -e flag
If you install package, i.e. without -e flag, all files from build/lib... are copied to /usr/lib/python...
And the local dir is clean because because building were made in the system /tmp dir
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:22
also, this seems to create the files src/_SHTOOLS-f2pywrappers.f and src/_SHTOOLSmodule.c.
Ok. I'll just add these file to make clean, with the understanding that you should also do a pip uninstall as well. The other option would be to just add these files to gitignore.
Ilya Oshchepkov
@ioshchepkov
Aug 05 2016 13:29
Also build dir by default is ignored from sdist command which creates distribution archive
MMesch
@MMesch
Aug 05 2016 13:32
@MarkWieczorek the coupling matrix is incorrect again. Did you undo the change?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:33
I though we agreed that it was right the way it was before?
MMesch
@MMesch
Aug 05 2016 13:33
no. It is wrong
The coupling matrix should map input power -> output power , like a convolution does in the ordinary Fourier transform
that's why the output power has size lmax + lwin + 1 and the input power has size lmax
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:35
but you said I was right the other day...
MMesch
@MMesch
Aug 05 2016 13:35
yeah. I mean it is just the transpose
I don't know. I'll check your paper again
The question is, do we want the coupling matrix to map: unbiased -> biased or biased-> unbiased ?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:40
I think that the matrix is correct, and should not be transposed, but that at the same time, you could modify its dimensions.
The way I look at this is S_PHI(L) depends upon S_phi(L+Lwin). Since S_phi(L+Lwin+1) = UNKNOWN, S_PHI(L+1) is undefined.
MMesch
@MMesch
Aug 05 2016 13:41
\partial test
Slbiased=MllSlunbiasedS_l^{biased} = M_{ll'} S_{l'}^{unbiased}
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:42
If you assume S_phi is zero beyond L+Lwin, we could modify the dimensions
MMesch
@MMesch
Aug 05 2016 13:45
To me the result of the convolution should be longer than the original spectrum. It makes more sense
If you don't include the boundaries, you get what you mean
but it is less general
if it is longer, the user can just cut it to the length he wants.
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:46
so, lets say that Lwin =10, and that Lmax of the data is 100. L=101 of the data is thus NaN. What is the largest degree of Sl_biased that is not NaN?
MMesch
@MMesch
Aug 05 2016 13:47
yes you are right, but in the other case it is just assumed zero
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:48
Sometimes you can't make that assumption though. If you do, I agree with you and furthermore, that you can then invert for S_unbiassed.
MMesch
@MMesch
Aug 05 2016 13:48
Take this example: I generate data with l max= 100 , then I window and compute the spectrum again, I will have degrees > lmax - lwin .
does the coupling_matrix predict them correctly if we assume that it is zero beyond the grid resolution !?
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:50
I would say that those degrees only make sense if the power of the original data was exactly zero between 101 and 101+lwin.
MMesch
@MMesch
Aug 05 2016 13:50
ok so I am really not sure how we should define it. As I said, I just thought about it as equivalent to a convolution, which often returns a vector which is longer than the original one
assuming zeros outside
mode : {‘full’, ‘valid’, ‘same’}, optional

    ‘full’:

        By default, mode is ‘full’. This returns the convolution at each point of overlap, with an output shape of (N+M-1,). At the end-points of the convolution, the signals do not overlap completely, and boundary effects may be seen.
    ‘same’:

        Mode ‘same’ returns output of length max(M, N). Boundary effects are still visible.
    ‘valid’:

        Mode ‘valid’ returns output of length max(M, N) - min(M, N) + 1. The convolution product is only given for points where the signals overlap completely. Values outside the signal boundary have no effect.
we are discussing 'full' against 'valid'
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:53
In that equation, they have the limits of the sum going to +-infinity!
MMesch
@MMesch
Aug 05 2016 13:53
yes. I agree. It is just a matter of definition
but we should choose an intelligent one. If we return 'full', we can just cut the matrix in python to 'valid'
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:54
we could add the same options.
MMesch
@MMesch
Aug 05 2016 13:54
make the Fortran routine return 'full' and we cut in python. That's easy
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 13:55
For me (until GRAIL), I never really worked with a field with ldata > 130. I've thus always needed to make sure that I was only using the valid part.
One thing to consider is that SHBias only returns the valid coefficients. Also, SHBias does not directly call SHCouplingMatrix. Thus, its not too big of a deal to me if we modify this routine.
MMesch
@MMesch
Aug 05 2016 13:58
coupling.png
# compute the biased spectrum using two different functions
from pyshtools.localizedspectralanalysis import SHBiasK
power_local2 = np.dot(power, coupling_matrix)
power_local3 = SHBiasK(tapers, power)

# plot the observed and predicted spectra
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(ls, power_global, label='measured global')
ax.plot(ls, power_local, label='measured local')
ax.plot(power_local2, label='predicted local (CouplingMatrix)')
ax.plot(power_local3, 'x', label='predicted local (SHBiasK)')
ax.legend(loc=3)
ax.set_xscale('log', basex=2)
ax.set_yscale('log', basey=2)
SHBiasK doesn't return only valid points
oh and SHBiasK and coupling matrix match exactly if we change it back
I vote heavily for 'full' and 'valid' as options and to change the fortran routine to 'full'
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:01
Actually, I'm wrong. According to the man page " It is implicitly assumed that the power spectrum of inspectrum is zero beyond degree ldata."
MMesch
@MMesch
Aug 05 2016 14:02
ok. I'll change it back then if you agree. We add the 'valid' option to the python code
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:02
Ok. but I don't think that this is a transpose. I think that we just need to extend the dimensions of the matrix in one direction.
MMesch
@MMesch
Aug 05 2016 14:03
this works as well. But it is the transpose also. Property of the Wigner symbols I think.
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:04
I don't think its symmetric, but I could be wrong.
MMesch
@MMesch
Aug 05 2016 14:04
size(inpower) = lmax + 1, size(outpower) = lmax + lwin + 1
Mmt.shape = (lmax + lwin + 1, lmax + 1)
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:05
But I don't think that means you just take the transpose and get the correct matrix.
MMesch
@MMesch
Aug 05 2016 14:05
it is surprising, but yes, it is the same
but for clarity, we can change the dimensions. It is better for the multiplication afterwards anyway
can you change it? I can do it as well if you want ...
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:09
Given that SHBias assumes zero padding, I think that we should do the same for the coupling matrix. We will have to update the documentation to make this explicit. As for the transpose, I'd like to see of a proof before doing this, otherwise, I would just resize it.
MMesch
@MMesch
Aug 05 2016 14:10
I'll resize it. I have saved the 'before' image. I'll show you the 'after' image and you'll see (hopefully I am right, haha)
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:14
The coupling matrix is not equal to its traspose. I just plotted it and it is obvious for the low degrees.
MMesch
@MMesch
Aug 05 2016 14:15
no it is not equal to its' transpose, but it is equal to its transpose multiplied with (2l'+1)/(2l+1)
the (2l+1) factor makes it asymmetric. You can have it exactly symmetric as well if you consider sqrt(2l + 1)
maybe I'll write a notebook about it ...
It can help for the inversion
that's why it was equal to the transpose when I exchanged the (2 i + 1) by (2 j + 1)
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:22
Let me think about this 15 minutes. I can change this, but I need to update the documentation files, describe what the assumptions are in the fortran file, and update the wrapper files.
MMesch
@MMesch
Aug 05 2016 14:30
I changed it already. Just the documentation needs to be adapted
If you want I'll push
It works, I tested it with all notebooks
This message was deleted
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 14:37
Can you put this on a separate branch? I don't want to push this to develop until all the docs are updated. I also want to review this and make sure that it gives the same result as SHBias.
MMesch
@MMesch
Aug 05 2016 14:41
I did the check. It is in the notebooks
even several
if I push to develop, you can run all the notebooks.
I am pretty sure it is correct, but otherwise we can revert the commit.
wait I'll make a branch with the test
I need it for the notebooks, which is why I am pushing, sorry :) pushing not in the github sense :)
here is the benchmark
MMesch
@MMesch
Aug 05 2016 14:46
coupling.png
MMesch
@MMesch
Aug 05 2016 15:01
Screenshot_2016-08-05_17-01-28.png
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 15:17
I need to go. Could you put this all on a separate branch, and I'll look at it tonight or tomorrow and merge it? I just want to make sure that all the documentation is updated, because as you saw, the only way I can remember how the code works is by the documentation. I'd like to put a little effort into making this as clear as possible, because there is a real risk that an idiot will use this incorrectly. I should probably update the documentation for SHBias as well.
MMesch
@MMesch
Aug 05 2016 15:20
ok I'll do so.
be aware that develop is not correct though
I'll make a PR
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 15:21
Otherwise, let me know what you want to do concerning a release. I would like to put something out by the end of next week when I leave on vacation. I suspect that we will have another release that will follow rather quickly when we update the notebooks, pip installer (PyPi), fortran STOP. These are loose that bother me a bit, but it is in a better state than the previous version.
MMesch
@MMesch
Aug 05 2016 15:22
the notebooks are ok once the coupling matrix is updated. Someone needs to read through them and check the texts
that's why I wonder if I don't just update develop a last time to bring it on a working level
and then continue with PR's
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 15:23
I still think that develop is "correct", but that the dimensions of the matrix aren't consistent with what is done is SHBias. We are just adding extra rows/columns....
MMesch
@MMesch
Aug 05 2016 15:24
the plot_couplingmatrix axis labels are inversed
so it is not completely correct.
Mark Wieczorek
@MarkWieczorek
Aug 05 2016 15:25
Just update the notebooks on that branch and I'll merge them all together. I'm done for the day, need to go.
MMesch
@MMesch
Aug 05 2016 15:25
ok, see you !
MMesch
@MMesch
Aug 05 2016 17:51
I am done for now. The PR contains all of my additions for the next release. Everything should hopefully work but needs to be checked. Especially the notebooks, that test most of the class interface. Fortran STOP is the last major issue for the python library in my opinion. Otherwise I think we are definitely ready for a real python release on all channels ...