PyCS : Python Curve Shifting

About

Warning

This package is no longer maintained, please consider using the new version for python3 PyCS3.

PyCS is a software toolbox to estimate time delays between multiple images of strongly lensed quasars, from resolved light curves such as obtained by the COSMOGRAIL monitoring program. It comes in the form of a python package, and heavily depends on numpy, scipy, and matplotlib for its core functionality. The repository is on GitHub.

To measure time delays with pycs, you’ll typically write a script calling some high-level functions provided by the package. PyCS allows you to compare different point estimators (including your own), without much code integration.

If you have already read our papers, you might want to proceed with Download & Installation, or the Tutorial. To get a quick first impression of how PyCS looks in practice, install PyCS and then go directly to the Demonstration scripts. These show you how to reproduce the figures from the method paper.

Warning

Please read this important warning about using PyCS.

Questions ?

Feel free to post an issue on GitHub, or to contact the code authors Malte Tewes and Vivien Bonvin.

Contents

Warning & Disclaimer

By itself, the use of PyCS does not guarantee realistic time-delay uncertainty estimates. All we have (blindly) demonstrated regarding the quality of PyCS uncertainty estimates assumes a careful generation of the simulated light curves. Keep in mind that the PyCS approach involves setting values or ranges for several parameters. In principle, the impact of these settings has to be assessed on each light curve set, as described in the PyCS paper (Tewes et al. 2013).

Warning

In particular, obtaining uncertainty estimates by running on simulated light curves that all have a same true delay (i.e., setting truetsr = 0, in PyCS terminology) can lead to vastly underestimated error bars.

It is of high importance to use simulations with a range of plausible true time delays to tune and/or verify the accuracy and precision of time-delay estimators. Tests on simulations with only a single true time delay do not probe reliably the quality of a time-delay estimation, as many time-delay estimators are prone to responding unsteadily to the true delay. Do not hesitate to contact us (see front page of this documentation) in case of doubts!

Note

Please try to avoid publishing time delays with overly optimistic uncertainty estimates. Seemingly accurate measurements might be propagated into unrealistic cosmological inferences, and harm the community. Be critical, and aim for the most comprehensive instead of the smallest error bar.

Thanks!

Intra-dependency chart

This chart lists the intra-dependencies of the various modules and files in PyCS.

Last updated on 2018-09-28. Open the image in a new tab for a better experience.

pycs package

Subpackages

pycs.disp package

Submodules

pycs.disp.disps module

Defines dispersion functions. No optimization is done here ! We just calculate the dispersion.

