pycs.gen package

Submodules

pycs.gen.lc module

Module containing all we need to manipulate lightcurves. Could ideally be used by all our curve-shifting algorithms.

class pycs.gen.lc.lightcurve(telescopename='Telescope', object='Object', plotcolour='red')

A class to handle light curves of lensed QSO images.

The actual information is stored into independent “vectors” (1D numpy arrays) : one for time, one for magnitude… It would seem nicer to have some kind of iterable lightcurve object, using perhaps python lists/dicts, but this would be slower to manipulate with numpy I guess.

An exception to this is the properties attribute : a list of dicts, one for each data point, usually all containing the same fields.

Note that a lightcurve MUST always be chronologically ordered, this allows faster algorithms. See methods sort() and validate(). Importation functions should automatically sort().

There is only one remaining philosophy about handling shifts of lightcurves :

  • use shifttime() and cie, this is fast and you do not change the actual data (recommended)
__init__(telescopename='Telescope', object='Object', plotcolour='red')

The plain constructor, binding all the instance variables. It gives a small default lightcurve with 5 points, labels, and a mask to play with.

@todo: nothing to do, that’s just a test of the “todo” field :-)

@type telescopename: string @param telescopename: The name of the telescope that took this curve @type object: string @param object: The name of the object @type plotcolour: “colour” @param plotcolour: A colour to use when plotting this lightcurve. In principle this is a string, you can use matplotlib names like “red”, “blue” etc, or “#FF0000” …

@note: Usually you will use I{factory-functions} like “importfromvirginie”, and not explicitely call this constructor.

telescopename = None

@type: string @ivar: Choose a name for the telescope

object = None

@type: string @ivar: A name for the object

plotcolour = None

@type: “colour” @ivar: “red”, “blue”, “#00FF00”…

jds = None

@type: 1D float array @ivar: The “dates”, as floats, in HJD, MHJD, … these jds should always be in chronological order please (and of course always correspond to the other arrays).

mags = None

@type: 1D float array @ivar: Magnitudes

magerrs = None

@type: 1D float array @ivar: Errors of magnitudes, as floats

mask = None

@type: 1D boolean array @ivar: A boolean mask, allows to “switch off” some points. True means “included”, False means “skip this one”. This is nice, it directly allows you to use the mask as indices; example (suppose C{mylc} is a lightcurve) :

>>> mylc.mask[34] = False
>>> mylc.mags[mylc.mask].std()

masks the 35th point and gives you the standard deviation of the unmasked mags. Masked points show up on graphs as circled in black, and are naturally skipped by dispersion methods etc.

labels = None

@type: list of strings @ivar: Each point has a string label. Plots can show these labels. Setting these is simple :

>>> mylc.labels[42] = "This point is funny"

We can use them for all kind of things (debugging…).

commentlist = None

@type: list of strings @ivar: This is a list of string comments for each full lightcurve, and NOT one comment for each point. You can “append” comments at any time; they can be shown by various functions, and will trace the history of the lightcurve.

properties = None
Variables:properties – This is a list of dicts (exactly one per data point) in which you can put whatever you want to. You can use them to carry image names, airmasses, sky levels, telescope names etc. Of course this slows down a bit things like merging etc.
ploterrorbars = None

@type: boolean @ivar: A flag to show or hide errorbars in plots of this lightcurve.

showlabels = None

@type: boolean @ivar: A flag to show or hide labels in plots of this lightcurve.

timeshift = None

@type: float @ivar: A float giving the shift in time (same unit as jds) that was applied to the lightcurve. This is updated at every shift, so that a lightcurve is always “aware” of its eventual shift.

magshift = None

@type: float @ivar: A float giving the shift in magnitude (similar to timeshift).

fluxshift = None

@type: float @ivar: A float giving the shift in flux (!). As we work with magnitudes, a shift in flux is a bit special. Note by the way that a flux shift does not commute with a magnitude shift !

ml = None

@type: microlensings object @ivar: The optional microlensings object allows to take microlensing into account during the optimizations and time-delay quests.

__str__()

We explicitly define how str(mylightcurve) will look. “[Euler/A]”, or “[Euler/A](0.00,0.10)” in case the lightcurve is shifted (here by 0.1 mag). For more info, look at longstr(). This will show up in plot legends etc. Plus, at any time, it allows you to :

>>> print mylc

Microlensing “constrains” are shown like |12|, meaning in this case two seasons : first one with 1 parameter, second one with two.

__len__()

We define what len(mylightcurve) should return : the number of points. I do not exclude masked points – this is the full length.

samplingstats(seasongap=60.0)

Calculates some statistics about the temporal sampling of the lightcurve object. I do not take into account the mask !

Parameters:seasongap – Minimal interval in days to consider a gap to be a season gap.

I return a dict with the following keys :

  • nseas : number of seasons
  • meansg : mean length (in days, of course) of the season gaps
  • minsg : min season gap
  • maxsg : max season gap
  • stdsg : standard deviation of the season gaps
  • med : the median sampling interval inside seasons
  • mean : the mean sampling interval inside seasons
  • min : …
  • max
  • std
commonproperties(notonlycommon=False)

Returns a list of strings containing the property field names that are common to all data points. If notonlycommon is True, then I return all fields, not only the ones common to all my data points.

longinfo(title='None')

Returns a multi-line string (small paragraph) with information about the current state of a lightcurve.

printinfo(title='None')

Prints the above longinfo about the lightcurve.

shifttime(days)

Update the value of timeshift, relative to the “current” shift. The jds array is not shifted; we just keep track of a timeshift float. So this is fast.

shiftmag(deltamag)

Shifts the curve in mag, with respect to the “current” shift. Same remark as for shifttime(). @todo: think about changing the errors as well… – no, the error in magnitude is a relative error.

shiftflux(flux)

Idem, but a shift in flux. @todo: Hmm, here we should change the errors as well – but flux shifts are usuall very small, no ?

setfluxshift(flux, consmag=False)

An absolute setter for the fluxshift, that directly tries to correct for the equivalent magnitude shift with respect to the mean magnitude of the curve. So be careful : I change magshift as well !

getjds()

A getter method returning the jds + timeshift. This one does the actual shift. Because of the addition, we return a copy.

getminfluxshift()

Returns the minimal flux shift (i.e. largest negative fluxshift) so that the flux of all points in this curve remain positive. -> returns - the minimum flux of the curve. This does not depend on timeshifts, magshifts, or microlensing. When optimizing fluxshifts, you should check by yourself that you respect this minimum value.

getrawfluxes()

Returns the raw fluxes corresponding to self.mags

calcfluxshiftmags(inverse=False)

Returns an array of mags that you have to add to the lc.mags if you want to take into account for the fluxshift. Be careful with the order of things ! Always take into account this fluxshift before the magshift / microlensing.

Inverse = True is a special feature, it returns the mags to subtract from fluxshiftedmags to get back the original mags. Perfectly clear, no ? Just changes signs in the fomula. You could also use - calcfluxshiftmags(self, inverse=True) and reverse the sign of the fluxshift, but this is more convenient.

addfluxes(fluxes)

Give me fluxes as a numpy array as long as self.mags, and I’ll add those to self.mags. I work directly on self.mags, thus I do not take into account magshifts or ml. As this changes the actual self.mags, it cannot be undone !

getmags(noml=False)
A getter method returning magnitudes. Now this is non trivial, as magnitudes are influenced by :
  • fluxshift
  • magshift
  • microlensing (= like a variable magshift)

We “add” this stuff in this order, i.e. fluxshift always comes first. We always return a copy of the array, to prevent you from changing the actual lightcurve.

lc.mags “+” fluxshift + magshift + microlensings if ml is present, and just lc.mags “+” fluxshift + magshift if not. You can overwrite this : if noml = True, I don’t include the ml even if a microlensing is present !

A bit ugly all these ifs, but I think it makes it faster (the fluxshift calculation is slow).

getmagerrs()

Hmm, how does fluxshift and stuff influence the magerrs ? Copy is important here as we don’t do any calculation.

addml(microlensing, verbose=False)

Adds a microlensing object (by add we mean it replaces one if it was already present) – we are not just setting the parameters here ! Note that we COPY the microlensing object, so that the one added is independent from yours.

rmml()

Simply deletes the microlensing object

