ElMaestro
★★★

Denmark,
2018-03-15 10:24
(1408 d 11:53 ago)

Posting: # 18549
Views: 6,209
 

 Nervous ticks in C - a reminder of the cruel logic of logic [Software]

Hi all,

sometimes logic totally sucks.

I am making code in C for plotting various graphs, and am making a little routine for inserting ticks on an axis. Ticks have to be "nice" (yes, that's also the word used in the R documentation).

So, let us say the maximum observed value is x=0.234567 and I want ticks evenly dispersed at a distance of dx=0.05 (say from 0.0 and upwards). So, I thought I was clever and calculated the "maximum" tick as
xmax = dx*ceil(x/dx)

You can see one definition of the ceil function here.

So, one way to test it is like this, which will print the actual ticks:

int Test10001(double x, double dx)
{
    double xmax;
    xmax=dx*ceil(x/dx);
    x=0.0;
    while (x<=xmax)
    {
        printf("x=%f\n",x);
        x=x+dx;
    }
    return(0);
}


When I call that function with arguments 0.234567 and 0.05 it prints:
x=0.000000
x=0.050000
x=0.100000
x=0.150000
x=0.200000


but it doesn't print 0.25 ?!???:-D:-D:-D:-D

If I add a line like
printf("xmax=%f\n", xmax);

then I get
xmax=0.250000

This is the beauty of C. Sometimes it is all-out warfare with this language and little tasks that should take 2 minutes take 2 days and end up with theory around little endians and data types and ANSI specifications and God knows what else :crying::crying:.


Edit: Category changed; see also this post #1. [Helmut]

Pass or fail!
ElMaestro
mittyri
★★  

Russia,
2018-03-15 11:55
(1408 d 10:23 ago)

@ ElMaestro
Posting: # 18552
Views: 5,637
 

 Nervous ticks in C - a reminder of the cruel logic of logic

Hi ElMaestro,

very strange
using your code
#include <stdio.h>
int Test10001(double x, double dx)
{
    double xmax;
    xmax=dx*ceil(x/dx);
    x=0.0;
    while (x<=xmax)
    {
        printf("x=%f\n",x);
        x=x+dx;
    }
    return(0);
}

int main()
{
    double x = 0.234567;
    double dx = 0.05;
    Test10001(x, dx);
    return(0);
}


in
https://www.onlinegdb.com/online_c_compiler

I got
x=0.000000
x=0.050000
x=0.100000
x=0.150000
x=0.200000
x=0.250000
as it should be
:confused:

Kind regards,
Mittyri
ElMaestro
★★★

Denmark,
2018-03-16 13:13
(1407 d 09:04 ago)

@ mittyri
Posting: # 18561
Views: 5,571
 

 Binary representation

Hi Mittyri,

» x=0.250000
» as it should be
» :confused:

How nice, then it is also compiler-dependent...
The issue is that 0.05 cannot be represented accurately in a binary brain.
So, the compiler will need to deal with such a fraction and either round it up or round it down. Whatever it does, it can have downstream implications. In my case it implied that my simple naïve decimal logic (0.25 <=0.25) never evaluated to TRUE because it was not tested as such.
The solution could be to add e.g. 1e-8 to the righthand side. Dumb, but that's logic.
Computers suck.
:-D:-D:-D

Pass or fail!
ElMaestro
mittyri
★★  

Russia,
2018-03-17 18:12
(1406 d 04:05 ago)

@ ElMaestro
Posting: # 18564
Views: 5,566
 

 bypassing eqaulity conditions

Hi ElMaestro,

never trust in '==' or '!=' for floating point numbers comparison
See for example here

Kind regards,
Mittyri
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-03-18 00:08
(1405 d 22:10 ago)

@ mittyri
Posting: # 18565
Views: 5,531
 

 And R?

Hi mittyri,

» never trust in '==' or '!=' for floating point numbers comparison
» See for example here

Great essay!

  I have stopped reading Stephen King novels.
Now I just read C code instead.
    Richard A. O’Keefe
  Being really good at C++ is like being really good
at using rocks to sharpen sticks.
    Thant Tessman


In R-lingo:

test.while <- function(x0, x1, dx) {
  x    <- double()
  x[1] <- x0
  j    <- 1L
  while (max(x) < x1) {
    j <- j + 1L
    x[j] <- x[j-1] + dx
  }
  return(x)
}

test.repeat1 <- function(x0, x1, dx) {
  x <- x0
  repeat {
    if (max(x) >= x1) break
    x <- c(x, max(x) + dx)
  }
  return(x)
}

test.repeat2 <- function(x0, x1, dx) {
  x    <- double(round((x1-x0)/dx, 0) + 1)
  x[1] <- x0
  j    <- 1L
  repeat {
    j    <- j + 1L
    if (x[j-1] >= x1) break
    x[j] <- x[j-1] + dx
  }
  return(x)
}

test.for <- function(x0, x1, dx) {
  x <- double(ceiling((x1-x0)/dx) + 1)
  for (j in seq_along(x)) {
    x[j] <- x0 + (j-1)*dx
  }
  return(x)
}

test.seq <- function(x0, x1, dx) {
  x <- seq(from=x0, to=dx*ceiling((x1-x0)/dx), by=dx)
  return(x)
}

test.pretty <- function(x0, x1, dx) {
  n <- ceiling((x1-x0)/dx) + 1
  x <- pretty(x=c(x0, x1), n=n)
  return(x)
}

x0 <- 0
x1 <- 0.234567
dx <- 0.05

test.while(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.repeat1(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.repeat2(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.for(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.seq(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.pretty(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25


Speed?

library(microbenchmark)
options(microbenchmark.unit="relative")
res <- microbenchmark(test.while(x0, x1, dx), test.repeat1(x0, x1, dx),
                      test.repeat2(x0, x1, dx), test.for(x0, x1, dx),
                      test.seq(x0, x1, dx), test.pretty(x0, x1, dx),
                      times=1000L, control=list("random", warmup=10))
print(res, digits=4)

Unit: relative
                     expr    min    lq  mean median  uq    max neval    cld
   test.while(x0, x1, dx)  2.333 2.249 2.233  2.222 2.1 0.9143  1000    d 
 test.repeat1(x0, x1, dx)  2.166 2.000 1.992  2.000 1.9 0.6714  1000   c   
 test.repeat2(x0, x1, dx)  1.333 1.375 1.342  1.333 1.3 0.3572  1000  b   
     test.for(x0, x1, dx)  1.000 1.000 1.000  1.000 1.0 1.0000  1000 a     
     test.seq(x0, x1, dx)  8.663 7.372 7.101  6.889 6.5 2.4857  1000     e
  test.pretty(x0, x1, dx) 10.662 8.746 8.393  8.111 7.6 3.5286  1000      f


Yes, we have a winner.
IMHO, seq() is most comfortable. One-liner: seq(0, 0.05*ceiling(0.234567/0.05), 0.05)
Good luck if you want nice ticks for a log-axis…

Dif-tor heh smusma 🖖
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
ElMaestro
★★★

Denmark,
2018-03-18 08:29
(1405 d 13:48 ago)

@ Helmut
Posting: # 18568
Views: 5,502
 

 And R?

Hi both,

in C the syntactic options are much more limited than in R. Fortunately :-D The solution I switched to is much like Hötzi's test.for case. This type of solution is in my opinon the "right" one in the sense that it doesn't do double comparisons and it can be also adapted to the case where the axis start is not at 0.0; for the purpose of getting some axis ticks where you want them this philosophy delivers and is safe.
There is e.g. no direct equivalent of the test.seq equivalent in C; it would require custom built helper functions or similar. Note also you actually need to do something with the ticks, so looping or applying would be necessary on top. Generating the ticks and using them in one and the same loop is likely fastest (but the latter is a postulate and not something I looked into in R). In C counting an integer down to 0 is often a bit faster, but not necessarily a lot. Highly settings- and compiler-dependent.
Speed is to me of no importance to me here but readibility is. This is part of an executable that does its job in 0.078 seconds and quits.:-D

Pass or fail!
ElMaestro
nobody
nothing

2018-03-18 13:13
(1405 d 09:04 ago)

@ ElMaestro
Posting: # 18569
Views: 5,485
 

 Nervous ticks in C - a reminder of the cruel logic of logic

» Hi all,
»
» sometimes logic totally sucks.
»
» I am making code in C for plotting various graphs,...

Problem identified. Simply use Graphpad Prism. The older versions came without copyright protection, more recent versions are PITA, connecting home once a month.

Kindest regards, nobody
ElMaestro
★★★

Denmark,
2018-03-19 09:44
(1404 d 12:33 ago)

@ nobody
Posting: # 18572
Views: 5,399
 

 Nervous ticks in C - a reminder of the cruel logic of logic

Hi Niemand,

» Simply use Graphpad Prism. The older versions came without copyright protection, more recent versions are PITA, connecting home once a month.

The executable produces a bunch of graphs in a splitsecond without intervention from the user.
Basically, it is a matter of just double-clicking an executable, waiting less than 0.1 seconds, and being happy about a document being generated and which contains about 50 graphs and various other info. All automatic. I bet Prism can't do that :-D
R can do it, too, but it is v-e-e-e-e-e-e-e-e-e-e-r-y slow in comparison and is a clumsy giant on comparison. My executable takes up just around 44kB of memory. R is more than 1000x that amount.

Pass or fail!
ElMaestro
nobody
nothing

2018-03-22 08:56
(1401 d 13:21 ago)

@ ElMaestro
Posting: # 18579
Views: 5,295
 

 Nervous ticks in C - a reminder of the cruel logic of logic

...but how many days have you already wasted for this millisecond code to do what you want it to do? Have you succeeded yet? :-D Just saying...

Prism you learn in 30 minutes and if you have an appropriate template your done in a few minutes to fine-tune your output.

If you have a hammer the world is full of nails. Unlucky if you have a sledgehammer and your nail is just 10 mm long.... ;-)

Kindest regards, nobody
ElMaestro
★★★

Denmark,
2018-03-22 17:16
(1401 d 05:01 ago)

(edited by ElMaestro on 2018-03-22 17:59)
@ nobody
Posting: # 18580
Views: 5,350
 

 Check it out for yourself - Dissolution Bootstrap:-)

Hi Nobody and everybody else,

to illustrate it please have a look at this archive. Unzip to its own folder, then drag and drop the csv file onto the executable. You may need to give green light if your security package protests or something. There is no virus in it (at least there wasn't when the programmer uploaded it), and it does not phone home.

It demonstrates what the programmer is doing with a bootstrap solution for dissolution data. It provides 1.000.000 bootstraps of two datasets and illustrates the results in the html file.
The hot stuff, whose detail I will not discuss here (that's for now solely a matter between me and someone in a public institution who likes to look at dissolution and to reject dossiers in murderous ways :-D ... I hope said person isn't reading this post), is the automatic generation of a series of graphs and plots and with an objective function, gamma. The user looks at these graphs, considers the gamma, and takes a decision. The user could in principle have a look at 1.000.000 such graphs and gammas (how about doing that in Prism, Nobody?).
The selection of curves to present is done in a pseudo-clever (actually pretty retarded, presently) fashion.

I hope this illustrates why a manual solution isn't viable (if not then then trust me, a manual solution is not viable. Not at all.)

For those who are into bootstraps, the result is at the end as a table, including raw bootstrap confidence intervals, Bias-corrected bootstrap interval, and bias-corrected and accelerated interval.

Hötzi: If you are reading this then I know you will absolutely torture my poor innocent little executable which is guilty of no fault :-D:-D:-D:-D:-D:-D. Yes, you can make it crash. Yes, the programmer (some wacko from Denmark) did not implement error trapping for whatever garbage the users fill in. No, this is not for the masses. Yes, the programmer has restricted the intended use to this dataset only.
By the way: The upload link function didn't work 2 minutes ago (it added bebac.at in the target string), so I pasted the raw link. Please forgive me.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-03-22 18:50
(1401 d 03:27 ago)

@ ElMaestro
Posting: # 18581
Views: 5,231
 

 Check it out for yourself - Dissolution Bootstrap:-)

Hi ElMaestro,

» Hötzi: If you are reading this then I know you will absolutely torture my poor innocent little executable which is guilty of no fault.

Why should I? Anyway, I have no win-OS with me (sitting in a Hotel on the other side of el Jeffe’s office waiting to be grilled tomorrow).

» By the way: The upload link function didn't work 2 minutes ago (it added bebac.at in the target string), so I pasted the raw link. Please forgive me.

Strange. I made a cosmetical change to your post.

Dif-tor heh smusma 🖖
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
ElMaestro
★★★

Denmark,
2018-03-22 19:04
(1401 d 03:13 ago)

@ Helmut
Posting: # 18582
Views: 5,292
 

 Check it out for yourself - Dissolution Bootstrap:-)

Hi Hötzi,

» Why should I? Anyway, I have no win-OS with me (sitting in a Hotel on the other side of el Jeffe's office waiting to be grilled tomorrow).

:-D:-D:-D I still recall your response to one of the previous executables I shared. :-D:-D:-D Forgot which one (Was it the Deep Sink project or something) - apart from disassembly it I think it was subjected to all sorts of challenges and weird settings and extremities.
Being away from office without windows....Og h the horror! I can't imagine what that would be like...
Good luck at the agency. I imagine you are well prepared and it will be a success.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-03-22 22:22
(1400 d 23:56 ago)

@ ElMaestro
Posting: # 18585
Views: 5,310
 

 Check it out for yourself - Dissolution Bootstrap:-)

Hi ElMaestro,

» I still recall your response to one of the previous executables I shared. [...] I think it was subjected to all sorts of challenges and weird settings and extremities.

Can’t remember but sounds familiar.

» Being away from office without windows....Og h the horror! I can't imagine what that would be like...

Simply great!

» Good luck at the agency. I imagine you are well prepared and it will be a success.

It’s a TSD (I’ll refer to your famous paper), parallel, n1 / group, GMR 0.90, high variability. Validated for extremely unequal group sizes and variance ratios of 0.25-4.0.

Dif-tor heh smusma 🖖
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Activity
 Admin contact
21,835 posts in 4,569 threads, 1,554 registered users;
online 6 (0 registered, 6 guests [including 4 identified bots]).
Forum time: Friday 22:18 CET (Europe/Vienna)

No problem can stand the assault of sustained thinking.    Voltaire

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5