• linintnp is the “default” one
• pelt95
• others should be ported into pycs from my f90 code :-(

They all share a common structure : give me 2 light curves as input, and you get a float describing the match. By construction the dispersion is often not symmetric, i.e. if you switch lc1 and lc2 you will get a different result. To symmetrize the techniques, use the provided symmetrize function. I return a dictionnary, containing a field “d2” that gives the dispersion (among possible other fields). Any microlensing is naturally taken into account, if present, simply using the getmags() method.

pycs.disp.disps.linintnp(lc1, lc2, interpdist=30.0, weights=True, usemask=True, plot=False)

This is the default one, fast implementation without loops (25 times faster then maltef90). If usemask == True, it is a bit slower… If you do not change the mask, it’s better to cutmask() it, then use usemask = False.

pycs.disp.disps.pelt95(lc1, lc2, decorlength=3.0, usemask=True, verbose=False, plot=False)

Dispersion method as described in Pelt Kayser Refsdal 93 & 95 For now this corresponds to D3, the third form, i.e. using only stricly neigbooring pairs, and no weighting. The fourth form uses all pairs within a given decorlength, non only neighboring.

pycs.disp.disps.linint90(lc1, lc2, interpdist=30.0, weights=True, verbose=False, plot=False, debug=False)

Do not use, historical only. Terribly slow variant of linintnp, with plenty of loops, mimics my f90 implementation of 2007…

Between two ORDERED light curves lc1 and lc2 It does skip masked points.

IT IS SLOW !!! THIS IS FOR DEMONSTRATION PURPOSES ONLY !

About the maths : interpolate linearly between the points of lc2, and sum the square of the differences between the points of lc1 and this interpolation. If you choose to weight this sum according to errors, then for each point-pair (i.e. “element”) of the sum, the weight is calculated using both the error on lc1 and the interpolated (again !) error on lc2

We try to do a minimum of loops, using the fact that the light curves are ordered.

lc1 will be the “reference”, and we will interpolate between the points of lc2 This algorithm absolutely NEEDS ordered light curves !!!!

The way we check for the mask looks a bit un-natural, but this way should be quite fast in fact. Another think to note is that there is for now no special procedure for “close-by” points, and subtraction of very close numbers is not good…

@return: a dict, containing for now :

• “n” : the number of “pairs” that could be used
• “d2” : the chi2 ( i.e. sum(squares)/(n-1) )

@type lc1: lightcurve @param lc1: The reference light curve. It can of course be shifted and or partially masked. @type lc2: lightcurve @param lc2: The light curve that will be interpolated. Same remark.

@type interpdist: float @param interpdist: a maximum timespan (in days) over which interpolations are done. See the code for exact details. If you are too conservative (low value), you might perhaps not use some isolated points at all… @type weights: boolean @param weights: True if you want to use the errorbars of the lightcurves @type verbose: boolean @param verbose: True if you want to see details about the curve matching @type plot: boolean @param plot: True if you want to see the actual interpolation as a matplotlib plot

@type debug: boolean @param debug: Turn this on to make it even slower…

pycs.disp.disps.symmetrize(lc1, lc2, dispersionmethod)

Calls your dispersion method on (A,B) and on (B,A) and returns the “average”.

pycs.disp.multiopt module

Functions to optimize timeshifts, microlensing, or all this at the same time between an arbitrary set of curves. These are building blocks, to be assembled to make a general optimizer.

pycs.disp.multiopt.opt_magshift(lcs, rawdispersionmethod, verbose=True, trace=False)

I optimize the magnitude shifts between lightcurves. This is usually a first thing to do. First curve does not get shifted.

You might think that this can be done by just using the median mag of a curve, but here we use the dispersion method, to get it right even in case of curves with non-overlapping regions.

pycs.disp.multiopt.opt_ml(lcs, rawdispersionmethod, verbose=True, maxit=3, stoprel=0.01, maxpowellit=5, splflat=True, trace=False)

Optimize the ml params of the lcs using the Powell method. Generic, give me as many lcs as you want, with either polynomial or spline microlensings.

Fluxshift optimization could perhaps be included in the loop later, for now we skip this.

If some of the lightcurves have no microlensing, I will “aim” at them for my first iteration. So the first iteration is special.

pycs.disp.multiopt.opt_ts_mix(lcs, rawdispersionmethod, movefirst=False, verbose=True)

Optimize the time shifts (nothing else) with a subtle mix of optimization methods.

Steps: * Run a simulated annealing, with a lower temperature bound and no other stop condition, i.e. a more or less constant number of iterations. * Starting from the best point found, explore a 5 x 5 cube * idem, but smaller * idem, but smaller * idem, but smaller * from this last point, simplex minimization

Expect about 600 calls of the rawdispersion on the couples, i.e. 600*12 = 7200 for 4 lightcurves.

Todo

Implement trace !

pycs.disp.multiopt.bruteranges(step, radius, center)

Auxiliary function for brute force exploration. Prepares the “ranges” parameter to be passed to brute force optimizer In other words, we draw a cube … radius is an int saying how many steps to go left and right of center. center is a list or array of the centers, it can be of any lenght.

You make 2*radius + 1 steps in each direction !, so radius=2 means 5 steps thus 125 calls for 4 curves.

pycs.disp.topopt module

Global optimizers that make use of the building blocks from multiopt.py They return a float, the dispersion value. They can be used with pycs.sim.run.multirun() We do not redefine the ML, but we keep it.

pycs.disp.topopt.opt_full(lcs, rawdispersionmethod, nit=5, verbose=True)

A first optimization of ml and ts, alternatively doing one after the other. Works great, seems to be a keeper. Note that I do keep the ML.

You can put a rawdispersionmethod (i.e. non symmetric) into these guys, as they will anyway apply it on AB and BA for every pair AB.

Module contents

This subpackage defines the “dispersion-inspired” family of curve shifting methods. It contains functions that optimize lightcurves (shifts, ML) so to minimize a given dispersion between them. Dispersion functions are defined in the module pycs.disp.disps.

• multiopt defines the building blocks of these optimizers.
• topopt assembles thoes blocks into wrapper functions.

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.

pycs.regdiff package

Submodules

pycs.regdiff.multiopt module

pycs.regdiff.pymcgp module

pycs.regdiff.rslc module

Module contents

pycs.regdiff2 package

Submodules

pycs.regdiff2.gpr module

pycs.regdiff2.multiopt module

pycs.regdiff2.rslc module

Module contents

pycs.sim package

Submodules

pycs.sim.draw module

This module generates fake curves by “drawing” them from common sourcesplines, adding noise, tweaking microlensing.

pycs.sim.draw.sample(lc, spline)

You give me a lightcurve, and I will set the lc.mags according to the spline. I DO modify the lc object in place ! I do NOT modify the spline, of course.

If your lightcurve has mircolensing, magshifts, or fluxshift, and of course also time shifts, I will set its mags so that it matches the spline with all of these applied !

If you think about it, this is WONDERFUL !!! To do this, I have to do something like a reverse getmags().

Note that I do not add any random noise here, this is sampling only.

pycs.sim.draw.saveresiduals(lcs, spline)

You give me some optimized lcs, and the fitting spline. This function saves the residuals of the optimized lightcurves (as an attribute), in terms of non-fluxshifted magnitudes. This means that you will have to ADD these residuals (in terms of magnitude) before any fluxshifts to the spline mags in order to recover the lc. This is made so that you can reuse this particular shotnoise once fake curves are drawn.

We have to keep this function separate from draw(), as you might want to shift the curves after saving the residuals…

Warning

Call this function if your spline matches to the lcs !

pycs.sim.draw.transfershifts(lcs, reflcs, transferml=True)

I will put (copies) of all the shifts and ML from the reflcs onto the lcs. I will check that you didn’t mix up any curves.

pycs.sim.draw.draw(lcs, spline, shotnoise=None, shotnoisefrac=1.0, tweakml=None, scaletweakresi=True, tweakspl=None, keepshifts=True, keeptweakedml=False, keeporiginalml=True, trace=False, tracedir='draw', inprint_fake_shifts=None)

Wrapper to produce one set of fake lightcurves that are similar to your lcs. Give me some lightcurves, that are optimized in fluxshift, magshift, microlensing, timeshift, and the resulting spline fit representing the intrinsic variations. I will tweak the source spline and the microlensing etc, and return to you a list of fake lightcurves. These lightcurves are just blank data points, “as observed”. I will build them from scratch. But note that by default the drawn light curve will be equiped with the same shifts and ML as your lcs.

Note

I do NOT modify lcs or spline !

Note

It is perfectly ok to shift your input curves in time after the optimization, to make me build fake lightcurves with other delays ! So if you want me to draw curves that look different from yours, it is perfectly ok to give me lcs that do not match to the spline ! Same logic applies for all other shifts and microlensing.

Parameters:

• lcs – reference lightcurves to “immitate”. I will use their epochs, their time/mag/flux shifts, their microlensing, their errorbars.
• spline – reference spline from which I will draw my magnitudes.
• shotnoise –

  Select among [None, “magerrs”, “res”, “mcres”, “sigma”]

  It tells what kind of noise to add to the fake mags. This noise will, as required, be added to the “observed” data (i.e. not fluxshifted).

• shotnoisefrac – a multiplier of the shotnoise added to the curve. Set to 0.5 and I’ll add only “half” of the noise …
• tweakml – either a function, or a list of functions, that takes a list of lightcurves and tweaks their ML in place. I will use this tweaked ML to draw them. If you give a single function, I will apply it to all the curves If you give a list of functions, I will apply them to the respective curves of your lcs (i.e., give the functions in the order corresponding to your lcs !).
• scaletweakresi – scales the “residuals” obtained by tweakml according to saved residuals
• tweakspl – a function that takes a spline and returns a tweaked spline. I will use this on the sourcespline you pass me, before drawing from it.
• keepshifts – by default I will set the time/flux/mag/ shifts from your lcs also to the fake curves.
• keeptweakedml – if keepshifts is True, and keeptweakedml is True, I will keep the tweaked ML, not the input ML, on the output curves. It makes no sens to keep the ML if not keeping the shift
• keeporiginalml – if keepshifts is True, and keeporiginalml is True, I will keep the the input ML, on the output curves.
• inprint_fake_shifts – give an array of shifts corresponding to your lcs that you want to inprint in the mock curves.

Note

I will tweak the ML only of those curves that have spline ML. You probably want me to tweak the ML of all your curves ! So be sure to add some spline ML to all your curves before calling me !

pycs.sim.draw.multidraw(lcs, spline=None, optfctnots=None, onlycopy=False, n=20, npkl=5, simset='draw', simdir=None, shotnoise=None, shotnoisefrac=1.0, truetsr=8.0, tweakml=None, scaletweakresi=True, tweakspl=None, shuffle=True, verbose=True, trace=False, destpath='./')

Even higher wrapper to produce mock + tweaked lightcurves, and save them into a directory (as pickle files), in preparation for analysing them with pycs.sim.run.multirun()

The curves I return are “like observed” : they have no shifts, no ML. Just datapoints.

Parameters:

• lcs (list) – The starting-point lightcurves from which I will draw the mock data. They must match (!), i.e., have appropriate microlensing and be optimized somehow.
• spline (spline) – The source spline used to draw the new curves. Is not used (-> None) if onlycopy=True, or if optfct is specified.
• optfctnots (function) – A function to fit a spline and the ML + all shifts execpt for timeshifts. It is called after setting the true time shifts. Put at None if you want to use always the same current ML and spline.
• onlycopy – If True, I will simply save copies of the input lcs, not drawing anything.
• n – number of lightcurve-lists to simulate per pickle
• npkl – number of pickle files
• simset – give a name to your simulation !
• simdir – where should I put these simulations ?
• shotnoise (string) – Select among None, “magerrs”, “mcres”, “res”. See definitions in pycs.sim.draw.draw().
• truetsr – radius of exploration for unifomly slectrec true time shifts (passed to multidraw) Not used if draw = False, of course.
• shuffle – Shuffle the curves before feeding them into optfctnots ? If you use this, be sure to pycs.gen.lc.objsort the curves before !

Todo

Possibility to add shot noise even with draw = False. Like in the good old days. -> Separate addshotnoise from multidraw ?

Note

I will tweak the ML only of those curves that have spline ML. You probably want me to tweak the ML of all your curves ! So be sure to add some spline ML to all your curves before calling me !

pycs.sim.draw.shareflux(lc1, lc2, frac=0.01)

I add “noise” to lc1 and lc2 by randomly sharing flux between the two sources.

Parameters:frac – The stddev of the gaussian “noise” in flux, with respect to the minimum flux in the curves.

pycs.sim.plot module

Subpackage with functions to plot all kind of results from runs.

pycs.sim.plot.mad(xs)

Return the median absolute deviation. Write it myself here instead of importing it from astropy, since it will add another depenency. Work with 1d array only

@todo: for PyCS 3, will use astropy as a default module (good) and use their functions

Parameters:xs – list of values Returns:median absolute deviation

class pycs.sim.plot.delaycontainer(data, name='No name', objects=None, plotcolour='black', marker=None)

Stores the delay or error bar measurement(s) (one for each curve pair). This object is usually produced by the plot-functions hists or meanvstrue below.

markers : [ 7 | 4 | 5 | 6 | ‘o’ | ‘D’ | ‘h’ | ‘H’ | ‘_’ | ‘’ | ‘None’ | ‘ ‘ | None | ‘8’ | ‘p’ | ‘,’ | ‘+’ | ‘.’ | ‘s’ | ‘*’ | ‘d’ | 3 | 0 | 1 | 2 | ‘1’ | ‘3’ | ‘4’ | ‘2’ | ‘v’ | ‘<’ | ‘>’ | ‘^’ | ‘|’ | ‘x’ | ‘$…$’ | tuple | Nx2 array ]

__init__(data, name='No name', objects=None, plotcolour='black', marker=None)

self.data is a list of dicts, the fields depend on if it’s delays or errorbars

• delays : label, mean, med, std
• errorbars : label, tot, sys, ran, bias

__module__ = 'pycs.sim.plot'

pycs.sim.plot.newdelayplot(plotlist, rplot=7.0, displaytext=True, hidedetails=False, showbias=True, showran=True, showerr=True, showlegend=True, text=None, figsize=(10, 6), left=0.06, right=0.97, top=0.99, bottom=0.08, wspace=0.15, hspace=0.3, txtstep=0.04, majorticksstep=2, filename=None, refshifts=None, refdelays=None, legendfromrefdelays=False, hatches=None, centershifts=None, ymin=0.2, hlines=None, tweakeddisplay=False, blindness=False, horizontaldisplay=False, showxlabelhd=True)

Plots delay measurements from different methods, telescopes, sub-curves, etc in one single plot. For this I use only delaycontainer objects, i.e. I don’t do any “computation” myself.

Parameters:plotlist – Give me a list of tuples (delays, errorbars), where delays and errorbars are delaycontainer objects as written into pkl files by hists and meanvstrue.

NEW : plotlist delaycont can handle asymmetric errors (e.g. to seamlessly compare pycs with other papers’ results). Instead of the “tot” key, new “plus” and “minus” keys are used. :type plotlist: list

Parameters:

• rplot – radius of delay axis, in days.
• displaytext (boolean) – Show labels with technique names and values of delays
• hidedetails (boolean) – Do not show (ran, sys) in labels
• refshifts (list) – This is a list of dicts like {“colour”:”gray”, “shifts”:(0, 0, 0, 90)}. Will be plotted as dashed vertical lines.
• refdelays (list) – a list of tuples (delays, errorbars) to be plotted as shaded vertical zones.
• legendfromrefdelays (boolean) – if you want to display the refdelays name in the legend panel
• hatches (list) – list of hatch keyword for the refdelays plotting
• showbias (boolean) – draws a little cross at the position of the delay “corrected” for the bias.
• showran (boolean) – draws “minor” error bar ticks using the random error only.
• text (list) – Text that you want to display, 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).
• blindness (boolean) – Shift the measurements by their mean, so the displayed value are centered around 0
• horizontaldisplay (boolean) – display the delay panels on a single line. Works only for three-delay containers.
• showxlabelhd (boolean) – display or not the x label when horizontal display is True

