Calculate total force of muscle contraction

[vellamike]

It is not obvious how to calculate force extended to contract one of the C elegans muscles in the Sibernetic simulation. This is in part due to the fact that it is not clear what the mass and acceleration of a particle is at any given time. We need a good way to calculate this.

[a-palyanov]

Test scene is ready

[a-palyanov]

It seems that we need significanlty higher dissipation within a muscle.

[vellamike]

agreed - those oscillations are quite extreme. What parameter controls damping?

[a-palyanov]

I expect that there is no dumping in such a system except inter-particle interaction forces (like viscosity and pressure forces) which are present in elastic matter because elastic matter particles also possess the same properties which are inherent to liquid particles. Possibly we should additionally introduce some actual dumping force?

[a-palyanov]

Relaxation in case of short muscle-activating impulse

[vellamike]

Andrey, what happens if the impulse is maintained? I expect there will be no oscillation?

On 12 January 2014 05:49, a-palyanov notifications@github.com wrote:

Relaxation in [image: 4]https://f.cloud.github.com/assets/1803774/1895665/38827eaa-7b4d-11e3-8d51-06b262e109ea.png case of short muscle-activating impulse

— Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-32115997 .

[vellamike]

Sorry I just saw the post with the long impulse, I'm confused about what the graphs mean exactly: If the muscle is contracted, why isn't the position of the centre particle moving in one net direction?

On 12 January 2014 05:50, Mike Vella vellamike@gmail.com wrote:

Andrey, what happens if the impulse is maintained? I expect there will be no oscillation?

On 12 January 2014 05:49, a-palyanov notifications@github.com wrote:

Relaxation in [image: 4]https://f.cloud.github.com/assets/1803774/1895665/38827eaa-7b4d-11e3-8d51-06b262e109ea.png case of short muscle-activating impulse

— Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-32115997 .

[vellamike]

Nevermind, I see it moves from 0 to -5. Could you upload a video of this?

On 12 January 2014 05:51, Mike Vella vellamike@gmail.com wrote:

Sorry I just saw the post with the long impulse, I'm confused about what the graphs mean exactly: If the muscle is contracted, why isn't the position of the centre particle moving in one net direction?

On 12 January 2014 05:50, Mike Vella vellamike@gmail.com wrote:

Andrey, what happens if the impulse is maintained? I expect there will be no oscillation?

On 12 January 2014 05:49, a-palyanov notifications@github.com wrote:

Relaxation in [image: 4]https://f.cloud.github.com/assets/1803774/1895665/38827eaa-7b4d-11e3-8d51-06b262e109ea.png case of short muscle-activating impulse

— Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-32115997 .

[vellamike]

I also notice that the oscillations are very fast: ~ 200Hz.

On 12 January 2014 05:52, Mike Vella vellamike@gmail.com wrote:

Nevermind, I see it moves from 0 to -5. Could you upload a video of this?

On 12 January 2014 05:51, Mike Vella vellamike@gmail.com wrote:

Sorry I just saw the post with the long impulse, I'm confused about what the graphs mean exactly: If the muscle is contracted, why isn't the position of the centre particle moving in one net direction?

On 12 January 2014 05:50, Mike Vella vellamike@gmail.com wrote:

Andrey, what happens if the impulse is maintained? I expect there will be no oscillation?

On 12 January 2014 05:49, a-palyanov notifications@github.com wrote:

Relaxation in [image: 4]https://f.cloud.github.com/assets/1803774/1895665/38827eaa-7b4d-11e3-8d51-06b262e109ea.png case of short muscle-activating impulse

— Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-32115997 .

[a-palyanov]

[a-palyanov]

Here is the video of muscle contraction test: www.youtube.com/watch?v=s1Di8qgltRo

[Neurophile]

A few questions: 1) How many "particles" make up the muscle shown? 2) What is the cross section area of a) the thin muscle, b) the thick muscle (the impulse energy of the thick muscle looks to be about 2x the thin muscle just based on the decay curve) 3) What is the viscosity parameter for elastic particles? Is this different (higher) than fluid particles?

Sorry if these questions are too naive, I am still trying to understand the operation of the SPH model.

[Neurophile]

