Changeset 25256 for trunk

Timestamp:

Sep 2, 2009, 2:36:52 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap_mops into trunk. ppMops now merges multiple skycells, and these get published for MOPS as a single file per exposure. Tested and works.

Location:

trunk

Files:

• . (modified) (1 prop)
• dbconfig/changes.txt (modified) (1 diff)
• dbconfig/publish.md (modified) (1 diff)
• doc/misc/docgen.pl (modified) (1 prop)
• ippScripts/scripts/magic_destreak_revert.pl (modified) (1 prop)
• ippScripts/scripts/publish_file.pl (modified) (8 diffs)
• ippTools/share/difftool_skyfile.sql (modified) (3 diffs)
• ippTools/share/pxadmin_create_tables.sql (modified) (1 diff)
• ippTools/src/pubtool.c (modified) (2 diffs)
• ippTools/src/pubtoolConfig.c (modified) (1 diff)
• ppMops/ICDlite.txt (modified) (2 diffs)
• ppMops/src/Makefile.am (modified) (1 diff)
• ppMops/src/ppMops.c (modified) (2 diffs)
• ppMops/src/ppMops.h (modified) (1 diff)
• ppMops/src/ppMopsArguments.c (copied) (copied from branches/pap_mops/ppMops/src/ppMopsArguments.c )
• ppMops/src/ppMopsData.c (deleted)
• ppMops/src/ppMopsDetections.c (copied) (copied from branches/pap_mops/ppMops/src/ppMopsDetections.c )
• ppMops/src/ppMopsMerge.c (copied) (copied from branches/pap_mops/ppMops/src/ppMopsMerge.c )
• ppMops/src/ppMopsRead.c (copied) (copied from branches/pap_mops/ppMops/src/ppMopsRead.c )
• ppMops/src/ppMopsWrite.c (copied) (copied from branches/pap_mops/ppMops/src/ppMopsWrite.c )
• ppStack/src/ppStackMatch.c (modified) (1 diff)
• psLib/src/sys/psType.h (modified) (4 diffs)
• psLib/src/types/psTree.c (modified) (21 diffs)
• psLib/src/types/psTree.h (modified) (5 diffs)
• psLib/test/optime (modified) (1 prop)
• psLib/test/types/tap_psTree.c (modified) (4 diffs)
• psModules/src/camera/pmReadoutFake.c (modified) (1 prop)
• psModules/src/objects/pmSourceMatch.c (modified) (1 diff)
• psphot/doc/efficiency.txt (modified) (1 prop)
• psphot/src/psphotFake.c (modified) (1 prop)

trunk/dbconfig/changes.txt

@@ -1206 +1206 @@
 ALTER TABLE magicDSRun ADD COLUMN inv_magic_id BIGINT AFTER magic_id;
 ALTER TABLE magicDSRun ADD FOREIGN KEY (inv_magic_id) REFERENCES magicRun(magic_id);
+
+-- Adding columns for debugging publishing problems
+
+ALTER TABLE publishDone ADD COLUMN hostname VARCHAR(64) AFTER path_base;
+ALTER TABLE publishDone ADD COLUMN dtime_script FLOAT AFTER hostname;

trunk/dbconfig/publish.md

@@ -1 +1 @@
 # Tables for publishing data to a science client
 
-publishClient   METADATA 
-    client_id   S64         0       # Primary Key AUTO_INCREMENT
-    product     STR         64
-    stage       STR         64
-    workdir     STR         255
-    comment     STR         255
+publishClient    METADATA 
+    client_id    S64         0       # Primary Key AUTO_INCREMENT
+    product      STR         64
+    stage        STR         64
+    workdir      STR         255
+    comment      STR         255
+END              
+                 
+publishRun       METADATA
+    pub_id       S64         0       # Primary Key AUTO_INCREMENT
+    client_id    S64         0
+    stage_id     S64         0
+    label        STR         64
+    state        STR         64
+END              
+                 
+publishDone      METADATA
+    pub_id       S64         0       # Primary Key
+    path_base    STR         255
+    hostname     STR         64
+    dtime_script F32         0.0
+    fault        S16         0
 END
-
-publishRun      METADATA
-    pub_id      S64         0       # Primary Key AUTO_INCREMENT
-    client_id   S64         0
-    stage_id    S64         0
-    label       STR         64
-    state       STR         64
-END
-
-publishDone     METADATA
-    pub_id      S64         0       # Primary Key
-    path_base   STR         255
-    fault       S16         0
-END

trunk/ippScripts/scripts/publish_file.pl

@@ -39 +39 @@
 
 # Parse the command-line arguments