Warning

Altough the code says I’m plotting mean and std for the measured delays, I might be using median and mad instead! This depends on how hists was called! Be careful with this…

pycs.sim.plot.newdelayplot2(plotlist, rplot=7.0, displaytext=True, hidedetails=False, showbias=True, showran=True, showlegend=True, text=None, figsize=(10, 6), left=0.06, right=0.97, top=0.99, bottom=0.08, wspace=0.15, hspace=0.3, txtstep=0.04, majorticksstep=2, filename='screen', refshifts=None, refdelays=None, legendfromrefdelays=False, hatches=None, centershifts=None, ymin=0.2, hlines=None)

Plots delay measurements from different methods, telescopes, sub-curves, etc in one single plot. For this I use only delaycontainer objects, i.e. I don’t do any “computation” myself.

Difference from newdelayplot is that the previously hatched/shaded regions are plotted as smaller points, without infos on the time-delay

Parameters:

• plotlist (list) – Give me a list of tuples (delays, errorbars), where delays and errorbars are delaycontainer objects as written into pkl files by hists and meanvstrue.
• rplot – radius of delay axis, in days.
• displaytext (boolean) – Show labels with technique names and values of delays
• hidedetails (boolean) – Do not show (ran, sys) in labels
• refshifts (list) – This is a list of dicts like {“colour”:”gray”, “shifts”:(0, 0, 0, 90)}. Will be plotted as dashed vertical lines.
• refdelays (list) – a list of tuples (delays, errorbars) to be plotted as shaded vertical zones.
• legendfromrefdelays – if you want to display the refdelays name in the legend panel
• hatches – list of hatch keyword for the refdelays plotting
• showbias (boolean) – draws a little cross at the position of the delay “corrected” for the bias.
• showran (boolean) – draws “minor” error bar ticks using the random error only.
• text (list) – Text that you want to display, 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).

pycs.sim.plot.normal(x, mu, sigma)

Plain normal distribution. You can directly apply me on numpy arrays x, mu, sigma.

pycs.sim.plot.hists(rrlist, r=10.0, nbins=100, showqs=True, showallqs=False, qsrange=None, title=None, xtitle=0.5, ytitle=0.95, titlesize=18, niceplot=False, displaytext=True, figsize=(16, 9), left=0.06, right=0.95, bottom=0.065, top=0.95, wspace=0.2, hspace=0.2, txtstep=0.04, majorticksstep=2, hideyaxis=True, trueshifts=None, filename=None, dataout=False, blindness=False, usemedian=False, outdir='./')

Comparing the delay distributions from different run result objects.

Parameters:

• rrlist – a list of runresults object.
• r – a range radius for the hists
• showqs – If True, I overplot the qs as scatter points.
• dataout – True means that I’ll write the pkl file needed to make the delayplot.
• removeoutliers – True means I remove estimates that are the farthest from the median. Use this with CAUTION !!!
• usemedian – if True, use the median and median absolute deviation instead of mean and std.

Warning

To avoid rewriting newdelayplot, if usemedian is True then I write the median and mad in the mean and std fields of the pickles. This is dangerous (and a bit stupid and lazy), but since hists() and newdelayplot() are usually called one after the other it should not create too much confusion.

Note

Actually, using median and mad as default estimators might be smarter…? To meditate for PyCS 3.0…

pycs.sim.plot.newcovplot(rrlist, r=6, rerr=3, nbins=10, nbins2d=3, binclip=True, binclipr=10.0, figsize=(13, 13), left=0.06, right=0.97, top=0.97, bottom=0.04, wspace=0.3, hspace=0.3, method='indepbin', minsamples=10, showplots=True, printdetails=True, printcovmat=True, detailplots=False, filepath=None, verbose=True)

pycs.sim.plot.measvstrue(rrlist, r=10.0, nbins=10, plotpoints=True, alphapoints=1.0, plotrods=True, alpharods=0.2, ploterrorbars=True, sidebyside=True, errorrange=None, binclip=False, binclipr=10.0, title=None, xtitle=0.75, ytitle=0.95, titlesize=30, figsize=(10, 6), left=0.06, right=0.97, top=0.99, bottom=0.08, wspace=0.15, hspace=0.3, txtstep=0.04, majorticksstep=2, displayn=True, filename=None, dataout=False, tweakeddisplay=False, blindness=False, outdir='./')

Plots measured delays versus true delays

Parameters:

• r – radius of simulation input delays to plot (x axis range)
• nbins – number of bins for the bar plot within this range.
• plotpoints – should I plot the points (scatter plot) ?
• plotrods – should I plot the avg within each bin ?
• ploterrorbars – should I add errorbars upon the bar plot ?
• sidebyside – should I plot bars side by side, or overplot them ?
• errorrange – radius of measurement errors to plot (y axis range). You can also give a tuple (low, high), to make asymetric plots.
• binclip – should I clip errors larger than binclipr days (catastrophic failures of methods) ?
• binclipr – see binclip …

pycs.sim.plot.covplot(rrlist, showpoints=False, showcontour=True, showdensity=False, fractionalresiduals=False, bins=50, smoothing=0.0, figsize=(12, 12), left=0.02, right=0.98, bottom=0.02, top=0.98, wspace=0.05, hspace=0.05, r=5.0, title=None, txtstep=0.04, filename=None)

Covariance scatter of all measurement errors. Give me a single runresults object (from a sim, with known true delays).

pycs.sim.run module

Functions to run curve shifting techniques on lightcurves produced by sim.multidraw. We define the class runresults that “holds” results obtained by curve shifting techniques (saved into pickles).

pycs.sim.run.applyopt(optfct, lcslist, **kwargs)

Applies optfct (an optimizing function that takes a list of lightcurves as single argument) to all the elements (list of lightcurves) of lcset (a list of lists of lightcurves).

Optimizes the lightcuves themselves, in place, and returns a list of the ouputs of the optimizers, corresponding to the lcslist. For instance, if the optfct output is a spline, it also contains the final r2s, that I will later save into the pkls !