I believe this paper has some of the proper parameters we should be attempting to match: "Biomechanical analysis of gait adaptation in the nematode Caenorhabditis elegans" http://www.pnas.org/content/107/47/20323.full

It also has interesting implications for the neural locomotion circuits. Again, apologies if this has all been covered before.

[a-palyanov]

1) How many "particles" make up the muscle shown? 2) What is the cross section area of a) the thin muscle, b) the thick muscle (the >impulse energy of the thick muscle looks to be about 2x the thin muscle just >based on the decay curve)

Thin muscle (red half) is 3_5_30 = 450 particles (as well as another, ordinary elastic matter, green half) Thick variant muscle is 10_11_30 = 3300 particles.

3) What is the viscosity parameter for elastic particles? Is this different (higher) >than fluid particles?

Viscosity parameter (=0.00005) is the same for all particles of the simulation.

[a-palyanov]

Liquid damper

[a-palyanov]

Same spring, less liquid around ==> less damping (no gravity); By the way, frequency is lower than for the same system without liquid. By the way, I don't know yet why small amplitude oscillation remains. In case of our latest worm body muscles there shouldn't be long oscillations - worm body should work as a big damper.

[Neurophile]

"Viscosity parameter (=0.00005) is the same for all particles of the simulation."

After researching SPH and Navier-Stokes (NS) methods in general, I still don't understand how this value for viscosity was chosen. My very shaky understanding is that the viscosity term in NS is generally responsible for energy dissipation. The above graphs roughly match my intuition about how the equations should behave given the energy input and the damping term due to viscosity. More particles = more viscous forces = more damping. The viscous term is also relative to the Laplacian of the particle velocity vector field, so that seems to explain the behavior where the oscillations decay to some low value and then remain stable.