resetml()

Resets the microlensing to “zero”. * for polynomial ML the coefficients are put to zero * for spline ML we reset the coeffs to zero

resetshifts(keeptimeshift=False)

Removes all shifts and microlensing, getting back to the observed data. Does keep the mask.

applyfluxshift()

It adds the fluxshift-float to the present flux, then puts this fluxshift-float to 0. So that “nothing” changes as seen from the outside. This is used for instance when you want to merge different lightcurves with different fluxshifts. Needless to say, use this carefully, only if it really makes sense for what you want to do.

Note that we do not touch microlensing here, it remains in place and does not change its meaning in any way.

applymagshift()

It adds the magshift-float to the present mags, then puts this magshift-float to 0. So that “nothing” changes as seen from the outside. This is used for instance when you want to merge different lightcurves with different magshifts. Needless to say, use this carefully, only if it really makes sense for what you want to do.

Note that we do not touch microlensing here, it remains in place and does not change its meaning in any way.

applyml()

We “add” the microlensing to the actual mags, then remove the microlensing object. The advantage is that getmags() gets faster, as it does not have to evaluate the microlensing anymore. It also allows you to do some other tricks like tweaking a curve, shifting seasons etc if you want. So as a conclusion : use this when you know what you are doing.

We do not touch magshift, it remains perfectly valid afterwards. But, if there is a fluxshift, we will have to make sure that it was “applied” before !

setindexlabels()

First point gets “0”, second point gets “1”, etc; handy to identify troublemakers on a plot !

setjdlabels()

Points are labeled by their mjd rounded to .1 days. Use these labels to write a skiplist.

maskskiplist(filepath, searchrange=0.2, accept_multiple_matches=False, verbose=True)

I mask points according to a skiplist. I do not modify the mask for points that are not on the skiplist, i.e. I do not reset the mask in any way. The structure of the skiplist is one “point” per line, where the first element in the line is the MJD (for instance as identified on a plot using setjdlabels() !). For each point from the skiplist, I will search for the corresponding point in the lightcurve within searchrange days. I will warn you in case of anything fishy.

Parameters:
  • filepath (string or path) – file to read
  • searchrange (float) – range in which I’ll look for points (in days).
maskinfo()

Returns a description of the masked points and available properties of them. Note that the output format can be directly used as a skiplist.

clearlabels()

Sets label = “” for all points. You could use this before setting the labels of just some of the points “by hand”, for instance.

clearproperties()

Removes all properties

copy()

Returns a “deep copy” of the lightcurve object. Try to avoid this within loops etc … it is slow ! Typically if you want to optmize time and mag shifts, think about using local backups of lc.getmags() and lc.getjds() etc …

We use the copy module, imported as “pythoncopy” to avoid confusion with this method.

@rtype: lightcurve @return: A copy of the lightcurve.

cutmask()

Erases all masked points from the lightcurve.

This is one way of handling the mask, but perhaps not the best/fastest one, depending on what you want do ! If you write a curve shifting algorithm, probably your code will be more elegant if you teach him to simply skip masked points, instead of making a copy and removing them using this method !

Warning

Cutmask will change the meaning of seasons, that’s why we are forced to delete any present microlensing. If you want to “keep” it, you’ll have to apply it first.

hasmask()

Returns True if some points are masked (i.e., if there is a False), False otherwise.

clearmask()

Sets all the mask to True.

validate(verbose=False)

Checks the “health” and coherency of a lightcurve. Are there as many mags as jds, etc ? No return value, but raises RuntimeError in case of failure.

Note

Here we also check that the curve is “ordered”, i.e. the jds are monotonously increasing. At this stage it is OK to have several points with the same julian date. In fact such points are a problem for splines, but this is addressed once you use splines.

remove_epochs(index)

Delete epochs in your LightCurve

Parameters:index – array or element with the position of the epoch to remove
Returns:
sort()

We sort the lightcurve points according to jds. Of course all the arrays and lists have to be sorted.

The meaning of seasons is changed, that’s why we have to delete any microlensing.

montecarlomags(f=1.0, seed=None)

We add gaussian noise to all mags according to their errors. We do not care about any microlensing of shifts here, but directly tweak the self.mags. The order of the points is not touched.

I modify all points, even masked ones.

Parameters:
  • seed (int or None) – if None, the clock is used, if not, the given seed.
  • f – a multiplier
montecarlojds(amplitude=0.5, seed=None, keepml=True)

We add a UNIFORM random variable from -amplitude to +amplitude to each jd. This is not trivial, as we have to take care about the order.

We need sorting, so normally the microlensing would be removed. But if you know that, using only a small amplitude, the seasons will remain of the same length after the bootstrapping, you can use keepml=True (default), and we will “keep” the microlensing. Note that the microlensing functions are defined with respect to the jds : so as the jds are bootstrapped, the mircolensing is slightly affected, but this is fine as typically you just want to keep a rough estimate of the microlensing, and run an ML optimization anyway ! Hence this obscure option.

I modify all points, even masked ones.

__module__ = 'pycs.gen.lc'
pseudobootstrap()

We randomly mask some points, but without duplicating them. Real bootstap would be to draw N points from a lightcurve of N points with replacements. But we do not replace here, instead we do as if we would replace, but then skip any “duplicates”.

I respect mask : masked points will stay masked, I do as if they were not there.

bootstrap()

Real bootstap = we draw the points with replacements. You will get points with identical jds in the curve !

I respect mask : masked points will stay masked, I do as if they were not there.

Hmm, a bit slow… would need to rearrange all the properties etc

In fact no, should be ok, do similar then above

merge(otherlc)

Merges a lightcurve into the current one. Warning : it takes the other lightcurve as it comes, with all shifts or eventual microlensings – “applied” or just specified ! In fact it behaves as if every shift or microlensing was applied before the merging.

It’s up to you to check that you don’t merge lightcurves with timeshifts, which would probably be nonsense. We delete an eventual current microlensing object -> You will have to define new seasons and parameters if you want.

Parameters:otherlc (lightcurve) – A lightcurve to merge into the current one.

The origin of each point of the otherlc is appended to the labels. Settings and colour etc from the current curve are not changed. After the merging process, the lightcurve is sorted and validated.

rdbexport(filename=None, separator='\t', writeheader=True, rdbunderline=True, properties=None)

Writes the lightcurve into an “rdb” file, that is tab-separeted columns and a header. Note that any shifts/ML are taken into account. So it’s a bit like if you would apply the shifts before writing the file.

Includes mask column only if required (if there is a mask)

Parameters:
  • filename (string or path) – where to write the file
  • separator (string) – how to separate the collumns
  • writeheader (boolean) – include rdb header ?
  • properties (list of strings, e.g. ["fwhm", "ellipticity"]) – properties of the lightcurves to be include in the file.
pycs.gen.lc.factory(jds, mags, magerrs=None, telescopename='Unknown', object='Unknown', verbose=False)

Returns a valid lightcurve object from the provided arrays. The numpy arrays jds and mags are mandatory. If you do not specify a third array containing the magerrs, we will calculate them “automatically” (all the same value), to avoid having 0.0 errors.

@type jds: 1D numpy array @param jds: julian dates @type mags: 1D numpy array @param mags: magnitudes @type magerrs: 1D numpy array @param magerrs: optional magnitude errors

@todo: improve it and use this in file importing functions

pycs.gen.lc.flexibleimport(filepath, jdcol=1, magcol=2, errcol=3, startline=1, flagcol=None, propertycols=None, telescopename='Unknown', object='Unknown', plotcolour='red', verbose=True, absmagerrs=False)

The general form of file reading. We read only one lightcurve object. To comment a line in the input file, start the line with # like in python.

Parameters:
  • jdcol – The column number containing the MHJDs. First column is number 1, not 0 !
  • magcol – The column number containing the magnitudes.
  • errcol – The column number containing the magnitude errorbars.
  • flagcol – (default : None) The colum containing a mask (if available). This one should contain False (= will be masked) or True (not masked).
  • propertycols (dictionary) – (default : None) is a dict : propertycols = {"fwhm":7, "skylevel":8} means that col 7 should be read in as property “fwhm” and 8 as “skylevel”.
