diff --git a/docs/configuration.md b/docs/configuration.md index 160d355..69a60ff 100644 --- a/docs/configuration.md +++ b/docs/configuration.md @@ -22,14 +22,14 @@ lens.x = unif 48 52 lens.y = unif 48 52 lens.r = unif 15 20 lens.q = unif 0 1 -lens.pa = unif 0 180 -source.x = unif 45 55 -source.y = unif 45 55 +lens.pa = wrap unif 0 180 +source.x = image unif 80 85 +source.y = image unif 35 40 source.r = unif 1 5 source.mag = unif -5 0 source.n = unif 0.5 8.0 source.q = unif 0 1 -source.pa = unif 0 180 +source.pa = wrap unif 0 180 [labels] lens.x = x_L @@ -148,7 +148,7 @@ file. They are given in the format ```ini [priors] -obj.param = [wrap] ... +obj.param = [wrap] [image] ... ``` An exception is the pseudo-prior that fixes the value of a parameter, which is @@ -162,11 +162,21 @@ of 170 ± 20 degrees falls out of the natural [0, 180] degree bounds. The following priors are known: Prior | Description ----------------|---------------------------------------------------------------- +---------------|-------------------------------------------------------------- `` | pseudo-prior that fixes the parameter to the given *value* `unif ` | uniform prior on the interval [*a*, *b*] `norm ` | normal prior with mean *m* and standard deviation *s* +Additionally, the following optional keywords can be specified: + +Keyword | Description +---------------|-------------------------------------------------------------- +`wrap` | prior is cyclic and wraps around at the boundaries +`image` | prior is specified on the image plane + +The `image` keyword can be used to specify a prior on the observed quantities +on the image plane. This is limited to position parameters. + Example: ```ini @@ -183,6 +193,10 @@ lens.r = unif 10 20 ; the orientation angle is between 0 and 180 degrees and wraps around lens.pa = wrap unif 0 180 + +; the source position is given by one of its images +source.x = image unif 80 85 +source.y = image unif 35 40 ``` diff --git a/examples/full_mock_nopsf.ini b/examples/full_mock_nopsf.ini index 67a143f..4b7734d 100644 --- a/examples/full_mock_nopsf.ini +++ b/examples/full_mock_nopsf.ini @@ -38,8 +38,8 @@ lens.y = unif 48 52 lens.r = unif 15 20 lens.q = unif 0 1 lens.pa = wrap unif 0 180 -source.x = unif 45 55 -source.y = unif 45 55 +source.x = image unif 34 36 +source.y = image unif 63 65 source.r = unif 1 5 source.mag = unif -5 0 source.n = unif 0.5 8.0 diff --git a/examples/full_mock_psf.ini b/examples/full_mock_psf.ini index 93d4336..dda5758 100644 --- a/examples/full_mock_psf.ini +++ b/examples/full_mock_psf.ini @@ -40,8 +40,8 @@ lens.y = unif 48 52 lens.r = unif 15 20 lens.q = unif 0 1 lens.pa = wrap unif 0 180 -source.x = unif 45 55 -source.y = unif 45 55 +source.x = image unif 34 36 +source.y = image unif 63 65 source.r = unif 1 5 source.mag = unif -5 0 source.n = unif 0.5 8.0 diff --git a/examples/test_sersic_bulge.ini b/examples/test_sersic_bulge.ini index c03f499..33549ba 100644 --- a/examples/test_sersic_bulge.ini +++ b/examples/test_sersic_bulge.ini @@ -16,20 +16,20 @@ lens = sie source = sersic [priors] -host.x = unif 55 65 -host.y = unif 55 65 -host.r = unif 1 50 -host.mag = unif -10 0 -host.n = 4. -host.q = unif 0.5 1 -host.pa = wrap unif 0 180 +host.x = unif 55 65 +host.y = unif 55 65 +host.r = unif 1 50 +host.mag = unif -10 0 +host.n = 4. +host.q = unif 0.5 1 +host.pa = wrap unif 0 180 lens.x = unif 55 65 lens.y = unif 55 65 lens.r = unif 20 35 lens.q = unif 0.1 1 lens.pa = wrap unif 0 180 -source.x = unif 40 80 -source.y = unif 40 80 +source.x = image unif 80 85 +source.y = image unif 35 40 source.r = unif 0.1 20 source.mag = unif -5 5 source.n = unif 0.5 8.0 @@ -37,13 +37,13 @@ source.q = unif 0.1 1 source.pa = wrap unif 0 180 [labels] -host.x = x_H -host.y = y_H -host.r = r_H -host.mag = mag_H -host.n = n_H -host.q = f_H -host.pa = \theta_H +host.x = x_H +host.y = y_H +host.r = r_H +host.mag = mag_H +host.n = n_H +host.q = f_H +host.pa = \theta_H lens.x = x_L lens.y = y_L lens.r = r_L diff --git a/src/input.c b/src/input.c index 6fcb903..04a4f47 100644 --- a/src/input.c +++ b/src/input.c @@ -275,7 +275,7 @@ void print_input(const input* inp) print_prior(inp->objs[i].pars[j].pri, buf, 99); // collect tags - snprintf(tag, 100, " [%s%s%s%s%s%s%s%s", + snprintf(tag, 100, " [%s%s%s%s%s%s%s%s%s", inp->objs[i].pars[j].type == PAR_POSITION_X ? "position x, " : "", inp->objs[i].pars[j].type == PAR_POSITION_Y ? @@ -291,7 +291,9 @@ void print_input(const input* inp) inp->objs[i].pars[j].bounded ? "bounded, " : "", inp->objs[i].pars[j].wrap ? - "wrap, " : "" + "wrap, " : "", + inp->objs[i].pars[j].ipp ? + "IPP, " : "" ); // check if tags were set diff --git a/src/input.h b/src/input.h index 6d8d3b5..ef4cd3d 100644 --- a/src/input.h +++ b/src/input.h @@ -85,6 +85,9 @@ typedef struct // flag for wrap-around parameters int wrap; + // flag for image plane priors + int ipp; + // label, used for output const char* label; } param; diff --git a/src/input/objects.c b/src/input/objects.c index bf6e453..48fcdca 100644 --- a/src/input/objects.c +++ b/src/input/objects.c @@ -219,6 +219,7 @@ void add_object(input* inp, const char* id, const char* name) obj->pars[i].upper = param_bounds[i].s[1]; obj->pars[i].pri = NULL; obj->pars[i].wrap = 0; + obj->pars[i].ipp = 0; obj->pars[i].label = NULL; } @@ -316,16 +317,21 @@ void set_param_prior(param* par, const char* str) // free existing prior free_prior(par->pri); - // get length of first word - len = strcspn(str, WS); - - // check for wrap keyword - if(strncmp(str, "wrap", len) == 0) + // parse possible keywords + while(1) { - // set wrap flag - par->wrap = 1; + // get length of first word + len = strcspn(str, WS); + + // check for keywords or break loop + if(strncmp(str, "wrap", len) == 0) + par->wrap = 1; + else if(strncmp(str, "image", len) == 0) + par->ipp = 1; + else + break; - // skip word + // skip word and whitespace str += len; str += strspn(str, WS); } diff --git a/src/kernel.c b/src/kernel.c index 6a116a6..faba6c1 100644 --- a/src/kernel.c +++ b/src/kernel.c @@ -114,6 +114,12 @@ static const char SETPHEAD[] = "kernel void set_params(ulong dsiz, global int* gdata, local int* ldata,\n" " constant float* params)\n" "{\n" + " // image plane priors\n" + " float2 x;\n" + " float2 a;\n" + " x = 0;\n" + " a = 0;\n" + " \n" " // load parameters to local memory\n" " for(size_t i = get_local_id(0); i < dsiz; i += get_local_size(0))\n" " ldata[i] = gdata[i];\n" @@ -121,6 +127,7 @@ static const char SETPHEAD[] = ; static const char SETPLEFT[] = " set_%s((local void*)(ldata + %zu)"; static const char SETPARGS[] = ", params[%zu]"; +static const char SETPIPPA[] = ", %s"; static const char SETPRGHT[] = ");\n"; static const char SETPFOOT[] = " \n" @@ -130,6 +137,16 @@ static const char SETPFOOT[] = " \n" "}\n" ; +static const char SETPIPP_POSINIT[] = + " x = (float2)(params[%zu], params[%zu]);\n" +; +static const char SETPIPP_POSLENS[] = + " a += deflection_%s((local void*)(ldata + %zu), x);\n" +; +static const char SETPIPP_POSDEFL[] = + " x -= dot(a, a) < HUGE_VALF ? a : (float2)(1E10f, 1E10f);\n" + " a = 0;\n" +; // object kernel static const char OBJHEAD[] = @@ -371,8 +388,14 @@ static const char* compute_kernel(size_t nobjs, object objs[]) static const char* set_params_kernel(size_t nobjs, object objs[]) { + // trigger for changing lens planes + int trigger; + + // index of object starting current plane + size_t plane; + // buffer for kernel - size_t buf_size; + size_t siz, len; char* buf; // current output position @@ -387,88 +410,235 @@ static const char* set_params_kernel(size_t nobjs, object objs[]) // parameter offset size_t p; - // calculate buffer size - d = p = 0; - buf_size = sizeof(FILEHEAD) + strlen("set_params"); - buf_size += sizeof(SETPHEAD); - for(size_t i = 0; i < nobjs; ++i) - { - buf_size += sizeof(SETPLEFT); - buf_size += strlen(objs[i].name); - buf_size += log10(1+d); - for(size_t j = 0; j < objs[i].npars; ++j) - buf_size += sizeof(SETPARGS) + log10(1+p+j); - buf_size += sizeof(SETPRGHT); - d += objs[i].size; - p += objs[i].npars; - } - buf_size += sizeof(SETPFOOT); - buf_size += sizeof(FILEFOOT); - - // allocate buffer - buf = malloc(buf_size); - if(!buf) - errori(NULL); - - // start at beginning of data and parameters - p = d = 0; + // start empty and with 0 length to prevent writing + buf = NULL; + out = NULL; + siz = 0; + len = 0; - // output tracks current writing position on buffer - out = buf; - - // write file header - wri = sprintf(out, FILEHEAD, "", "set_params"); - if(wri < 0) - errori(NULL); - out += wri; - - // write header - wri = sprintf(out, SETPHEAD); - if(wri < 0) - errori(NULL); - out += wri; - - // write body - for(size_t i = 0; i < nobjs; ++i) + // two-pass: calculate buffer size and allocate, then fill + for(int pass = 0; pass < 2; ++pass) { - // write left side of line - wri = sprintf(out, SETPLEFT, objs[i].name, d); + // allocate buffer after first pass + if(pass > 0) + { + // allocate + buf = malloc(siz + 1); + if(!buf) + errori(NULL); + + // output tracks writing + out = buf; + + // maximum length is now huge + len = -1; + } + + // start with invalid trigger + trigger = 0; + + // keep track of where current plane starts + plane = 0; + + // start at beginning of data and parameters + p = d = 0; + + // write file header + wri = snprintf(out, len, FILEHEAD, "", "set_params"); if(wri < 0) errori(NULL); - out += wri; + if(pass > 0) + out += wri; + else + siz += wri; + + // write header + wri = snprintf(out, len, SETPHEAD); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; - // write arguments - for(size_t j = 0; j < objs[i].npars; ++j) + // write body + for(size_t i = 0; i < nobjs; ++i) { - wri = sprintf(out, SETPARGS, p+j); + // check if lens plane change is triggered + if(objs[i].type != trigger && objs[i].type != OBJ_FOREGROUND) + { + // when triggering from lenses to sources, change planes + if(trigger == OBJ_LENS && objs[i].type == OBJ_SOURCE) + plane = i; + + // new trigger + trigger = objs[i].type; + } + + // search for image plane priors in this object + for(size_t j = 0; j < objs[i].npars; ++j) + { + if(objs[i].pars[j].ipp) + { + // image plane prior depends on parameter type + switch(objs[i].pars[j].type) + { + // position: lens through all previous planes + case PAR_POSITION_X: + { + // inner trigger for changing lens planes + int trigger2 = 0; + + // inner data offset + size_t d2 = 0; + + // initialise position IPP + wri = snprintf(out, len, SETPIPP_POSINIT, p + j, p + j + 1); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + + // deflect IPP through previous planes + for(size_t k = 0; k < plane; ++k) + { + if(objs[k].type != trigger2 && objs[k].type != OBJ_FOREGROUND) + { + // deflect when triggering from lens + if(trigger2 == OBJ_LENS) + { + wri = snprintf(out, len, SETPIPP_POSDEFL); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + } + + // reset trigger + trigger2 = objs[k].type; + } + + // compute deflection for lens + if(objs[k].type == OBJ_LENS) + { + wri = snprintf(out, len, SETPIPP_POSLENS, objs[k].name, d2); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + } + + // increase data offset + d2 += objs[k].size; + } + + // apply final deflection when stopped with lens + if(trigger2 == OBJ_LENS) + { + wri = snprintf(out, len, SETPIPP_POSDEFL); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + } + } + + // handled together with X + case PAR_POSITION_Y: + break; + + // should not happen + default: + break; + } + } + } + + // write left side of line + wri = snprintf(out, len, SETPLEFT, objs[i].name, d); if(wri < 0) errori(NULL); - out += wri; + if(pass > 0) + out += wri; + else + siz += wri; + + // write arguments + for(size_t j = 0; j < objs[i].npars; ++j) + { + // special argument for image plane priors + if(objs[i].pars[j].ipp) + { + const char* arg; + + switch(objs[i].pars[j].type) + { + case PAR_POSITION_X: arg = "x.x"; break; + case PAR_POSITION_Y: arg = "x.y"; break; + default: arg = "0"; break; + } + + wri = snprintf(out, len, SETPIPPA, arg); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + } + else + { + wri = snprintf(out, len, SETPARGS, p + j); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + } + } + + // write right side of line + wri = snprintf(out, len, SETPRGHT); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; + + // increase offsets + d += objs[i].size; + p += objs[i].npars; } - // write left side of line - wri = sprintf(out, SETPRGHT); + // write footer + wri = snprintf(out, len, SETPFOOT); if(wri < 0) errori(NULL); - out += wri; + if(pass > 0) + out += wri; + else + siz += wri; - // increase offsets - d += objs[i].size; - p += objs[i].npars; + // write file footer + wri = snprintf(out, len, FILEFOOT); + if(wri < 0) + errori(NULL); + if(pass > 0) + out += wri; + else + siz += wri; } - // write footer - wri = sprintf(out, SETPFOOT); - if(wri < 0) - errori(NULL); - out += wri; - - // write file footer - wri = sprintf(out, FILEFOOT); - if(wri < 0) - errori(NULL); - out += wri; - // this is our code return buf; } diff --git a/src/lensed.c b/src/lensed.c index 8dd6519..8e70d09 100644 --- a/src/lensed.c +++ b/src/lensed.c @@ -108,68 +108,102 @@ int main(int argc, char* argv[]) for(size_t i = 0; i < inp->nobjs; ++i) lensed->npars += inp->objs[i].npars; - // get all parameters + // list of all parameters lensed->pars = malloc(lensed->npars*sizeof(param*)); if(!lensed->pars) errori(NULL); - for(size_t i = 0, p = 0; i < inp->nobjs; ++i) - for(size_t j = 0; j < inp->objs[i].npars; ++j, ++p) - lensed->pars[p] = &inp->objs[i].pars[j]; - // process parameters - for(size_t i = 0; i < lensed->npars; ++i) + // collect and check all parameters + for(size_t i = 0, p = 0; i < inp->nobjs; ++i) { - // current parameter - param* par = lensed->pars[i]; + // current object + object* obj = &inp->objs[i]; - // apply default bounds if none provided - if(!par->lower && !par->upper) + for(size_t j = 0; j < inp->objs[i].npars; ++j, ++p) { - switch(par->type) + // current parameter + param* par = &obj->pars[j]; + + // apply default bounds if none provided + if(!par->lower && !par->upper) { - case PAR_RADIUS: - par->lower = 0; - par->upper = HUGE_VAL; - break; + switch(par->type) + { + case PAR_RADIUS: + par->lower = 0; + par->upper = HUGE_VAL; + break; - case PAR_AXIS_RATIO: - par->lower = 0; - par->upper = 1; - break; + case PAR_AXIS_RATIO: + par->lower = 0; + par->upper = 1; + break; - default: - break; + default: + break; + } } - } - - // mark if parameter is bounded - par->bounded = (par->lower || par->upper); - - // check that prior is compatible with bounds - if(par->bounded && (prior_lower(par->pri) < par->lower || prior_upper(par->pri) > par->upper)) - { - // check if there is overlap between the bounds - if(prior_lower(par->pri) >= par->upper || prior_upper(par->pri) <= par->lower) + + // mark if parameter is bounded + par->bounded = (par->lower || par->upper); + + // check that prior is compatible with bounds + if(par->bounded && (prior_lower(par->pri) < par->lower || prior_upper(par->pri) > par->upper)) { - // prior falls completely outside of allowed range - error("%s: prior does not include parameter bounds [%g, %g]\n" - "The prior does not include the allowed range of values " - "for the parameter. It is impossible to draw a valid " - "parameter value. Please fix the prior to include the " - "range of values indicated above.", - par->id, par->lower, par->upper); + // check if there is overlap between the bounds + if(prior_lower(par->pri) >= par->upper || prior_upper(par->pri) <= par->lower) + { + // prior falls completely outside of allowed range + error("%s: prior does not include parameter bounds [%g, %g]\n" + "The prior does not include the allowed range of values " + "for the parameter. It is impossible to draw a valid " + "parameter value. Please fix the prior to include the " + "range of values indicated above.", + par->id, par->lower, par->upper); + } + else + { + // prior overlaps parameter bounds, issue a warning + warn("%s: prior partially outside parameter bounds [%g, %g]\n" + "Part of the prior lies outside of the allowed range of " + "parameter values. Values will be drawn from the prior " + "until a valid one is found. If the bounds of prior and " + "parameter do not overlap significantly, this can be " + "slow. You might want to change the prior accordingly.", + par->id, par->lower, par->upper); + } } - else + + // check image plane priors + if(par->ipp) { - // prior overlaps parameter bounds, issue a warning - warn("%s: prior partially outside parameter bounds [%g, %g]\n" - "Part of the prior lies outside of the allowed range of " - "parameter values. Values will be drawn from the prior " - "until a valid one is found. If the bounds of prior and " - "parameter do not overlap significantly, this can be " - "slow. You might want to change the prior accordingly.", - par->id, par->lower, par->upper); + // not all parameters can have image plane priors + switch(par->type) + { + // check that X is followed by IPP Y + case PAR_POSITION_X: + if(j + 1 == obj->npars || obj->pars[j + 1].type != PAR_POSITION_Y) + error("object `%s`: image plane prior requires pair (X,Y) of parameters", obj->id); + if(!obj->pars[j + 1].ipp) + error("object `%s`: pair (X,Y) must have image plane priors", obj->id); + break; + + // check that Y is preceded by IPP X + case PAR_POSITION_Y: + if(j == 0 || obj->pars[j - 1].type != PAR_POSITION_X) + error("object `%s`: image plane prior requires pair (X,Y) of parameters", obj->id); + if(!obj->pars[j - 1].ipp) + error("object `%s`: pair (X,Y) must have image plane priors", obj->id); + break; + + // all other parameters cannot have image plane priors + default: + error("prior `%s`: cannot give image plane prior for this parameter", par->id); + } } + + // store parameter + lensed->pars[p] = par; } }