parameter bounds

glenco · Jun 22, 2015 · 31bb54d · 31bb54d
1 parent e153556
commit 31bb54d
Show file tree

Hide file tree

Showing 12 changed files with 102 additions and 28 deletions.
diff --git a/kernel/object.cl b/kernel/object.cl
@@ -21,6 +21,7 @@ enum
 // structure that holds parameter definition
 struct param
 {
-    char name[32];
-    int  type;
+    char   name[32];
+    int    type;
+    float2 bounds;
 };
diff --git a/objects/sersic.cl b/objects/sersic.cl
@@ -2,13 +2,13 @@ type = SOURCE;
 
 params
 {
-    { "x",      POSITION_X  },
-    { "y",      POSITION_Y  },
-    { "r",      RADIUS      },
-    { "mag",    MAGNITUDE   },
-    { "n",      PARAMETER   },
-    { "q",      AXIS_RATIO  },
-    { "pa",     POS_ANGLE   }
+    { "x",   POSITION_X },
+    { "y",   POSITION_Y },
+    { "r",   RADIUS     },
+    { "mag", MAGNITUDE  },
+    { "n",   PARAMETER, { 0.5f, 8.0f } },
+    { "q",   AXIS_RATIO },
+    { "pa",  POS_ANGLE  }
 };
 
 data

diff --git a/src/input.c b/src/input.c
@@ -269,7 +269,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",
+                snprintf(tag, 100, " [%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 ?
@@ -282,6 +282,8 @@ void print_input(const input* inp)
                         "axis ratio, " : "",
                     inp->objs[i].pars[j].type == PAR_POS_ANGLE ?
                         "pos. angle, " : "",
+                    inp->objs[i].pars[j].bounded ?
+                        "bounded, " : "",
                     inp->objs[i].pars[j].wrap ?
                         "wrap, " : ""
                 );

diff --git a/src/input.h b/src/input.h
@@ -72,6 +72,11 @@ typedef struct
     // type of parameter
     int type;
 
+    // hard parameter bounds
+    double lower;
+    double upper;
+    int bounded;
+
     // prior for parameter
     prior* pri;
 

diff --git a/src/input/objects.c b/src/input/objects.c
@@ -46,6 +46,8 @@ void add_object(input* inp, const char* id, const char* name)
     cl_char*    param_names;
     cl_mem      param_types_mem;
     cl_int*     param_types;
+    cl_mem      param_bounds_mem;
+    cl_float2*  param_bounds;
 
     // realloc space for one more object
     inp->nobjs += 1;
@@ -146,7 +148,8 @@ void add_object(input* inp, const char* id, const char* name)
     // buffers for kernel parameters
     param_names_mem  = clCreateBuffer(lcl->context, CL_MEM_WRITE_ONLY, obj->npars*meta_nnam*sizeof(cl_char), NULL, NULL);
     param_types_mem  = clCreateBuffer(lcl->context, CL_MEM_WRITE_ONLY, obj->npars*sizeof(cl_int),            NULL, NULL);
-    if(!param_types_mem || !param_names_mem)
+    param_bounds_mem = clCreateBuffer(lcl->context, CL_MEM_WRITE_ONLY, obj->npars*sizeof(cl_float2),         NULL, NULL);
+    if(!param_types_mem || !param_names_mem || !param_bounds_mem)
         error("object %s: failed to create buffer for parameters", id);
 
     // setup and run kernel to get parameters
@@ -156,6 +159,7 @@ void add_object(input* inp, const char* id, const char* name)
         error("object %s: failed to create kernel for parameters", id);
     err |= clSetKernelArg(param_kernel, 0, sizeof(cl_mem), &param_names_mem );
     err |= clSetKernelArg(param_kernel, 1, sizeof(cl_mem), &param_types_mem );
+    err |= clSetKernelArg(param_kernel, 2, sizeof(cl_mem), &param_bounds_mem);
     if(err != CL_SUCCESS)
         error("object %s: failed to set kernel arguments for parameters", id);
     err = clEnqueueTask(queue, param_kernel, 0, NULL, NULL);
@@ -165,12 +169,14 @@ void add_object(input* inp, const char* id, const char* name)
     // arrays for kernel parameters
     param_names  = malloc(obj->npars*meta_nnam*sizeof(cl_char));
     param_types  = malloc(obj->npars*sizeof(cl_int));
-    if(!param_types || !param_names)
+    param_bounds = malloc(obj->npars*sizeof(cl_float2));
+    if(!param_types || !param_names || !param_bounds)
         errori("object %s", id);
 
     // get kernel parameters from buffer
     err |= clEnqueueReadBuffer(queue, param_names_mem,  CL_TRUE, 0, obj->npars*meta_nnam*sizeof(cl_char), param_names,  0, NULL, NULL);
     err |= clEnqueueReadBuffer(queue, param_types_mem,  CL_TRUE, 0, obj->npars*sizeof(cl_int),            param_types,  0, NULL, NULL);