pycs.gen.lc.rdbimport(filepath, object='Unknown', magcolname='mag', magerrcolname='magerr', telescopename='Unknown', plotcolour='red', mhjdcolname='mhjd', flagcolname=None, propertycolnames='lcmanip', verbose=True, absmagerrs=False)

The relaxed way to import lightcurves, especially those from cosmouline, provided they come as rdb files. Don’t care about column indices, just give the column headers that you want to read.

Propertycolnames is a list of column names that I will add as properties. Possible settings : “lcmanip” : (default) I will import the standard cols from lcmanip / cosmouline, if those are available. None : I will not import any properties [“property1”, “property2”] : just give a list of properties to import.

The default column names are those used by the method pycs.gen.lc.lightcurve.rdbexport() : “mhjd”, “mag”, “magerr”, “mask”.

We use flexibleimport under the hood.

pycs.gen.lc.multidisplay(setlist=[], title=None, style=None, showlegend=True, showlogo=False, logopos='left', showdates=False, showdelays=False, nicefont=False, text=None, jdrange=None, magrange=None, initfigsize=(12, 8), plotsize=(0.08, 0.96, 0.09, 0.95), showgrid=False, markersize=6, showerrorbars=True, errorbarcolour='#BBBBBB', capsize=3, knotsize=0.015, legendloc='best', showspldp=False, colourprop=None, hidecolourbar=False, transparent=False, collapseref=False, jdmintickstep=100, magmintickstep=0.2, filename='screen', verbose=False)

Function that uses matplotlib to plot a list of lightcurves/splines/GPRs, either on screen or into a file. It uses lightcurve attributes such as lc.plotcolour and lc.showlabels, and displays masked points as black circles. It’s certainly a key function of pycs. You can also put tuples (lightcurve, listofseasons) in the lclist, and the seasons will be drawn. This function is intended both for interactive exploration and for producing publication plots.

Parameters:
  • lclist (list) – A list of lightcurve objects [lc1, lc2, …] you want to plot.
  • splist (list) – A list of spline or rslc (e.g., GPR) objects to display on top of the data points.
  • title (string) – Adds a title to the plot, center top of the axes, usually used for lens names. To nicely write a lens name, remember to use raw strings and LaTeX mathrm, e.g. : title = r"$\mathrm{RX\,J1131-1231}$"
  • style (string) –

    A shortcut to produce specific kinds of stylings for the plots. Available styles:

    • homepagepdf : for cosmograil homepage, ok also with long magnitude labels (like -13.2)
  • showlegend (boolean) – Automatic legend (too technical/ugly for publication plots, uses str(lightcurve))
  • showlogo (boolean) – Adds an EPFL logo + www.cosmograil.org on the plot.
  • logopos (string) – Where to put it, “left” or “right”
  • showdates (boolean) – If True, the upper x axis will show years and 12 month minor ticks.
  • showdelays (boolean) – If True, the relative delays between the curves are written on the plot.
  • nicefont (boolean) – Sets default to serif fonts (terrible implementation, but works)
  • text (list) – Generic text that you want to display on top of the plot, in the form : [line1, line2, line3 …] where line_i is (x, y, text, kwargs) where kwargs is e.g. {“fontsize”:18} and x and y are relative positions (from 0 to 1).
  • jdrange (tuple) – Range of jds to plot, e.g. (53000, 54000).
  • magrange (tuple) – Range of magnitudes to plot, like magrange = [-11, -13]. If you give only a float, like magrange=1.0, I’ll plot +/- this number around the mean curve level Default is None -> automatic mag range.
  • figsize (tuple) – Figure size (width, height) in inches, for interactive display or savefig.
  • plotsize (tuple) – Position of the axes in the figure (left, right, bottom, top), default is (0.065, 0.96, 0.09, 0.95).
  • showgrid (boolean) – Show grid, that is vertical lines only, one for every year.
  • markersize (float) – Size of the data points, default is 6
  • showerrorbars (boolean) – If False, the ploterrorbar settings of the lightcurves are disregared and no error bars are shown.
  • errorbarcolour (string) – Color for the error bars
  • capsize (float) – Size of the error bar “ticks”
  • knotsize (float) – For splines, the length of the knot ticks, in magnitudes.
  • legendloc (string or int) –

    Position of the legend. It is passed to matplotlib legend. It can be useful to specify this if you want to make animations etc.

    • string int
    • upper right 1
    • upper left 2
    • lower left 3
    • lower right 4
    • right 5
    • center left 6
    • center right 7
    • lower center 8
    • upper center 9
    • center 10
  • showspldp (boolean) – Show the acutal data points of spline objects (for debugging etc)
  • colourprop (tuple) – If this is set I will make a scatter plot with points coloured according to the given property, disregarding the lightcurve attribute plotcolour. Format : (property_name, display_name, minval, maxval), where display_name is a “nice” version of property_name, like “FWHM [arcsec]” instead of “seeing”. Note that this property will be used in terms of a float. So you cannot use this for properties that are not floats.
  • hidecolourbar (boolean) – Set to True to hide the colourbar for the colourprop
  • transparent (boolean) – Set transparency for the plot, if saved using filename
  • collapseref (boolean) – Plot one single dashed line as the reference for the microlensing. Use this if you would otherwise get ugly overplotted dashed lines nearly at the same level … This option is a bit ugly, as it does not correct the actual microlensing curves for the collapse.
  • jdmintickstep (float) – Minor tick step for jd axis
  • magmintickstep (float) – Minor tick step for mag axis
  • filename (string) – If this is not “screen”, I will save the plot to a file instead of displaying it. Try e.g. “test.png” or “test.pdf”. Success depends on your matplotlib backend.
  • verbose (boolean) – Set to True if you want me to print some details while I process the curves
pycs.gen.lc.display(lclist=[], splist=[], title=None, titlexpos=None, style=None, showlegend=True, showlogo=False, logopos='left', showdates=False, showdelays=False, nicefont=False, text=None, keeponlygrid=False, jdrange=None, magrange=None, figsize=None, plotsize=(0.08, 0.96, 0.09, 0.95), showgrid=False, markersize=6, showerrorbars=True, showdatapoints=True, errorbarcolour='#BBBBBB', capsize=3, knotsize=0.015, legendloc='best', showspldp=False, colourprop=None, hidecolourbar=False, transparent=False, show_ylabel=True, labelfontsize=14, collapseref=False, hidecollapseref=False, jdmintickstep=100, magmintickstep=0.2, filename='screen', showinsert=None, insertname=None, verbose=False, ax=None)

Function that uses matplotlib to plot a list of lightcurves/splines/GPRs, either on screen or into a file. It uses lightcurve attributes such as lc.plotcolour and lc.showlabels, and displays masked points as black circles. It’s certainly a key function of pycs. You can also put tuples (lightcurve, listofseasons) in the lclist, and the seasons will be drawn. This function is intended both for interactive exploration and for producing publication plots.

