cw-metric-template-banks (OctApps Reference Manual)

2.7 `cw-metric-template-banks`

Function File: AnalyseSegmentList ( segment_list_2col )

Function File: AnalyseSegmentList ( segment_list_4col, Ndet, Tsft=1800 )

Function File: props = AnalyseSegmentList ( … )

Analyse a segment list and print/return various properties.

Arguments

segment_list_2col: 2-column segment list: [start_times, end_times]
segment_list_4col: 4-column segment list: [start_times, end_times, unused, num_SFTs]
Ndet: number of detectors
Tspan: SFT time-span in seconds

Examples

props = AnalyseSegmentList([0:2:10; 1:2:11]' * 86400);
sert(props.num_segments, 6)
sert(props.coh_mean_Tspan, 86400)
sert(props.inc_Tobs, 6 * 86400)
sert(props.inc_Tspan, 11 * 86400)
sert(props.inc_duty, 6 / 11, 1e-3)

Function File: iseq = CompareSegmentLists ( seglist, seglist_1, seglist_2, … )

Function File: iseq = CompareSegmentLists ( seglist, {seglist_1, seglist_2, ...} )

Return true if the segment list ’seglist’ equals any of the segment lists seglist_1, seglist_2, etc., and false otherwise.

Examples

seglists = CreateSegmentList(123, [2, 3, 5, 9, 17], [1.2, 4.5, 9.7], [], [0.7, 0.9]);
for n = 1:numel(seglists)
  assert(find(CompareSegmentLists(seglists{n}, seglists)) == n);
endfor

Function File: [ metric, coordIDs ] = ComputeDopplerMetric ( opt, val, … )

Create various Doppler phase- and Fstat metrics.

Arguments

metric

Doppler metric, with struct-elements (depending on metric_type):

metric.g_ij: phase metric
metric.gF_ij: full F-stat metric
metric.gFav_ij: averaged F-stat metric

coordIDs

coordinate IDs of the chosen coordinates

Options

coords

comma-separated list of coordinates:

alpha_delta: physical sky coordinates
ssky_equ: super-sky in equatorial coordinates
ssky_ecl: super-sky in ecliptic coordinates
spin_equ: spin sky in x-y equatorial coordinates
orbit_ecl: orbit sky in x-y-z ecliptic coordinates
freq: frequency in SI units
fdots: frequency spindowns in SI units
gct_nu: GCT frequency/spindowns in SI units
gct_nx_ny_equ: GCT constrained equatorial sky coordinates

spindowns

number of spindown coordinates: 0=none, 1=1st spindown, 2=1st+2nd spindown, etc.

segment_list

list of segments [start_time, end_time; start_time, end_time; …] in GPS seconds

ref_time

reference time in GPS seconds [default: mean of segment list start/end times]

detectors

comma-separated list of detector names

ephemerides

Earth/Sun ephemerides from loadEphemerides()

fiducial_freq

fiducial frequency for sky-position coordinates

det_motion

which detector motion to use (default: spin+orbit)

alpha

for physical sky coordinates, right ascension to compute metric at

delta

for physical sky coordinates, declination to compute metric at

metric_type

compute either phase metric (phase), F-stat metric (Fstat), or both (all)

cosi

cos(iota) signal parameter for full F-stat metric (gF_ij)

psi

polarization angle signal parameter for full F-stat metric (gF_ij)

Examples

See the tutorial on ComputeDopplerMetric().

Function File: metrics = ComputeSuperskyMetrics ( opt, val, … )

Create the supersky parameter-space metrics.

Arguments

metrics

a struct containing the following fields:

coh_rssky_metric: coherent reduced supersky metric for each segment
coh_rssky_transf: coherent reduced supersky metric coordinate transform data for each segment
semi_rssky_metric: semicoherent reduced supersky metric
semi_rssky_transf: semicoherent reduced supersky metric coordinate transform data

Options

spindowns: number of spindown coordinates: 0=none, 1=1st spindown, 2=1st+2nd spindown, etc.
segment_list: list of segments [start_time, end_time; start_time, end_time; ...] in GPS seconds
ref_time: reference time in GPS seconds [default: mean of segment list start/end times]
fiducial_freq: fiducial frequency for sky-position coordinates [required]
detectors: comma-separated list of detector names [required]
detector_weights: vector of weights used to combine single-detector metrics [default: unit weights]
detector_motion: which detector motion to use [default: spin+orbit]
ephemerides: Earth/Sun ephemerides [default: load from loadEphemerides()]
cache_size: size of cache of previously-computed metrics to use [default: 100; 0=no cache]

Examples

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
metrics = ComputeSuperskyMetrics( ...
                                 "spindowns", 1, ...
                                 "segment_list", CreateSegmentList(1e9, 5, 86400, [], 0.75), ...
                                 "fiducial_freq", 123.4, ...
                                 "detectors", "H1,L1", ...
                                 "cache_size", 0 ...
                               );
metrics = ComputeSuperskyMetrics( ...
                                 "spindowns", 1, ...
                                 "segment_list", CreateSegmentList(1e9, 5, 86400, [], 0.75), ...
                                 "fiducial_freq", 123.4, ...
                                 "detectors", "H1,L1", ...
                                 "cache_size", 100 ...
                               );

Function File: segment_lists = CreateSegmentList ( mean_time, Nseg, Tseg, Tspan, duty, [ always_cell=false ] )

Create segment list(s) with the given mean timestamp mean_time, and exactly 3 of the following 4 input parameters:

number of segments Nseg
segment time-span Tseg
total time span Tspan
incoherent duty cycle duty

Aside from mean_time, each parameter may be vectorised to return multiple segment lists with all possible parameter combinations.

Examples

assert( iscell(CreateSegmentList(123, 10, [3, 4], [], 1.0)) )

assert( numel(CreateSegmentList(123, 10, [3, 4], [], 1.0)) == 2 )

assert( !iscell(CreateSegmentList(123, 10, 3, [], 1.0)) )

assert( iscell(CreateSegmentList(123, 10, 3, [], 1.0, true)) )

assert( AnalyseSegmentList(CreateSegmentList(456, [], 7.3, 140.9, 0.87)).mean_time, 456, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(456, [], 7.3, 140.9, 0.87)).num_segments, 16, 1e-5 )