-my ( $pub_id, $camera, $stage, $stage_id, $format, $product, $workdir );
-my ( $dbname, $verbose, $no_update, $save_temps, $redirect );
+my ( $pub_id, $camera, $stage, $stage_id, $fileset, $format, $product, $workdir );
+my ( $dbname, $verbose, $no_update, $no_op, $save_temps, $redirect );
 
 GetOptions(
@@ -48 +48 @@
     'stage_id=s'        => \$stage_id,    # Stage identifier
     'product=s'         => \$product,     # Datastore product name
+    'fileset=s'         => \$fileset,     # Fileset name
     'workdir=s'         => \$workdir,     # Working directory
     'dbname=s'          => \$dbname,    # Database name
     'verbose'           => \$verbose,   # Print to stdout
     'no-update'         => \$no_update, # Don't update the database?
+    'no-op'             => \$no_op, # Don't do any operations
     'save-temps'        => \$save_temps, # Save temporary files?
     'redirect-output'   => \$redirect,   # Redirect output to log file?
@@ -77 +79 @@
 my $mdcParser = PS::IPP::Metadata::Config->new;
 
-
-# Retrieve file name of interest
-my %files;                      # Input filenames
-my %zp;                         # Zero points
-my %exp_id;                     # Exposure identifiers
-my %exp_name;                   # Exposure names
-my %direction;                  # Direction of subtraction
-{
-    my $command;                # Command to run
-
-    if ($stage eq 'diff') {
-        $command =  "difftool -diffskyfile -diff_id $stage_id";
-    } elsif ($stage eq 'camera') {
-        $command =  "camtool -processedexp -cam_id $stage_id";
-    } else {
-        &my_die( "Unrecognised stage: $stage", $pub_id, $PS_EXIT_CONFIG_ERROR );
-    }
+my ($dsFile, $dsFileName) = tempfile("/tmp/publish.$pub_id.ds.XXXX", UNLINK => !$save_temps );
+
+if ($stage eq 'camera') {
+    my $command =  "camtool -processedexp -cam_id $stage_id";
     $command .= " -dbname $dbname" if defined $dbname;
-
     my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
         run(command => $command, verbose => $verbose);
@@ -106 +94 @@
         &my_die("Unable to parse metadata list", $pub_id, $PS_EXIT_PROG_ERROR);
 
+    &my_die("More than one entry for cam_id $stage_id", $pub_id, $PS_EXIT_PROG_ERROR) if scalar @$components > 1;
+
+    my $comp = $$components[0]; # Component of interest
+    my $path_base = $comp->{path_base}; # Base name for file
+    my $file = $ipprc->filename( "PSASTRO.OUTPUT", $path_base );
+    $file = $ipprc->file_resolve($file);
+    my $cam_id = $comp->{cam_id};
+    my $name = "cam_$cam_id";
+    print $dsFile "$file|||$product|$name|\n";
+} elsif ($stage eq 'diff') {
+    my $command =  "difftool -diffskyfile -diff_id $stage_id";
+    $command .= " -dbname $dbname" if defined $dbname;
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+        run(command => $command, verbose => $verbose);
+    &my_die( "Unable to retrieve filename", $pub_id, $PS_EXIT_SYS_ERROR) unless $success;
+
+    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or
+        &my_die("Unable to parse metadata config", $pub_id, $PS_EXIT_PROG_ERROR);
+
+    my $components = parse_md_list($metadata) or
+        &my_die("Unable to parse metadata list", $pub_id, $PS_EXIT_PROG_ERROR);
+
+    my ($mopsPositiveFile, $mopsPositiveFileName) = tempfile("/tmp/publish.$pub_id.mops.pos.XXXX", UNLINK => !$save_temps ) if $product eq 'IPP-MOPS';
+    my ($mopsNegativeFile, $mopsNegativeFileName) = tempfile("/tmp/publish.$pub_id.mops.neg.XXXX", UNLINK => !$save_temps ) if $product eq 'IPP-MOPS';
+
+    my %positive;               # Data for positive diff detections
+    my %negative;               # Data for negative diff detections
     foreach my $comp ( @$components ) {
         my $path_base = $comp->{path_base}; # Base name for file
@@ -112 +127 @@
         (carp "Bad zpt_obs or exp_time for component" and next) if not defined $comp->{zpt_obs} or not defined $comp->{exp_time};
         my $zp = $comp->{zpt_obs} + 2.5 * log($comp->{exp_time}) / log(10);
-
-        if ($stage eq 'diff') {
-            my $skycell_id = $comp->{skycell_id};
-
-            # Positive direction
-            $files{"$skycell_id.pos"} = $ipprc->filename( "PPSUB.OUTPUT.SOURCES", $path_base );
-            $zp{"$skycell_id.pos"} = $zp;
-            $exp_id{"$skycell_id.pos"} = $comp->{exp_id_1};
-            $exp_name{"$skycell_id.pos"} = $comp->{exp_name_1};
-            $direction{"$skycell_id.pos"} = 1;
-
-            # Negative direction
-            if (defined $comp->{bothways} and $comp->{bothways}) {
-                $files{"$skycell_id.neg"} = $ipprc->filename( "PPSUB.INVERSE.SOURCES", $path_base );
-                $zp{"$skycell_id.neg"} = $zp;
-                $exp_id{"$skycell_id.neg"} = $comp->{exp_id_2};
-                $direction{"$skycell_id.neg"} = 0;
-                $exp_name{"$skycell_id.neg"} = $comp->{exp_name_2};
+        my $astrom = sqrt($comp->{sigma_ra_1}**2 + $comp->{sigma_dec_1}**2);
+
+        my $skycell_id = $comp->{skycell_id};
+        my $filename = $ipprc->filename( "PPSUB.OUTPUT.SOURCES", $path_base );
+        $filename = $ipprc->file_resolve($filename);
+
+        my $data = { zp => $zp,
+                     zp_err => $comp->{zpt_stdev},
+                     astrom => sqrt($comp->{sigma_ra_1}**2 + $comp->{sigma_dec_1}**2),
+                     exp_name => $comp->{exp_name_1},
+                     exp_id => $comp->{exp_id_1},
+                     chip_id => $comp->{chip_id_1},
+                     cam_id => $comp->{cam_id_1},
+                     fake_id => $comp->{fake_id_1},
+                     warp_id => $comp->{warp1},
+                     diff_id => $comp->{diff_id},
+                     direction => 1,
+        };
+
+        diff_check(\%positive, $data, "positive");
+
+        if ($product eq 'IPP-MOPS') {
+            print $mopsPositiveFile "$filename\n";
+        } else {
+            print $dsFile "$filename|||$product|${skycell_id}.pos|\n";
+        }
+
+        # Negative direction
+        if (defined $comp->{bothways} and $comp->{bothways}) {
+            my $filename = $ipprc->filename( "PPSUB.INVERSE.SOURCES", $path_base );
+            $filename = $ipprc->file_resolve($filename);
+
+            my $data = { zp => $zp,
+                         zp_err => $comp->{zpt_stdev},
+                         astrom => sqrt($comp->{sigma_ra_2}**2 + $comp->{sigma_dec_2}**2),
+                         exp_name => $comp->{exp_name_2},
+                         exp_id => $comp->{exp_id_2},
+                         chip_id => $comp->{chip_id_2},
+                         cam_id => $comp->{cam_id_2},
+                         fake_id => $comp->{fake_id_2},
+                         warp_id => $comp->{warp2},
+                         diff_id => $comp->{diff_id},
+                         direction => 0,
+            };
+
+            diff_check(\%negative, $data, "negative");
+
+            if ($product eq 'IPP-MOPS') {
+                print $mopsNegativeFile "$filename\n";
+            } else {
+                print $dsFile "$filename|||$product|${skycell_id}.neg|\n";
             }
-
-        } elsif ($stage eq 'camera') {
-            my $cam_id = $comp->{cam_id};
-            $files{$cam_id} = $ipprc->filename( "PSASTRO.OUTPUT", $path_base );
-            $zp{$cam_id} = $zp;
-            $exp_id{$cam_id} = $comp->{exp_id};
-            $exp_name{$cam_id} = $comp->{exp_name};
-            $direction{$cam_id} = 1;
         }
     }
-}
-
-# Prepare for data store input
-my ($listFile, $listFileName) = tempfile("/tmp/publish.$pub_id.list.XXXX", UNLINK => !$save_temps );
-
-# Process each file
-foreach my $comp ( keys %files ) {
-    my $file = $ipprc->file_resolve( $files{$comp} ) or
-        &my_die("Unable to resolve file $files{$comp}", $pub_id, $PS_EXIT_SYS_ERROR);
-    my_die("Unable to find file $file", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $file );
-
-    my $zp = $zp{$comp};
-    my $exp_id = $exp_id{$comp};
-    my $exp_name = $exp_name{$comp};
-    my $direction = $direction{$comp};
-    if ($product eq "IPP-MOPS") {
-        my $outuri = "$outroot.$comp.fits";
-        my $out = $ipprc->file_resolve( $outuri, 'create' ) or
-            &my_die( "Unable to resolve output file $outuri", $pub_id, $PS_EXIT_SYS_ERROR);
-
-        my $command = "$ppMops $file $zp $exp_id $exp_name $direction $out";
-        my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
-            run(command => $command, verbose => $verbose);
-        &my_die( "Unable to translate $file", $pub_id, $PS_EXIT_SYS_ERROR) unless $success;
-        &my_die( "Unable to find translated file $out", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $out );
-        $file = $out;
+
+    close $mopsPositiveFile;
+    close $mopsNegativeFile;
+
+    if ($product eq 'IPP-MOPS') {
+        if (scalar keys %positive > 0) {
+            my $output = mops_combine(\%positive, "$outroot.pos.mops", $mopsPositiveFileName);
+            print $dsFile "$output|||$product|positive|\n";
+        }
+        if (scalar keys %negative > 0) {
+            my $output = mops_combine(\%negative, "$outroot.neg.mops", $mopsNegativeFileName);
+            print $dsFile "$output|||$product|negative|\n";
+        }
     }
-
-    # format: filename|filesize|md5sum|filetype|
-    # note: since we omit filesize and md5sum, dsreg will calculate them
-    print $listFile "$file|||$product|$comp|\n";
-}
-
+}
+
+close $dsFile;
 
 unless ($no_update) {
-    my $command = "$dsreg --add $stage.$stage_id --copy --abspath --product $product --type $product --list $listFileName";
+    my $command = "$dsreg --add $stage.$stage_id --copy --abspath --product $product --type $product --list $dsFileName";
 
     my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
@@ -185 +209 @@
 unless ($no_update) {
     my $command = "$pubtool -add -pub_id $pub_id -path_base $outroot";
+    $command .= (" -dtime_script " . ((DateTime->now->mjd - $mjd_start) * 86400));
+    $command .= " -hostname $host" if defined $host;
+
     $command .= " -dbname $dbname" if defined $dbname;
 
@@ -195 +222 @@
 
 ### Pau.
+
+# Check inputs for a diff
+sub diff_check
+{
+    my $data = shift;           # Data hash
+    my $comp = shift;           # Component data
+    my $name = shift;           # Name of component
+
+    $data->{zp} = $comp->{zp} unless defined $data->{zp};
+    $data->{zp_err} = $comp->{zp_err} unless defined $data->{zp_err};
+    $data->{astrom} = $comp->{astrom} unless defined $data->{astrom};
+    $data->{exp_name} = $comp->{exp_name} unless defined $data->{exp_name};
+    $data->{exp_id} = $comp->{exp_id} unless defined $data->{exp_id};
+    $data->{chip_id} = $comp->{chip_id} unless defined $data->{chip_id};
+    $data->{cam_id} = $comp->{cam_id} unless defined $data->{cam_id};
+    $data->{fake_id} = $comp->{fake_id} unless defined $data->{fake_id};
+    $data->{warp_id} = $comp->{warp_id} unless defined $data->{warp_id};
+    $data->{diff_id} = $comp->{diff_id} unless defined $data->{diff_id};
+    $data->{direction} = $comp->{direction} unless defined $data->{direction};
+
+    &my_die("zp value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{zp} and $comp->{zp} != $data->{zp};
+    &my_die("zp_err value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{zp_err} and $comp->{zp_err} != $data->{zp_err};
+    &my_die("astrom value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{astrom} and $comp->{astrom} != $data->{astrom};
+    &my_die("exp_name value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{exp_name} and $comp->{exp_name} ne $data->{exp_name};
+    &my_die("exp_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{exp_id} and $comp->{exp_id} != $data->{exp_id};
+    &my_die("chip_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{chip_id} and $comp->{chip_id} != $data->{chip_id};
+    &my_die("cam_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{cam_id} and $comp->{cam_id} != $data->{cam_id};
+    &my_die("fake_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{fake_id} and $comp->{fake_id} != $data->{fake_id};
+    &my_die("warp_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{warp_id} and $comp->{warp_id} != $data->{warp_id};
+    &my_die("diff_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{diff_id} and $comp->{diff_id} != $data->{diff_id};
+    &my_die("direction value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{direction} and $comp->{direction} != $data->{direction};
+
+    return 1;
+}
+
+# Combine multiple files for MOPS
+sub mops_combine
+{
+    my $data = shift;           # Data
+    my $output = shift;         # Output name
+    my $input = shift;          # Input name
+
+    $output = $ipprc->file_resolve( $output, 'create' ) or
+        &my_die( "Unable to resolve output file $output", $pub_id, $PS_EXIT_SYS_ERROR);
+
+    my $command = "$ppMops $input $output";
+    $command .= " -exp_name " . $data->{exp_name} if defined $data->{exp_name};
+    $command .= " -exp_id " . $data->{exp_id} if defined $data->{exp_id};
+    $command .= " -chip_id " . $data->{chip_id} if defined $data->{chip_id};
+    $command .= " -cam_id " . $data->{cam_id} if defined $data->{cam_id};
+    $command .= " -fake_id " . $data->{fake_id} if defined $data->{fake_id};
+    $command .= " -warp_id " . $data->{warp_id} if defined $data->{warp_id};
+    $command .= " -diff_id " . $data->{diff_id} if defined $data->{diff_id};
+    $command .= " -inverse" if defined $data->{direction} and $data->{direction} == 0;
+    $command .= " -zp " . $data->{zp} if defined $data->{zp};
+    $command .= " -zp_error " . $data->{zp_err} if defined $data->{zp_err};
+    $command .= " -astrom_rms " . $data->{astrom} if defined $data->{astrom};
+
+    unless ($no_op) {
+        my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+            run(command => $command, verbose => $verbose);
+        &my_die( "Unable to translate", $pub_id, $PS_EXIT_SYS_ERROR) unless $success;
+        &my_die( "Unable to find translated file $output", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $output );
+    } else {
+        print "Not running: $command\n";
+    }
+
+    return $output;
+}
 
 
@@ -210 +306 @@
         $command .= " -pub_id $pub_id";
         $command .= " -path_base $outroot";
+        $command .= (" -dtime_script " . ((DateTime->now->mjd - $mjd_start) * 86400));
+        $command .= " -hostname $host" if defined $host;
         $command .= " -fault $fault";
         $command .= " -dbname $dbname" if defined $dbname;

trunk/ippTools/share/difftool_skyfile.sql

@@ -12 +12 @@
     -- The following are only valid for warps
     -- XXX This needs to be more clever to handle diffs between stacks
+    -- Zero points are appropriate for both forward and backward diffs
     camProcessedInput.zpt_obs,
     camProcessedInput.zpt_stdev,
@@ -19 +20 @@
     rawInput.camera,
     rawInput.exp_name AS exp_name_1,
+    rawInput.exp_id AS exp_id_1,
+    chipInput.chip_id AS chip_id_1,
+    camInput.cam_id AS cam_id_1,
+    fakeInput.fake_id AS fake_id_1,
+    camProcessedInput.sigma_ra AS sigma_ra_1,
+    camProcessedInput.sigma_dec AS sigma_dec_1,
     rawTemplate.exp_name AS exp_name_2,
-    rawInput.exp_id AS exp_id_1,
-    rawTemplate.exp_id AS exp_id_2
+    rawTemplate.exp_id AS exp_id_2,
+    chipTemplate.chip_id AS chip_id_2,
+    camTemplate.cam_id AS cam_id_2,
+    fakeTemplate.fake_id AS fake_id_2,
+    camProcessedTemplate.sigma_ra AS sigma_ra_2,
+    camProcessedTemplate.sigma_dec AS sigma_dec_2
 FROM diffRun
 JOIN diffSkyfile USING(diff_id)
@@ -43 +54 @@
 LEFT JOIN camRun AS camTemplate
     ON camTemplate.cam_id = fakeTemplate.cam_id
+LEFT JOIN camProcessedExp AS camProcessedTemplate
+    ON camProcessedTemplate.cam_id = camTemplate.cam_id
 LEFT JOIN chipRun AS chipTemplate
     ON chipTemplate.chip_id = camTemplate.chip_id

trunk/ippTools/share/pxadmin_create_tables.sql

@@ -1407 +1407 @@
     pub_id BIGINT AUTO_INCREMENT, -- link to publishRun
     path_base VARCHAR(255),     -- base path of output
+    hostname VARCHAR(64),       -- name of host
+    dtime_script FLOAT,         -- run time for script
     fault SMALLINT NOT NULL DEFAULT 0, -- Fault code
     PRIMARY KEY(pub_id),

trunk/ippTools/src/pubtool.c

@@ -259 +259 @@
     // required
     PXOPT_LOOKUP_S64(pub_id, config->args, "-pub_id", true, false);
-    PXOPT_LOOKUP_STR(path_base, config->args, "-path_base",  true, false);
+    PXOPT_LOOKUP_STR(path_base, config->args, "-path_base", true, false);
 
     // optional
+    PXOPT_LOOKUP_STR(hostname, config->args, "-hostname", false, false);
+    PXOPT_LOOKUP_F32(dtime_script, config->args, "-dtime_script", false, false);
     PXOPT_LOOKUP_S32(fault, config->args, "-fault", false, false);
 
@@ -269 +271 @@
     }
 
-    if (!publishDoneInsert(config->dbh, pub_id, path_base, fault)) {
+    if (!publishDoneInsert(config->dbh, pub_id, path_base, hostname, dtime_script, fault)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to add file");
         if (!psDBRollback(config->dbh)) {

trunk/ippTools/src/pubtoolConfig.c

@@ -68 +68 @@
     psMetadataAddS64(addArgs, PS_LIST_TAIL, "-pub_id", 0, "define pub_id (required)", 0);
     psMetadataAddStr(addArgs, PS_LIST_TAIL, "-path_base", 0, "define path_base (required)", NULL);
+    psMetadataAddStr(addArgs, PS_LIST_TAIL, "-hostname", 0, "define hostname", NULL);
+    psMetadataAddF32(addArgs, PS_LIST_TAIL, "-dtime_script", 0, "define time for script", NAN);
     psMetadataAddS32(addArgs, PS_LIST_TAIL, "-fault", 0, "define fault code", 0);
 

trunk/ppMops/ICDlite.txt

@@ -79 +79 @@
  * ASTRORMS added to make absolute (i.e., across exposures) astrometry errors more accurate
  * DE_MAGnn and DE_EFFnn replace DE1 through DE10, which were never well defined
+ * Removed LIMITMAG, which was never well defined, and is unnecessary given the DE_MAGnn and DE_EFFnn
 
 === Table ===
@@ -119 +120 @@
 
 
+=== Example ===
 
-
-                            Table 3: IPP-MOPS Transient Detection Keywords
-FITS Keyword      Precision                     Description
-FPA ID            char(20)                      IPP-assigned identifier that can be used to trace back to this FPA
-MJD-OBS           F64                           midpoint time of the exposure as a MJD and day fraction
-RA                HH:MM:SS.SSS                 field center RA at exposure midpoint, string
-DEC               sDD:MM:SS (s is + or -) field center declination at exposure midpoint, string
-EXPTIME           F64                          exposure time in seconds
-ROTANGLE          F64                          image rotation of the y-axis in degrees, from local N toward E
-FILTER            char(3)                       effective filter used for the exposure, one of g, r, i, z, y, B, V, w
-AIRMASS           F64                          airmass at exposure midpoint
-LIMITMAG          F64                           limiting magnitude of detections for this FPA
-DE1 through DE10  F64                           detection efficiency coefficients
-OBSCODE           char(5)                       MPC three-character observatory code
-TEL ALT           F64                          telescope altitude at exposure midpoint, in degrees
-TEL AZ            F64                          telescope azimuth at exposure midpoint, in degrees
-                         Table 4: IPP-MOPS Transient Detection Table Columns
-             Column       Precision Description
-             RA DEG       F64          detection center coordinates in degrees
-             RA SIG       F64          error in the detection center, in degrees
-             DEC DEG      F64          detection center coordinates in degrees
-             DEC SIG      F64          error in the detection center, in degrees
-             FLUX         F64          flux
-             FLUX SIG     F64          error in the flux value
-             STARPSF      F64          probability that the PSF matches a starlike PSF
-             ANG          F64          detection’s elongation angle in degrees, local N toward E
-             ANG SIG      F64          error in the angle, in degrees
-             LEN          F64          elongation length of the detection in degrees
-             LEN SIG      F64          error in the length, in degrees
+{{{
+XTENSION= 'BINTABLE'           / binary table extension
+BITPIX  =                    8 / 8-bit bytes
+NAXIS   =                    2 / 2-dimensional binary table
+NAXIS1  =                   72 / width of table in bytes
+NAXIS2  =                42032 / number of rows in table
+PCOUNT  =                    0 / size of special data area
+GCOUNT  =                    1 / one data group (required keyword)
+TFIELDS =                   13 / number of fields in each row
+TTYPE1  = 'RA      '           / label for field   1
+TFORM1  = '1D      '           / data format of field: 8-byte DOUBLE
+TTYPE2  = 'RA_ERR  '           / label for field   2
+TFORM2  = '1D      '           / data format of field: 8-byte DOUBLE
+TTYPE3  = 'DEC     '           / label for field   3
+TFORM3  = '1D      '           / data format of field: 8-byte DOUBLE
+TTYPE4  = 'DEC_ERR '           / label for field   4
+TFORM4  = '1D      '           / data format of field: 8-byte DOUBLE
+TTYPE5  = 'MAG     '           / label for field   5
+TFORM5  = '1E      '           / data format of field: 4-byte REAL
+TTYPE6  = 'MAG_ERR '           / label for field   6
+TFORM6  = '1E      '           / data format of field: 4-byte REAL
+TTYPE7  = 'STARPSF '           / label for field   7
+TFORM7  = '1E      '           / data format of field: 4-byte REAL
+TTYPE8  = 'ANGLE   '           / label for field   8
+TFORM8  = '1E      '           / data format of field: 4-byte REAL
+TTYPE9  = 'ANGLE_ERR'          / label for field   9
+TFORM9  = '1E      '           / data format of field: 4-byte REAL
+TTYPE10 = 'LENGTH  '           / label for field  10
+TFORM10 = '1E      '           / data format of field: 4-byte REAL
+TTYPE11 = 'LENGTH_ERR'         / label for field  11
+TFORM11 = '1E      '           / data format of field: 4-byte REAL
+TTYPE12 = 'FLAGS   '           / label for field  12
+TFORM12 = '1J      '           / data format of field: 4-byte INTEGER
+TZERO12 =           2147483648 / offset for unsigned integers
+TSCAL12 =                    1 / data are not scaled
+TTYPE13 = 'DIFF_SKYFILE_ID'    / label for field  13
+TFORM13 = '1K      '           / data format of field: 8-byte INTEGER
+SWSOURCE= '60eb6cdc-a59c-4636-a4e0-dba66a9721fd' / Software source
+SWVERSN = 'branches/pap_mops/ppMops@25227' / Software version
+HISTORY ppMops at 2009-09-02T03:48:46.695783
+HISTORY psLib version: branches/pap_mops/psLib@25227
+HISTORY psLib source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd
+HISTORY ppMops version: branches/pap_mops/ppMops@25227
+HISTORY ppMops source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd
+EXP_NAME= 'o4995g0129o'        / Exposure name
+EXP_ID  =                77164 / Exposure identifier
+CHIP_ID =                24019 / Chip stage identifier
+CAM_ID  =                17726 / Cam stage identifier
+FAKE_ID =                10227 / Fake stage identifier
+WARP_ID =                 8842 / Warp stage identifier
+DIFF_ID =                    0 / Diff stage identifier
+DIFF_POS=                    F / Positive subtraction?
+MJD-OBS =     54995.4740598313 / MJD of exposure midpoint
+RA      = '18:25:01.988'       / Right Ascension of boresight
+DEC     = '-17:20:40.069'      / Declination of boresight
+TEL_ALT =            51.951873 / Telescope altitude
+TEL_AZ  =           179.483883 / Telescope azimuth
+EXPTIME =                  38. / Exposure time (sec)
+ROTANGLE=             333.1039 / Rotator position angle
+FILTER  = 'r.00000 '           / Filter name
+AIRMASS =                1.269 / Airmass of exposure
+OBSCODE = 'F51     '           / IAU Observatory code
+SEEING  =             1.678401 / Mean seeing
+MAGZP   =             28.65226 / Magnitude zero point
+MAGZPERR=             0.353063 / Error in magnitude zero point
+ASTRORMS=            0.3111496 / RMS of astrometric fit
+EXTNAME = 'MOPS_TRANSIENT_DETECTIONS'
+END
+}}}

trunk/ppMops/src/Makefile.am

@@ -28 +28 @@
         ppMops.c                \
         ppMopsVersion.c         \
-        ppMopsData.c                    
+        ppMopsArguments.c       \
+        ppMopsDetections.c      \
+        ppMopsRead.c            \
+        ppMopsWrite.c           \
+        ppMopsMerge.c
 
 noinst_HEADERS = \

trunk/ppMops/src/ppMops.c

@@ -6 +6 @@
 int main(int argc, char *argv[])
 {
-    if (argc != 7) {
-        fprintf(stderr, "Insufficient arguments.\n");
-        fprintf(stderr, "Usage: %s DETECTIONS ZP EXP_ID EXP_NAME DIRECTION OUTPUT\n", argv[0]);
+    psLibInit(NULL);
+
+    ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments
+    if (!args) {
+        psErrorStackPrint(stderr, "Error parsing arguments");
         exit(PS_EXIT_CONFIG_ERROR);
     }
 
-    ppMopsData *data = ppMopsDataAlloc(); // Configuration data
-    data->detections = psStringCopy(argv[1]);
-    data->zp = atof(argv[2]);
-    data->exp_id = atoll(argv[3]);
-    data->exp_name = psStringCopy(argv[4]);
-    data->direction = atoi(argv[5]);
-    data->output = psStringCopy(argv[6]);
-
-    if (!isfinite(data->zp)) {
-        psErrorStackPrint(stderr, "Zero point is unknown\n");
-        exit(PS_EXIT_CONFIG_ERROR);
-    }
+    psArray *detections = ppMopsRead(args); // Detections from each input
+    if (!detections) {
+        psErrorStackPrint(stderr, "Unable to read detections");
+        exit(PS_EXIT_SYS_ERROR);
+    }
+
+    ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections
+    psFree(detections);
+    if (!merged) {
+        psErrorStackPrint(stderr, "Unable to merge detections");
+        exit(PS_EXIT_SYS_ERROR);
+    }
+
+    if (!ppMopsWrite(merged, args)) {
+        psErrorStackPrint(stderr, "Unable to write detections");
+        exit(PS_EXIT_SYS_ERROR);
+    }
+
+    psFree(merged);
+    psFree(args);
+
+    psLibFinalize();
+
+    return PS_EXIT_SUCCESS;
+}
+
+
+#if 0
+    ps
+
 
     psArray *detections = NULL;         // Detections
@@ -211 +231 @@
     psFree(data);
 
-    return PS_EXIT_SUCCESS;
-}
+#endif
+

trunk/ppMops/src/ppMops.h

@@ -11 +11 @@
                      PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SKY_FAILURE) // Flags to exclude
 
-
 // Configuration data
 typedef struct {
-    psString detections;                // Detections filename
-    float zp;                           // Magnitude zero point
+    psArray *input;                     // Input filenames
+    psString exp_name;                  // Exposure name
     psS64 exp_id;                       // Exposure identifier
-    psString exp_name;                  // Exposure name
-    bool direction;                     // Direction of subtraction, 1=positive, 0=negative
+    psS64 chip_id;                      // Chip stage identifier
+    psS64 cam_id;                       // Camera stage identifier
+    psS64 fake_id;                      // Fake stage identifier
+    psS64 warp_id;                      // Warp stage identifier
+    psS64 diff_id;                      // Diff stage identifier
+    bool positive;                      // Sense of subtraction, T=positive, F=negative
+    float zp, zpErr;                    // Magnitude zero point and error
+    float rmsAstrom;                    // Astrometric solution RMS
     psString output;                    // Output filename
-} ppMopsData;
+} ppMopsArguments;
 
-// Allocator
-ppMopsData *ppMopsDataAlloc(void);
+/// Parse arguments
+ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]);
+
+typedef struct {
+    psString raBoresight, decBoresight; // RA,Dec of telescope boresight
+    psString filter;                    // Filter for exposure
+    float airmass;                      // Airmass of exposure
+    float exptime;                      // Exposure time
+    double posangle;                    // Position angle
+    double alt, az;                     // Telescope altitude and azimuth
+    double mjd;                         // Modified Julian Date
+    float seeing;                       // Seeing of exposure
+    long num;                           // Number of detections
+    psVector *x, *y;                    // Image coordinates
+    psVector *ra, *dec;                 // Sky coordinates
+    psVector *raErr, *decErr;           // Error in sky coordinates
+    psVector *mag, *magErr;             // Magnitude and associated error
+    psVector *extended;                 // Measure of extendedness
+    psVector *angle, *angleErr;         // Angle of trail and associated error
+    psVector *length, *lengthErr;       // Length of trail and associated error
+    psVector *flags;                    // psphot flags
+    psVector *diffSkyfileId;            // Identifier for source image
+    psVector *naxis1, *naxis2;          // Size of image
+    psVector *mask;                     // Mask for detections
+} ppMopsDetections;
+
+ppMopsDetections *ppMopsDetectionsAlloc(long num);
+
+/// Copy a detection
+bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index);
+
+/// Purge the detections list of masked detections
+bool ppMopsDetectionsPurge(ppMopsDetections *detections);
+
+
+/// Read detections
+psArray *ppMopsRead(const ppMopsArguments *args);
+
+/// Merge detections
+ppMopsDetections *ppMopsMerge(const psArray *detections);
+
+/// Write detections
+bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args);
 
 

trunk/ppStack/src/ppStackMatch.c

@@ -87 +87 @@
     x->n = y->n = numGood;
 
-    psTree *tree = psTreePlant(2, 2, x, y); // kd tree
+    psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree
 
     psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources

trunk/psLib/src/sys/psType.h

@@ -141 +141 @@
 } psDataType;
 
-// macros to abstract the generic mask type : these values must be consistent 
+// macros to abstract the generic mask type : these values must be consistent
 #define PS_TYPE_MASK PS_TYPE_U8        /**< the psElemType to use for mask image */
 #define PS_TYPE_MASK_DATA U8           /**< the data member to use for mask image */
@@ -152 +152 @@
 // alternate versions if needed
 // #define PS_NOT_MASK(A)(UINT16_MAX-(A))
-// #define PS_NOT_MASK(A)(UINT32_MAX-(A)) 
+// #define PS_NOT_MASK(A)(UINT32_MAX-(A))
 // #define PS_NOT_MASK(A)(UINT64_MAX-(A))
 
-// macros to abstract the vector mask type : these values must be consistent 
+// macros to abstract the vector mask type : these values must be consistent
 #define PS_TYPE_VECTOR_MASK PS_TYPE_U8        /**< the psElemType to use for mask image */
 #define PS_TYPE_VECTOR_MASK_DATA U8           /**< the data member to use for mask image */
@@ -161 +161 @@
 #define PS_MIN_VECTOR_MASK_TYPE 0             /**< minimum valid Vector Mask value */
 #define PS_MAX_VECTOR_MASK_TYPE UINT8_MAX     /**< maximum valid Vector Mask value */
-typedef psU8 psVectorMaskType;                    ///< the C datatype for a mask image
+typedef psU8 psVectorMaskType;                    ///< the C datatype for a mask image
 #define PS_NOT_VECTOR_MASK(A)(UINT8_MAX-(A))
 
-// macros to abstract the image mask type : these values must be consistent 
-#define PS_TYPE_IMAGE_MASK PS_TYPE_U16       /**< the psElemType to use for mask image */
-#define PS_TYPE_IMAGE_MASK_DATA U16          /**< the data member to use for mask image */
-#define PS_TYPE_IMAGE_MASK_NAME "psU16"      /**< the data type for mask as a string */
+// macros to abstract the image mask type : these values must be consistent
+#define PS_TYPE_IMAGE_MASK PS_TYPE_U16       /**< the psElemType to use for mask image */
+#define PS_TYPE_IMAGE_MASK_DATA U16          /**< the data member to use for mask image */
+#define PS_TYPE_IMAGE_MASK_NAME "psU16"      /**< the data type for mask as a string */
 #define PS_MIN_IMAGE_MASK_TYPE 0             /**< minimum valid Image Mask value */
 #define PS_MAX_IMAGE_MASK_TYPE UINT16_MAX    /**< maximum valid Image Mask value */
@@ -246 +246 @@
 };
 
-/// Macro to get the bad pixel reason code (stored as part of mask value)
-#define PS_BADPIXEL_BITMASK 0x0f
-#define PS_GET_BADPIXEL(maskValue) (maskValue & PS_BADPIXEL_BITMASK)
-
-#define PS_IS_BADPIXEL(maskValue) (PS_GET_BADPIXEL(maskValue) != 0)
-
-/// Macro to apply a bad pixel reason code to mask image
-#define PS_SET_BADPIXEL(maskValue, reasonCode) \
-{ \
-    maskValue = (psMaskType)((reasonCode & PS_BADPIXEL_BITMASK) | (maskValue & ~PS_BADPIXEL_BITMASK)); \
-}
-
 /// Macro to determine if the psElemType is an integer.
 #define PS_IS_PSELEMTYPE_INT(x) ((x & 0x100) == 0x100)

trunk/psLib/src/types/psTree.c

@@ -14 +14 @@
 #include "psTree.h"
 
+//#define INPUT_CHECK                   // Check inputs for functions that may be in a tight loop?
+
 
 // XXX Upgrades:
@@ -84 +86 @@
 }
 
-psTree *psTreeAlloc(int dim, int maxLeafContents, long numNodes)
+psTree *psTreeAlloc(int dim, int maxLeafContents, psTreeType type, long numNodes)
 {
     psAssert(dim > 0, "Dimensionality (%d) must be positive", dim);
@@ -95 +97 @@
     tree->dim = dim;
     tree->maxLeafContents = maxLeafContents;
+    tree->type = type;
 
     tree->numNodes = numNodes;
@@ -115 +118 @@
 
 
-psTree *psTreePlant(int dim, int maxLeafContents, ...)
+psTree *psTreePlant(int dim, int maxLeafContents, psTreeType type, ...)
 {
     PS_ASSERT_INT_POSITIVE(dim, NULL);
@@ -122 +125 @@
     // Parse coordinate list
     va_list args;                       // Variable argument list
-    va_start(args, maxLeafContents);
+    va_start(args, type);
     psArray *coords = psArrayAlloc(dim); // Array of coordinates
     long numData = 0;                   // Number of data points
@@ -145 +148 @@
         long pow2;                      // Smallest power of two >= numData
         for (pow2 = 1; pow2 < numData; pow2 <<= 1);
-        numNodes = PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1);
-    }
-
-    psTree *tree = psTreeAlloc(dim, maxLeafContents, numNodes);
+        numNodes = PS_MAX(PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1), 1);
+    }
+
+    psTree *tree = psTreeAlloc(dim, maxLeafContents, type, numNodes);
     tree->data = psTreeCoordArrayAlloc(numData, dim);
     tree->numData = numData;
@@ -199 +202 @@
     }
 
+    if (numData <= maxLeafContents) {
+        // Don't need to do any more work
+        return tree;
+    }
+
     psArray *work = psArrayAlloc(numNodes); // Work queue
     work->data[0] = root;
@@ -365 +373 @@
 psTreeNode *psTreeLeaf(const psTree *tree, const psVector *coords)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, NULL);
     PS_ASSERT_VECTOR_NON_NULL(coords, NULL);
@@ -403 +411 @@
 {
     int dim = tree->dim;                // Dimensionality
-    switch (dim) {
-      case 2:
-        return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
-            PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]);
-      case 3:
-        return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
-            PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) +
-            PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]);
-      default: {
-          double distance2 = 0.0;             // Distance of interest
-          for (int i = 0; i < dim; i++) {
-              distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]);
+    switch (tree->type) {
+      case PS_TREE_EUCLIDEAN:
+        switch (dim) {
+          case 2:
+            return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
+                PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]);
+          case 3:
+            return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) +
+                PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) +
+                PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]);
+          default: {
+              double distance2 = 0.0;             // Distance of interest
+              for (int i = 0; i < dim; i++) {
+                  distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]);
+              }
+              return distance2;
           }