Parameters:
  • lclist (list) – A list of lightcurve objects [lc1, lc2, …] you want to plot.
  • splist (list) – A list of spline or rslc (e.g., GPR) objects to display on top of the data points.
  • title (string) – Adds a title to the plot, center top of the axes, usually used for lens names. To nicely write a lens name, remember to use raw strings and LaTeX mathrm, e.g. : title = r"$\mathrm{RX\,J1131-1231}$"
  • titlexpos (float) – Specify where you want your to anchor the center of your title in the x axis. the value is in the x-axis fraction. Default is center. (x = 0.5)
  • style (string) –

    A shortcut to produce specific kinds of stylings for the plots. Available styles:

    • homepagepdf and “homepagepdfnologo” : for cosmograil homepage, ok also with long magnitude labels (like -13.2)
    • 2m2 : for 2m2 data
    • posterpdf : for 2m2 data
    • internal : for 2m2 data
  • showlegend (boolean) – Automatic legend (too technical/ugly for publication plots, uses str(lightcurve))
  • showlogo (boolean) – Adds an EPFL logo + www.cosmograil.org on the plot.
  • logopos (string) – Where to put it, “left” or “right”
  • showdates (boolean) – If True, the upper x axis will show years and 12 month minor ticks.
  • showdelays (boolean) – If True, the relative delays between the curves are written on the plot.
  • nicefont (boolean) – Sets default to serif fonts (terrible implementation, but works)
  • text (list) – Generic text that you want to display on top of the plot, in the form : [line1, line2, line3 …] where line_i is (x, y, text, kwargs) where kwargs is e.g. {“fontsize”:18} and x and y are relative positions (from 0 to 1).
  • jdrange (tuple) – Range of jds to plot, e.g. (53000, 54000).
  • magrange (tuple) – Range of magnitudes to plot, like magrange = [-11, -13]. If you give only a float, like magrange=1.0, I’ll plot +/- this number around the mean curve level Default is None -> automatic mag range.
  • figsize (tuple) – Figure size (width, height) in inches, for interactive display or savefig.
  • plotsize (tuple) – Position of the axes in the figure (left, right, bottom, top), default is (0.065, 0.96, 0.09, 0.95).
  • showgrid (boolean) – Show grid, that is vertical lines only, one for every year.
  • markersize (float) – Size of the data points, default is 6
  • showerrorbars (boolean) – If False, the ploterrorbar settings of the lightcurves are disregared and no error bars are shown.
  • showdatapoints – If False, no data points are shown. Useful if you want e.g. to plot only the microlensing
  • keeponlygrid (boolean) – If True, keeps the yearly grid from showdates but do not display the dates above the plot.
  • showinsert (boolean) – If True, display the insertname image in the top-right corner of the main image
  • insertname – path to the image you want to insert
  • errorbarcolour (string) – Color for the error bars
  • capsize (float) – Size of the error bar “ticks”
  • knotsize (float) – For splines, the length of the knot ticks, in magnitudes.
  • legendloc (string or int) –

    Position of the legend. It is passed to matplotlib legend. It can be useful to specify this if you want to make animations etc.

    • string int
    • upper right 1
    • upper left 2
    • lower left 3
    • lower right 4
    • right 5
    • center left 6
    • center right 7
    • lower center 8
    • upper center 9
    • center 10
  • showspldp (boolean) – Show the acutal data points of spline objects (for debugging etc)
  • colourprop (tuple) – If this is set I will make a scatter plot with points coloured according to the given property, disregarding the lightcurve attribute plotcolour. Format : (property_name, display_name, minval, maxval), where display_name is a “nice” version of property_name, like “FWHM [arcsec]” instead of “seeing”. Note that this property will be used in terms of a float. So you cannot use this for properties that are not floats.
  • hidecolourbar (boolean) – Set to True to hide the colourbar for the colourprop
  • transparent (boolean) – Set transparency for the plot, if saved using filename
  • collapseref (boolean) – Plot one single dashed line as the reference for the microlensing. Use this if you would otherwise get ugly overplotted dashed lines nearly at the same level … This option is a bit ugly, as it does not correct the actual microlensing curves for the collapse.
  • jdmintickstep (float) – Minor tick step for jd axis
  • magmintickstep (float) – Minor tick step for mag axis
  • filename (string) – If this is not “screen”, I will save the plot to a file instead of displaying it. Try e.g. “test.png” or “test.pdf”. Success depends on your matplotlib backend.
  • ax (matplotlib axes) – if not None, I will return what I plot in the given matplotlib axe you provide me with instead of plotting it.
  • verbose (boolean) – Set to True if you want me to print some details while I process the curves
pycs.gen.lc.displayrange(lcs, margin=0.05)

returns a plausible range of mags and hjds to plot, so that you can keep this fixed in your plots

pycs.gen.lc.getnicetimedelays(lcs, separator='\n', sorted=False)

Returns a formatted piece of text of the time delays between a list of lc objects.

Parameters:
  • separator (string) – try " | " to get all the delays in one line.
  • sorted (boolean) – If True, I will sort my output according to l.object. But of course I will not modify the order of your lcs ! By default this is False : if the curves are B, A, C and D in this order, I would return BA, BC, BD, AC, AD, CD

Warning

This is the function that defines the concept of “delay” (and its sign) used by pycs. The delay AB corresponds to timeshift(B) - timeshift(A) Other places where this is hard-coded : pycs.sim.plot.hists()

pycs.gen.lc.getnicetimeshifts(lcs, separator='\n')

I return the timeshifts as a text string, use mainly for checks and debugging.

pycs.gen.lc.gettimeshifts(lcs, includefirst=True)

I simply return the absolute timeshifts of the input curves as a numpy array. This is used by time shift optimizers, and usually not by the user.

Parameters:includefirst (boolean) – If False, I skip the first of the shifts. For timeshift optimizers, it’s indeed often ok not to move the first curve.
pycs.gen.lc.settimeshifts(lcs, shifts, includefirst=False)

I set the timeshifts like those returned from gettimeshifts(). An do a little check.

pycs.gen.lc.shuffle(lcs)

I scramble the order of the objects in your list, in place.

pycs.gen.lc.objsort(lcs, ret=False, verbose=True)

I sort the lightcurve objects in your list (in place) according to their object names.

Parameters:ret (boolean) – If True, I do not sort them in place, but return a new sorted list (containing the same objects).

pycs.gen.mrg module

Functions to match and merge lightcurves usually of different telescopes, prior to exploring time delays. We typically adjust the magshifts and fluxshifts of curves from telescope X so that they match the same curves from telescope Y, without touching the time axis. Of course we need overlapping curves to do this.

pycs.gen.mrg.matchtels(lcsref, lcsmatch, dispersionmethod, fluxshifts=True)

Shifts the curves of lcsmatch so that they fit to lcsref. lcsmatch and lcsref are lists of corresponding lightcurves. We do not modify lscref.

To help me a bit, set fluxshift and magshifts of the lcsmatch to some good starting values before calling me !

For now, I do :
  • one single magshift for all curves
  • one fluxshift for each curve (if fluxshifts = True)
pycs.gen.mrg.merge(lcss)

Function to merge (across the first index) a list of lists of lightcurves. example :

lcss = [[EulerC2_A, EulerC2_B], [Mercator_A, Mercator_B]]

# will *return* [merged_A, merged_B]
# where merged_i is EulerC2_i.merge(Mercator_i)

Note

I do not modify the input lightcurves !

pycs.gen.mrg.colourise(lcs)

Attributes different colours to the lightcurves of the list lcs.

http://www.w3schools.com/html/html_colornames.asp

pycs.gen.mrg.colorize(lcs)

pycs.gen.polyml module

One season – one microlensing This potentially allows to optimize one season after the other, thus reducing the number of simultaneous parameters, but of course you can still choose to optimize evrything at once (when using dispersion). When using a spline fit, optimizing a polynom corresponds to weighted to linear least squares – fast !

Idea for the future : “mutable” objects should be used here, with pointers in tables to pass the positions of each ml parameter very quickly. The present implementation is not that elegant… we make to many copies.

We still need a more sophisticated factory function, so that you can set masks on individual params. This is not possible for now, only because of the lack of such a factory. (edit : not sure if really needed, I never use these…) Let’s forget about these “free params” and masks.

pycs.gen.polyml.polyfit(jds, mags, magerrs, nparams)

I define my own polyfit, as numpy polyfit does not accept weights until numpy version 1.5 This is linear least squares – can only be used when optimizing polyml so to fit a “model” like a spline.

What we want to do : fitp = np.polyfit(jds, mags, deg = nparams-1)

Parameters:
  • jds – x values
  • mags – y values
  • magerrs – y errors
  • nparams – number of parameters of polynom to fit.

Code to test that the “heteroskedasticity correction” works :

import matplotlib.pyplot as plt

x = np.linspace(0, 10, 20)
ystd = 1.0 * np.ones(len(x))
yerrs = ystd * np.random.randn(len(x))
yerrs[-3] = 40.0
ystd[-3] = 40.0 
print yerrs

y = 5.0 + 2.0*x - 1.0*x*x + yerrs

finex = np.linspace(0, 10, 100)

plt.errorbar(x, y, ystd)

polyp = polyfit(x, y, ystd, 3)
#polyp = np.polyfit(x, y, 3)

finey = np.polyval(polyp, finex)

plt.plot(finex, finey)
plt.show()
class pycs.gen.polyml.seasonfct(season, mltype='poly', params=array([0.]), mask=array([ True]))

A seasonfct object is one season + a certain polynom (simple poly or legendre for now – from outside you do not see what it is) + some stuff to play with that. By itself, it makes no sense. It needs to be linked to a lightcurve : this is done by adding it to lc.ml, wich is a microlensing object, i.e. essentially a list of seasonfct objects.

