@acbecker can you please edit this issue description with a link to the code for function we need to import?

Even better might be a tiny C++ standalone source code that instantiates the class and has reproducible input/output.

from kbmod.

Here are the places I found:

• skyToPixel: does what we need, calling wcss2p
• wcss2p: this probably does the bulk of the math, haven't found the source code for this yet

from kbmod.

Wcss2p can be found at ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib.tar.bz2 and then inside wcslib-4.24/C/wcs.c

from kbmod.

And the use of the wcss2p code comes from here:

https://dev.lsstcorp.org/cgit/LSST/DMS/afw.git/tree/src/image/TanWcs.cc#n272

I will try and figure out the minimal amount of information to copy the contents of a Wcs instance, that will allow us to instantiate a new one in a UDF.

from kbmod.

At the risk of ticket clutter (vs. having half our discussions on email and half on this ticket) I should have a ~1% database available for ingest by early next week. Would be useful to understand if we want to move to EC2 by then.

from kbmod.

Also, for ref, from the comments in the TanWcs.h definition file:

/**
  *  @brief Implementation of the WCS standard for the special case of the Gnomonic
  *  (tangent plane) projection.
  * 
  *  This class treats the special case of tangent plane projection. It extends the Wcs standard by
  *  optionally accounting for distortion in the image plane using the Simple Imaging Polynomial (SIP)
  *  convention.
  *  This convention is described in Shupe et al. (2005) (Astronomical Data Analysis Software and Systems
  *  XIV, Asp Conf. Series Vol XXX, Ed: Shopbell et al.), and descibed in some more detail in
  *  http://web.ipac.caltech.edu/staff/fmasci/home/wise/codeVdist.html
  * 
  *  To convert from pixel coordintates to radec ("intermediate world coordinates"), first use the matrices
  *  _sipA and _sipB to calculate undistorted coorinates (i.e where on the chip the image would lie if
  *  the optics gave undistorted images), then pass these undistorted coorinates wcsp2s() to calculate radec.
  * 
  *  For the reverse, radec -> pixels, convert the radec to undistorted coords, and then use the _sipAp and
  *  _sipBp matrices to add in the distortion terms.

from kbmod.

And with a little bit of gymnastics I can extract the below information from an instance. It looks like we'll want to store the AP_* and BP_0* terms (42 columns) as well as the 8 nominal terms CRPIX,CRVAL,CD_*. So 50 columns of double precision floats, and ~1e6 rows.

 NAXIS = 2
 EQUINOX = 2000.0000000000
 RADESYS = "FK5"
 CRPIX1 = 0.99999937976196
 CRPIX2 = 0.99999999917061
 CD1_1 = -6.0998100447528e-09
 CD1_2 = 0.00010999609605758
 CD2_1 = 0.00010997941628578
 CD2_2 = -1.6012342099554e-08
 CRVAL1 = 318.61941480048
 CRVAL2 = -0.21434393244613
 CUNIT1 = "deg"
 CUNIT2 = "deg"
 A_ORDER = 4
 A_0_2 = -3.5896369352650e-09
 A_0_3 = -9.8193522641534e-17
 A_0_4 = -3.2395053794349e-21
 A_1_1 = 6.0811031828094e-12
 A_1_2 = 1.8423040866085e-12
 A_1_3 = 1.0999158271664e-20
 A_2_0 = -5.6483937909758e-07
 A_2_1 = -2.3067576899243e-15
 A_2_2 = -2.4761199575771e-19
 A_3_0 = 1.8848400025438e-10
 A_3_1 = 6.9195454224547e-19
 A_4_0 = 5.1600503906739e-19
 B_ORDER = 4
 B_0_2 = 2.4544356030624e-12
 B_0_3 = 1.2285110191299e-12
 B_0_4 = 8.8596091343505e-21
 B_1_1 = 7.1809238938035e-09
 B_1_2 = -1.0989603613886e-16
 B_1_3 = 2.2394496582017e-20
 B_2_0 = -7.1087610143063e-09
 B_2_1 = -4.0830284671424e-15
 B_2_2 = -1.2863491244552e-20
 B_3_0 = 1.5702746516556e-12
 B_3_1 = 1.3655274519653e-18
 B_4_0 = 1.0827786190003e-20
 AP_ORDER = 5
 AP_0_0 = 3.7250546440694e-09
 AP_0_1 = -1.1412528166471e-11
 AP_0_2 = 3.5896513217540e-09
 AP_0_3 = 9.0185436003823e-17
 AP_0_4 = -1.0282141088766e-20
 AP_0_5 = -8.0932044215343e-27
 AP_1_0 = 1.7519653075116e-11
 AP_1_1 = -6.0500454000118e-12
 AP_1_2 = -1.8383256630685e-12
 AP_1_3 = -1.6702922974262e-21
 AP_1_4 = 6.9145914248588e-24
 AP_2_0 = 5.6483904480165e-07
 AP_2_1 = 2.3253193165425e-15
 AP_2_2 = -4.8640124822731e-18
 AP_2_3 = -2.7269329555650e-24
 AP_3_0 = -1.8784469840417e-10
 AP_3_1 = -7.2093132938818e-19
 AP_3_2 = 1.3862327754001e-21
 AP_4_0 = -5.3384586880658e-16
 AP_4_1 = 6.9598270418341e-24
 AP_5_0 = 1.0680958500467e-19
 BP_ORDER = 5
 BP_0_0 = 3.8413237939587e-11
 BP_0_1 = -1.1597618127228e-13
 BP_0_2 = -2.4542882038954e-12
 BP_0_3 = -1.2285368791464e-12
 BP_0_4 = -8.8241402845071e-21
 BP_0_5 = 4.5276129853359e-24
 BP_1_0 = 2.3648654999189e-13
 BP_1_1 = -7.1809235913996e-09
 BP_1_2 = 1.6078142232942e-16
 BP_1_3 = 2.6217226073356e-20
 BP_1_4 = -1.1550193291221e-26
 BP_2_0 = 7.1087571371124e-09
 BP_2_1 = 7.8220891526374e-17
 BP_2_2 = -5.6267493315872e-20
 BP_2_3 = -3.1149801367721e-26
 BP_3_0 = -1.5622815651210e-12
 BP_3_1 = -1.1995802346489e-20
 BP_3_2 = 1.4423882402653e-23
 BP_4_0 = -5.3493477312959e-18
 BP_4_1 = 2.4359310266642e-26
 BP_5_0 = 8.8947447722507e-22
 CTYPE1 = "RA---TAN-SIP"
 CTYPE2 = "DEC--TAN-SIP"

from kbmod.

OK, so I think I have the structure of a UDF put together in Python. Have a look at

https://github.com/uwescience/kbmod/blob/master/python/parseWcs.py

Everything that I access via a md.get() will instead come from a column in the database. We can run everything in our own methods except the following code block, which we'll have to rewrite to use wcslib instead of pywcs and then use the equivalent internal call to "s2p":

import pywcs 
prm = pywcs.Wcsprm()
prm.crpix   = np.array((md.get("CRPIX1"), md.get("CRPIX2")))
prm.crval   = np.array((md.get("CRVAL1"), md.get("CRVAL2")))
prm.cunit   = ["deg", "deg"]
prm.ctype   = ["RA---TAN", "DEC--TAN"]
prm.equinox = md.get("EQUINOX")
prm.radesys = md.get("RADESYS")
prm.cd      = np.array(((md.get("CD1_1"), md.get("CD1_2")), (md.get("CD2_1"), md.get("CD2_2"))))
xlin, ylin  = prm.s2p(((raref,decref),), origin=0)["pixcrd"][0] 

Once we get the linear distortion run as above, we have to apply additional non-linear terms as in distortPixel() which are entirely straightforward.

from kbmod.

I have an initial UDF that I am working on. However, I am having a difficult time debugging things because I cannot DROP the FUNCTION within a session, and re-CREATE it with a debugged version. I have found that if I DROP the function, log out, and then log in, my CREATE statement picks up on the newly compiled shared library. But if I don't log out, I get the old stale library.

So something is not committing, or flushing, or whatnot (COMMIT does not work since there is no on-going transaction in progress). Any ideas?

Use case below:

kbmod=#  CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
     AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
     LANGUAGE C STRICT;
kbmod=# SELECT hello( name ) FROM test;
    hello     
---------------
 Hello, Xavier
 Hello, Yari
 Hello, Zack
### Modify the function to print "Hello2", rebuild the shared library.  DROP function just to be explicit
kbmod=# DROP FUNCTION hello(TEXT);
kbmod=#  CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
     AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
     LANGUAGE C STRICT;
kbmod=# SELECT hello( name ) FROM test;
    hello     
---------------
 Hello, Xavier
 Hello, Yari
 Hello, Zack
#### But if I \q and then log back in
kbmod=#  CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
     AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
     LANGUAGE C STRICT;
kbmod=# SELECT hello( name ) FROM test;
     hello      
----------------
 Hello2, Xavier
 Hello2, Yari
 Hello2, Zack

from kbmod.

A bit more info:

http://www.postgresql.org/docs/8.1/static/sql-createfunction.html

"When repeated CREATE FUNCTION calls refer to the same object file, the file is only loaded once. To unload and reload the file (perhaps during development), use the LOAD command."

But LOAD does not fix anything:

kbmod=# SELECT hello( name ) FROM test;
     hello      
----------------
 Hello4, Xavier
 Hello4, Yari
 Hello4, Zack
(3 rows)
kbmod=# LOAD '/Users/acbecker/src/github/kbmod/src/hello.so';
LOAD
kbmod=# SELECT hello( name ) FROM test;
     hello      
----------------
 Hello4, Xavier
 Hello4, Yari
 Hello4, Zack
(3 rows)
kbmod=# \q
acbecker7:/opt/local/lib/postgresql93/bin/psql -U postgres -d kbmod
psql (9.3.5)
Type "help" for help.
kbmod=# CREATE TEMP TABLE test( name ) AS VALUES ('Xavier'), ('Yari'), ('Zack');
SELECT 3
kbmod=# SELECT hello( name ) FROM test;
     hello      
----------------
 Hello5, Xavier
 Hello5, Yari
 Hello5, Zack
(3 rows)

from kbmod.

Andy,

Two questions:

• If you drop the function (without replacing it), can you still invoke it?
• What happens if you try using a different file name: for example,
  compiling to hello2.so?

On Fri, Oct 31, 2014 at 1:40 PM, Andy Becker [email protected]
wrote:

I have an initial UDF that I am working on. However, I am having a
difficult time debugging things because I cannot DROP the FUNCTION within a
session, and re-CREATE it with a debugged version. I have found that if I
DROP the function, log out, and then log in, my CREATE statement picks up
on the newly compiled shared library. But if I don't log out, I get the old
stale library.

So something is not committing, or flushing, or whatnot (COMMIT does not
work since there is no on-going transaction in progress). Any ideas?

Use case below:

kbmod=# CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
LANGUAGE C STRICT;

kbmod=# SELECT hello( name ) FROM test;

hello

Hello, Xavier
Hello, Yari
Hello, Zack

Modify the function to print "Hello2", rebuild the shared library. DROP function just to be explicit

kbmod=# DROP FUNCTION hello(TEXT);

kbmod=# CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
LANGUAGE C STRICT;

kbmod=# SELECT hello( name ) FROM test;

hello

Hello, Xavier
Hello, Yari
Hello, Zack

But if I \q and then log back in

kbmod=# CREATE OR REPLACE FUNCTION hello(TEXT) RETURNS TEXT
AS '/Users/acbecker/src/github/kbmod/src/hello.so', 'hello'
LANGUAGE C STRICT;

kbmod=# SELECT hello( name ) FROM test;

hello

Hello2, Xavier
Hello2, Yari
Hello2, Zack

—
Reply to this email directly or view it on GitHub
#1 (comment).

from kbmod.

kbmod=# DROP FUNCTION hello(TEXT);
DROP FUNCTION
kbmod=# SELECT hello( name ) FROM test;
ERROR:  function hello(text) does not exist
LINE 1: SELECT hello( name ) FROM test;
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I can save it is hello2.so and load it (once) but changes to the shared library do not get reloaded

from kbmod.

This is not what you want, but my workflow is to put all the code in a
script and re-run the script on each test as opposed to using the client
directly. At the top of the script, I drop and recreate the functions,
views, etc. every time.

This also allows you to use a real editor (say, with syntax highlighting)
as opposed to the psql line-oriented interpreter editor.

On Fri, Oct 31, 2014 at 1:55 PM, Andy Becker [email protected]
wrote:

kbmod=# DROP FUNCTION hello(TEXT);
DROP FUNCTION

kbmod=# SELECT hello( name ) FROM test;
ERROR: function hello(text) does not exist
LINE 1: SELECT hello( name ) FROM test;
^
HINT: No function matches the given name and argument types. You might need to add explicit type casts.

I can save it is hello2.so and load it (once) but changes to the shared
library do not get reloaded

—
Reply to this email directly or view it on GitHub
#1 (comment).

from kbmod.

Mailed the pgsql-general mailing list and got the following response back:

Hi - I seem to be unable to reLOAD a shared library within the session
that I LOADed it.

Nope, you can't, there's no such functionality.

Hints as to what is going wrong here? I would certainly expect to be
able to re-load a shared library while debugging my UDF.

That would require being able to unload it, which is an operation fraught
with hazards. We used to allow that, but gave it up after observing that
practically every extant extension could be made to crash on unload.
There is for instance no safe way to get out of a function hook --- the
code pattern you may have seen of restoring the prior value is wrong and
unsafe, because it doesn't account for some other extension having plugged
into the hook after you. You'd leave that other extension kneecapped,
with some of its hook callbacks disabled but others perhaps not.

from kbmod.

I have checked in and (somewhat) tested a UDF. For a whopping single image, and the mapping of all of its 4 corners from sky->pixel, I get the following comparison between the "original" mapping that comes from the SDSS data and the wcs table + UDF mapping that we now have available to us. It all seems pretty consistent. We can chat about closing out this issue tomorrow.

LLC

orig : 3.63807384218e-09 1.72617435013e-07
UDF : 0.000000 -0.000000

LRC

orig : 1999.99388942 2.49034561473e-05
UDF : 2000.000001 -0.000001

URC

orig : 1999.99388912 1000.00002493
UDF : 2000.000001 999.999999

ULC

orig : 2.13230300083e-09 1000.00000019
UDF : 0.000000 1000.000000

from kbmod.

Very cool!

On Mon, Nov 3, 2014 at 9:19 PM, Andy Becker [email protected]
wrote:

I have checked in and (somewhat) tested a UDF. For a whopping single
image, and the mapping of all of its 4 corners from sky->pixel, I get the
following comparison between the "original" mapping that comes from the
SDSS data and the wcs table + UDF mapping that we now have available to us.
It all seems pretty consistent. We can chat about closing out this issue
tomorrow.

LLC

orig : 3.63807384218e-09 1.72617435013e-07
UDF : 0.000000 -0.000000

LRC

orig : 1999.99388942 2.49034561473e-05
UDF : 2000.000001 -0.000001

URC

orig : 1999.99388912 1000.00002493
UDF : 2000.000001 999.999999

ULC

orig : 2.13230300083e-09 1000.00000019
UDF : 0.000000 1000.000000

—
Reply to this email directly or view it on GitHub
#1 (comment).

Bill Howe
Associate Director and Senior Data Science Fellow, UW eScience Institute
Affiliate Faculty, Computer Science & Engineering
University of Washington
To acknowledge eScience: "Supported in part by the University of Washington
eScience Institute"

from kbmod.

For reference, to build the UDF on OS X:

rm skyToIdx.o skyToIdx.so ; cc -fpic -c skyToIdx.c -I/opt/local/include -I/opt/local/include/postgresql93/server/ ; cc -bundle -flat_namespace -undefined suppress -lwcs -L/opt/local/lib -o skyToIdx.so skyToIdx.o

And then to install it in a pgsql session:

kbmod=# CREATE OR REPLACE FUNCTION c_skyToIdx(wcs, double precision, double precision) RETURNS integer
     AS '/Users/acbecker/src/github/kbmod/src/skyToIdx', 'c_skyToIdx'
     LANGUAGE C STRICT;

And to call it:

select c_skyToIdx(wcs, -41.0748270035, 0.833286131146) as idx from wcs;

And to interpret the results, the returned integer is the x,y coordinate of the pixel stored using the mapping (x + 2048*y) where 2048 is the width of the image in x. We have this value in the pixels table as column "pidx", and an index built on it, so we should be able to start testing this all tomorrow.

from kbmod.

I was able to build the UDF on the EC2 instance[*] and have run the command at scale. The nice thing is that is appears to work! And its speed things up by a factor of x5-7 over the next best query option. I'm going to close this issue and work towards optimizing the query. FWIW here is the EXPLAIN:

explain SELECT p.pixelId, p.pidx, p.ra, p.decl, p.fval
FROM
  ST_SetSRID(ST_MakePoint(-41.1, 0.9),3786) as traj,
  pixels as p,
  fields as f,
  wcs as w,
  c_skyToIdx(w, -41.1, 0.9) as idx
WHERE
  TIMESTAMP WITH TIME ZONE '1999-10-14 03:49:01.772609z' <@ f.trange
AND
  ST_INTERSECTS(traj, f.bbox)
AND
 f.fieldId = w.fieldId
AND
 f.fieldId = p.fieldId
AND
 p.pidx = idx;


                                                 QUERY PLAN                                                  
-------------------------------------------------------------------------------------------------------------
 Nested Loop  (cost=25141.23..3897860.81 rows=38 width=32)
   Join Filter: (f.fieldid = p.fieldid)
   ->  Nested Loop  (cost=0.57..16.91 rows=1 width=20)
         ->  Nested Loop  (cost=0.57..16.89 rows=1 width=452)
               ->  Nested Loop  (cost=0.28..8.57 rows=1 width=8)
                     Join Filter: ((traj.traj && f.bbox) AND _st_intersects(traj.traj, f.bbox))
                     ->  Function Scan on traj  (cost=0.00..0.01 rows=1 width=32)
                     ->  Index Scan using trangeidx on fields f  (cost=0.28..8.29 rows=1 width=128)
                           Index Cond: ('1999-10-14 03:49:01.772609+00'::timestamp with time zone <@ trange)
               ->  Index Scan using fieldwcsidx on wcs w  (cost=0.29..8.31 rows=1 width=444)
                     Index Cond: (fieldid = f.fieldid)
         ->  Function Scan on c_skytoidx idx  (cost=0.00..0.01 rows=1 width=4)
   ->  Bitmap Heap Scan on pixels p  (cost=25140.66..3881056.66 rows=1342979 width=40)
         Recheck Cond: (pidx = idx.idx)
         ->  Bitmap Index Scan on pidxidx  (cost=0.00..24804.91 rows=1342979 width=0)
               Index Cond: (pidx = idx.idx)

[*] : Note on building this. It turns out I had to "yum reinstall" the postgres packages because when I first "yum install"ed them I had not added /usr/local/lib to /etc/ld.so.conf and re-ran ldconfig. The failure mode was that postgres complained:

 libwcs.so.4: cannot open shared object file: No such file or directory

which is in /usr/local/lib

from kbmod.