About multi cpu : First try using the multiprocessing module -> failed, as the optfct cannot be pickled. Perhaps need to rethink the strategy. Second try using good old forkmap… it works !

ncpu : None = I will use all CPUs, -1 = I will use all - 1 CPUs, and otherwise I will use ncpu CPUs.

class pycs.sim.run.runresults(lcslist, qs=None, name='None', plotcolour='#008800', success_dic=None)

Summarizes the huge list of list of lightcurves as a numpy array of timeshifts and some further info, to serve as input for plots, and actual time delay determinations. This replaces the old “boot.pkl” files … The TRUE shifts are also saved (if available)

All this is not related to a particular optimization technique. Please also provide the success_dic to remove the curves where the optimiser failed.

Note the mask funcitonallity.

__init__(lcslist, qs=None, name='None', plotcolour='#008800', success_dic=None)

lcslist may or may not have “truetimeshifts”. If not, I will put 0.0 as true shifts.

qs should be a numpy array (as long as lcslist) that contains some chi2, r2, or d2 stuff to quantify how good the fit was.

All the lcs in lcslist should be in the same order (I will check this, even if slow). I will not sort these lcs (as you might want them in an unsorted order).

__len__()

The number of runs

nimages()

The number of images (4 for a quad, 2 for a double) …

__str__()

copy()

check()

applymask(mask)

Removes some of the runresults according to your mask.

gettruets()

Returns some summary stats about the true delays. Used to find histogram ranges for plots, etc.

getts()

A bit similar to gettruets, we return the median of the measured ts… Used for plots etc, not for calculations.

get_delays_from_ts()

Return the time delays, from the timeshifts. I do not account for the true timeshift. :return: dictionary containing the median, max, and min delays + delay labels

__module__ = 'pycs.sim.run'

pycs.sim.run.joinresults(rrlist)

Give me a list of runresults objects, I join those into a single one an return the latter.

pycs.sim.run.collect(directory='./test', plotcolour='#008800', name=None)

Collects the runresult objects from a directory (typically from a multirun run), and returns the joined runresults.

pycs.sim.run.multirun(simset, lcs, optfct, kwargs_optim=None, optset='multirun', tsrand=10.0, analyse=True, shuffle=True, keepopt=False, trace=False, verbose=True, destpath='./')

Top level wrapper to get delay “histograms” : I will apply the optfct to optimize the shifts between curves that you got from pycs.sim.draw.multidraw(), and save the results in form of runresult pickles.

It is perfectly ok to launch several instances of myself on the same simset, to go faster. I will process every pkl of the simset only once, and prevent other instances from processing the same files.

You can use me for a lot of different tasks. (note from VB : not to make coffee apparently)

Parameters:

• simset – The name of the simulations to run on. Those are in a directory called sims_name.
• lcs – Lightcurves that define the initial shifts and microlensings you want to use. I will take the lightcurves from the simset, and put these shifts and ML on them.
• kwargs_optim – kwargs to be passed to your optfct
• optset – A new name for the optimisation.
• optfct (function) – The optimizing function that takes lcs as single argument, fully optimizes the curves, and returns a spline, or a d2 value. Can be None if argument analyse is False (used for tests).
• tsrand – I will randomly shift the simulated curves before running the optfct This randomizes the initial conditions. (uniform distrib from -tsrand to tsrand)
• shuffle – if True, I will shuffle the curves before running optc on them, and then sort them immediatly afterwards.
• keepopt – a bit similar to Trace, but simpler : we write the optimized lightcurves as well as the output of the optimizers into one pickle file per input pickle file. {“optfctoutlist”:optfctouts, “optlcslist”:simlcslist}

pycs.sim.run.clean_simlist(simlcslist, success_dic)

pycs.sim.src module

Stuff to represent a light curve as points on a regular grid, and power spectrum plots. This allows e.g. to tweak its power spectrum by adding “correlated noise”, resample it, etc.

class pycs.sim.src.Source(spline=None, name='Source', range=(0, 10000), sampling=0.2)

A class representing a “source” lightcurve, i.e. any artificial signal whose powerspectrum we can modify, and still evaluate the curve at given irregular jds.

To do this, we work with “internal regular arrays” (and do FFTs and so on on these), and then interpolate these fine arrays to get back values for given jds.

We store the data in magnitudes.

Ideas : use the mask carried around by self.inispline.datapoints to define regions that really matter for the power specturm and stuff.

__init__(spline=None, name='Source', range=(0, 10000), sampling=0.2)

range are in jds, only used if you don’t give me a sourcespline. sampling is in days.

I will store the spline as inispline, but I will not modify it, so no need to pass me a copy.

setup()

Sets the internal regular arrays by sampling the spline. We update self.sampling, so to get an even number of points (for FFTs).

__str__()

copy()

setmean(mean=-12.0)

Shifts the magnitudes so that their mean value are at the specified level. Don’t do this, just for testing purposes ! When doing flux PS, this scales the power, that’s all.

addgn(sigma=0.1, seed=None)

Don’t do this, just for testing purposes ! There is no reason to add white noise to the source !

addrw(sigma=0.1, seed=None)

Add a random walk, also for experiement (power law param is -2)

addplaw2(beta=-3.0, sigma=0.01, flux=False, fmin=None, fmax=None, hann=False, seed=None)

Next version, better Adds noise according to a power law PSD. See Timmer & Koenig 1995

power law would be -2

if hann, we soften the window (Hann window instead of tophat).

eval(jds)

I interpolate my ijds/imags to give you an array of mags corresponding to your jds This could in principle be done using the spline object made by the function below, but this is safer and faster.

spline()

I return a new pycs.gen.spl.Spline object corresponding to the source. So this is a bit the inverse of the constructor. You can then put this spline object as ML of a lightcurve, as source spline, or whatever.

Note that my output spline has LOTs of knots… it is an interpolating spline, not a regression spline !

..note:: This spline is a priori for display purposes only. To do an interpolation, it might be safer (and faster) to use

the above linear interpolation eval() function. But making a plot, you’ll see that this spline seems well accurate.

__module__ = 'pycs.sim.src'

pycs.sim.src.sourceplot(sourcelist, filename=None, figsize=(12, 8), showlegend=True, showspline=True, marker=None)

I show you a plot of a list of Source objects.

pycs.sim.src.window_hanning(x)

return x times the hanning window of len(x)

pycs.sim.src.detrend_mean(x)

Return x minus the mean(x)

class pycs.sim.src.PS(source, flux=False)

A class representing a power spectrum of a Source object.

__init__(source, flux=False)

Constructor, simply takes a Source object and calculates the power spectrum

The Source is always expressed in magnitudes, but you might want to get a powerspectrum in terms of flux or “counts”. Put flux = True, and I will convert your magnitudes into fluxes before calculating the power spectrum.

__module__ = 'pycs.sim.src'

__str__()

copy()

calcslope(fmin=0.001, fmax=0.5)

Measures the slope of the PS, between fmin and fmax. All info about this slope is sored into the dictionnary self.slope.

Should this depend on self.flux ?? No, only the constructor depends on this ! This is just fitting the powerspectrum as given by the constructor.

pycs.sim.src.psplot(pslist, nbins=0, filename=None, figsize=(12, 8), showlegend=True)

Plots a list of PS objects. If the PS has a slope, it is plotted as well.

if nbins > 0, I bin the spectra.

add option for linear plot ?

pycs.sim.twk module

Top level functions that use the src module to tweak microlensing and source splines. These are the function to pass them to draw.draw or draw.multidraw.

pycs.sim.twk.tweakml(lcs, spline, beta=-2.0, sigma=0.05, fmin=0.002, fmax=None, psplot=False, sampling=0.1)

I tweak the SplineML of your curves by adding small scale structure. I DO modify your lcs inplace.

pycs.sim.twk.tweakspl(spline, beta=-2.5, sigma=0.03, fmin=0.03333333333333333, fmax=0.2, hann=False, psplot=False)

Give me a spline, I return a tweaked version with added small scale structure. Note that the spline I return will have a LOT of knots.

I DO NOT modify your spline, but return a new one.

pycs.sim.twk.addspl(spl1, spl2, op='add')

Give me two splines, I return the addition/substraction of the two: spl1 +/- spl2