-          return distance2;
-      }
-    }
+        }
+        break;
+      case PS_TREE_SPHERICAL:
+        switch (dim) {
+          case 2: {
+              // Haversine formula
+              double dphi = coords->data.F64[1] - tree->data->F64[index][1];
+              double sindphi = sin(dphi / 2.0);
+              double dlambda = coords->data.F64[0] - tree->data->F64[index][0];
+              double sindlambda = sin(dlambda / 2.0);
+              return PS_SQR(sindphi) +
+                  cos(coords->data.F64[1]) * cos(tree->data->F64[index][1]) * PS_SQR(sindlambda);
+          }
+          default:
+            psAbort("Spherical distances not supported for more than 2 dimensions");
+        }
+      default:
+        psAbort("Unrecognised type: %x", tree->type);
+    }
+
     return NAN;
 }
@@ -430 +459 @@
         double minDiff = tree->min->F64[index][i] - coords->data.F64[i];
         if (minDiff > 0) {
-            distance += PS_SQR(minDiff);
+            switch (tree->type) {
+              case PS_TREE_EUCLIDEAN:
+                distance += PS_SQR(minDiff);
+                break;
+              case PS_TREE_SPHERICAL: {
+                  double sinDiff = sin(minDiff / 2.0);
+                  distance += PS_SQR(sinDiff);
+                  break;
+              }
+              default:
+                psAbort("Unrecognised type: %x", tree->type);
+            }
             continue;
         }
         double maxDiff = coords->data.F64[i] - tree->max->F64[index][i];
         if (maxDiff > 0) {
-            distance += PS_SQR(maxDiff);
+            switch (tree->type) {
+              case PS_TREE_EUCLIDEAN:
+                distance += PS_SQR(maxDiff);
+                break;
+              case PS_TREE_SPHERICAL: {
+                  double sinDiff = sin(maxDiff / 2.0);
+                  distance += PS_SQR(sinDiff);
+                  break;
+              }
+              default:
+                psAbort("Unrecognised type: %x", tree->type);
+            }
             continue;
         }
