Speed improvement: Use `logit` transformed transition kernel

Introduction
--

This is a pretty exciting change.  I'll go slow so this is understandable.

To illustrate the change, we'll solve a simple Bayesian inverse problem in three dimensions, the seed of the chain is near the boundary of the domain and the proposal covariance is about three times the size of the domain.  It's a little convoluted but this is just an illustration.  For the more detail-oriented readers, I've included the example code (and input file) at the bottom of this post.

The posterior we will sample is Gaussian with mean at the mid-point of the domain.

Result
--

Here's the main timing result before this change:

```
Total run time = 236 seconds
```

And here are the samples from the first component of the chain (before the change):

![samples-before](https://cloud.githubusercontent.com/assets/328518/4153533/dbc7b038-3459-11e4-931b-14fe909192b2.png)

For comparison, here's the timing result after this change:

```
Total run time = 57 seconds
```

And here are the samples from the first component of the chain (after the change):

![samples-after](https://cloud.githubusercontent.com/assets/328518/4153648/edcfc3c8-345a-11e4-9e21-73ab44cb6191.png)

Explanation
--

So, in summary, we're getting a nice speedup with this change.  Why?  Well, `queso` uses a Gaussian proposal distribution, so if the problem we're solving is posed on a bounded domain, it can be the case that `queso` will propose a state outside of the domain of the problem, which is a probability zero event.  In fact, as the dimension of the problem increases, the likelihood of generating a proposal outside of the domain grows (and grows very quickly).

Instead, this PR applies a coordinate transformation to the Gaussian transition kernel so that all proposals generated are within the domain of the problem.  This way, we can still draw a state X from a Gaussian transition kernel, which is easy to do, and simply propose f(X), where

f(x) = (b exp(x) + a) / (1 + exp(x))

This transformation ensures that f(X) is contained in [a, b].

What is supported?
--

Of course, for some problems, the domain is not bounded.  Or there might be directions in the domain that are bounded and some that are not.  Or there might be directions of the domain that are only bounded on one side.  This code handles all of those cases.

For directions of the domain that are unbounded, the transformation above is just the identity transformation and the transition kernel falls back to the usual Gaussian distribution.  For directions that are only bounded on one side, a half-transform is done.

Code design
--

The PR is designed such that the user should not have to change their code.

The only thing I'm not certain on is where, exactly, should all the transform logic sit.  Currently, there's transform logic in the PDF of the random variable f(X), which is unavoidable.  There's also transform logic in the transition kernel, and this is needed because when the Markov chain wants to update the mean of the transition kernel, the transition kernel can map between coordinate spaces and set the mean of the *underlying* Gaussian distribution, not the mean of the transformed distribution.

It would be nice to just have all the transform logic sitting in one place, rather than in several, but I can't think of how to do that.  Maybe it is not possible.

What needs to be done?
--

I've had to update a couple of the regression test reference solutions, because different samples are drawn after this change due to the change in transition kernel.

I need to write some tests for the `InvLogit*` classes, but that's not hard to do.

Appendix
--

Example code:

```cpp
#include <queso/GenericScalarFunction.h>
#include <queso/GslVector.h>
#include <queso/GslMatrix.h>
#include <queso/UniformVectorRV.h>
#include <queso/StatisticalInverseProblem.h>
#include <queso/ScalarFunction.h>
#include <queso/VectorSet.h>

template<class V, class M>
class Likelihood : public QUESO::BaseScalarFunction<V, M>
{
public:

  Likelihood(const char * prefix, const QUESO::VectorSet<V, M> & domain)
    : QUESO::BaseScalarFunction<V, M>(prefix, domain)
  {
  }

  virtual ~Likelihood()
  {
  }

  virtual double lnValue(const V & domainVector, const V * domainDirection,
      V * gradVector, M * hessianMatrix, V * hessianEffect) const
  {
    // Just evaluate a Gaussian likelihood
    double sum = 0.0;

    for (unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
      sum += (domainVector[i] - 4.5) * (domainVector[i] - 4.5);
    }

    return -50.0 * sum;
  }

  virtual double actualValue(const V & domainVector, const V * domainDirection,
      V * gradVector, M * hessianMatrix, V * hessianEffect) const
  {
    return std::exp(this->lnValue(domainVector, domainDirection, gradVector,
          hessianMatrix, hessianEffect));
  }
};

int main(int argc, char ** argv) {
  MPI_Init(&argc, &argv);

  QUESO::FullEnvironment env(MPI_COMM_WORLD, argv[1], "", NULL);

  QUESO::VectorSpace<QUESO::GslVector, QUESO::GslMatrix> paramSpace(env,
      "param_", 3, NULL);

  double min_val = 4.0;
  double max_val = 5.0;

  QUESO::GslVector paramMins(paramSpace.zeroVector());
  paramMins.cwSet(min_val);
  QUESO::GslVector paramMaxs(paramSpace.zeroVector());
  paramMaxs.cwSet(max_val);

  QUESO::BoxSubset<QUESO::GslVector, QUESO::GslMatrix> paramDomain("param_",
      paramSpace, paramMins, paramMaxs);

  QUESO::UniformVectorRV<QUESO::GslVector, QUESO::GslMatrix> priorRv("prior_",
      paramDomain);

  Likelihood<QUESO::GslVector, QUESO::GslMatrix> lhood("llhd_", paramDomain);
  
  QUESO::GenericVectorRV<QUESO::GslVector, QUESO::GslMatrix>
    postRv("post_", paramSpace);

  QUESO::StatisticalInverseProblem<QUESO::GslVector, QUESO::GslMatrix>
    ip("", NULL, priorRv, lhood, postRv);

  QUESO::GslVector paramInitials(paramSpace.zeroVector());

  paramInitials[0] = 4.99;
  paramInitials[1] = 4.99;
  paramInitials[2] = 4.99;
  
  QUESO::GslMatrix proposalCovMatrix(paramSpace.zeroVector());

  for (unsigned int i = 0; i < 3; i++) {
    proposalCovMatrix(i, i) = 3.0;
  }
  
  ip.solveWithBayesMetropolisHastings(NULL, paramInitials, &proposalCovMatrix);

  MPI_Finalize();
 
  return 0;
}
```

Example input file:

```
###############################################
# UQ Environment
###############################################
#env_help                = anything
env_numSubEnvironments   = 1
env_subDisplayFileName   = outputData/display
env_subDisplayAllowAll   = 0
env_subDisplayAllowedSet = 0
env_displayVerbosity     = 2
env_syncVerbosity        = 0
env_seed                 = 0

###############################################
# Statistical inverse problem (ip)
###############################################
#ip_help                 = anything
ip_computeSolution      = 1
ip_dataOutputFileName   = outputData/sipOutput
ip_dataOutputAllowedSet = 0

###############################################
# 'ip_': information for Metropolis-Hastings algorithm
###############################################
#ip_mh_help                 = anything
ip_mh_dataOutputFileName   = outputData/sipOutput
ip_mh_dataOutputAllowedSet = 0 1

ip_mh_rawChain_dataInputFileName    = .
ip_mh_rawChain_size                 = 1000000
ip_mh_rawChain_generateExtra        = 0
ip_mh_rawChain_displayPeriod        = 50000
ip_mh_rawChain_measureRunTimes      = 1
ip_mh_rawChain_dataOutputFileName   = outputData/ip_raw_chain
ip_mh_rawChain_dataOutputAllowedSet = 0 1
ip_mh_rawChain_computeStats         = 1

ip_mh_displayCandidates             = 0
ip_mh_putOutOfBoundsInChain         = 0
ip_mh_tk_useLocalHessian            = 0
ip_mh_tk_useNewtonComponent         = 1
ip_mh_dr_maxNumExtraStages          = 1
ip_mh_dr_listOfScalesForExtraStages = 5.
ip_mh_am_initialNonAdaptInterval    = 0
ip_mh_am_adaptInterval              = 10
ip_mh_am_eta                        = 1.92
ip_mh_am_epsilon                    = 1.e-5

ip_mh_filteredChain_generate             = 0
```

Ok, here's the best explanation I can come up with.

When the user specifies the covariance of the proposal distribution, there's not enough information here to know exactly what is happening.  Why?  Well, because now we have introduced the concept of a transformed coordinate space.  Before, there was no question the space the covariance was being set in was always the same as the space the user's parameters lived in.  In this space, the proposal distribution was always Gaussian.  Now, there is a coordinate transformation that sits between the space the user's parameters live in and the space the proposals live in.  The space the proposals live in is always the real line, and in this space the proposal distribution is Gaussian.  This is the space the proposal covariance is being set by the user (and what is being set in this example).  Once that proposed state is mapped back to the user space (which may be bounded), the covariance doesn't quite mean the same thing.  It's a stretched covariance.  This is why I had to 'unstretch' it here.

All that said, I think what might be useful is if I sit down and document what happens to the covariance of the proposal distribution under the logit coordinate transformation.  I think it's probably unfair to ask the user to specify the covariance of a distribution in a space they are otherwise completely shielded from.