Important: I assume that the two splines have the correct jds ! The resulting spline

Parameters:

• spl1 – first spline
• spl2 – second spline
• op – ‘add’ for addition, ‘sub’ for subtraction

Returns:

spline

Module contents

Simulating “real” light curves, started from scratch in summer 2011

pycs.spl package

Submodules

pycs.spl.multiopt module

Functions to optimize time shifts and microlensing between lcs, using spline fits.

By default the functions don’t touch the ML and sourcespline knots, it’s up to you to enable BOK iterations. And check the results by eye …

Typical order of operations :

Put smooth microlensing on the curves Fit a source to the curve without microlensing opt_magshift (does not use the source) opt_ml (does not change the source) Fit a new source, to all the curves Put a less smooth microlensing on them opt_ml_source

pycs.spl.multiopt.opt_magshift(lcs, sourcespline=None, verbose=False, trace=False)

If you don’t give any sourcespline, this is a dirty rough magshift optimization, using the median mag level (without microlensing), once for all. We don’t touch the magshift of the first curve.

If you do give a sourcespline, I’ll optimize the magshift of each curve one by one to match the spline. New : I even touch the magshift of the first curve.

We use the l1-norm of the residues, not a usual r2, to avoid local minima. Important ! This is done with the nosquare=True option to the call of r2 !

pycs.spl.multiopt.opt_source(lcs, sourcespline, dpmethod='extadj', bokit=0, bokmethod='BF', verbose=True, trace=False)

Just the source spline, without touching the ML of the lcs. At each call, I update the sourcespline with the merged lcs. The internal knots of the sourcespline stay where they are, only the external ones are ajusted.

pycs.spl.multiopt.opt_fluxshift(lcs, sourcespline, verbose=True, trace=False)

Optimizes the flux shift and the magshift of the lcs (not the first one) to get the best fit to the “sourcespline”. Does not touch the microlensing, nor the spline. So this is a building block to be used iteratively with the other optimizers. Especially of the sourcespline, as we fit here even on regions not well constrained by the spline !

The spline should typically well fit to the first curve.

pycs.spl.multiopt.opt_ml(lcs, sourcespline, bokit=0, bokmethod='BF', splflat=False, verbose=True, trace=False)

Optimizes the microlensing of the lcs (one curve after the other) so that they fit to the spline. I work with both polynomial and spline microlensing. NO YOU DONT !!! For spline micorlensing, I can do BOK iterations to move the knots.

Note

Does not touch the sourcespline at all !

But this it what makes the problem linear (except for the BOK iterations) for both splines and polynomial ML, and thus fast !

Parameters for spline ML :

Parameters:

• bokit –
• bokeps –
• boktests –
• bokwindow –
• splflat –

Parameters for poly ML :

None ! We just to a linear weighted least squares on each season ! So for poly ML, all the params above are not used at all.

We do not return anything. Returning a r2 would make no sense, as we do not touch the sourcepline !

pycs.spl.multiopt.redistribflux(lc1, lc2, sourcespline, verbose=True, maxfrac=0.2)

Redistributes flux between lc1 and lc2 (assuming these curves suffer form flux sharing), so to minimize the r2 with respect to the sourcespline. I do not touch the sourcespline, but I do modify your curves in an irreversible way !

Parameters:

• lc1 – a lightcurve
• lc2 – another lightcurve
• sourcespline – the spline that the curves should try to fit to

pycs.spl.multiopt.opt_ts_powell(lcs, sourcespline, optml=True, movefirst=False, verbose=True, trace=False)

If you use this again, implement mlsplflat for optml, might be important !

Optimize the timeshifts of all four lightcurves AND the spline coeffs (but not the knots) so to get the best possible fit.

If optml is True, I will optimize the ml for every trial delays (and then optimize the source again). For this I will always start from the initial settings, I don’t cumulate the ML - source optimizations. This is required as ML - source optimizations are degenarate. If optml is false, this would not be needed, as the spline fitting is unique anyway.

Note that even if I don’t modify the sourcespline knots, I do modify the sourcespline datapoints as the time shifts move !

Improvement ideas : see if better to not optimize the first curves time shift ? Don’t try to optimize the timeshifts at the same time ? Do one after the other, with respect to the first curve ? Special function that takes a spline fit to the first curve as argument, to do this ?

-> Yes, see opt_ts_indi below !

pycs.spl.multiopt.comb(*sequences)

http://code.activestate.com/recipes/502199/ combinations of multiple sequences so you don’t have to write nested for loops

>>> from pprint import pprint as pp
>>> pp(comb(['Guido','Larry'], ['knows','loves'], ['Phyton','Purl']))
[['Guido', 'knows', 'Phyton'],
['Guido', 'knows', 'Purl'],
['Guido', 'loves', 'Phyton'],
['Guido', 'loves', 'Purl'],
['Larry', 'knows', 'Phyton'],
['Larry', 'knows', 'Purl'],
['Larry', 'loves', 'Phyton'],
['Larry', 'loves', 'Purl']]
>>> 

pycs.spl.multiopt.opt_ts_brute(lcs, sourcespline, movefirst=True, optml=False, r=2, step=1.0, verbose=True, trace=False)

If you use this again, implement mlsplflat, might be important.

Given the current delays, I will explore a hypercube (r steps in each direction on each axis) of possible time shift combinations. And choose the shifts that gave the smallest chi2.

For each trial shift, I optimize the sourcespline to fit all curves, and optionally also the ML of every curve.

This is very slow, not really used.

pycs.spl.multiopt.opt_ts_indi(lcs, sourcespline, method='fmin', crit='r2', optml=False, mlsplflat=False, brutestep=1.0, bruter=5, verbose=True, trace=False)

We shift the curves one by one so that they match to the spline, using fmin or brute force. A bit special : I do not touch the spline at all ! Hence I return no r2. Idea is to have a fast ts optimization building block.

The spline should be a shape common to the joined lcs. No need to work on copies, as we do not change the ML or spline iteratively, but only the ML -> nothing can go wrong.

Parma brutestep:  step size, in days Parameters:bruter – radius in number of steps

pycs.spl.topopt module

Higher level wrappers around the multiopt optimizers. These attempt to be a as multi-puropose as possible, but please understand that there is no general “optimal” way to optimize a spline fit. Depening on your data, custom optimizers might be needed. That’s what we did for the TDC.

In principle we do not define the microlensing inside of the optimizer, we really just optimize it.

The functions must return a final optimal spline. This spline object also contains its final r2, saved into the pkl, and used for plots etc.

pycs.spl.topopt.opt_rough(lcs, nit=5, shifttime=True, crit='r2', knotstep=100, stabext=300.0, stabgap=20.0, stabstep=4.0, stabmagerr=-2.0, verbose=True)

Getting close to the good delays, as fast as possible : no BOK (i.e. knot positions are not free), only brute force without optml. Indeed with optml this tends to be unstable.

Note

This function is here to correct for shift errors in the order of 20 days, hence its name.

Parameters:

• nit – a number of iterations (ok to leave default : 5)
• splstep – the step (in days) between the knots of the spline. 100 is default. Lower this to eg 50 if your curve has very short term intrinsic variations, increase it if your curve is very noisy.

pycs.spl.topopt.opt_fine(lcs, spline=None, nit=10, shifttime=True, crit='r2', knotstep=20, stabext=300.0, stabgap=20.0, stabstep=4.0, stabmagerr=-2.0, bokeps=10, boktests=10, bokwindow=None, redistribflux=False, splflat=True, verbose=True, trace=False, tracedir='opt_fine')

Fine approach, we assume that the timeshifts are within 2 days, and ML is optimized.

Watch out for splflat, a bit dangerous all this … Default is True for the iterations, and we release the splines only at the end. If you put it to False, spline are left free at any stage. This should be fine if you leave one curve without ML.

Parameters:spline – If this is None, I will fit my own spline, but then I expect that your lcs already overlap (i.e., ML is set) ! If you give me a spline, I will use it (in this case I disregard your splstep setting !)

pycs.spl.topopt.opt_fine2(lcs, knotstepfact=1.0, maxit=10, minchange=1.0, verbose=True)

Version 2 Generalisation of pycs.tdc.splotp.spl3, but for a random number of lcs. Assumes a really good initial time shift, does not care about inital magshift. Also, contrary to what was done for tdc.splopt, we do NOT add ML in here, but use whatever is there!