@@ -479 +530 @@
 }
 
-// Return the index of the nearest neighbour to given coordinates, within some radius
+// Return the index of the nearest neighbour to given coordinates, within some distance measure
 // This is the engine for psTreeNearest() and psTreeNearestWithin()
 static inline long treeNearestWithin(const psTree *tree, // Tree
                                      const psVector *coordinates, // Coordinates of interest
-                                     double bestDistance // Distance (radius-squared) to best point
+                                     double bestDistance // Distance measure to best point
     )
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, -1);
     PS_ASSERT_VECTOR_NON_NULL(coordinates, -1);
@@ -545 +596 @@
 
 
+// Convert a radius to our internal "distance measure"
+// Often, we're given a search radius, but for efficiency reasons, we don't use that internally.
+static double treeRadiusToDistance(const psTree *tree, double radius)
+{
+    switch (tree->type) {
+      case PS_TREE_EUCLIDEAN:
+        // Using the square of the distance as the distance measure
+        return PS_SQR(radius);
+      case PS_TREE_SPHERICAL: {
+          // Using a rearrangement of the Haversine formula
+          double sindist = sin(radius / 2.0);
+          return PS_SQR(sindist);
+      }
+      default:
+        psAbort("Unrecognised type: %x", tree->type);
+    }
+}
+
+
 long psTreeNearestWithin(const psTree *tree, const psVector *coords, double radius)
 {
-    return treeNearestWithin(tree, coords, PS_SQR(radius));
-}
-
-
-// Search a leaf node for points within radius squared
-static inline long treeLeafSearchWithin(double radius2, // Radius squared to search
+    return treeNearestWithin(tree, coords, treeRadiusToDistance(tree, radius));
+}
+
+
+// Search a leaf node for points within distance
+static inline long treeLeafSearchWithin(double distance, // Distance to search
                                         const psTree *tree, // Tree of interest
                                         const psTreeNode *leaf, // Leaf to search
@@ -561 +631 @@
     for (int i = 0; i < leaf->num; i++) {
         long index = leaf->contents[i]; // Index of point
-        double distance = treeContentDistance(tree, index, coords); // Distance to point
-        if (distance < radius2) {
+        if (treeContentDistance(tree, index, coords) < distance) {
             num++;
         }
@@ -572 +641 @@
 long psTreeWithin(const psTree *tree, const psVector *coordinates, double radius)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, -1);
     PS_ASSERT_VECTOR_NON_NULL(coordinates, -1);
@@ -581 +650 @@
                         psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
 
-    radius *= radius;                   // We work with the radius-squared
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
     long num = 0;                       // Number of points in circle
 
@@ -588 +657 @@
     // Find the closest point in the leaf that contains the point of interest
     psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
-    num += treeLeafSearchWithin(radius, tree, leaf, coords);
+    num += treeLeafSearchWithin(distance, tree, leaf, coords);
 
     psArray *work = psArrayAlloc(tree->numNodes); // Work queue
@@ -605 +674 @@
             }
             // Leaf node
-            num += treeLeafSearchWithin(radius, tree, node, coords);
+            num += treeLeafSearchWithin(distance, tree, node, coords);
             work->data[workIndex] = NULL;
             workIndex--;
@@ -618 +687 @@
 }
 