assert( mean(AnalyseSegmentList(CreateSegmentList(456, [], 7.3, 140.9, 0.87)).coh_Tspan), 7.375, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(456, [], 7.3, 140.9, 0.87)).inc_Tspan, 140, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(456, [], 7.3, 140.9, 0.87)).inc_duty, 0.84286, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(345, 50, [], 256.7, 0.77)).mean_time, 345, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(345, 50, [], 256.7, 0.77)).num_segments, 50, 1e-5 )

assert( mean(AnalyseSegmentList(CreateSegmentList(345, 50, [], 256.7, 0.77)).coh_Tspan), 4, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(345, 50, [], 256.7, 0.77)).inc_Tspan, 256, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(345, 50, [], 256.7, 0.77)).inc_duty, 0.78125, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(234, 200, 3.7, [], 0.77)).mean_time, 234, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(234, 200, 3.7, [], 0.77)).num_segments, 200, 1e-5 )

assert( mean(AnalyseSegmentList(CreateSegmentList(234, 200, 3.7, [], 0.77)).coh_Tspan), 3.7, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(234, 200, 3.7, [], 0.77)).inc_Tspan, 962, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(234, 200, 3.7, [], 0.77)).inc_duty, 0.76923, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(123, 100, 2.3, 365.25, [])).mean_time, 123, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(123, 100, 2.3, 365.25, [])).num_segments, 100, 1e-5 )

assert( mean(AnalyseSegmentList(CreateSegmentList(123, 100, 2.3, 365.25, [])).coh_Tspan), 2.34, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(123, 100, 2.3, 365.25, [])).inc_Tspan, 366, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(123, 100, 2.3, 365.25, [])).inc_duty, 0.63934, 1e-5 )

assert( AnalyseSegmentList(CreateSegmentList(123, 1, 2.3, 365.25, [])).inc_duty, 1.0, 1e-5 )

Function File: [ metric, dnorm, invdnorm ] = DiagonalNormaliseMetric ( metric, [ tolerant ] )

Diagonally normalise a metric.

Arguments

metric: metric
dnorm: diagonal normalisation coefficients
invdnorm: inverse of diagonal normalisation coefficients

If tolerant option is given, tolerate zero or negative diagonal elements.

Examples