The idea is that the number of params and the mask remain fixed for one given seasonfct object. This makes the code simpler and faster. So once the instance is made, only the values of the parameters can be updated. If you want to change something else, simply make a new object.

One idea of the relatively simple get and set methods is that one day we might want to use something else then polynoms, and at that time all this might be handy.

for the mltype argument, see discussion of ml.factory function.

__init__(season, mltype='poly', params=array([0.]), mask=array([ True]))

Also see the factory functions below … they will be much easier to use for your task. Default values shown above : 1 parameter, i.e. a simple additive constant for the given season, which is free to be optimized.

season = None

@type: season @ivar: The season for which this microlensing if for.

mltype = None

@type: string @ivar: The type of microlensing

params = None

@type: ndarray @ivar: A 1D numpy array of params (floats)

mask = None

@type: ndarray @ivar: boolean array : True if the parameter is “free”.

nfree = None

@type: int @ivar: How many free parameters

copy()
__str__()

This is just the number of params of the polynom (i.e. degree + 1) (Not the number of free params !)

longinfo()

Returns a longer description of the object.

printinfo()

Prints the longer description of the object.

setparams(p)

Remember that you are not allowed to change the number of parameters !

getparams()
setfreeparams(p)

Here we should really do some testing … be careful when using this !!!

getfreeparams()
validate()

@todo: code this in case of ideas…

checkcompatibility(lightcurve)

Not sure if needed.

calcmlmags(lightcurve)

Returns a “lc.mags”-like array made using the ml-parameters. It has the same size as lc.mags, and contains the microlensing to be added to them. The lightcurve object is not changed !

For normal use, call getmags() from the lightcurve.

Idea : think about only returning the seasons mags to speed it up ? Not sure if reasonable, as no seasons defined outside ?

smooth(lightcurve)

Only for plotting purposes : returns jds, mlmagshifts, and refmags with a tight and regular sampling, over the range given by the season. Note that this time we are interested in the acutal shifted jds, so that it looks right when plotted ! We return arrays that can directly be plotted to illustrate the microlensing.

TODO : return refmag, not refmags !

__module__ = 'pycs.gen.polyml'
class pycs.gen.polyml.microlensing(mllist)

Contains a list of seasonfct objects OF ONE SAME LIGHTCURVE OBJECT, and some methods to treat them.

You probably do not want your seasonfct seasons to overlap. But no problem if they do.

Again : do not change the contents of such a microlensing object otherwise then retrieving and updating the params through the provided methods ! Do not change the list of seasonfct objects in any other way ! 1) build the seasonfct objects 2) put them into a microlensing object 3) do not touch the seasonfct objects anymore if you do not know what you are doing.

Instead, make a brand new object if you need, using factory functions provided below or to be written.

__init__(mllist)
mllist = None

@type: list @ivar: The list of seasonfct objects to be applied.

nfree = None

@type: int @ivar: The number of free parameters.

copy()
__str__()
longinfo()
printinfo()
getfreeparams()

Straightforward.

setfreeparams(p)

This method distributes the params on the microlensing objects as fast as I could … it’s a bit delicate and unelegant. As we want to be fast – do not mess with this and give a p with exactly the right size.

reset()

Puts all coefs back to 0.0 Allows to start from scratch without having to rebuild a new ML.

__module__ = 'pycs.gen.polyml'
checkcompatibility(lightcurve)
calcmlmags(lightcurve)

Returns one a “lc.mags”-like array made using the parameters of all seasonfct objects. This array has the same size as lc.mags, and contains the microlensing to be added to lc.mags.

Idea : perhaps we can make this faster by not calling calcmags of the seasonfct objects ?

stats(lightcurve)

Calculates some statistics on the microlensing deformation evaluated at the same sampling than the lightcurve. The idea is to get the flux ratios, in case of calm microlensing.

pycs.gen.polyml.multigetfreeparams(lclist)

Give me a list of lightcurves (with or without ml !), I give you a single flat array of parameters of all MLs concatenated. Note that the order of lightcurves in lclist is important ! You will have to give the same order for multisetfreeparams() !

For now a bit simple … but it’s not that slow.

pycs.gen.polyml.multisetfreeparams(lclist, params)

Be careful to respect the order of lcs in lclist … otherwise everything gets messed up. Again this seems a bit slower then it really is – prefer simplicity.

pycs.gen.polyml.factory(seasons, nparams, mltype='poly')

A factory function to create a microlensings object filled by seasonfct objects. seasons is a list of season objects nparams is an array or list of ints. “default” = one constant per season. All parameters will be set to 0.0, and free to be optimized.

mltype = “poly” : simple polynomial microlensing, very stupid but fast, ok for degree <= 3 default type.

mltype = “leg” : legendre polynomials, very clever but slow :-) These are fine for degree <= 25 Above deg 25, some numerical errors set in. Could perhaps be rewritten to make this faster. Or implement in C !!!!

pycs.gen.polyml.addtolc(l, seasons=None, nparams=1, autoseasonsgap=60.0)

Adds polynomial ML to the lightcurve l. Top level function, to make it really easy ! We just wrap the factory above.

If seasons = None, I will make some autoseasons. You can specify autoseasonsgap. Else, seasons should be a list of seasons (see factory function above)

If nparams is an int, each season will get a polynom of nparams parameters. 1 = contant, 2 = slope, … Else nparams should be a list of ints (see factory function above)

If you want one single polynom for the full curve, just set autoseasonsgap to 10000.0 …

pycs.gen.polyml.addtolc(l, nparams=1) # One constant on each season.

pycs.gen.sea module

What we call a season is in fact a way of specifying any “connected” range of points in a lightcurve. Seasons are needed for microlensing, polynomial methods, etc. Functions that want or return “seasons” usually handle lists of these season objects. And as such lists of seasons are so common, there are a few functions below that directly handle lists of seasons.

class pycs.gen.sea.season(indices, name='')

A season is essentially a numpy integer array containing the indexes of the lightcurve to consider. I made this class to allow for a bit more explicit error checking, and display capabilities etc.

__init__(indices, name='')

Prefer to use the factory functions below.

@type indices: int array @param indices: The included indices of the season

indices = None

@type: array of ints @ivar: An array of indices, that codes a season. You can use it like lc.jds[season.indices] to get the jds of that season.

__str__()
getjdlims(lc)

Returns a tuple of two jds giving the extremities of the season.

validate()

This checks some inherent properties of a season, that must be true independent of any lightcurve.

checkcompatibility(lc)

This checks if the season is compatible with the given lightcurve.

__module__ = 'pycs.gen.sea'
pycs.gen.sea.autofactory(l, seasongap=60, tpe='seas')

A factory function that automatically create season objects. Formerly known as the lightcurve method “autoseasons”. Returns a LIST of season objects, describing the “shape” of the seasons of lc. In this case, these seasons are determined simply by looking for gaps larger then seasonsgap in the lightcurve.

mode option added 11/2013:

mode=’seasons’ (default option) , do create the season objects descibed above mode=’interseasons’ , return the extremity days between seasons, useful for cutting seasons in a synthetic lightcurve.

Small example : suppose you have a lightcurve called C{mylc}, then :

>>> myseasons = sea.autofactory(mylc)
>>> mylc.mags[myseasons[0].indices] += 0.2

… increases the mags of the first season by 0.2 ! Or perhaps you want to mask the second season :

>>> mylc.mask[myseasons[1].indices] = False

@rtype: list @return: list of season objects.

@todo: rewrite the example above.

pycs.gen.sea.manfactory(lc, jdranges)

A factory function that returns a LIST of season objects. As input, you give a lightcurve object as well as a list of JD ranges. Note that if the lightcurve is shifted, shifted jds are used !

Example:
>>> seas1 = sea.manfactory(lc1, [[2800, 3500], [3500, 4260]])

@rtype: list @return: list of season objects.

pycs.gen.sea.validateseasons(seasons)

Checks if your seasons (i.e. a list of season objects) can really be considered as seasons … This is typically called by the microlensing constructor. Seasons can well be “incomplete”, i.e. you can give only one season for your 4 year lightcurve, but we do not allow for overlapping seasons, as this would -up microlensing.