Parameters:

• maxit – maximum number of iterations
• minchange – minimum decrease in percent of the r2. I stop if the decrease gets smaller.

Module contents

Simultaneous spline fit technique. This subpackage contains functions to optimize lightcurves (their shifts and ML) so that they fit to a common spline.

• multiopt contains “builing blocks” for such optimizers, and
• topopt assembles these building blocks into big wrappers.

pycs.spldiff package

Submodules

pycs.spldiff.multiopt module

Optimization using regdiff

Would be great to have an opt_flux in here !

pycs.spldiff.multiopt.opt_ts(lcs, method='weights', pd=5, plotcolour=None, knotstep=20.0, n=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, magshift=True)

Give me lightcurves (with more or less good initial time shifts) I run a spline regression on them, optimize the difference, and set their delays to the optimal values. I also optimise the magnitude shifts for display purpose, but it does not affect the time shifts at all !

Parameters:pd – the point density of the regularly sampled lightcurves, in points per days.

pycs.spldiff.rslc module

Defines a class that represents a regularly sampled lightcurve

class pycs.spldiff.rslc.rslc(jds, mags, magerrs, pad, pd, timeshift=0.0, name='Name', plotcolour='black')

A regularly sampled lightcurve, typically obtained by regression. To make such a rslc from a usual lightcurve object, look at the factory function below. One idea is that we want to be able to add and subtract those, propagating errors. There is no “microlensing” or similar stuff – only time shifts.

__init__(jds, mags, magerrs, pad, pd, timeshift=0.0, name='Name', plotcolour='black')

__str__()

shifttime(timeshift)

copy()

getjds()

getmags()

getmagerrs()

mask(maxmagerr=0.1, target=20.0)

wtv(method='weights')

Return some weighted average variation WAV. Usuall called on a “difference” lightcurve.

__module__ = 'pycs.spldiff.rslc'

pycs.spldiff.rslc.factory(l, pad=300, pd=2, plotcolour=None, knotstep=20.0, n=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)

Give me a lightcurve, I return a regularly sampled light curve, by performing some spline regression. !!! New: I also return the spline used for the regression :param pad: the padding, in days :param pd: the point density, in points per days.

The points live on a regular grid in julian days, 0.0, 0.1, 0.2, 0.3 …

pycs.spldiff.rslc.subtract(rs1, rs2)

I subtract rs2 from rs1. This means I keep the jds and timeshift of rs1, and only change the mags and magerrs, interpolating rs2. I return a brand new rslc object, that has no timeshift (as we do not care about a timeshift, for a difference).

Parameters:

• rs1 (rslc object) –
• rs2 (rslc object) –

pycs.spldiff.rslc.wtvdiff(rs1, rs2, method)

Returns the wtv (weighted TV) of the difference between 2 curves.

This is symmetric (no change if you invert rs1 and rs2), up to some small numerical errors.

pycs.spldiff.rslc.bruteranges(step, radius, center)

Auxiliary function for brute force exploration. Prepares the “ranges” parameter to be passed to brute force optimizer In other words, we draw a cube … radius is an int saying how many steps to go left and right of center. center is an array of the centers, it can be of any lenght.

You make 2*radius + 1 steps in each direction !, so radius=2 means 5 steps thus 125 calls for 4 curves.

pycs.spldiff.rslc.opt_ts(rslcs, method='weights', verbose=True)

I optimize the timeshifts between the rslcs to minimize the wtv between them. Note that even if the wtvdiff is only about two curves, we cannot split this into optimizing AB AC AD in a row, as this would never calculate BC, and BC is not contained into AB + AC. !!! New : I also return a spline to optimise the magshifts

Parameters:rslcs – a list of rslc objects

pycs.spldiff.splreg module

spline regression module…

pycs.spldiff.splreg.splreg(jds, mags, magerrs, knotstep=20.0, n=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)

Give me datapoints x y, with errorbars on y I return a function : you pass an array of new x, the func returns (newy, newyerr)

I return a spline as well, so I can use it for magshifts later performed with a BOK spline fitting on the datas and the errors. The errors on y take into account the seasons in y

Module contents

Spline difference curve shifting technique, using spline enveloppes.

pycs.tdc package

Submodules

pycs.tdc.combiconf module

pycs.tdc.combiconf.combiconf1(estimates)

Give me a group of estimates of the same TDC pair, I return the combined confidence. 1 is for doubtless 2 for plausible 3 for multimodal 4 for uninformative

Simple version

@param estimates: list of Estimate objects

@return: dictionary that contains the combined confidence code and the estimate set, rung and pair.

pycs.tdc.combiconf.combiconf2(estimates)

Give me a group of estimates of the same TDC pair, I return the combined confidence. 1 is for doubtless 2 for plausible 3 for multimodal 4 for uninformative

Updated version of combiconf1, now with accordance criterias between the estimates, and weight given to Vivien/Malte estimates

@param estimates: list of Estimate objects

@return: dictionary that contains the combined confidence code and the estimate set, rung and pair.

pycs.tdc.combiconf.reroll(estimates)

Give me a list of estimates. I remove all estimates not from Malte or Vivien, then recompute the combiconfcode Return the new combiconfcode, along with the id of the estimate

@param estimates: list of Estimate objects

@return: dictionnary containing the new combiconfcode computed by combiconf2, the estimate set, rung and pair and the corresponding estimates from Malte and Vivien only

pycs.tdc.combiconf.combiquads(estimates)

This function is based on a exploit discovered during the rungs generation: the quads with the same pair indices are the same amongst all rungs, thus have the same delay.

Not used for TDC1 official submissions, just for fun.

Give me a list of already combined estimates (i.e. one per pair only). I compute for each quad pair the best td and tderr among the corresponding rungs (which is an exploit !) I return the modified estimates list, with all the quads modified, and with only the doubtless double

# BIG WARNING !!! This function (pycs.tdc.combiconf.combiquads) actually modify estimates, we DO NOT want that !

@param estimates: list of Estimate objects. Each id must be present no more than once.

@return: list of modified Estimate objects.

pycs.tdc.est module

Stuff to manipulate time delay estimations from different techniques, specific to the TDC.

class pycs.tdc.est.Estimate(set='tdc0', rung=0, pair=1, method='None', methodpar='None', td=0.0, tderr=0.0, ms=0.0, confidence=0, timetaken=0.0)

Class to hold delay estimates for TDC curves, as obtained by various delay estimators.

__init__(set='tdc0', rung=0, pair=1, method='None', methodpar='None', td=0.0, tderr=0.0, ms=0.0, confidence=0, timetaken=0.0)

Parameters:

• set – What set (e.g., “tdc0”) the curves are from
• rung – rung
• pair – pair
• method – Name of the method
• methodpar – String containing parameters of this method, or username…
• td – time delay point estimate
• tderr – 1sigma error estimate on the time delay
• ms – magnitude shift estimate
• confidence – confidence level. 0 = not estimated, 1 = doubtless, 2 = plausible, …
• timetaken – seconds it took to get this estimate

__str__()

valstr()

aslist()

check()

copy()

getcolor()

setid()

applytolcpair(lcb)

__module__ = 'pycs.tdc.est'

pycs.tdc.est.readcsv(filepath)

Read a CSV file of estimates. You can specify a directory instead of a file, and I will read in all .csv files form there.

pycs.tdc.est.printestimates(estimates)

pycs.tdc.est.writecsv(estimates, filepath, append=False)

Write a CSV file of estimates. If append = True, I do NOT overwrite an existing file, but append to it !

pycs.tdc.est.importfromd3cs(filepath, set='tdc1')

Reads a d3cs log file and returns the list of estimates

pycs.tdc.est.select(estimates, sets=None, rungs=None, pairs=None, idlist=None)

Returns a sublist of the estimates selected according to the specified arguments. It can be sets/rungs/pairs values, or a list of ids

pycs.tdc.est.group(estimates, verbose=True)

Groups estimates by “quasar” In other words : takes a messy list of mixed estimates and returns a list of lists of estimates for a pair.

pycs.tdc.est.checkunique(estimates)

Checks that there is only one estimate per pair

pycs.tdc.est.sort(estimates)

