# LSC-Virgo Burst Analysis Working Group

## Stand-Alone Tests

Many of the BlockNormal pipeline functions can be tested in isolation without use of LIGO frame data and data-conditioning. The following details some test methods.

### Setup

The BlockNormal pipeline is written in MATLAB, so the easiest test method is using MATLAB. The code should work on any version from R14SP3 onward.

#### Get source code

The latest code is always in the 'lscsoft' CVS repository under 'matapps/src/searches/burst/BlockNormal'. Tar-balls are also prepared for tagged releases of the code. The version for the August, 2007 result now under review is BN_2007-07-09, stored as file BN_2007-07-09.tar.gz.

#### Set MATLAB path

You need the directories for the BlockNormal package in your MATLAB path. This can be done at the MATLAB prompt. Assuming you are in the directory you expanded the tar-ball at:
>> addpath src/searches/burst/BlockNormal/BNETG;
>> bnetgpkgpath;


### Create time-series data

BlockNormal works on time-series data. Typical tests use white-noise data. The function wnseries in 'matapps' can be used to create such a series. A copy is also available here at wnseries.m. We will call it at the MATLAB prompt to create a white-noise time-series of unit variance.
>> dfs = 256;
>> duration = 1800;
>> d = wnseries(duration,dfs,1)


### Testing main change-point routine 'cp'

The main change-point routine 'cp.m' has the following interface:
  CP - Generate change-points
cps = cp(d,dfs,r2t,cpcache,cpparms[,bkgdvar])

d         data
dfs       data sample rate
r2t       threshold for recognizing a change-point
cpcache   cache for change-point constants
cpparms    parameters for controlling cp generation
.maxintl - maximum iteration segment (sec)
bkgdvar  [OPTIONAL] variance of background
-

cps(:)     struct array of change-points, sorted by field ndx
.ndx       index of first element new block
.dur       duration of block in samples
.r2        rho2 values corresponding to change-point
.mu        mean of block starting at change-point
.nu        variance of block starting at change-point
.gamma     0th-order centroid param
.lambda    1st-order centroid param
.omega     2nd-order centroid param

So you need a time-series vector, the sample-rate of the series, a ρ2 threshold, a change-point constants cache, and parameter structure

#### Defining change-point parameters

This can be defined at the MATLAB prompt. To set to the default of 10 seconds, use
>> cpparms = struct('maxintl',10);


#### Getting cache of change-point constants

These save time in running 'cp'. They only depend on the maximum number of samples in the time series. This can be done at the MATLAB prompt
>> cpcache = mkr2cache(numel(d));

assuming 'd' is the time-series.

#### Calling 'cp'

Now we can test cp at the MATLAB prompt. The rho2 threshold is logarithmic. Values from 0 to 10 (4 being typical) can be used
r2t = 4;
cps = cp(d,dfs,r2t,cpcache,cpparms);


### Testing ρ2 statistic routine 'rho2'

The ρ2 statistic routine 'rho2.m' at the heart of the change-point determination has the following interface:
 RHO2 - evaluate block-normal statistic for two normal
distributions

[rslt,fom,const] = rho2(data,cache)

data      data to analyze
cache     (optional) Compute-once constants

rslt
.r2     log ratio describing likelihood that two normal
distributions explain data, as opposed to one
.rp     log of the peak contribution to r (value at ndx)
.ndx    index of likeliest location of changepoint
(i.e., ndx is the location of the
first point of the second segment)
fom       vector of log of values contributing in sum to \rho_2
const     cache data for use in subsequent call to rho2

So (again) you need a time-series vector, and a change-point constants cache.

#### Getting cache of change-point constants

These save time in running 'rho2'. They only depend on the maximum number of samples in the time series. This can be done at the MATLAB prompt
>> cpcache = mkr2cache(numel(d));

assuming 'd' is the time-series.

#### Calling 'rho2'

Now we can test rho2 at the MATLAB prompt.
[rslt,fom] = rho2(d,cpcache);

Note that 'fom' will have all the per-sample values of log(r2). The value in rslt.rp is from the sample 'ndex' with the largest ρ2
$Id: Testing_BlockNormal_Code.html,v 1.3 2008/04/09 21:22:43 kathorne Exp$