But why 5E-5 ? The measured value for water is about 1E-3 (http://www.engineeringtoolbox.com/water-dynamic-kinematic-viscosity-d_596.html) Is this due to scaling or some other property of the SPH method?

[a-palyanov]

Here is an illustration of how the system relaxation goes during a long muscle activating impulse .

[a-palyanov]

Muscle fiber contraction force is defined by the following code in Sibernetic (line 781-782 in sphFluid.cl code at contraction_force_test branch; might be a bit other lines at other branches, but corresponding code is the same): if(muscle_activation_signal[i]>0.f) acceleration[ id ] += -(vect_r_ij/r_ij) * muscle_activation_signal[i] * 800.f;

I've tried to change the value of 800 to others, in the interval from 500 to 4500, and built a graph of central point displacement (showed by orange double-arrow at the previous image) depending from this value:

[a-palyanov]

Possibly we need to implement something like this in Sibernetic for more realistic muscle fibers simulation: There is one more way, which might be even better - do no changes in 3D geometry of already existing worm muscles, but change a way of contraction force calculation - make it dependent from current shrinkage of muscle (difference between current and relaxed state) and reproduce the dependence similar to real one.

[Neurophile]

How much biological realism are we trying to achieve? Worm Atlas has a very nice primer on worm muscle anatomy: http://www.wormatlas.org/hermaphrodite/musclesomatic/MusSomaticframeset.html

Some diagrams, created mostly for my own understanding.

[a-palyanov]

I used browser.openworm.org as a sourse of worm's musles 3D geometry data. According to this, muscles has not so regular structure. Due to quite low 'resolution' of our worm body model I had to approximate them like this:

[a-palyanov]

Anyway, my recent picture with relaxed and contracted muscle states is about the following: we need to derive all necessary macroscopic features of c. elegans muscle cell contraction and reproduce them in our simulation. The most simple example - when we know, how real muscle contraction force ('muscle power') denepds on current muscle contraction, we may reproduce it and use it to calculate force for each muscle fiber of the simulation depending on its current state. Currently there is no any realistic depencence, it is linear.

At my recent picture with relaxed and contracted muscle states I've tried to offer a simple model of actin-myosin elementary structure represented in terms of Sibernetic primitives. But now I think that it will lead to growth of the model size (which is highly undesirable), and we can solve the problem by just defining the dependence between 'muscle power' and current muscle degree of contraction, inherent to natural process of muscle contraction.

[a-palyanov]

Neurophile>But why 5E-5 ? The measured value for water is about 1E-3 I don't know, but possibly because the size of 'elementary real water particle' is significanty less that the size of a particle in Sibernetic. It we'll use 1e-3, simulated water will be visually too viscous to be realistic. I think it can be useful to use the experience of other researchers who also faced the same problem.

[Neurophile]

@a-palyanov "we can solve the problem by just defining the dependence between 'muscle power' and current muscle degree of contraction, inherent to natural process of muscle contraction."

The real worm modulates muscle power over several orders of magnitude depending on the external viscous resistance. Now we are always using water, but even changing the water depth will affect the power needed to maintain a normal swimming gait. The paper I linked earlier shows some values for peak power. The muscles should be strong enough to produce that level, and then be modulated by nervous inputs. For the current model without a nervous system, the muscles can be activated with some reduced duty cycle.

[a-palyanov]

While trying to solve a problem of muscle contraction force measurement, I've revealed one more problem. A very simple scene with a single particle falling down the flat ground, dependence of its coordinate from time behaves like g=-2.45, but its value defined in Sibernetic is g=-9.8. So, something obviously works wrong. Maybe integration. Currently I don't know when this problem appeared - before or after optimization, or it already existed a long time ago.

[a-palyanov]

(-4.9*2)/2 = -9.8/2, so we lose multiplier of 2 somewhere. So pcisph_integrate() is highly probable.

[a-palyanov]

I was right - the problem was connected with integration method. We used some which existed in the code from very old times. When I replaced it with a quite good one (so-called 'Velocity Verlet', a second order method), free falling single particle test worked perfectly. This method needs as input position(t), velocity(t), acceleration(t) and acceleration(t+1). Previous one (old) method used only position(t), velocity(t), acceleration(t), so most probably it was an Euler's method (first order, low precision) with an error in realization.

[a-palyanov]

Free falling particle test was working perfectly for new integration method. After that I've performed another one - also a very simple configuration, which included a particle hanging on the end of the vertical spring at its resting length at time=0, with its upper end attached to the floor. Here are y(t) dependencies at different timeStep values Unfortunately, at our usual timeStep = 5e-06 this test shows quite noticeable error, leading to amplitude grouth, which is inconvenient and physically incorrect. Old method error was larger, but it caused damping, which 'looked' natural. I've found the following opinion on 'Velocity Verlet' integrator: "All the dynamics objects need a time step. Choosing it too small will waste computer time, choosing it too large will make the dynamics unstable, typically the energy increases dramatically (the system “blows up”). If the time step is only a little to large, the lack of energy conservation is most obvious in Velocity Verlet dynamics, where energy should otherwise be conserved." So, my next step is to find the most appropriate integrator for our case.

[a-palyanov]

There is one more widely-used second order explicit Runge-Kutta method which is also called the improved Euler's method, or modified trapezoidal method, or Heun's method. Two function evaluations are required per time step, which in the case of multibody systems implies the solution of the equations of motion to obtain the accelerations twice at the given time step. It will require quite a lot of changes in the code and it will result in still 2-nd order precision. Fourth order Runge-Kutta integration method is even more complex. But what if we will modify Velocity Verlet method by introducing damping into the system. We can compensate grouth of the oscillations amplitude or use any other necessary damping value:

[a-palyanov]

I'm currently on the final stage of choosing the best integration method for equations of motion in Sibernetic - a number of different methods were implemented and compared at different integration time steps. It is very important, because different methods show noticeably different maximal adequate integration time steps, and choosing the correct one can give significant impact to overall computational performance of the system. Also I've fixed the problem of correct values of rigidity coefficient and muscle contraction force values - now they correspond to real values and can be described in the paper without problems. Worm muscle in actual model is more difficult than just a single 'muscle connection' between two adjacent particles - there are many other elastic connections within a muscle structure - transversal and diagonal to the direction of contraction - which will result in another relation between muscle rigidity and adequate contraction force, so I'll need to make additional calculations to get this values. So, things are getting better with muscular system implementation and its description. Returning to integration methods: I've implemented a set of them: Explicit Euler, Implicit Euler, Improved Euler, Velocity Verlet, Runge-Kutta 2nd order and LeapFrog along with analytical solution for a very simple system - single mass point attached to a vertical spring (its another end is attached to the ceiling).

[a-palyanov]

Comparing LeapFrog integration (2nd oder of precision; high stability and energy conservation) results with analytical solution (blue curve) I've noticed that the precision is quite poor, impossibly poor (green curve). But all the formulas are correct. Then I've reproduced the same math within Excel sheet (red curve), which gave the best precision I've ever seen during all recent calculation experiments. So, probably something is wrong with the work of C++ or OpenCL code.

[tarelli]

@borismarin do you have any idea about what Andrey is seeing?

[borismarin]

My two cents: Andrey's previous analysis makes total sense given the theory: for systems without dissipation (conservative) -- as the harmonic oscillator in the ceiling -- it is important to use symplectic integrators, such as Euler-Cromer or Verlet, unlike Euler or usual RK schemes. Notice that the increase in energy seen above is expected for nonsymplectic stepping. But, as Andrey noticed, this problem is mitigated in dissipative systems. Try adding a -bv viscous friction to the example and comparing the results...

One (possibly) easy way to verify the correctness of your leapfrog implementation is testing it for time reversibility -- going backwards in time should bring you to the starting point.

It may as well be that, for the complicated cases, you are experiencing numerical instabilities due to stiffness. In that case, only implicit methods will help -- such as Backwards (or implicit) Euler, trapezoid rule. Unfortunately, implicit methods mean that there is an extra rootfinding step made at each time step. Using Newton-like iteration is usually ok, but can slow things down if derivatives are expensive.

[a-palyanov]

Here is a code which I use as a test (LeapFrog integration, symplectic, 2nd order): void do_test() { float x_t1=33.4f,v_t1=0.f,a_t1=-9.8f; float x_t2,v_t2,a_t2; float dt = 0.00003; float t;

for(t=0.f;t<16.f*dt;t+=dt)
{
    printf("\n===[ x_t1= %e, v_t1= %e, a_t1= %e ]===",x_t1,v_t1,a_t1);

    x_t2 = x_t1 + v_t1*dt + a_t1*dt*dt/2.f;
    a_t2 = -1e-7*x_t2/3.25e-14 + (-9.8f); // 1e-7 - coefficient of rigidity, 3.25e-14 - particle mass
    v_t2 = v_t1 + (a_t1 + a_t2)*dt/2.f;

    x_t1 = x_t2;
    v_t1 = v_t2;
    a_t1 = a_t2;
}

return;

}

If I start with x_t1= 0.0f, I get completely correct solution (red curve on previous image): ===[ x_t1= 0.000000e+000, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= -4.410000e-009, v_t1= -2.937965e-004, a_t1= -9.786431e+000 ]=== ===[ x_t1= -1.762779e-008, v_t1= -5.867793e-004, a_t1= -9.745761e+000 ]=== ===[ x_t1= -3.961676e-008, v_t1= -8.781372e-004, a_t1= -9.678102e+000 ]=== ===[ x_t1= -7.031602e-008, v_t1= -1.167063e-003, a_t1= -9.583643e+000 ]=== ...

but if, for example, x_t1=33.4f (the system is just located in other point of the space, everything else is completely the same, so velocities and accelerations should keep the same, and coordinates should only shift for a value of 33.4), everything changes (green curve on previous image): ===[ x_t1= 3.340000e+001, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= 3.340000e+001, v_t1= -1.541539e+003, a_t1= -1.027692e+008 ]=== ===[ x_t1= 3.330751e+001, v_t1= -4.620347e+003, a_t1= -1.024847e+008 ]=== ===[ x_t1= 3.312278e+001, v_t1= -7.686361e+003, a_t1= -1.019163e+008 ]=== ===[ x_t1= 3.284633e+001, v_t1= -1.073109e+004, a_t1= -1.010656e+008 ]=== ...

Currenlty I don't know how to fix this. And yes, it looks like numerical instabilities.

[Neurophile]

Excel uses 64 bits for floating point numbers. Try declaring all variables as double and compare to the excel results. I need more time to think about the other issue.

On Wed, Feb 19, 2014 at 12:17 PM, a-palyanov notifications@github.comwrote:

Here is a code which I use as a test: void do_test() { float x_t1=33.4f,v_t1=0.f,a_t1=-9.8f; float x_t2,v_t2,a_t2; float dt = 0.00003; float t;

for(t=0.f;t<16.f*dt;t+=dt) { printf("\n===[ x_t1= %e, v_t1= %e, a_t1= %e ]===",x_t1,v_t1,a_t1);

x_t2 = x_t1 + v_t1*dt + a_t1*dt*dt/2.f;
a_t2 = -1e-7*x_t2/3.25e-14 + (-9.8f);
v_t2 = v_t1 + (a_t1 + a_t2)*dt/2.f;

x_t1 = x_t2;
v_t1 = v_t2;
a_t1 = a_t2;

}

return;

}

If I start with x_t1= 0.0f, I get completely correct solution (red curve on previous image): ===[ x_t1= 0.000000e+000, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= -4.410000e-009, v_t1= -2.937965e-004, a_t1= -9.786431e+000 ]=== ===[ x_t1= -1.762779e-008, v_t1= -5.867793e-004, a_t1= -9.745761e+000 ]=== ===[ x_t1= -3.961676e-008, v_t1= -8.781372e-004, a_t1= -9.678102e+000 ]=== ===[ x_t1= -7.031602e-008, v_t1= -1.167063e-003, a_t1= -9.583643e+000 ]=== ...

but if, for example, x_t1=33.4f (the system is just located in other point of the space, everything else is completely the same, so velocities and accelerations should keep the same, and coordinates should only shift for a value of 33.4), everything changes (green curve on previous image): ===[ x_t1= 3.340000e+001, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= 3.340000e+001, v_t1= -1.541539e+003, a_t1= -1.027692e+008 ]=== ===[ x_t1= 3.330751e+001, v_t1= -4.620347e+003, a_t1= -1.024847e+008 ]=== ===[ x_t1= 3.312278e+001, v_t1= -7.686361e+003, a_t1= -1.019163e+008 ]=== ===[ x_t1= 3.284633e+001, v_t1= -1.073109e+004, a_t1= -1.010656e+008 ]=== ...

Currenlty I don't know how to fix this.

Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-35529899 .

[Neurophile]

Actually there is an error in the math.

a_t2 = -1e-7*x_t2/3.25e-14 + (-9.8f);

the x value needs to be normalized vs x0, otherwise moving the spring in space will give you crazy values.

On Wed, Feb 19, 2014 at 12:40 PM, Chris Jensen < chrisandmichele.jensen@gmail.com> wrote:

Excel uses 64 bits for floating point numbers. Try declaring all variables as double and compare to the excel results. I need more time to think about the other issue.

On Wed, Feb 19, 2014 at 12:17 PM, a-palyanov notifications@github.comwrote:

Here is a code which I use as a test: void do_test() { float x_t1=33.4f,v_t1=0.f,a_t1=-9.8f; float x_t2,v_t2,a_t2; float dt = 0.00003; float t;

for(t=0.f;t<16.f*dt;t+=dt) { printf("\n===[ x_t1= %e, v_t1= %e, a_t1= %e ]===",x_t1,v_t1,a_t1);

x_t2 = x_t1 + v_t1*dt + a_t1*dt*dt/2.f;
a_t2 = -1e-7*x_t2/3.25e-14 + (-9.8f);
v_t2 = v_t1 + (a_t1 + a_t2)*dt/2.f;

x_t1 = x_t2;
v_t1 = v_t2;
a_t1 = a_t2;

}

return;

}

If I start with x_t1= 0.0f, I get completely correct solution (red curve on previous image): ===[ x_t1= 0.000000e+000, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= -4.410000e-009, v_t1= -2.937965e-004, a_t1= -9.786431e+000 ]=== ===[ x_t1= -1.762779e-008, v_t1= -5.867793e-004, a_t1= -9.745761e+000 ]=== ===[ x_t1= -3.961676e-008, v_t1= -8.781372e-004, a_t1= -9.678102e+000 ]=== ===[ x_t1= -7.031602e-008, v_t1= -1.167063e-003, a_t1= -9.583643e+000 ]=== ...

but if, for example, x_t1=33.4f (the system is just located in other point of the space, everything else is completely the same, so velocities and accelerations should keep the same, and coordinates should only shift for a value of 33.4), everything changes (green curve on previous image): ===[ x_t1= 3.340000e+001, v_t1= 0.000000e+000, a_t1= -9.800000e+000 ]=== ===[ x_t1= 3.340000e+001, v_t1= -1.541539e+003, a_t1= -1.027692e+008 ]=== ===[ x_t1= 3.330751e+001, v_t1= -4.620347e+003, a_t1= -1.024847e+008 ]=== ===[ x_t1= 3.312278e+001, v_t1= -7.686361e+003, a_t1= -1.019163e+008 ]=== ===[ x_t1= 3.284633e+001, v_t1= -1.073109e+004, a_t1= -1.010656e+008 ]=== ...

Currenlty I don't know how to fix this.

Reply to this email directly or view it on GitHubhttps://github.com/openworm/OpenWorm/issues/158#issuecomment-35529899 .

[a-palyanov]

Thanks, @borismarin and @Neurophile ! Improved code is here: void do_test() { double x0=33.4f; double x_t1=x0,v_t1=0.f,a_t1=-9.8f; double x_t2,v_t2,a_t2; double dt = 0.00003; double t;

for(t=0.f;t<16.f*dt;t+=dt)
{
    printf("\nt=%f, x(t)= %e, v(t)= %e, a(t)= %e",t,x_t1-x0,v_t1,a_t1);

    x_t2 = x_t1 + v_t1*dt + a_t1*dt*dt/2.0;
    a_t2 = -1e-7*(x_t2-x0)/3.25e-14 + (-9.8);// 1e-7 - coefficient of rigidity, 3.25e-14 - particle mass
    v_t2 = v_t1 + (a_t1 + a_t2)*dt/2.0;

    x_t1 = x_t2;
    v_t1 = v_t2;
    a_t1 = a_t2;
}

return;

}

now works fine when all values are double, both for x0 = 0 and x0 = 33.4 (I substract x0 from x(t) at output for easier comparison) t=0.000000, x(t)= 0.000000e+000, v(t)= 0.000000e+000, a(t)= -9.800000e+000 t=0.000030, x(t)= -4.409998e-009, v(t)= -2.937965e-004, a(t)= -9.786431e+000 t=0.000060, x(t)= -1.762778e-008, v(t)= -5.867793e-004, a(t)= -9.745761e+000 t=0.000090, x(t)= -3.961675e-008, v(t)= -8.781373e-004, a(t)= -9.678102e+000 t=0.000120, x(t)= -7.031601e-008, v(t)= -1.167063e-003, a(t)= -9.583643e+000

and works wrong for float: t=0.000000, x(t)= 0.000000e+000, v(t)= 0.000000e+000, a(t)= -9.800000e+000 t=0.000030, x(t)= 0.000000e+000, v(t)= -2.940000e-004, a(t)= -9.800000e+000 t=0.000060, x(t)= 0.000000e+000, v(t)= -5.880000e-004, a(t)= -9.800000e+000 t=0.000090, x(t)= 0.000000e+000, v(t)= -8.820000e-004, a(t)= -9.800000e+000 t=0.000120, x(t)= 0.000000e+000, v(t)= -1.176000e-003, a(t)= -9.800000e+000

[Neurophile]

This is an unavoidable problem with numerical methods using finite precision. If it is not possible to periodically calculate an error term against an analytical form, the best we can do is to minimize the sources of accumulated error and pick a method that is stable in the domain of our problem space. Note that C/C++ is slightly different from OpenCL in how it performs math on single precision float variables.

Tricks for lowering error:

• Only add/subtract terms that are relatively close in magnitude.
• Avoid subtracting terms that are very close in value
• Combine constant terms
• Only use as many math operations as absolutely necessary
• Take control of the rounding method

More details of binary floating point math and it's error sources can be found here: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

[tarelli]

gotta love this project :)

[slarson]

I believe this is now settled.