Sorts by set, rung, pair… by id !

pycs.tdc.est.checkallsame(estimates)

Checks that the estimates are all about the same quasar If yes, return the set, rung and pair

pycs.tdc.est.match(candestimates, refestimates)

candestimates and refestimates are two lists of estimates. I return a tuple of two lists :

• those estimates of candestimates which are about quasars that are present in refestimates.
• the remaining candestimates, about quasars not in refestimates

pycs.tdc.est.removebad(estimates, verbose=True, crit='conflevel')

Remove the estimates with bad confidence level (typically 4)

pycs.tdc.est.combine(estimates, method='meanstd', methodcl=None)

Combine estimates according the method of your choice. Return an estimate object

list of methods (add yours): - meanstd - initialestimation - d3cscombi1 - …

Methods should define : td, tderr, confidence, ms, timetaken

pycs.tdc.est.multicombine(estimates, method='meanstd')

Wrapper around combine. Multiple calls to the combine function. Eats a list of messy estimates, groups and combine them by quasar according to method, and return the list of combinated estimates. Can be used to feed writesubmission.

pycs.tdc.est.show(estimates)

Show the estimates

pycs.tdc.est.d3cs(set, rung, pair)

Open d3cs with the rung/pair curve in your default browser

pycs.tdc.est.interactivebigplot(estimates, shadedestimates=None, plotpath=None, interactive=True, minibox=False, groupests=False, minradius=100)

WARNING !! No more TDC0-proof !!

Large graphical representation of your estimates.

Parameters:

• minradius – Minimal half-width of the time delay axes, in day.
• shadedestimates – Here you can give me a list of estimates that I will show as shaded bars instead of errorbars. However, it’s estimates that determines which panels I will draw and with what range, so that you can get panels shown without any shadedestimate.
• plotpath – if you want to save your figures instead of displaying them. Work only when interactive is set to False.
• interactive – Add a button to show the curve with shifts corresponding to the displayed estimates, and another button redirecting to d3cs
• minibox –

  add a minibox to the right of each panel, displaying the mean value of each set of estimate TODO : do the same implementation as shadedestimates, but with a list of different estimates

  also add captions

• groupests – interactive must be set to True. If groupests is True, display a command center allowing the user to navigate between the different curves selected, displaying them 10 by 10.

pycs.tdc.est.bigplot(estimates, shadedestimates=None, plotpath=None, minradius=100)

Large graphical representation of your estimates.

Parameters:

• minradius – Minimal half-width of the time delay axes, in day.
• shadedestimates – Here you can give me a list of estimates that I will show as shaded bars instead of errorbars. However, it’s estimates that determines which panels I will draw and with what range, so that you can get panels shown without any shadedestimate.

pycs.tdc.metrics module

Functions related to the TDC metrics

pycs.tdc.metrics.fN(estimates)

@param estimates: list of Estimate objects @return: length of the estimates list

pycs.tdc.metrics.f(estimates, N)

@param estimates: list of Estimate objects @param N: maximum number of Estimates expected @return: fractional length of the estimates list with respect to N

pycs.tdc.metrics.P(estimates)

@param estimates: list of Estimate objects @return: Approximation of the P metric…

pycs.tdc.metrics.sortbyP(estimates)

I sort your estimates according to their claimed precision lowest “precision” (== highest tderr/td) first -> select from end !

@param estimates: list of Estimate objects @return: sorted list of Estimate objects.

pycs.tdc.metrics.sortbyabstd(estimates)

I sort your estimates according to the absolute value of their time delay. lowest “precision” (== lowest delays) first -> select from end !

@param estimates: list of Estimate objects @return: sorted list of Estimate objects.

pycs.tdc.metrics.maxPplot(estslist, N, filepath=None)

Give me a list of estimate-lists, I plot P as a function of f

@param estslist: list of lists of Estimate objects @param N: total number of curves (e.g. 56 for tdc0) (defines f=1) @param filepath: if not None, it defines the path where I will save the plot

pycs.tdc.metrics.getP(db, method, median=False)

Compute P and Perr - std(P)/sqrt(len(P)) - for a given method stored in the database

pycs.tdc.metrics.getA(db, method, median=False)

Compute A and Aerr - std(A)/sqrt(len(A)) - for a given method stored in the database

pycs.tdc.metrics.getAmod(db, method, median=False)

Compute Amod and Amoderr - std(Amod)/sqrt(len(Amod)) - for a given method stored in the database

pycs.tdc.metrics.getchi2(db, method, median=False)

Compute chi2 and chi2err - std(chi2)/sqrt(len(chi2)) - for a given method stored in the database

pycs.tdc.metrics.getf(db, method, N)

Compute f for a given method stored in the database N is the total number of curves

pycs.tdc.metrics.Pplotall(db, methods, N, zoomchi2=False, zoomA=False, errorbar=False, median=False)

give me the db and a list of methods [method1,method2,…], I plot P vs f for each method and chi2 vs f, A vs f… for the same arrangement according to P

If as a method_i you give me a tuple (method_i,chi2max), I will remove the values wich chi2 > chi2max from the db

I can do the same exercice with median values instead of mean, just specify it in the kwargs

pycs.tdc.metrics.Pplotallperrung(wholedb, methods, errorbar=False)

Ugly beast ahead…

give me the db and a list of methods [method1,method2,…], I plot P vs f for each method and chi2 vs f, A vs f… for the same arrangement according to P

If as a method_i you give me a tuple (method_i,chi2max), I will remove the values wich chi2 > chi2max from the db

pycs.tdc.metrics.combigauss(subtds, subtderrs, truetds, lensmodelsigma=0.0)

Give me submission delays and error bars, as well as the corresponding true delays, in form of numpy arrays. I compute the mean and sigma of the combined posterior on the fractional time delay distance error.

pycs.tdc.metrics.Pplotcombi(db, methods, N, lensmodelsigma=0.0)

give me the db and a list of methods, I plot combigauss params (center, sigma, zerovalue) vs f f is selected according to bestP

pycs.tdc.optfct module

Here we collect the optimizers specifically designed for the TDC These optimizers are to be used on lcs that have a spline already ! Typically, they are called by the runsim and runobs functions of run.py

pycs.tdc.optfct.spldiff(lcs, verbose=True, magshift=False)

Custom spldiff optimizer for TDC

@param lcs: list of light curves @param magshift: boolean, do you want to allow magnitude shifts in the pycs.spldiff.multiopt.opt_ts function called here

pycs.tdc.optfct.regdiff2(lcs)

Regdiff usign scikit-learn Let’s see if this works without any initial shift.

@param lcs: list of light curves

pycs.tdc.run module

Wrapper stuff to run PyCS on TDC data

We change the philosophy here: we make a clear separation between functions that draw lcs, and functions that analyse them.

The goal is to create “once and for all” a huge amount of simulated curves, and then run the optimizer on a limited number of simulated curves, randomly chosen.

The copied and simulated lcs are stored each into an individual pkl (one per lcs)

pycs.tdc.run.createdir(estimate, path)

pycs.tdc.run.drawcopy(estimate, path, n=1, maxrandomshift=None, datadir='')

Draw n times one single copy of lcs (the ones your estimate is about) into the path directory. NEW: set can be the name of a COSMOGRAIL lens.

@param estimates: list of estimate objects, with id, td and tderr @param path: where your copy will be written @param addmlfct: If you want to add different microlensing to your obslcs @param n: Number of time you will run drawobs (each run produce new copycurves) @param maxrandomshift: Maximum range for the random time shift added to the copycurves

pycs.tdc.run.drawsim(estimate, path, sploptfct, n=1, maxrandomshift=None, datadir='')

Draw n times one single sim curves of lcs (the ones your esimate is about) into the path directory.

@param estimates: list of estimate objects, with id, td and tderr @param path: where your copy will be written @param addmlfct: If you want to add different microlensing to your obslcs @param n: Number of time you will run drawsim (each run produce new simcurves) @param maxrandomshift: Maximum range for the random “true” and “guessed” time shift added to the simcurves

pycs.tdc.run.runcopy(estimate, path, optfct, n=1, clist=None)

Run the optimizer (optfct) on n different copycurves Return the optimised timeshift and magnitude for each copycurves

The n copycurves are chosen randomly in the copydir, unless you give a clist of numbers i.e. if clist = [1,3] runobs will run on c001.pkl and c003.pkl