+    err |= clEnqueueReadBuffer(queue, param_bounds_mem, CL_TRUE, 0, obj->npars*sizeof(cl_float2),         param_bounds, 0, NULL, NULL);
     if(err != CL_SUCCESS)
         error("object %s: failed to get parameters", id);
 
@@ -207,6 +213,8 @@ void add_object(input* inp, const char* id, const char* name)
         obj->pars[i].name   = name;
         obj->pars[i].id     = id;
         obj->pars[i].type   = param_types[i];
+        obj->pars[i].lower  = param_bounds[i].s[0];
+        obj->pars[i].upper  = param_bounds[i].s[1];
         obj->pars[i].pri    = NULL;
         obj->pars[i].wrap   = 0;
         obj->pars[i].label  = NULL;
@@ -219,6 +227,8 @@ void add_object(input* inp, const char* id, const char* name)
     free(param_names);
     clReleaseMemObject(param_types_mem);
     free(param_types);
+    clReleaseMemObject(param_bounds_mem);
+    free(param_bounds);
 
     free(param_kernam);
     clReleaseKernel(param_kernel);

diff --git a/src/kernel.c b/src/kernel.c
@@ -51,15 +51,17 @@ static const char METAKERN[] =
 
 // kernel to get parameters for object
 static const char PARSKERN[] = 
-    "kernel void params_<name>(global char* names, global int* types)\n"
+    "kernel void params_<name>(global char*   names, global int* types,\n"
+    "                          global float2* bounds)\n"
     "{\n"
     "    const size_t n = sizeof(parlst_<name>)/sizeof(struct param);\n"
     "    const size_t l = sizeof(((struct param*)0)->name);\n"
     "    for(size_t i = 0; i < n; ++i)\n"
     "    {\n"
     "        for(size_t j = 0; j < l; ++j)\n"
     "            names[i*l + j] = parlst_<name>[i].name[j];\n"
-    "        types[i] = parlst_<name>[i].type;\n"
+    "        types[i]  = parlst_<name>[i].type;\n"
+    "        bounds[i] = parlst_<name>[i].bounds;\n"
     "    }\n"
     "}\n"
 ;

diff --git a/src/lensed.c b/src/lensed.c
@@ -92,13 +92,62 @@ int main(int argc, char* argv[])
         for(size_t j = 0; j < inp->objs[i].npars; ++j, ++p)
             lensed.pars[p] = &inp->objs[i].pars[j];
 
-    // get all priors that will be needed when running
-    lensed.pris = malloc(lensed.npars*sizeof(prior*));
-    if(!lensed.pris)
-        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.pris[p] = inp->objs[i].pars[j].pri;
+    // process parameters
+    for(size_t i = 0; i < lensed.npars; ++i)
+    {
+        // current parameter
+        param* par = lensed.pars[i];
+
+        // apply default bounds if none provided
+        if(!par->lower && !par->upper)
+        {
+            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;
+
+            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)
+            {
+                // 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);
+            }
+        }
+    }
 
 
     /*****************

diff --git a/src/lensed.h b/src/lensed.h
@@ -12,7 +12,6 @@ struct lensed
     // parameter space
     size_t npars;
     param** pars;
-    prior** pris;
 
     // results
     const char* fits;

diff --git a/src/nested.c b/src/nested.c
@@ -18,7 +18,13 @@ void loglike(double cube[], int* ndim, int* npar, double* lnew, void* lensed_)
 
     // transform from unit cube to physical
     for(size_t i = 0; i < lensed->npars; ++i)
-        apply_prior(lensed->pris[i], &cube[i]);
+    {
+        double phys;
+        do
+            phys = apply_prior(lensed->pars[i]->pri, cube[i]);
+        while(lensed->pars[i]->bounded && (phys < lensed->pars[i]->lower || phys > lensed->pars[i]->upper));
+        cube[i] = phys;
+    }
 
     // error flag
     cl_int err = 0;

diff --git a/src/prior.c b/src/prior.c
@@ -164,9 +164,9 @@ void print_prior(const prior* pri, char* buf, size_t n)
     pri->print(pri->data, buf, n);
 }
 
-void apply_prior(const prior* pri, double* u)
+double apply_prior(const prior* pri, double u)
 {
-    *u = pri->prior(pri->data, *u);
+    return pri->prior(pri->data, u);
 }
 
 double prior_lower(const prior* pri)

diff --git a/src/prior.h b/src/prior.h
@@ -10,7 +10,7 @@ void free_prior(prior* pri);
 void print_prior(const prior* pri, char* buf, size_t n);
 
 // apply prior to unit variate
-void apply_prior(const prior* pri, double* u);
+double apply_prior(const prior* pri, double u);
 
 // get lower bound for prior
 double prior_lower(const prior* pri);

diff --git a/src/prior/norm.c b/src/prior/norm.c
@@ -83,10 +83,10 @@ double prior_norm(const void* data, double u)
 
 double prior_lower_norm(const void* data)
 {
-    return -INFINITY;
+    return -HUGE_VAL;
 }
 
 double prior_upper_norm(const void* data)
 {
-    return INFINITY;
+    return +HUGE_VAL;
 }