-// Search a leaf node for any points within radius squared
-static inline bool treeLeafSearchWithinAny(double radius2, // Radius squared to search
-                                           const psTree *tree, // Tree of interest
-                                           const psTreeNode *leaf, // Leaf to search
-                                           const psVector *coords // Coordinates of interest
+// Search a leaf node for points within distance
+static inline void treeLeafSearchAllWithin(psVector *result,       // Result vector
+                                          double distance, // Distance to search
+                                          const psTree *tree, // Tree of interest
+                                          const psTreeNode *leaf, // Leaf to search
+                                          const psVector *coords // Coordinates of interest
     )
 {
     for (int i = 0; i < leaf->num; i++) {
         long index = leaf->contents[i]; // Index of point
-        double distance = treeContentDistance(tree, index, coords); // Distance to point
-        if (distance < radius2) {
-            return true;
-        }
-    }
-    return false;
-}
-
-// Given an arbitrary point and a matching radius, return whether there are any points within that radius
-bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius)
-{
-#if 1 // Might be in a tight loop
-    PS_ASSERT_TREE_NON_NULL(tree, false);
-    PS_ASSERT_VECTOR_NON_NULL(coordinates, false);
-    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false);
+        if (treeContentDistance(tree, index, coords) < distance) {
+            psVectorAppend(result, index);
+        }
+    }
+    return;
+}
+
+// Given an arbitrary point and a matching radius, return the index of all points within that radius
+psVector *psTreeAllWithin(const psTree *tree, const psVector *coordinates, double radius)
+{
+#ifdef INPUT_CHECK // Might be in a tight loop
+    PS_ASSERT_TREE_NON_NULL(tree, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(coordinates, NULL);
+    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, NULL);
 #endif
 
@@ -647 +716 @@
                         psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
 
-    radius *= radius;                   // We work with the radius-squared
-
-    // This is essentially the same as psTreeWithin, except we can bail as soon as we find something
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
+
+    psVector *result = psVectorAllocEmpty(4, PS_TYPE_S64); // Indices of points within match radius
+
+    // This is essentially the same as psTreeNearest, except we're not allowed to prune
 
     // Find the closest point in the leaf that contains the point of interest
     psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
-    if (treeLeafSearchWithinAny(radius, tree, leaf, coords)) {
-        psFree(coords);
-        return true;
-    }
+    treeLeafSearchAllWithin(result, distance, tree, leaf, coords);
 
     psArray *work = psArrayAlloc(tree->numNodes); // Work queue
@@ -673 +741 @@
             }
             // Leaf node