g = [7,3,5; 3,6,2; 5,2,5];
sert(DiagonalNormaliseMetric(g), g ./ sqrt(diag(g) * diag(g)'), 1e-3)

Function File: mtwoF = EmpiricalFstatMismatch ( Tcoh, Tsemi, mcoh, msemi )

Function File: [ mtwoF, stwoF ] = EmpiricalFstatMismatch ( Tcoh, Tsemi, mcoh, msemi, scoh, ssemi )

Compute an empirical fit, derived from numerical simulations, to the F-statistic mismatch as a function of the coherent and semicoherent time-spans and metric mismatches.

Arguments

mtwoF: mean F-statistic mismatch
stwoF: standard deviation of F-statistic mismatch
Tcoh: coherent segment time-span, in seconds
Tsemi: semicoherent search time-span, in seconds
mcoh: mean coherent metric mismatch
msemi: mean semicoherent metric mismatch
scoh: standard deviation of coherent metric mismatch
ssemi: standard deviation of semicoherent metric mismatch

Input variables may also be vectors. Note that a non-interpolating search implies mcoh = 0.

Examples

assert(EmpiricalFstatMismatch(0, 0, 0, 0) == 0)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0, 0) == 0)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0, 0.5) > 0)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0, 0.5) < 1)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0.1, 0.5) > 0)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0.1, 0.5) < 1)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0.1, 50) > 0)

assert(EmpiricalFstatMismatch(86400, 10*86400, 0.1, 50) < 1)

assert(EmpiricalFstatMismatch(86400, 10*86400, 50, inf) == 1)

assert(EmpiricalFstatMismatch(86400, 10*86400, inf, inf) == 1)

Function File: closest = FindClosestTemplate ( points, metric, max_mismatch, lattice )

Given a set of points, find the closest template to each in a ’virtual’ lattice template bank, constructed using a given metric and maximum mismatch.

Arguments

points: set of points (in columns) to find closest template
metric: parameter-space metric associated with points’
max_mismatch: maximum mismatch of lattice template bank
lattice: type of lattice to use; see LatticeFindClosestPoint()
closest: closest template to each point

Examples

x = octforge_unifrnd(-10, 10, [3,100]);
metric = [7,3,5; 3,6,2; 5,2,5];
max_mismatch = 0.4;
dx = x - FindClosestTemplate(x, metric, max_mismatch, "Zn");
assert(dot(dx, metric * dx) <= max_mismatch);

Function File: g = GCTCoherentTaylorMetric ( opt, val, … )

Computes the coherent global correlation transform metric, as given in Pletsch, PRD 82 042002 (2010), which Taylor-expands the Earth’s orbital motion.

Arguments

g: coherent GCT metric using Taylor-expanded phase model

Options

smax: number of spindowns (up to second spindown)
tj: value of tj, the mid-point of the coherent time span
t0: value of t0, an overall reference time
T: value of T, the coherent time span
Omega: value of Omega, the Earth’s angular rotation frequency (default: 2*pi / (sidereal day in seconds)

Examples

[ ...
  0.33333333333333,-0.66666666666666,1.2,0,-0.31830988618379; ...
  -0.66666666666666,1.422222222222222,-2.666666666666666,0.20264236728467,0.63661977236758; ...
  1.2,-2.666666666666666,5.142857142857143,-0.60792710185402,-1.079730338135965; ...
  0,0.20264236728467,-0.60792710185402,0.5,0; ...
  -0.31830988618379,0.63661977236758,-1.079730338135965,0,0.5; ...
]);

[ ...
  0.33333333333333,0,0.2,0,0.31830988618379; ...
  0,0.088888888888888,0,-0.20264236728467,0; ...
  0.2,0,0.14285714285714,0,0.12480067958459; ...
  0,-0.20264236728467,0,0.5,0; ...
  0.31830988618379,0,0.12480067958459,0,0.5; ...
]);

[ ...
  0.33333333333333,-0.088888888888888,0.21777777777777,-0.014039475036123,-0.0081056946913869; ...
  -0.088888888888888,0.11259259259259,-0.093629629629629,-0.036633359944481,0.072096915013279; ...
  0.21777777777777,-0.093629629629629,0.17936084656084,-0.025350712749409,-0.051935118925438; ...
  -0.014039475036123,-0.036633359944481,-0.025350712749409,0.4959471526543,0.0070197375180618; ...
  -0.0081056946913869,0.072096915013279,-0.051935118925438,0.0070197375180618,0.48784145796291; ...
]);