@param estimate: give me the estimate you want me to run on. (for the id) @param path: where the copycurves are written @param optfct: which optimisation function do you want me to use @param n: on how many copies do you want me to run the optimiser @param clist: ids of specific copycurves you want me to run on.

pycs.tdc.run.runsim(estimate, path, optfct, n=1, slist=None)

Run the optimizer (optfct) on n different simcurves Return the optimised timeshift and magnitude for each copycurves Also return the true delays of each simcurve

The n simcurves are chosen randomly in the simdir, unless you give a slist of numbers i.e. if slist = [1,3] runobs will run on s001.pkl and s003.pkl

@param estimate: give me the estimate you want me to run on. (for the id) @param path: where the simcurves are written @param optfct: which optimisation function do you want me to use @param n: on how many sims do you want me to run the optimiser @param slist: ids of specific simcurves you want me to run on.

pycs.tdc.run.multirun(estimate, path, optfct, ncopy, nsim, clist=None, slist=None)

Wrapper around runsim and runobs Run the optimizer ncopy and nsim times on the copy and sim Return the results in a list of estimates objects

@param estimate: give me the estimate you want me to run on. (for the id) @param path: where the copycurves and simcurves are written @param optfct: which optimisation function do you want me to use @param ncopy: on how many copies do you want me to run the optimiser @param ncopy: on how many copies do you want me to run the optimiser @param slist: ids of specific copycurves you want me to run on. @param slist: ids of specific simcurves you want me to run on.

pycs.tdc.run.viz(estimate, path, datadir)

Look at some light curves. change in place…

pycs.tdc.run.summarize(estimate, path, makefig=False, skipdone=True)

Create an Estimate instance of a pair, that contains the results from the PyCS optimisation.

I will silently skip non-yet-ready pairs. If skipdone, I will silently skip those pairs that have already been summarized.

@param estimate: give me the estimate you want me to run on. (for the id) @param path: where the copycurves and simcurves are written @param makefig: do you want me to create a summary figure @param skipdone: do you want me to skip the pairs that are already summarized

pycs.tdc.run.collect(estimates, path)

Gather the output estimates created by summarize and return them in a single tuple

@param estimates: list of estimates you want me to collect @param path: where the summarized results are written

pycs.tdc.run2 module

pycs.tdc.run3 module

pycs.tdc.splopt module

Here we collect some spline optimizers specifically designed for the TDC

pycs.tdc.splopt.calcknotstep(varios)

Give me some outputs of vario, I try to return a good knotstep, based on the highest vratio between the curves and the sampling.

pycs.tdc.splopt.spl1(lcs, verbose=True)

Custom spline optimizer for TDC Assumes reasonable initial time shift, but no magshift.

pycs.tdc.splopt.spl2(lcs, maxit=7, minchange=1.0, verbose=True)

Custom spline optimizer for TDC Assumes good initial time shift, but does not care about inital magshift. Optimizes any ML that is present.

Parameters:

• maxit – maximum number of iteartions
• minchange – minimum decrease in percent of the r2. I stop if the decrease gets smaller.

pycs.tdc.splopt.spl3(lcs, knotstepfact=1.0, mlknotstep=365, maxit=7, minchange=1.0, verbose=True)

Version 3 Assumes a really good initial time shift, does not care about inital magshift.

ML is added inside this function to a random lc. The ML signal should not get strong if it is not required.

Parameters:

• maxit – maximum number of iteartions
• minchange – minimum decrease in percent of the r2. I stop if the decrease gets smaller.

pycs.tdc.splopt.splml1(lcs)

Some rather simple spline ML model. Randomly added on one curve of lcs

pycs.tdc.splopt.splml2(lcs)

Trying to get this better.

pycs.tdc.stats module

Compute some statistics about TDC1

pycs.tdc.stats.progress(estimates, htmlcode=False, path='leaderboard.html')

Print the overall progress of TDC1 estimations through D3CS

@param estimates: list of Estimate objects. @param htmlcode: boolean. If True, I will write the results in an html file @param path: path of the html output file.

pycs.tdc.util module

General purpose functions related to the TDC.

pycs.tdc.util.listtdc1v2pairs()

A list of all pairs that are available in TDC1 v2 (useful, as some pairs have been removed…)

@return: a list with all the existing pairs

pycs.tdc.util.tdcfilepath(set, rung, pair, skipset=False)

Return the path to a TDC light curve in the form of “set/rung/set_rung_pair.txt”

@param set: “tdc0” or “tdc1” @param rung: in [0, 4] for tdc1 @param pair: in [1, 1036] for tdc1 (some are missing) @param skipset: boolean. If true, skip the first “set” in the path

@return: string. Path “set/rung/set_rung_pair.txt”

pycs.tdc.util.pogmag(flux, fluxerr, m0=22.5)

Computes a “normal” Pogson magnitude from TDC info. Used for TDC0

@param flux: numpy array of floats @param fluxerr: numpy array of floats @param m0: float. Zero-point of the magnitude system. @return: tuple containing two lists of fluxes converted into magnitudes.

pycs.tdc.util.asinhmag(flux, fluxerr, m0=22.5, f0=1.0, b=0.01)

Implements asinh magnitudes, following http://ssg.astro.washington.edu/elsst/opsim.shtml?lightcurve_mags

@param flux: numpy array of floats @param fluxerr: numpy array of floats @param m0: float. Zero-point of the magnitude system. @param f0: float. asinh mag normalisation constant. @param b: float. Asinh mag normalisation constant. @return: tuple containing two lists of fluxes converted into magnitudes

pycs.tdc.util.read(filepath, mag='asinh', verbose=True, shortlabel=True)

Imports TDC light curves as a PyCS Lightcuve object. So far we expect exactly 2 curves in each file, as TDC simulates only doubles.

@param filepath: string. Path to the tdc lightcurves, typically the output of tdcfilepath @param mag: can be “asinh” of “pog”. Defines the type of magnitude considered @param verbose: boolean. Defines verbosity @param shortlabel: boolean. If True, gives a shorter object name to each light curve, omitting its TDC id.

@return: Return two PyCS Lightcurve objects, one per tdc light curve

pycs.tdc.util.writesubmission(estimates, filepath, commentlist=None, theseonly=False)

Writes a list of estimates into a TDC1 submission

@param estimates: list of Estimate objects, for which you want to write the time delay and error estimations. @param filepath: path of the file to be written. Accordign to the doc, filepath should be something like pycs_tdc1_algorithm.dt @param commentlist: list of srings. Additionnal comments to be written in the header of the submission file. @param theseonly: boolean. If set to True, I will not write -99 flags for those estimates that are missing in your list !

pycs.tdc.util.setnicemagshift(lcs)

Sets a “nice” magshift to the n-1 last curves of lcs, so that they appear nicely below each other when displayed. Also works if you curves contain ML models or are alredy shifted.

@param lcs: list of Lightcurve objects

pycs.tdc.util.goingon()

Asks the user if he wants to proceed. If not, exits python.

pycs.tdc.util.readsubmission(filepath)

Read a TDC1 submission file, and return its values in PyCS standards.

@param filepath: string. Path to the TDC1 submission

@return: list of tuples. Each tuple are (estimate id, estimated delay, estimated error)

pycs.tdc.util.godtweak(estimates)

I return a list of estimates in which I have modified some by divine inspiration. This should be systematically done before writing a TDC1-D3CS submission.

@param estimates: list of Estimate objects @return: tweaked list of estimates

pycs.tdc.vario module

Variability analysis stuff

pycs.tdc.vario.vario(l, plot=False, filepath=None, nsamp=1000000, verbose=False)

A simple stochatic variogram-like plot, and return the ratio delta_mag around 50-75 over delta_mag for the smallest step

@param l: light curve Object. @param nsamp: integer. The number of random samples to take @param plot: boolean. Do you want to plot the vario analysis ? @param filepath: if not None, it is the path to where the plot will be saved

Module contents

This subpackage contains code more or less specific to the Time Delay Challenge (tdc)

Module contents

This is PyCS, a python toolbox to estimate time delays between COSMOGRAIL light curves.

License:GPLv3 (see README.md) Website:www.cosmograil.org

Last build of this documentation : Sep 01, 2020.