-            if (treeLeafSearchWithinAny(radius, tree, node, coords)) {
+            treeLeafSearchAllWithin(result, distance, tree, node, coords);
+            work->data[workIndex] = NULL;
+            workIndex--;
+        }
+
+        leaf = up;
+    }
+    psFree(work);
+    psFree(coords);
+
+    return result;
+}
+
+// Search a leaf node for any points within distance measure
+static inline bool treeLeafSearchWithinAny(double distance, // Distance to search
+                                           const psTree *tree, // Tree of interest
+                                           const psTreeNode *leaf, // Leaf to search
+                                           const psVector *coords // Coordinates of interest
+    )
+{
+    for (int i = 0; i < leaf->num; i++) {
+        long index = leaf->contents[i]; // Index of point
+        if (treeContentDistance(tree, index, coords) < distance) {
+            return true;
+        }
+    }
+    return false;
+}
+
+// Given an arbitrary point and a matching radius, return whether there are any points within that radius
+bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius)
+{
+#ifdef INPUT_CHECK // Might be in a tight loop
+    PS_ASSERT_TREE_NON_NULL(tree, false);
+    PS_ASSERT_VECTOR_NON_NULL(coordinates, false);
+    PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false);
+#endif
+
+    psVector *coords = (coordinates->type.type == PS_TYPE_F64 ? psMemIncrRefCounter((psVector*)coordinates) :
+                        psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates
+
+    double distance = treeRadiusToDistance(tree, radius); // Distance measure
+
+    // This is essentially the same as psTreeWithin, except we can bail as soon as we find something
+
+    // Find the closest point in the leaf that contains the point of interest
+    psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest
+    if (treeLeafSearchWithinAny(distance, tree, leaf, coords)) {
+        psFree(coords);
+        return true;
+    }
+
+    psArray *work = psArrayAlloc(tree->numNodes); // Work queue
+    while (leaf->up) {
+        psTreeNode *up = leaf->up;      // Parent node
+
+        long workIndex = 0;
+        work->data[workIndex] = (leaf == up->left ? up->right : up->left); // The other node
+        while (workIndex >= 0) {
+            psTreeNode *node = work->data[workIndex];
+            if (node->left) {
+                // Branch node
+                work->data[workIndex] = node->right;
+                work->data[++workIndex] = node->left;
+                continue;
+            }
+            // Leaf node
+            if (treeLeafSearchWithinAny(distance, tree, node, coords)) {
                 // Clear out the work queue
                 memset(work->data, 0, (workIndex + 1) * sizeof(void));
@@ -695 +830 @@
 psVector *psTreeCoords(psVector *out, const psTree *tree, long index)
 {
-#if 1 // Might be in a tight loop
+#ifdef INPUT_CHECK // Might be in a tight loop
     PS_ASSERT_TREE_NON_NULL(tree, NULL);
     PS_ASSERT_INT_NONNEGATIVE(index, NULL);

trunk/psLib/src/types/psTree.h

@@ -16 +16 @@
 } psTreeCoordArray;
 
+/// Type of tree
+///
+/// Specifies how distances are measured
+typedef enum {
+    PS_TREE_EUCLIDEAN,                  // d^2 = dx^2 + dy^2 + ...
+    PS_TREE_SPHERICAL,                  // sin(dist/2)^2 = sin(dphi/2)^2 + cos(phi1)cos(phi2)sin(dlambda/2)^2
+} psTreeType;
+
 /// A simple kd-tree implementation
 ///
@@ -23 +31 @@
     int dim;                            ///< Dimensionality
     int maxLeafContents;                ///< Maximum number of points on a leaf
+    psTreeType type;                    ///< Type of tree
     long numNodes;                      ///< Number of nodes
     long numData;                       ///< Number of data points
@@ -67 +76 @@
 psTree *psTreeAlloc(int dim,            ///< Dimensionality
                     int maxLeafContents,///< Maximum number of points on a leaf
+                    psTreeType type,    ///< Type of tree
                     long numNodes       ///< Number of nodes in tree
                     );
@@ -75 +85 @@
 psTree *psTreePlant(int dim,            ///< Dimensionality
                     int maxLeafContents,///< Maximum number of points on a leaf
+                    psTreeType type,    ///< Type of tree
                     ...                 ///< psVector for each coordinate
                     );
@@ -108 +119 @@
                      );
 
+/// Return the index of all points within some radius of given coordinates
+psVector *psTreeAllWithin(const psTree *tree,          ///< Tree
+                          const psVector *coordinates, ///< Coordinates of interest
+                          double radius                ///< Radius of interest
+                          );
+
 /// Return the coordinates of a point in the tree, specified by its index
 psVector *psTreeCoords(psVector *out,   ///< Output vector, or NULL

trunk/psLib/test/optime

@@ -2 +2 @@
 Makefile.in
 .deps
+tap_psStatsTiming

@@ -2 +2 @@
 Makefile.in
 .deps
+tap_psStatsTiming

trunk/psLib/test/types/tap_psTree.c

@@ -3 +3 @@
 #include "pstap.h"
 
-#define NUM 10000                       // Number of points
+#define NUM 1000000                      // Number of points
+#define SPHERICAL_DISTANCE 3.0          // Distance of interest for spherical test
 
 int main(int argc, char *argv[])
 {
     psLibInit(NULL);
-    plan_tests(6);
+    plan_tests(13);
 
+    // Euclidean geometry: 6 tests
     {
         psMemId id = psMemGetId();
@@ -23 +25 @@
         psFree(rng);
 
-        psTree *tree = psTreePlant(2, 2, x, y);
+        psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y);
 
         ok(tree, "Tree planted");
@@ -35 +37 @@
             long closeIndex = psTreeNearest(tree, coords);
             psFree(coords);
-            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found point: %ld", closeIndex);
+            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex);
 
             long bestIndex = -1;
@@ -68 +70 @@
     }
 
+    // Spherical geometry: 7 tests
+    {
+        psMemId id = psMemGetId();
+
+        psVector *ra = psVectorAlloc(NUM, PS_TYPE_F64);
+        psVector *dec = psVectorAlloc(NUM, PS_TYPE_F64);
+
+        psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 0);
+        for (int i = 0; i < NUM; i++) {
+            // Using http://mathworld.wolfram.com/SpherePointPicking.html
+            ra->data.F64[i] = psRandomUniform(rng);
+            dec->data.F64[i] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
+        }
+
+        psTree *tree = psTreePlant(2, 2, PS_TREE_SPHERICAL, ra, dec);
+
+        ok(tree, "Tree planted");
+        skip_start(!tree, 4, "tree died");
+        {
+            //            psTreePrint(stderr, tree);
+
+            psVector *coords = psVectorAlloc(2, PS_TYPE_F64);
+            coords->data.F64[0] = psRandomUniform(rng);
+            coords->data.F64[1] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
+
+            psVector *indices = psTreeAllWithin(tree, coords, DEG_TO_RAD(SPHERICAL_DISTANCE));
+            ok(indices && indices->type.type == PS_TYPE_S64, "got list of indices (%ld points)", indices->n);
+            long closeIndex = psTreeNearest(tree, coords);
+            ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex);
+
+            ok(psVectorSortInPlace(indices), "sorted indices");
+
+            bool allgood = true;        // All points in the appropriate place?
+            double bestDistance = INFINITY; // Distance to best point
+            long bestIndex = -1;        // Index of best point
+            for (long i = 0, j = 0; i < NUM; i++) {
+#if 1
+                // Traditional formula
+                double dist = acos(sin(coords->data.F64[1]) * sin(dec->data.F64[i]) +
+                                   cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
+                                   cos(coords->data.F64[0] - ra->data.F64[i]));
+#else
+                // Haversine formula: used in psTree
+                double dphi = coords->data.F64[1] - dec->data.F64[i];
+                double sindphi = sin(dphi/2.0);
+                double dlambda = coords->data.F64[0] - ra->data.F64[i];
+                double sindlambda = sin(dlambda/2.0);
+                double dist = 2.0 * asin(sqrt(PS_SQR(sindphi) +
+                                              cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
+                                              PS_SQR(sindlambda)));
+#endif
+                                              if (dist < bestDistance) {
+                    bestDistance = dist;
+                    bestIndex = i;
+                }
+                if (i == indices->data.S64[j]) {
+                    j++;
+                    if (dist > DEG_TO_RAD(SPHERICAL_DISTANCE)) {
+                        diag("%ld is in the list, but shouldn't be (%lf vs %lf)",
+                             i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                        allgood = false;
+                    }
+                } else if (dist <= DEG_TO_RAD(SPHERICAL_DISTANCE)) {
+                    diag("%ld is not in the list, but should be (%lf vs %lf)",
+                         i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                    allgood = false;
+                }
+            }
+            ok(allgood, "list is acurate");
+            ok(bestIndex == closeIndex, "correct point: %ld vs %ld", closeIndex, bestIndex);
+
+            psFree(coords);
+            psFree(indices);
+
+        }
+        skip_end();
+
+        psFree(rng);
+        psFree(tree);
+        psFree(ra);
+        psFree(dec);
+
+        ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks");
+    }
+
     psLibFinalize();
 }

trunk/psModules/src/objects/pmSourceMatch.c

@@ -221 +221 @@
         } else {
             // Match with the master list
-            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
+            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree
             long numMatch = 0;          // Number of matches
 

@@ -2 +2 @@
 Makefile.in
 .deps
+tap_psStatsTiming