T = 1:10;
GCT_det = arrayfun(@(T) det(GCTCoherentTaylorMetric("smax", 1, "tj", 0, "t0", 0, "T", T, "Omega", 2*pi)), T);
GCT_det_ref = arrayfun(@(T) detg1(2*pi*T/2), T);
assert(all(abs(GCT_det - GCT_det_ref) < 1e-10));

T = 1:10;
GCT_det = arrayfun(@(T) det(GCTCoherentTaylorMetric("smax", 2, "tj", 0, "t0", 0, "T", T, "Omega", 2*pi)), T);
GCT_det_ref = arrayfun(@(T) detg2(2*pi*T/2), T);
assert(all(abs(GCT_det - GCT_det_ref) < 1e-10));

Function File: gctco = GCTCoherentMetric ( togct, phyco, opt, val, … )

Function File: phyco = GCTCoherentMetric ( tophy, gctco, opt, val, … )

Computes the coherent global correlation coordinates, as given in Pletsch, PRD 82 042002 (2010)

Arguments

phyco: matrix of physical coordinates: [alpha; delta; freq; f1dot; ...]
gctco: matrix of GCT coordinates; order matches that of GCTCoherentMetric()

Options

t0: value of t0, an overall reference time
T: value of T, the coherent time span
detector: detector name, e.g. H1
fmax: maximum frequency to assume when converting sky coordinates
sgndelta: sign of declination when converting to physical coordinates [default: 1]
ephemerides: Earth/Sun ephemerides from loadEphemerides()
ptolemaic: use Ptolemaic orbital motion

Examples

alpha_ref = [1.1604229, 3.4494408, 0.3384569, 3.2528905, 4.2220307, 4.0216501, 4.4586725, 2.6643348, 4.5368049, 4.1968573, 1.0090024, 2.3518049, 6.1686511, 2.7418963, 5.7136839, 4.8999429, 3.5196677, 3.5462285];
delta_ref = [0.1932254, 0.1094716,-0.9947154, 0.0326421, 0.8608284,-0.5364847,-0.7566707, 0.3503306,-0.5262539,-0.7215617,-0.2617862,-0.4193604,-1.0175479, 0.4204785,-0.2233229,-0.2677687,-0.1032384, 1.4105888];
fndot_ref = [6.2948e-3, 1.3900e-3, 7.4630e-3, 3.4649e-3, 8.1318e-3, 4.6111e-3, 9.9765e-3, 2.4069e-3, 8.0775e-3, 6.0784e-3, 1.2569e-3, 8.3437e-3, 8.1557e-3, 4.2130e-3, 7.0748e-4, 4.0846e-3, 3.1204e-3, 3.8858e-3;
             0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000,-1.859e-10,-1.7734e-9,-6.7179e-9,-7.6231e-9,-5.9919e-9,-4.7396e-9,-8.0516e-9,-4.0123e-9,-9.1149e-9,-5.3355e-9,-2.0422e-9,-7.7313e-9;
             0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 2.596e-17, 9.455e-17, 1.913e-17, 6.869e-17, 2.286e-17, 1.743e-17];