pycs.gen.sea.printinfo(seasons)
pycs.gen.sea.easycut(lclist, keep=1, seasongap=60, mask=False, verbose=False)

Wrapper function to easily cut out some seasons from a list of lightcurves. Seasons are determined for each lightcurve individually

keep : a tuple of “indices” (human convention, starting at 1) of seasons to keep.

mask : if True, I will not cut the seasons, but only mask them out.

pycs.gen.spl module

Module defining the Spline class, something easy to wrap around SciPy splines. Includes BOK algorithms (Mollinari et al)

Some rules of splrep (k = 3)
  • do not put more then 2 knots between data points.
  • splrep wants inner knots only, do not give extremal knots, even only “once”.
class pycs.gen.spl.DataPoints(jds, mags, magerrs, splitup=True, deltat=1e-06, sort=True, stab=False, stabext=300.0, stabgap=30.0, stabstep=5.0, stabmagerr=-2.0, stabrampsize=0, stabrampfact=1.0)

An ultralight version of a lightcurve, made for fast computations. Can be “merged” from a list of lightcurves, see factory function below.

A Spline object has such a DataPoints object as attribute.

ATTENTION Datapoints are expected to be ALWAYS SORTED BY JDS, and no two datapoints have the same jd ! See the splitup option of the constructor. Note that this is not the case for lightcurves ! Hence the existence of datapoints. Should be enforced in every function that builds datapoints.

ABOUT STAB POINTS With scipy splines, we always get the last knots at the extrema of data points. So to get knots “outside” of the real datapoints, we have to insert fake points. And while we are at it, these fake points can also be used to stabilize the spline in gaps. The mask is used to differentiate between actual data points and “stabilization points” that are inserted to make the spline behave well at the extrema and in season gaps. It is modified by the two addgappts and addextpts.

The info about stabpoints is written into the object, so that they can be reconstrucuted from any new jds and mags.

__init__(jds, mags, magerrs, splitup=True, deltat=1e-06, sort=True, stab=False, stabext=300.0, stabgap=30.0, stabstep=5.0, stabmagerr=-2.0, stabrampsize=0, stabrampfact=1.0)

Constructor Always leave splitup and sort on True ! Only if you know that you are already sorted you can skip them. You cannot specify a mask, I do this myself. (could be done in principle).

stab : do you want stabilization points ? Don’t forget to run splitup, sort, and addstab again if you change the data !

splitup()

TO WRITE !!! We avoid that two points get the same jds… Note that this might change the order of the jds, but only of very close ones, so one day it would be ok to leave the mags as they are.

sort()

Absolutely mandatory, called in the constructor.

validate()

We check that the datapoint jds are increasing strictly :

rmstab()

Deletes all stabilization points

putstab()

Runs only if stab is True. I will : add datapoints (new jds, new mags, new magerrs) modify the mask = False for all those new datapoints.

calcstabmagerr()

Computes the mag err of the stabilisation points.

addgappts()

We add stabilization points with low weights into the season gaps to avoid those big excursions of the splines. This is done by a linear interpolation across the gaps.

addextpts()

We add stabilization points at both extrema of the lightcurves This is done by “repeating” the extremal points, and a ramp in the magerrs

getmaskbounds()

Returns the upper and lower bounds of the regions containing stabilization points. This is used when placing knots, so to put fewer knots in these regions. Crazy stuff…

ntrue()

Returns the number of real datapoints (skipping stabilization points)

__module__ = 'pycs.gen.spl'
pycs.gen.spl.merge(lcs, olddp=None, splitup=True, deltat=1e-06, sort=True, stab=False, stabext=300.0, stabgap=30.0, stabstep=5.0, stabmagerr=2.0, stabrampsize=0, stabrampfact=1.0)

Factory function for DataPoints objects, starting from lightcurves. Takes a list of lightcurves and quickly concatenate the jds, mags, and magerrs.

Instead of specifying all the stab point parameters, you can give me an old datapoints object, and I will reuse its settings… This is useful if you want to “update” the data points.

If overlap is True, I will keep only points that are “covered” by all four lightcurves ! This is useful when you want to build a first source spline, and your microlensing is messy at the borders. NOT YET IMPLEMENTED …

class pycs.gen.spl.Spline(datapoints, t=None, c=None, k=3, bokeps=2.0, boktests=5, bokwindow=None, plotcolour='black')

A class to represent a spline, that is essentially a set of knots and coefficients. As finding knots and coefficients requires access to some data points, these are included in the form of a DataPoints object.

Abount knots : Spline.t are all the knots, including extremas with multiplicity. But splrep wants internal knots only ! By internal we mean : not even the data extremas ! Spline.getintt() returns only these internal knots.

__init__(datapoints, t=None, c=None, k=3, bokeps=2.0, boktests=5, bokwindow=None, plotcolour='black')

t : all the knots (not only internal ones !) c : corresponding coeffs k : degree : default = cubic splines k=3 -> “order = 4” ??? whatever … 3 means that you can differentiate twice at the knots.

__str__()

Returns a string with: * degree * knot placement * number of intervals

copy()

Returns a “deep copy” of the spline.

shifttime(timeshift)

Hard-shifts your spline along the time axis. By “hard-shift”, I mean that unlike for a lightcurve, the spline will not know that it was shifted ! It’s up to you to be sure that you want to move it. We shift both the datapoints and the knots.

shiftmag(magshift)

Hard-shifts your spline along the mag axis. By “hard-shift”, I mean that unlike for a lightcurve, the spline will not know that it was shifted ! It’s up to you to be sure that you want to move it. We shift both the datapoints and the knots.

updatedp(newdatapoints, dpmethod='stretch')

Replaces the datapoints of the spline, and makes sure that the knots stay compatible.

If you tweaked your datapoints, I will have to tweak my knots to make sure that my external knots fit. Hence this method ! Due to the splitup, this is needed even if you just tweaked the mags ! And anyway in this case I have to rebuild the stab points.

Warning

IT’S UP TO YOU TO CHECK THAT YOU DON’T REPLACE DATATOINTS WITH DIFFERENT STAB SETTINGS Anyway it would work, just look ugly !

Replaces the datapoints (jds, mags, and magerrs) touching the knots and coeffs as less as possible. Note that we also have to deal with stab points here !

This is made for instance for time shifts that only very slightly change the datapoints, and you don’t want to optimize the knots all the time from scratch again. The current knots are “streched” (keeping their relative spacings) accross the new datapoints.

Options for “dpmethod” : - “stretch” : changes all the knots - “extadj” : does not touch the internal knots, but adjusts the external ones only, to fit the new datapoints. Probably the method to use when optimizing time shifts. - “leave” : does not touch the knots -> ok to evaluate the spline, but you will not be able to fit it anymore, as the external knots don’t correspond to datapoints.

Todo

In principle, why don’t we just update the real datapoints here, and leave the stab as they are ?

uniknots(nint, n=True)

Uniform distribution of internal knots across the datapoints (including any stab points). We don’t make a difference between stab and real points.

Parameters:
  • nint – The number of intervals, or the step
  • n (boolean) – If True, nint is the number of intervals (== piecewise polynoms) you want. If False : nint is a step in days you want between the knots (approximately).

Note

I also put all coeffs back to 0.0 !

resetc()

Sets all coeffs to 0.0 – if you want to start again your fit, keeping the knot positions.

reset()

Calls uniknots, i.e. resets both coeffs and knot positions, keeping the same number of knots.

buildbounds(verbose=True)

Build bounds for bok. By default I will make those bounds as wide as possible, still respecting epsilon. The parameter epsilon is the minimum distance two knots can have. If you give me a window size, I will not make the bounds as wide as possible, but only put them 0.5*window days around the current knots (still respecting all this epsilon stuff of course).

I look where your current knots are, and for each knots I build the bounds so that epsilon distance is respected between adjacent upper and lower bounds. But, there might already be knots only epsilon apart. So I’m a bit tricky, not so straightforward as my predecessors.

Knots at the extrema are not allowed to move.

Requires existing knots, puts lims in between them, and builds the bounds. @todo: Optimize me using numpy ! This is experiemental code for now.

bok(bokmethod='BF', verbose=True, trace=False)

We optimize the positions of knots by some various techniques. We use fixed bounds for the exploration, run buildbounds (with low epsilon) first. This means that I will not move my bounds.

