diff --git a/src/data/timing/epfold.sl b/src/data/timing/epfold.sl index 54c09491c551062044e0774e9effb2aceab05cf9..da47ce80c0df376ece4ecc992429db1eb7f7869e 100644 --- a/src/data/timing/epfold.sl +++ b/src/data/timing/epfold.sl @@ -1,31 +1,209 @@ -%%%%%%%%%%%%%%%%%%%%%%%% +% -*- mode: slang; mode: fold; -*- + +require("stats"); + +private define _epfold_event_avg_var (t, tstart, tstop) { + variable ttotal = sum(tstop - tstart); + variable avg = length(t)*1.0/ttotal; + + variable sd = 0.0; + variable i, w; + _for i (0, length(tstart)-1) + { + w = where(tstart[i] <= t <= tstop[i]); + if (length(w) > 0) + sd += sumsq(floor(interpol([tstart[i],t[w]]+1./avg, [tstart[i],t[w]], + [1:length(w)+1]))-[1:length(w)+1] - avg); + } + + return avg, sd*1./(length(t)+length(tstart))*avg, ttotal; +} + +private define _epfold_lc_avg_var (dt, r, fracexp) { + variable ttotal; + if (length(dt)==1 && length(fracexp)==1) + ttotal = length(r)*dt*fracexp; + else + ttotal = sum(dt*fracexp); + + variable avg = sum(r*dt*fracexp)/ttotal; + + variable var = sumsq(r - avg)/(length(r)-1); + + return avg, var, ttotal; +} + +private define _epfold_create_grid (pmin, pmax, ttotal, verbose) { + variable c, p; + + if (verbose > 0) + vmessage("Determining number of trial periods between %f and %f.", pmin, pmax); + + if (qualifier_exists("lineardp")) { + c = qualifier("lineardp"); + if (verbose > 1) + vmessage("Making linear search grid with delta period = %f.", c); + + p = [pmin:pmax:c]; + } else if (qualifier_exists("linearn")) { + c = qualifier("linearn"); + if (verbose > 1) + vmessage("Making linear search grid with %d periods.", c); + + p = [pmin:pmax:#c]; + } else { + % Not sure what the logic here is exactly, my experience is that + % this does not work so well. JJRS + if (verbose > 1) + vmessage("Making intelligent search grid!"); + + p = {pmin}; + while (p[-1] < pmax) + list_append(p, p[-1] + sqr(p[-1])/(ttotal*qualifier("sampling", 10))); + p = list_to_array(p); + } + + if (verbose > 0) + vmessage("number of periods = %d", length(p)); + + return p; +} + +private define _epfold_eval_stat (fun, avg, var, pp, arg) { + return @fun(pp, avg, var, __push_list(arg);; __qualifiers()); +} + +private define epfold_chisqrstat (pp, avg, var) { + % chi sqr test, see Davies 1990, eq (4) + return sum(pp.npts*(pp.value-avg)^2)/var; +} + +private define epfold_lstat (pp, avg, var) { + % lstat test, see Davies 1990, eq (11), (12) + variable snpts = sum(pp.npts); + variable f = (snpts-length(pp.value))/(length(pp.value)-1); + variable ssa = sum(pp.npts*(pp.value-avg)^2); + variable sse = sum(pp.error^2*_max(pp.npts-1, 0.)); + + return f*ssa/sse; +} + +private define epfold_zstat (pp, avg, var, n) { + % zstat test, see Bachetti 2021, eq (11), appendix b + variable f = 2./sum(pp.value); + variable s = 0.0; + variable phi = PI*(pp.bin_hi + pp.bin_lo); + + variable i; + _for i (1, n) { + s += sum(pp.value*cos(i*phi))^2 + sum(pp.value*sin(i*phi))^2; + } + + return f*s; +} + +private variable _epfold_stats = Assoc_Type[Struct_Type]; +_epfold_stats["chi"] = struct { statfun = & epfold_chisqrstat }; +_epfold_stats["l"] = struct { statfun = & epfold_lstat }; +_epfold_stats["z"] = struct { statfun = & epfold_zstat, statarg = {1} }; + +private define _epfold_run (i, ep, fun, avg, var, pp, arg) { + variable p0, wvalid; + + p0 = ep.p[i]; + + wvalid = wherenot(isnan(pp.value) or isinf(pp.value)); + + if (length(wvalid) != length(pp.value)) { + vmessage("WARNING: period %f has faulty phase bins! Ignoring those phase bins, but please check!", p0); + ep.badp[i] = 1; + struct_filter(pp, wvalid); + } + + ep.stat[i] = _epfold_eval_stat(fun, avg, var, pp, arg;; __qualifiers()); +} + +private define _epfold_search_grid (nbins, pgrid, avg, var, verbose) +{ + variable ep = struct { + p = pgrid, + stat = Double_Type[length(pgrid)], + nbins = nbins, + badp = Int_Type[length(pgrid)], + avg = avg, + var = var, + }; + + if (verbose > 0) { + variable tstart = _time(); + variable tlast = tstart; + variable titer = 0.0; + variable remain = 0.0; + } + + variable i; + _for i (0, length(pgrid)-1) { + _epfold_run(i, ep, qualifier("statfun", &epfold_chisqrstat), + avg, var, + pfold(__push_list(qualifier("runarg")), pgrid[i];; qualifier("runqual")), + qualifier("statarg", {});; + qualifier("statqual", struct { @NULL })); + + if(verbose>0) + { + if (_time()-tlast > 10) + { + tlast = _time(); + titer = (tlast-tstart)/double(i); + remain = (length(pgrid)-i)*titer; + + vmessage("%.2f%% done", i*100./length(pgrid)); + vmessage("its taking %.2fsec per iteration", titer); + + if (remain <= 600) + vmessage("%.2f sec remaining...", remain); + else if (remain <= 7200) + vmessage("%.2f min remaining...", remain/60.); + else + vmessage("%.2f h remaining...", remain/3600.); + } + } + } + + return ep; +} + +%%%%%%%%%%%%%%%% define epfold () -%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%% %!%+ %\function{epfold} %\synopsis{peforms epoch folding on a lightcurve in a given period %interval} -%\usage{(Struct_Type res) = epfold(Double_Type t, r, pstart, pstop); -%\altusage{(Struct_Type res) = epfold(Double_Type t, pstart, pstop) ; (event data)}} +%\usage{(Struct_Type ep) = epfold(Double_Type t, r, pstart, pstop); +%\altusage{(Struct_Type ep) = epfold(Double_Type t, pstart, pstop) ; (event data)}} % %\qualifiers{ %\qualifier{nbins}{number of bins for the pulse profile (default=20)} %\qualifier{exact}{calculate the pulse profile in a more exact way, see description of pfold (not recommended as it takes a very long time!).} %\qualifier{dt}{exposure of every lightcurve time bin, should be given to ensure correct results.} +%\qualifier{fracexp}{Fractional exposure time for light-curve data (default=1.0)} %\qualifier{sampling}{how many periods per peak to use (default=10)} %\qualifier{nsrch}{how many periods to search in a linear grid (default not set)} %\qualifier{dp}{delta period of linear period grid (default not set)} -%\qualifier{lstat}{use L-statistics instead of chi^2 statistics} +%\qualifier{lstat}{use L-statistic measure} +%\qualifier{chistat}{use chisqr-statistic measure} +%\qualifier{zstat}{use Z-statistic measure with given order (default=1)} %\qualifier{chatty}{set the verbosity, chatty-1 is piped to pfold (default=0)} -%\qualifier{gti}{GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}} +%\qualifier{gti}{GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}} %} +%#c%{{{ % % %\description % Performs epoch folding on a given lightcurve between the periods % pstart and pstop. The function was adopted from -% the IDL routine of the same name. GTI correction only implemented -% for event-data yet. +% the IDL routine of the same name. % % By default, periods are sampled according to the triangular rule % for estimating the period error, using "sampling" periods per peak. @@ -33,20 +211,20 @@ define epfold () % between two consecutive periods or "nsrch" for a given number of % periods can be given. These qualifiers are mutually exclusive. % -% The returned structure "res" contains four fields: "p" for the +% The returned structure "ep" contains six fields: "p" for the % evaluated period and "stat" for the value of the statistic used. % Additionally the field "nbins" contains the number of phase bins -% used and "badp" contains indices for res.p marking periods +% used and "badp" contains indices for ep.p marking periods % where the pulse profile showd empty phase bins. Values of -% res.stat[res.badp] should be taken with great care! -% -% Compared to the similar function sitar_epfold_rate, epfold.sl -% uses the chi^2 statistic as a default and is based on pfold.sl, -% which can take errors of the lightcurve into account. -% If the qualifier "lstat" is given, the statistic is switched to -% the l-stat statistic as in sitar_epfold_rate, but errors of the -% lightcurve are no longer taken into account. lstat is not -% available for event-data. +% ep.stat[ep.badp] should be taken with great care! The field "avg" +% is the average rate, and "var" the total variance. +% +% Compared to the similar function sitar_epfold_rate, epfold +% uses the chi^2 statistic as a default and is based on pfold. +% Two other statistics are available: L-statistic (with qualifier +% 'lstat') and Z^2-statistic (with qualifier 'zstat'). The later takes +% an integer value that gives the number of harmonics used to +% calculate the statistic. % % If the "exact" qualifier is given, the function takes the exposure % time of every time bin into account in the sense that, that a @@ -64,175 +242,101 @@ define epfold () % bugs or missing features to Felix % (felix.fuerst@sternwarte.uni-erlangen.de). % -% NOTE: Please read +% NOTE: Please read % Davies, S.R., 1990, MNRAS 244, 93 (L-statistics!) % Larsson, S., 1996, A&AS 117, 197 % Leahy, D.A., Darbro, W., Elsner, R.F., et al.,1983, ApJ 266, 160-170 % Schwarzenberg-Czerny, A., 1989, MNRAS 241, 153 -% -%!%- +% Bachetti, M., Pilia, M., Huppenkothen, D., et al., 2021, ApJ 909, 33 +% +%#c%}}} +%!%- { - variable t,r,pstart, pstop; + variable pmin, pmax, t, r = NULL; + variable verbose = qualifier("verbose", qualifier("chatty", _Isisscripts_Verbose)); - variable eventdata = 0; - % disable chatty if not set - variable chatty = qualifier("chatty",qualifier("verbose",_Isisscripts_Verbose)); - - switch(_NARGS) - { case 3: (t, pstart, pstop) = (); - eventdata = 1 ; - if (chatty>0) { vmessage("using event data"); } + switch (_NARGS) + { case 3: (t, pmin, pmax) = (); + if (verbose) + vmessage("using event data"); } - { case 4: (t,r,pstart, pstop) = (); } - { help(_function_name()); return; } - - variable t0 = qualifier ("t0", t[0]); - variable nbins = qualifier ("nbins", 20); - variable sampl = qualifier("sampling", 10); %periods per peak - variable dp = qualifier("dp", -10) ; %delta p for linear period sampling - variable nsrch = qualifier("nsrch", -10); %number of periods in linear period grid - - nsrch = int(nsrch); - - ifnot(qualifier_exists("dt") or eventdata==1) { - message("WARNING (epfold): no dt qualifier given, creating differential grid! Might lead to unexpected results!") ; + { case 4: (t, r, pmin, pmax) = (); + if (verbose) + vmessage("using light-curve data"); } - variable dt = qualifier("dt" , make_hi_grid(t)-t); %exposure times + { help(_function_name()); } - ifnot (qualifier_exists("gti") or eventdata==0) { - message("WARNING (epfold): using eventdata without gti qualifier, taking time of first and last event. This is most likely NOT what you want!") ; + variable flagc = 0; + variable gridq; + if (qualifier_exists("nsrch")) { + flagc++; + gridq = struct { linearn = qualifier("nsrch") }; + } else if (qualifier_exists("dp")) { + flagc++; + gridq = struct { lineardp = qualifier("dp") }; } - variable gti = qualifier("gti", struct{ start=min(t), stop=max(t) }); + else + gridq = struct { sampling = qualifier("sampling", 10) }; - variable tges = sum(dt) ; %effective observation time + if (flagc > 1) + throw UsageError, "Qualifiers 'nsrch' and 'dp' are mutually exclusive"; - if (tges<0) { - throw UsageError, "Please provide lightcurve with continously increasing time (tges was < 0)."; + flagc = 0; + variable statq = _epfold_stats["chi"]; + if (qualifier_exists("chistat")) + flagc++; + else if (qualifier_exists("lstat")) { + flagc++; + statq = _epfold_stats["l"]; + } else if (qualifier_exists("zstat")) { + flagc++; + statq = _epfold_stats["z"]; + statq.statarg = { qualifier("zstat") ? qualifier("zstat") : 1 }; } - if (chatty>0) { vmessage(">> Starting data analysis <<"); } + if (flagc>1) + throw UsageError, "Statistic qualifiers are mutually exclusive"; - if (chatty>0) { - vmessage("Determining number of trial periods between %f and %f.", pstart, pstop); + variable gti = NULL, dt = NULL; + if (NULL == r and not qualifier_exists("gti")) { + vmessage("***Warning (epfold): using eventdata without gti qualifier, taking time of first and last event. This is most likely NOT what you want!"); + gti = struct { start = min(t), stop = max(t) }; + } else if (NULL != r and not qualifier_exists("dt")) { + vmessage("***Warning (epfold): no dt qualifier given, creating differential grid! Might lead to unexpected results!"); + dt = mean(t[[1:]]-t[[:-2]]); } - if (eventdata and qualifier_exists("lstat")) { - vmessage("WARNING (epfold): using eventdata, lstat qualifier not sensible! Reverting to use of chi^2 statistic."); - } - - variable ergdim=1; - variable p; - variable i = 1; + variable t0 = qualifier("t0", 0.0); + variable nbins = qualifier("nbins", 20); + variable fracexp = qualifier("fracexp", 1.0); + dt = qualifier("dt", dt); + gti = qualifier("gti", gti); - if (dp > 0 && nsrch == -10) - { - if (chatty>0) { vmessage("Making linear period grid with delta P = %f.", dp); } - p = [pstart:pstop:dp]; - ergdim = length(p); - } - else if (dp == -10 && nsrch > 0) - { - if (chatty>0) { vmessage("Making linear period grid with %d entries.", nsrch); } - p = [pstart:pstop:#nsrch]; - ergdim = nsrch; - } - else if (dp > 0 && nsrch > 0) - { - throw UsageError, "Unable to construct search grid, provide either 'dp' or 'nsrch'"; - } - else { - if (chatty>0) { vmessage("Making intelligent grid!"); } - variable p_bin_hi = pstart; - variable jj = 1; - while (p_bin_hi <= pstop) { % This may look pedestrian, but I determine first how much periods the grid will contain, because doing the period calculation twice will still be WAY faster than copying the arrays - p_bin_hi = p_bin_hi + sqr(p_bin_hi)/(tges*sampl); - jj++; - }; - p = Double_Type[jj]; % initialize array long enough for all test periods - p[0] = pstart; - while (p[i-1] <= pstop) - { ergdim++ ; - p[i]= p[i-1]+p[i-1]*p[i-1]/(tges*sampl); - i++ ; - } - p = p[[0:i-1]]; % cut the empty array entries + variable avg, var, ttotal; + variable runarg, runqual; + if (NULL == r) { + (avg, var, ttotal) = _epfold_event_avg_var(t, gti.start, gti.stop); + runarg = {t}; + runqual = struct { nbins = nbins, t0 = t0, gti = gti }; + } else { + (avg, var, ttotal) = _epfold_lc_avg_var (dt, r, fracexp); + runarg = {t, r}; + runqual = struct { nbins = nbins, t0 = t0, dt = dt, fracexp = fracexp }; } - if (chatty>0) { vmessage("number of periods = %d .... done", ergdim); } - - if (chatty>0) { variable tstart = _time(); variable tlast = _time(); } - - variable ptest = pstart; - variable pp, ndx ; - variable parr = Double_Type[ergdim]; - variable chierg = Double_Type[ergdim]; - variable perc, titer, remain; - - variable factor = (1.* (length(t)-nbins)) / (nbins -1.); % (N-M)/(M-1) as in Davies Eq. 11 - variable finitepp, faulty = Integer_Type[0]; - - %subtract mean from rate, as we need only the - %difference for the statistics and it does not change - %the generality - ifnot(eventdata) { - variable rm = moment(r); - r -= rm.ave; - - if (chatty>0) { - vmessage("mean = %f", rm.ave); - vmessage("effective observation time: %f", tges ); - } + if (ttotal < 0) { + throw UsageError, "Total exposure time is negative, usually this is due to light-curve bins not order with increasing time."; } - _for i (0, length(parr)-1, 1) - { - ptest = p[i]; - if (eventdata) { pp = pfold(t,ptest ; nbins=nbins,chatty=chatty, t0=t0,gti=gti); } - else if(qualifier_exists("exact")) { pp = pfold(t,r,ptest ; nbins = nbins, dt = dt, chatty=chatty, t0=t0, exact); } - else { pp = pfold(t,r,ptest ; nbins = nbins, dt = dt, chatty=chatty, t0=t0 ); } + if (qualifier_exists("exact")) + runqual = struct { @runqual, exact }; - finitepp = wherenot(isnan(pp.value) or isinf(pp.value)); + variable pgrid = _epfold_create_grid (pmin, pmax, ttotal, verbose;; gridq); + variable ep = _epfold_search_grid (nbins, pgrid, avg, var, verbose;; + struct { runarg = runarg, runqual = runqual, @statq }); - if (length(finitepp) < nbins) { - vmessage("WARNING: period %f has faulty phase bins! Ignoring those phase bins, but please check!", ptest); - faulty = [faulty, i]; - struct_filter(pp, finitepp); - } - - if (eventdata) { - chierg[i] = sum((pp.value[*]-mean(pp.value))^2. / (mean(pp.value)/pp.ttot[*])); - } else { - if (qualifier_exists("lstat")) - { - % see Davies Eq. 11 - chierg[i] = factor * ( sum(pp.npts*pp.value^2))/(sum (pp.error*(pp.npts-1.))); - } - else - { - % see Davies Eq. 4 or Larsson Eq. 1 - chierg[i] = sum(pp.npts*pp.value^2.)/rm.var; - } - } - - if(chatty>0) - { - if (_time()-tlast > 10) - { - tlast=_time(); - titer = (tlast-tstart)/double(i); - remain=(ergdim-i)*titer; - - vmessage("%.2f%% done", i*100./ergdim ); - vmessage("its taking %.2fsec per iteration", titer); - - if (remain <= 600) vmessage("%.2f sec remaining...", remain); - else if (remain <= 7200) vmessage("%.2f min remaining...", remain/60.); - else vmessage("%.2f h remaining...", remain/3600.); - } - } - } + % convert bad phase bins to index array + ep.badp = where(ep.badp); - %structure containing searched periods, statistics value, - %phasebins used and indices of periods with gaps in the pulse profile - return struct{ p=p, stat=chierg, nbins=nbins, badp = faulty }; + return ep; } diff --git a/src/data/timing/pfold.sl b/src/data/timing/pfold.sl index bd234fc82bec8fe32deead71885301c12d3b5ef3..99b16adb63ff6c1fe8976f52571085eae77107e5 100644 --- a/src/data/timing/pfold.sl +++ b/src/data/timing/pfold.sl @@ -172,7 +172,7 @@ private define periodFoldEvents (time, p) %{{{ % working variables variable phaseLo, phaseHi, pulses, k, j, startj, n, pStart, pStop, nStart, - nStop, inStop, inStart, tTime, gtiStart, gtiStop, + nStop, inStop, inStart, tTime, gtiStart, gtiStop, addN, phase, eTime, bTime, gtiIdx, pulseStart, pulseStop, mapH; result = struct { @@ -181,7 +181,11 @@ private define periodFoldEvents (time, p) %{{{ value, % count rate error = Double_Type[nbins], % count rate error counts = UInt_Type[nbins], % counts - npts, % how often a bin was hit by start time (not sure what this is used for) + npts = Double_Type[nbins], % how much the data contributes to this bin + % As we construct it on the phase grid every + % this is just the number each phase bin was + % covered NOT the number of counts in each bin + % (fractional contribution due to gti windows). ttot = Double_Type[nbins], % per bin exposure }; phaselc = NULL; @@ -254,6 +258,7 @@ private define periodFoldEvents (time, p) %{{{ % gti window is completely inside bin tTime = gtiStop-gtiStart; result.ttot[inStart] += tTime; + result.npts[inStart] += abs((pStop-pStart)/(result.bin_hi[inStart]-result.bin_lo[inStart])); } else if ((nStop - nStart) > 5) { @@ -267,16 +272,21 @@ private define periodFoldEvents (time, p) %{{{ bTime = timeOfPhase([grid[[inStart:nbins-1]]+nStart, grid[[inStart+1:nbins]]+nStart], p, dp, ddp); bTime[0] = gtiStart; result.ttot[[inStart:nbins-1]] += bTime[[length(bTime)/2:length(bTime)-1]]-bTime[[0:length(bTime)/2-1]]; + result.npts[[inStart+1:nbins-1]] += 1.0; + result.npts[inStart] += abs((result.bin_hi[inStart]-(pStart-nStart))/(result.bin_hi[inStart]-result.bin_lo[inStart])); bTime = timeOfPhase([grid[[0:inStop]]+nStop, grid[[1:inStop+1]]+nStop], p, dp, ddp); bTime[-1] = gtiStop; result.ttot[[0:inStop]] += bTime[[length(bTime)/2:length(bTime)-1]]-bTime[[0:length(bTime)/2-1]]; + result.npts[[0:inStop-1]] += 1.0; + result.npts[inStop] += abs((pStop-nStop-result.bin_lo[inStop])/(result.bin_hi[inStop]-result.bin_lo[inStop])); _for j (0, nbins-1) { bTime = timeOfPhase([grid[j]+n, grid[j+1]+n], p, dp, ddp); result.ttot[j] += sum(bTime[[length(bTime)/2:length(bTime)-1]]-bTime[[0:length(bTime)/2-1]]); + result.npts[j] += 1.0; } } @@ -284,19 +294,28 @@ private define periodFoldEvents (time, p) %{{{ { % assume constant p over GTI such that there is a linear relation between dPhi and dT bTime = timeOfPhase([nStart+grid,nStop+grid], p, dp, ddp); tTime = (bTime[nbins+1]-bTime[nbins])/nbins; + result.ttot += tTime; result.ttot[[inStart:nbins-1]] += (bTime[nbins]-gtiStart)/(nbins-inStart); result.ttot[[0:inStop]] += (gtiStop-bTime[nbins+1])/(inStop+1); + + result.npts += (nStop - nStart - 2)*1.0; + result.npts[[inStart:nbins-1]] += 1.0; + result.npts[[0:inStop]] += 1.0; } } else { bTime = timeOfPhase(n+grid, p, dp, ddp); startj = j; + addN = abs((result.bin_hi[startj]-(pStart-nStart))/(result.bin_hi[startj]-result.bin_lo[startj])); do { if (j==nbins) { % next period? % phase grid change, so write it first result.ttot[[startj:nbins-1]] += _min(gtiStop,bTime[[startj+1:nbins]])-_max(gtiStart,bTime[[startj:nbins-1]]); + result.npts[[startj+1:nbins-1]] += 1.0; + result.npts[startj] += addN; + addN = 1.0; % only the first bin is != 1 n++; bTime = timeOfPhase(n+grid, p, dp, ddp); j = 0; startj = j; @@ -304,8 +323,11 @@ private define periodFoldEvents (time, p) %{{{ j++; } while ((j <= inStop) || (n < nStop)); % last partial period not covered yet, j==0 if last bin had been covered - if (j!=0) + if (j!=0) { result.ttot[[startj:j-1]] += _min(gtiStop,bTime[[startj+1:j]])-_max(gtiStart,bTime[[startj:j-1]]); + result.npts[[startj:j-2]] += 1.0; + result.npts[j-1] += abs((pStop-nStop-result.bin_lo[j-1])/(result.bin_hi[j-1]-result.bin_lo[j-1])); + } } if (length(phase)) @@ -314,7 +336,6 @@ private define periodFoldEvents (time, p) %{{{ result.value = result.counts*1./result.ttot; % count rate result.error = sqrt(result.counts)/result.ttot; % count rate error - result.npts = result.counts; % not sure what this is used for if ( NULL != phaselc ) phaselc = struct { phaselc = phaselc }; @@ -349,7 +370,9 @@ private define periodFoldLCfast (input) %{{{ value, % count rate error = Double_Type[length(fastIdx)], % count rate error counts = Double_Type[length(fastIdx)], % counts - npts = npts, % how often a bin was hit by start time (not sure what this is used for) + npts = 1.0*npts, % how much the data contributes to this bin + % For the fast algorithm this is simply the number hits + % in a phase bin. ttot = Double_Type[length(fastIdx)], % per bin exposure }; phaselc = NULL; @@ -359,8 +382,8 @@ private define periodFoldLCfast (input) %{{{ cDt = (length(dt) == 1) ? dt : dt[fastIdx[j]]; cFracExp = (length(fracexp) == 1) ? fracexp*cDt : fracexp[fastIdx[j]]*cDt; result.counts[j] = sum((rate[fastIdx[j]])*cFracExp); - result.ttot[j] = sum(cFracExp); - result.error[j] = (NULL != rateerr) ? sum(sqr((rateerr[fastIdx[j]])*cFracExp)) : result.counts[j]; + result.ttot[j] = (length(cFracExp) == 1) ? npts[j]*cFracExp : sum(cFracExp); + result.error[j] = (NULL != rateerr) ? sumsq((rateerr[fastIdx[j]])*cFracExp) : result.counts[j]; } if (qualifier_exists("phaselc")) @@ -424,7 +447,7 @@ private define periodFoldLCexact (input) %{{{ variable npts, tStart, tStop, pStart, pStop, nStart, nStop, inStart, inStop, cFracExp, n, j, k, tTime, tVal, result, - bTime, startj, phaselc, pulses, pulseStart, pulseStop; + bTime, startj, phaselc, pulses, pulseStart, pulseStop, fracnpts; result = struct { bin_lo = grid[[:-2]], @@ -432,7 +455,9 @@ private define periodFoldLCexact (input) %{{{ value, % count rate error = Double_Type[nbins], % count rate error counts = Double_Type[nbins], % counts - npts = Int_Type[nbins], % how often a bin was hit by start time (not sure what this is used for) + npts = Double_Type[nbins], % how much the data contributes to this bin + % For the fast algorithm this is simply the number hits + % in a phase bin. ttot = Double_Type[nbins], % per bin exposure }; phaselc = NULL; @@ -490,6 +515,7 @@ private define periodFoldLCexact (input) %{{{ result.counts[j] += tVal; result.ttot[j] += tTime; result.error[j] += (NULL != rateerr) ? sqr(rateerr[k]*tTime) : tVal; + result.npts[j] += 1.0; if (NULL != phaselc) { @@ -504,6 +530,10 @@ private define periodFoldLCexact (input) %{{{ result.counts += rate[k]*tTime; result.ttot += tTime; + % average over all bins + fracnpts = 1./((nStop - nStart - 2)*nbins + inStop + 1 + (nbins - inStart)); + result.npts += (nStop-nStart)*fracnpts; + pulseStop = pulseStart + nStop-nStart; if (NULL != phaselc) { @@ -514,6 +544,8 @@ private define periodFoldLCexact (input) %{{{ tTime = (bTime[nbins]-tStart)/(nbins-inStart)*cFracExp; result.ttot[[inStart:nbins-1]] += tTime; result.counts[[inStart:nbins-1]] += rate[k]*tTime; + % we overcounted those + result.npts[[0:inStart-1]] -= fracnpts; if (NULL != phaselc) { @@ -524,6 +556,8 @@ private define periodFoldLCexact (input) %{{{ tTime = (tStop-bTime[nbins+1])/(inStop+1)*cFracExp; result.ttot[[0:inStop]] += tTime; result.counts[[0:inStop]] += rate[k]*tTime; + % we overcounted those + result.npts[[inStop+1:nbins-1]] -= fracnpts; if (NULL != phaselc) { @@ -536,6 +570,7 @@ private define periodFoldLCexact (input) %{{{ % otherwise we loop through the phases until we hit the interval limit bTime = timeOfPhase(n+grid, input.P[0], input.P[1], input.P[2]); startj = j; + fracnpts = 1./((nStop - nStart - 2)*nbins + inStop + 1 + (nbins - inStart)); do { if (j == nbins) { @@ -544,6 +579,7 @@ private define periodFoldLCexact (input) %{{{ result.ttot[[startj:nbins-1]] += tTime; result.counts[[startj:nbins-1]] += tVal; result.error[[startj:nbins-1]] += (NULL != rateerr) ? rateerr[k]*tTime : tVal; + result.npts[[startj:nbins-1]] += fracnpts; if (NULL != phaselc) { @@ -566,6 +602,7 @@ private define periodFoldLCexact (input) %{{{ result.ttot[[startj:j-1]] += tTime; result.counts[[startj:j-1]] += tVal; result.error[[startj:j-1]] += (NULL != rateerr) ? sqr(rateerr[k]*tTime) : tVal; + result.npts[[startj:j-1]] += fracnpts; if (NULL != phaselc) { diff --git a/src/fitting/statistics/residual_runs.sl b/src/fitting/statistics/residual_runs.sl index b3a8e109bd5f2affaa086718bfbaf36c34b4d73f..9b6bcbca4392c396e8ed7e5a3753533d362ee54f 100644 --- a/src/fitting/statistics/residual_runs.sl +++ b/src/fitting/statistics/residual_runs.sl @@ -111,8 +111,8 @@ define residual_runs () foreach id (active) { noti = get_data_info(id).notice_list; - model = get_model(id).value; - data = get_data(id).value; + model = get_model(id).value[noti]; + data = get_data(id).value[noti]; s = data/model-1; w = where(s != 0); % ignore exact matches