gctco_ref = [1.7085e+3, 3.7728e+2, 2.0258e+3, 9.4045e+2, 2.2072e+3, 1.2517e+3, 2.7081e+3, 6.5316e+2, 2.1921e+3, 1.6492e+3, 3.4180e+2, 2.2644e+3, 2.2141e+3, 1.1432e+3, 1.9251e+2, 1.1085e+3, 8.4669e+2, 1.0548e+3;
             5.9246e-4,-1.4654e-4, 2.8904e-4,-3.4929e-4,-3.5437e-4,-4.7675e-4,-1.0908e+0,-1.0396e+1,-3.9393e+1,-4.4699e+1,-3.5128e+1,-2.7784e+1,-4.7212e+1,-2.3520e+1,-5.3450e+1,-3.1288e+1,-1.1973e+1,-4.5325e+1;
             1.1578e-6, 1.3190e-7,-1.0014e-6, 4.9297e-7,-7.7033e-8,-6.1422e-7,-1.6373e-6, 4.2624e-6, 2.0069e-5, 2.8420e-5,-2.3158e-5, 8.9556e-6, 2.1817e-3, 7.9903e-3, 1.5984e-3, 5.8097e-3, 1.9400e-3, 1.4695e-3;
            -1.1593e-3, 2.2515e-3,-1.2265e-3, 2.3005e-3, 8.7039e-4, 1.4456e-3, 6.2425e-4, 1.7776e-3, 5.9593e-4, 1.0377e-3,-1.4158e-3, 1.2768e-3,-1.1741e-3, 1.8152e-3,-1.7200e-3,-1.3115e-4, 2.2189e-3, 3.5332e-4;
            -1.9392e-3, 4.0833e-4,-2.6143e-4,-3.9448e-5, 1.2223e-3, 1.3511e-3, 1.5531e-3,-1.2310e-3, 1.8993e-3, 1.3822e-3,-1.7147e-3,-1.6706e-3, 2.9104e-4,-1.0591e-3, 1.4427e-3, 2.2162e-3, 5.6570e-4, 1.0014e-4];

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
gctco = GCTCoordinates("togct", [alpha_ref; delta_ref; fndot_ref], "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "ptolemaic", false);
assert(all(abs(gctco - gctco_ref) < 1e-4 * abs(gctco_ref)));

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
gctco = GCTCoordinates("togct", [alpha_ref(1:16); delta_ref(1:16); fndot_ref(1:2, 1:16)], "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "ptolemaic", false);
assert(all(abs(gctco - gctco_ref([1,2,4,5], 1:16)) < 1e-4 * abs(gctco_ref([1,2,4,5], 1:16))));

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
gctco = GCTCoordinates("togct", [alpha_ref(1:8); delta_ref(1:8); fndot_ref(1:1, 1:8)], "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "ptolemaic", false);
assert(all(abs(gctco - gctco_ref([1,4,5], 1:8)) < 1e-3 * abs(gctco_ref([1,4,5], 1:8))));

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
phyco = GCTCoordinates("tophy", gctco_ref, "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "sgndelta", delta_ref, "ptolemaic", false);
assert(all(abs(phyco(1, :) - alpha_ref) < 1e-4 * abs(alpha_ref)));
assert(all(abs(phyco(2, :) - delta_ref) < 2e-2 * abs(delta_ref)));
for i = 1:3
  assert(all(abs(phyco(2+i, :) - fndot_ref(i, :)) < max(1e-10, 1e-4 * abs(fndot_ref(i, :)))));
endfor

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
phyco = GCTCoordinates("tophy", gctco_ref(:, 1:16), "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "sgndelta", delta_ref(1:16), "ptolemaic", false);
assert(all(abs(phyco(1, 1:16) - alpha_ref(1:16)) < 1e-4 * abs(alpha_ref(1:16))));
assert(all(abs(phyco(2, 1:16) - delta_ref(1:16)) < 2e-2 * abs(delta_ref(1:16))));
for i = 1:2
  assert(all(abs(phyco(2+i, 1:16) - fndot_ref(i, 1:16)) < max(1e-10, 1e-4 * abs(fndot_ref(i, 1:16)))));
endfor

try
  lal; lalpulsar;
catch
  disp("skipping test: LALSuite bindings not available"); return;
end_try_catch
phyco = GCTCoordinates("tophy", gctco_ref(:, 1:8), "t0", 987654321, "T", 86400, "detector", "L1", "fmax", 0.02, "sgndelta", delta_ref(1:8), "ptolemaic", false);
assert(all(abs(phyco(1, 1:8) - alpha_ref(1:8)) < 1e-4 * abs(alpha_ref(1:8))));
assert(all(abs(phyco(2, 1:8) - delta_ref(1:8)) < 2e-2 * abs(delta_ref(1:8))));
for i = 1:1
  assert(all(abs(phyco(2+i, 1:8) - fndot_ref(i, 1:8)) < max(1e-10, 1e-4 * abs(fndot_ref(i, 1:8)))));
endfor

Function File: g = GCTSemicoherentTaylorMetric ( opt, val, … )

Computes the semicoherent global correlation transform metric, as given in Pletsch, PRD 82 042002 (2010), which Taylor-expands the Earth’s orbital motion.

Arguments

g: semicoherent GCT metric using Taylor-expanded phase model

Options

smax: number of spindowns (up to second spindown)
tj_list: list of values of tj, the mid-point of the coherent time span of each segment
t0: value of t0, an overall reference time
T: value of T, the coherent time span
Omega: value of Omega, the Earth’s angular rotation frequency [default: 2*pi / (sidereal day in seconds)]

Examples

[ ...
  0.33333333333333,0,8.2,0,0.31830988618379; ...
  0,10.75555555555555,0,-0.20264236728467,0.0; ...
  8.2,0,342.5428571428571,0.0,7.76423794799557; ...
  0,-0.20264236728467,0.0,0.5,0; ...
  0.31830988618379,0.0,7.76423794799557,0,0.5; ...
]);

[ ...
  0.33333333333333,0,32.2,0,0.31830988618379; ...
  0,42.75555555555555,0,-0.20264236728467,0.0; ...
  32.2,0,5286.542857142857,0.0,30.6825497532285; ...
  0,-0.20264236728467,0.0,0.5,0; ...
  0.31830988618379,0.0,30.6825497532285,0,0.5; ...
]);

[ ...
  0.33333333333333,0.53333333333334,32.32800000000034,8.881784197001253e-17,0.016211389382774; ...
  0.53333333333334,42.92622222222272,80.2005333333349,0.080754439908228,0.025938223012437; ...
  32.32800000000034,80.2005333333349,4626.441737142834,0.19381065577974,1.609575867543293; ...
  8.881784197001253e-17,0.080754439908228,0.19381065577974,0.48378861061722,8.639163721880971e-18; ...
  0.016211389382774,0.025938223012437,1.609575867543293,8.639163721880971e-18,0.5; ...
]);

T = 5;
N = 10;
t_j = ((1:N) - (N+1)/2)*T;
for s = 1:2
  GCT_semi = GCTSemicoherentTaylorMetric("smax", 1, "tj_list", t_j, "t0", mean(t_j), "T", T, "Omega", 2*pi);
  GCT_coh = GCTCoherentTaylorMetric("smax", 1, "tj", mean(t_j), "t0", mean(t_j), "T", T, "Omega", 2*pi);
  GCT_refine = sqrt(det(GCT_semi) / det(GCT_coh));
  GCT_refine_ref = refinement(1, N);
  assert(abs(GCT_refine - GCT_refine_ref) < 0.005*GCT_refine_ref);
endfor

Function File: bound_box = metricBoundingBox ( metric, max_mismatch )

Compute the metric ellipse bounding box, given a metric and max_mismatch.

Examples

assert(metricBoundingBox([6,3;3,2], 0.77), [1.4329;2.4819], 1e-3)

Function File: [ xcross, mucross ] = metricEllipseCrossSections ( x, mu, metric, cross )

Return the ellipse centers and mismatches for a template bank cross section.

Arguments

xcross: Centers of metric ellipses in the cross section.
mucross: Mismatches of metric ellipses in the cross section.
x: Template bank to find cross section of.
mu: Maximum mismatch of templates.
metric: Parameter-space metric.
cross: NAs indicate dimensions to cross section (only 2 allowed), otherwise give values of cross section in other dimensions

Examples

[xcross,mucross] = metricEllipseCrossSections([1;2;3], 0.5, [6,4,4;4,4,2;4,2,1], [NA;2;NA]);
assert(xcross, [1;3], 1e-3);
assert(mucross, 0.5, 1e-3);

Function File: function [ xx, yy, zz ] = metricEllipsoid ( gij, mismatch, Nsteps=20, method=1 )

return a metric iso-mismatch ellipsoid for given metric gij and mismatch, using ’Nsteps points per surface direction. The function can generate both 2D or 3D ellipses, depending on the input dimension of gij. using either method=1: Cholesky, or method=2: Eigenvectors

The output in the 3D case is a meshgrid, which can be plotted using mesh(), surface(), ...

The output in the 2D case lies in the x-y plane, zz is returned as zeros

Note

this function agrees with and supercedes the older calcMetric2DEllipse() and plotMetricEllipse()

Examples

[xx, yy, zz] = metricEllipsoid([4,1;1,1], 0.25, 4);
assert(xx, [0.25, -0.25, 0, 0.25], 1e-3);
assert(yy, [0.0, 0.5, -0.5, 0.0], 1e-3);
assert(zz, [0.0, 0.0, 0.0, 0.0], 1e-3);

Function File: N = NumberOfLatticeBankTemplates ( opt, val, … )

Estimate the number of templates in a lattice template bank.

Arguments

N: Number of templates

Options

lattice: Type of lattice; see LatticeNormalizedThickness()
metric: Parameter-space metric
max_mismatch: Maximum metric mismatch
param_vol: Parameter space volume

Examples

assert(round(NumberOfLatticeBankTemplates("lattice", "Ans", "metric", [5,1;1,1], "max_mismatch", 0.2, "param_vol", 100)), 385, 1e-3)

Function File: gOut_ij = projectMetric ( g_ij, c = 1 )

project out metric dimension c from the input n x n metric g_ij, by projecting onto the subspace orthogonal to the coordinate-axis of c, namely gOut_ij = g_ij - ( g_ic * g_jc / g_cc )

Returns a ’projected’ n x n metric, where the projected-out dimension c is replaced by zeros, consistent with the behavior of XLALProjectMetric()

Examples

assert(projectMetric(hadamard(4)), -[0,0,0,0;0,2,0,2;0,0,2,2;0,2,2,0], 1e-3)

Function File: gProj_ss = projectMetric2Subspace ( g_ij, sSpace )

function to compute the ’projection’ of a given symmetric metric g_ij’ onto the s-subspace defined by coordinates sSpace = [s1,s2,...], which is a vector of coordinate-indices

If we denotes the indices of the projection-subspace as s and the remaining coordinates as k, such that for an n-dimensional space, writing vn = {1,2,...n}, we have k = vn \ sSpace, then the projection operation is defined as

gProj_{ss’} = g_{ss’} - gkInv^{kk’} g_{ks} g_{k’s’}

where gkInv is the inverse matrix of g_{kk’}, ie the inverse of g restricted to the k-subspace excluding coordinates ’s’

Notes

this function generalizes projectMetric(g_ij,k1) to projecting out more than one dimension k = [k1,k2,...], contrary to projectMetric(), we only return the projected metric of the subspace, ie gProj_{ss’}
gProj_{ss’} gives the smallest possible distance mMin = gProj_{ss’} dl^{s} dl^{s’} of the general distance function m = g_{ij} dl^{i} dl^{j} if we fix the offset in the s-subspace, if for fixed dl^{s}, if one can freely vary the remaining coordinate offsets dl^{k}
gProj_{ss’} is just the Schur complement of g_{kk’}, see https://en.wikipedia.org/wiki/Schur_complement

Examples

rand("seed", 2); tmp = rand(4,4); gij = (tmp + tmp'); ## 4x4 pos-definite matrix, det~3.6, cond~9
tol = 10*eps;

sSpace1 = [2,3,4]; g1_old = projectMetric ( gij, 1 )(sSpace1,sSpace1); g1_new = projectMetric2Subspace(gij, sSpace1);
assert ( g1_old, g1_new, tol );

sSpace2 = [1,3,4]; g2_old = projectMetric ( gij, 2 )(sSpace2,sSpace2); g2_new = projectMetric2Subspace(gij, sSpace2);
assert ( g2_old, g2_new, tol );

sSpace12 = [3,4]; g12_old = projectMetric( projectMetric ( gij, 2 ), 1)(sSpace12,sSpace12); g12_new = projectMetric2Subspace(gij, sSpace12);
assert ( g12_old, g12_new, tol );

Function File: g_aa = projectSuperskyMetric2Sky ( g_nn, alpha0, delta0 )

returns metric in coordinates [a, ...] = [alpha, delta, ... ] at a given sky-position alpha0,delta0, for the given input ’supersky’ metric in coordinates [n, ...] = [ nx,ny,nz, ...].

Note

the supersky-coordinates n = [nx,ny,nz] must be the first 3 coordinates of theinput metric g_nn and the output metric has a = [alpha,delta] as the first 2 coordinates

Examples

## lalapps_FstatMetric_v2 -o gnn.dat --Alpha=1 --Delta=1 --coords="n3x_equ,n3y_equ,n3z_equ,freq,f1dot,f2dot"
## lalapps_FstatMetric_v2 -o gaa.dat --Alpha=1 --Delta=1 --coords="alpha,delta,freq,f1dot,f2dot"
gnn = [
 9.7162064407348633e+04,  1.6162401439666748e+05,  6.9920271533966064e+04,  2.0355045456878051e+07, -1.4588751710937500e+09,  6.5922495973127800e+14;
 1.6162401439666748e+05,  2.6886818181991577e+05,  1.1631187004089355e+05,  3.3860346089623548e+07, -1.2671848879843750e+09,  1.0971813972709620e+15;
 6.9920271533966064e+04,  1.1631187004089355e+05,  5.0317088542461395e+04,  1.4648181998421121e+07, -8.1788653728125000e+08,  4.7469871192193475e+14;
 2.0355045456878051e+07,  3.3860346089623548e+07,  1.4648181998421121e+07,  4.2643423299322696e+09, -2.4811145795203839e+11,  1.3819371887744971e+17;
-1.4588751710937500e+09, -1.2671848879843750e+09, -8.1788653728125000e+08, -2.4811145795203839e+11,  9.2138764319991376e+16, -1.3400518435962927e+19;
 6.5922495973127800e+14,  1.0971813972709620e+15,  4.7469871192193475e+14,  1.3819371887744971e+17, -1.3400518435962927e+19,  5.3316716910947309e+24;
];
gaa = [
 9.4392684805173118e+01, -3.7432745970430642e+03,  6.3033986158817098e+05,  2.9335076680322605e+08,  2.0580635047683145e+13;
-3.7432745970430642e+03,  1.5028899244227595e+05, -2.5315559717512231e+07,  1.1186296727757726e+09, -8.2011995357922125e+14;
 6.3033986158817098e+05, -2.5315559717512231e+07,  4.2643423299322696e+09, -2.4811145795203839e+11,  1.3819371887744971e+17;
 2.9335076680322605e+08,  1.1186296727757726e+09, -2.4811145795203839e+11,  9.2138764319991376e+16, -1.3400518435962927e+19;
 2.0580635047683145e+13, -8.2011995357922125e+14,  1.3819371887744971e+17, -1.3400518435962927e+19,  5.3316716910947309e+24;
];
alpha=1; delta=1;
gaa2 = projectSuperskyMetric2Sky ( gnn, alpha, delta );
tol = -1e-8;        ## relative tolerance passes as <0 !
assert ( gaa, gaa2, tol );

Function File: ret = RandMismatchPDF ( mt, dim, falseDismissal )

Return the mismatch probability density of "relative effective mismatch" mt = m/meff of a random-template bank with given false-dismissal probability

Note

can handle vector input in ’mis’

Examples

assert(RandMismatchPDF(0.5, 3, 0.01), 0.95877, 1e-3)

Function File: ret = RandNormalizedThickness ( dim, falseDismissal )

Return the normalized thickness of an ’dim’-dimensional random-template grid with given falseDismissal probability

Note

can handle vector input in dim

Examples

assert(RandNormalizedThickness(1, 0.05) > 0)

assert(RandNormalizedThickness(2, 0.05) > 0)

assert(RandNormalizedThickness(3, 0.05) > 0)

assert(RandNormalizedThickness(4, 0.05) > 0)

assert(RandNormalizedThickness(5, 0.05) > 0)

Function File: ret = RandThickness ( dim, falseDismissal )

Return the thickness of an dim-dimensional random-template grid with given falseDismissal probability

Note

can handle vector input in dim

Examples

assert(RandThickness(1, 0.05) > 0)

assert(RandThickness(2, 0.05) > 0)

assert(RandThickness(3, 0.05) > 0)

assert(RandThickness(4, 0.05) > 0)

assert(RandThickness(5, 0.05) > 0)

Function File: cache_dir = SuperskyMetricsCache ( )

Function File: SuperskyMetricsCache ‘action’

Manage cache of computed supersky metrics

Arguments

cache_dir

location of cache directory

‘action’

one of

install: Install the precomputed cache from the OctApps repository
clear: Clear the cache
copytorepo: Copy the current cache to the OctApps repository

Examples

cache_dir = SuperskyMetricsCache();

2.7 cw-metric-template-banks

Arguments

Examples

Examples

Arguments

Options

Examples

Arguments

Options

Examples

Examples

Arguments

Examples

Arguments

Examples

Arguments

Examples

Arguments

Options

Examples

Arguments

Options

Examples

Arguments

Options

Examples

Examples

Arguments

Examples

Note

Examples

Arguments

Options

Examples

Examples

Notes

Examples

Note

Examples

Note

Examples

Note

Examples

Note

Examples

Arguments

Examples

2.7 `cw-metric-template-banks`