For each knot, i will try ntestpos linearly spaced positions within its bounds. In this version, the bounds are included : I might put a knot on a bound ! The way the bounds are placed by buildbounds ensures that in any case the minimal distance of epsilon is respected.

Using this sheme, it is now possible to iteratively call mybok and buildbounds in a loop and still respect epsilon at any time.

bokmethods :
  • MCBF : Monte Carlo brute force with ntestpos trial positions for each knot
  • BF : brute force, deterministic. Call me twice
  • fminind : fminbound on one knot after the other.
  • fmin :global fminbound

Exit is automatic, if result does not improve anymore…

getintt()

Returns the internal knots (i.e., not even the datapoints extrema) This is what you need to feed into splrep ! There are nint - 1 such knots

getinttex()

Same as above, but we include the extremal points “once”.

knotstats()

Returns a string describing the knot spacing

setintt(intt)

Give me some internal knots (not even containing the datapoints extrema), and I build the correct total knot vector t for you. I add the extremas, with appropriate multiplicity.

@TODO: check consistency of intt with datapoints !

setinttex(inttex)

Including extremal knots

getnint()

Returns the number of intervals

getc(m=0)

Returns all active coefficients of the spline, the ones it makes sense to play with. The length of this guy is number of intervals - 2 !

setc(c, m=0)

Puts the coeffs from getc back into place.

getco(m=0)

Same as getc, but reorders the coeffs in a way more suited for nonlinear optimization

setco(c, m=0)

The inverse of getco.

setcflat(c)

Give me coeffs like those from getc(m=1), I will set the coeffs so that the spline extremas are flat (i.e. slope = 0).

setcoflat(c)

idem, but for reordered coeffs.

r2(nostab=True, nosquare=False)

Evaluates the spline, compares it with the data points and returns a weighted sum of residuals r2.

If nostab = False, stab points are included This is precisely the same r2 as is used by splrep for the fit, and thus the same value as returned by optc !

This method can set lastr2nostab, so be sure to end any optimization with it.

If nostab = True, we don’t count the stab points

tv()

Returns the total variation of the spline. Simple ! http://en.wikipedia.org/wiki/Total_variation

optc()

Optimize the coeffs, don’t touch the knots This is the fast guy, one reason to use splines :-) Returns the chi2 in case you want it (including stabilization points) !

Sets lastr2stab, but not lastr2nostab !

optcflat(verbose=False)

Optimizes only the “border coeffs” so to get zero slope at the extrema Run optc() first … This has to be done with an iterative optimizer

eval(jds=None, nostab=True)

Evaluates the spline at jds, and returns the corresponding mags-like vector. By default, we exclude the stabilization points ! If jds is not None, we use them instead of our own jds (in this case excludestab makes no sense)

display(showbounds=True, showdatapoints=True, showerrorbars=True, figsize=(16, 8))

A display of the spline object, with knots, jds, stab points, etc. For debugging and checks.

__module__ = 'pycs.gen.spl'
pycs.gen.spl.fit(lcs, knotstep=20.0, n=None, knots=None, stab=True, stabext=300.0, stabgap=20.0, stabstep=5.0, stabmagerr=-2.0, stabrampsize=0, stabrampfact=1.0, bokit=1, bokeps=2.0, boktests=5, bokwindow=None, k=3, verbose=True)

The highlevel function to make a spline fit.

lcs : a list of lightcurves (I will fit the spline through the merged curves)

Specify either knotstep : spacing of knots or n : how many knots to place or knots : give me actual initial knot locations, for instance prepared by seasonknots.

stab : do you want to insert stabilization points ? stabext : number of days to the left and right to fill with stabilization points stabgap : interval of days considered as a gap to fill with stab points. stabstep : step of stab points stabmagerr : if negative, absolte mag err of stab points. If positive, the error bar will be stabmagerr times the median error bar of the data points.

bokit : number of BOK iterations (put to 0 to not move knots) bokeps : epsilon of BOK boktests : number of test positions for each knot

pycs.gen.spl.seasonknots(lcs, knotstep, ingap, seasongap=60.0)

A little helper to get some knot locations inside of seasons only

knotstep is for inside seasons ingap is the number of knots inside gaps.

pycs.gen.spl.r2(lcs, spline, nosquare=False)

I do not modify the spline (not even its datapoints) ! Just evaluate the quality of the match, returning an r2 (without any stab points, of course).

This is used if you want to optimize something on the lightcurves without touching the spline.

Of course, I do not touch lastr2nostab or lastr2stab of the spline ! So this has really nothing to do with source spline optimization !

pycs.gen.spl.mltv(lcs, spline, weight=True)

Calculates the TV norm of the difference between a lightcurve (disregarding any microlensing !) and the spline. I return the sum over the curves in lcs.

Also returns a abs(chi) like distance between the lcs without ML and the spline

If weight is True, we weight the terms in sums according to their error bars.

Idea : weight the total variation somehow by the error bars ! Not sure if needed, the spline is already weighted.

pycs.gen.spl.optcmltv(lcs, spline, verbose=True)

I will optimize the coefficients of the spline so to minimize the mltv. I do not use the microlensing of the lcs at all !

Simple powell optimization, slow. A pity.

Add BOK and time shifts in there and it might be bingo !

Would be more efficient if we add knots on the fly

pycs.gen.splml module

Microlensing represented by splines. Rewamped in July 2011

class pycs.gen.splml.SplineML(spline)

A lightcurve can have such a microlensing object as attribute “ml”. For now we have only one spline per lightcurve (i.e., no seasons).

The splml.spline.datapoints are relative, first point is always at 0.0 That’s why we store a jdoffset EDIT : No, we don’t store it anymore. The spline just starts at 0.0, nobody here cares about this offset… use the timeshift of my lightcurve, if you want to know where I am located on an absolute time axis.

We dont use stabililzation points for microlensing splines (not needed I hope) Hence, the datapoints.jds are directly the jds - jdoffset from the associated lightcurve (They will not change) EDIT : we now use stab points, but still, the jds remain fixed.

The datapoint mags = what has to be added to the lightcurve so that it matches the sourcespline (or something else). This is the same convention set as for ML by polynoms.

__init__(spline)

EDIT : I do not store jdoffset anymore, makes no sense. jdoffset = Offset to add to spline.datapoints so to get back jds as they were at the moment of the factory call. Hmm, this is a bit stange, I think nobody needs to know this offset. If you want to plot the ML for instance, you should not use it, but take the current lc.timeshift, as the curve might have been shifted in time, and this jdoffset does not get updated !

__str__()
copy()
checkcompatibility(lightcurve)

It should be checked here that the self.datapoints are comptabile with the lightcurve (i.e. same “relative” jds)

settargetmags(lightcurve, sourcespline)

We update the self.spline.datapoints.mags so that a fit will results in a microlensing approximation.

Sourcespline is a spline object that is the reference where you want your lightcurve to go.

Only the mags are updated ! Indeed for this microlensing spline I don’t care if your lightcurve was shifted in time. But of course a shift in time with restpect to the sourcespline will give you differnent mags !

I hope you did not remove points or other tweaks.

@todo: treat magerrs with a getter method !

EDIT : Let’s try to add some stabilization points in the datapoints… The constructor of the SplineML now defines how to calculate these stabilisation points. These settings are permanent for one SplineML, hence like without any stabpoints, all we need to do is to update the mags. No need to update the stab point jds, or the datapoints.mask etc. But, to update the mags, we need to calculate the mags for these stab points.

replacespline(newspline)

If you want to change the spline of this SplineML object I should maybe check that you didn’t change stab attribute stuff here ?

reset()

Puts all coeffs back to zero, and redistributes the knots uniformly. The number of knots does not change of course !

calcmlmags(lightcurve=None)

Required by lc (for getmags, applyml, etc…) Returns a mags-like vector containing the mags to be added to the lightcurve.

We don’t need the lightcurve object ! Yes, as we do not care about time shifts, returning only mags. I leave it in the arguments as it is present in the calcmlmags of the polynomial splines.

smooth(lightcurve, n=1000)

Returns a dict of points to plot when displaying this microlensing. Just for plots.

Warning : we suppose here that lightcurve is the same as was used to build the ml !

Here we directly return stuff that can be plotted, so we do care about timeshifts.

n is the number of points you want. the more the smoother.

__module__ = 'pycs.gen.splml'
pycs.gen.splml.addtolc(lc, targetlc=None, n=5, knotstep=None, stab=True, stabgap=30.0, stabstep=3.0, stabmagerr=1.0, bokeps=10.0, boktests=10, bokwindow=None)

Adds a SplineML to the lightcurve. SplineML splines have NO external stabilization points (stabext = 0.0) ! We ignore any mask of the lc, cut it before putting this ML. This is just about putting things in place. Of course we cannot optimize anything here !

If targetlc is not None, then pass me another light curve that already have a microlensing spline. I will “copy” that microlensing to your lc and adjust the knots coefficients accordingly.

The stab stuff inserts stab points into gaps only.

We use uniknots, n is the number of uniform intervals you want. Or specify knotstep : if this is specified, I don’t use n.

pycs.gen.splml.addtolc(l, knotstep=200)

pycs.gen.stat module

Statistics related stuff.

pycs.gen.stat.normal(x, mu, sigma)
pycs.gen.stat.sf(l, binsize=200, ssf=False)

Structure function of a lightcurve

ssf gives a 2D density plot, otherwise binned.

For definition see for instance : De Vries, W.H. de, Becker, R., White, R., Loomis, C., 2005. Structure Function Analysis of Long-Term Quasar Variability. The Astronomical Journal 129, 615-615-629-629.

pycs.gen.stat.mad(data, axis=None)

Median absolute deviation :param data: array from which to compute the MAD :param axis: axis along to compute the MAD :return: float, MAD of the array

pycs.gen.stat.erf(x)

Error function. There is one in scipy, but this way we do it without scipy…

scipy.special.erf(z)

pycs.gen.stat.runstest(residuals, autolevel=False, verbose=True)

One-sample runs test of randomness as presented in Practical Statistics for Astronomers by J. V. Wall and C. R. Jenkins, paragraph 5.3.3 WARNING : ERROR IN THE BOOKS EXPECTATION FORMULA ! confirmed by author

residuals is a numpy array of floats

returns z (in units of sigmas) and p = assuming that the data are independent, the probability to get a result worse then what you got.

pycs.gen.stat.subtract(lcs, spline)

I return a list of residual light curves (“lcs - spline”). Technically, “residual” light curves are nothing but normal lightcurves objects. Of course, I take into account any shifts or microlensing of your lcs. I do not modify my input arguments.

pycs.gen.stat.resistats(rl)

Give me a residual lightcurve, I return a dict with some descriptive stats about its magnitudes.

pycs.gen.stat.mapresistats(rls)
pycs.gen.stat.anaoptdrawn(optoriglcs, optorigspline, simset='simset', optset='optset', npkl=1000, plots=True, nplots=3, r=0.11, plotjdrange=None, plotcurveindexes=None, showplot=False, directory='./', resihist_figsize=None)

Not flexible but very high level function to analyse the spline-fit-residuals of drawn curves and comparing them to the real observations. This can be used to tune the parameters of the “drawing”. .. warning:: The simset must have been optimized using spline fits, with option keepopt=True

Parameters:
  • optoriglcs – optimized original curves
  • optorigspline – spline that matches to these curves
  • simset – Put this e.g. to plotcurveindex=(0,2) if you want to plot only the first 2 curves …
  • plotcurveindexes – allows you to plot only a subset of lcs (smaller plots). Give a tuple like eg (0, 2, 3)
  • npkl – I read only the first npkl picke files.
pycs.gen.stat.plotresiduals(rlslist, jdrange=None, magrad=0.1, errorbarcolour='#BBBBBB', showerrorbars=True, showlegend=True, nicelabel=True, showsigmalines=True, filename=None, ax=None)

We plot the residual lightcurves in separate frames.

The arguement rlslist is a list of lists of lightcurve objects. Ths sublists should have the same length, I’ll choose my number of panels accordingly. The structure is : [[lca, lcb], [lca_sim1, lcb_sim1], …] If you have only one lightcurve object, you can of course pass [[l]] …

Parameters:rlslist –

I disregard the timeshift of the curves !

pycs.gen.util module

Various useful stuff. For now there are some wrappers to pickle objects.

pycs.gen.util.writepickle(obj, filepath, verbose=True, protocol=-1)

I write your python object obj into a pickle file at filepath. If filepath ends with .gz, I’ll use gzip to compress the pickle. Leave protocol = -1 : I’ll use the latest binary protocol of pickle.

pycs.gen.util.readpickle(filepath, verbose=True)

I read a pickle file and return whatever object it contains. If the filepath ends with .gz, I’ll unzip the pickle file.

pycs.gen.util.oldwritepickle(obj, filepath)

DO NOT USE ME ANYMORE Simplistic wrapper around pickle, to writes an object into a file, using cpickle.

@type obj: object @param obj: e.g. a dict of lightcurves @type filepath: string @param filepath: filename or path to write

pycs.gen.util.oldreadpickle(filepath)

DO NOT USE ME ANYMORE Reads a pickle and returns it.

@type filepath: string @param filepath: filename or path to read

@rtype: object @return: whatever was in that pickle

pycs.gen.util.readidlist(filepath, verbose=True)

Reads a textfile with “one point per line”, probably a skiplists. Accepts blank lines, and lines starting with # will not be read. Format of the lines is : id [comment] If this is a skiplist, id is a MJD.

pycs.gen.util.trace(lclist=[], splist=[], tracedir='trace', level='Full')

Function to save a “trace” of processes modifying lightcurves and splines, like optimizers do. Just call this from inside your loop etc, I will save your current lightcurves and splines into a pickle inside the tracedir. I increment the filenames.

The argument “level” is about what should be saved. level = “Full” : Everything you give me is saved in the pickle. Now this is large … level = “Light” : I try to reduce filesize of the pickle, by removing the splines datapoints etc. You can still plot these objects.

pycs.gen.util.plottrace(tracedir='trace', reset=False, showspl=True, **kwargs)

Turns a trace into plots … reset = True : I will remove all shifts/ML etc, just show the real “observations”. kwargs are passed to the display function.

pycs.gen.util.multilcsexport(lclist, filepath, separator='\t', rdbunderline=True, verbose=True, properties=None)

Writes the lightcurves as flat acscii files into one single file. Normally you should prefer writing each lightcurve into a single file, using pycs.gen.lc.lightcurve.rdbexport().

Note that only lightcurves of same length and sampling can be written with this function !

Parameters:
  • lclist (list) – A list of lightcurve objects to write
  • filepath (string) – where to write
  • separator (string) – how to separate the collumns
  • rdbunderline (boolean) – do you want the “=====” underlining ?
  • properties (list of strings) – properties of the lightcurves to include in the file.

Todo

We might need here an extra argument that specifies how to format the properties.

pycs.gen.util.datetimefromjd(JD)

Copy and past from cosmouline. Can be of use here to plot lightcurves with nice dates.

Returns the Gregorian calendar (i.e. our “normal” calendar) Based on wikipedia:de and the interweb :-)

Parameters:JD (float) – julian date
Return type:datetime object
Returns:corresponding datetime
pycs.gen.util.flatten(x)

Source : http://kogs-www.informatik.uni-hamburg.de/~meine/python_tricks flatten(sequence) -> list

Returns a single, flat list which contains all elements retrieved from the sequence and all recursively contained sub-sequences (iterables).

Examples:

>>> [1, 2, [3,4], (5,6)]
[1, 2, [3, 4], (5, 6)]
>>> flatten([[[1,2,3], (42,None)], [4,5], [6], 7, MyVector(8,9,10)])
[1, 2, 3, 42, None, 4, 5, 6, 7, 8, 9, 10]
pycs.gen.util.strtd(td)

To print out time differences … Could be improved a bit :-)

pycs.gen.util.zipdirs(pattern='rrs_*')

I will tgz all directories matching the pattern, except if the tgz already exists. (handy to transfert multirun output to another computer …)

Module contents

This subpackage contains general classes and functions to manipulate and display lightcurves and microlensing representations. “gen” stands for general : all this is as independent as possible from the actual curve shifting techniques.

The lc module is the central module of PyCS, it defines the lightcurve class.