R's char functions

Previously, I posted some functions for manipulating strings in R. I was surprised not to have found these in R. Obviously, I didn't look hard enough. RSiteSearch("strings") led me to a post mentioning the following:

> tolower("hI THerE, guYs!")
[1] "hi there, guys!"
> toupper("hI THerE, guYs!")
[1] "HI THERE, GUYS!"
> casefold("Hi")
[1] "hi"
> casefold("Hi", upper=TRUE)
[1] "HI"
> gsub("-", "", "what-what-what---?", fixed=TRUE)
[1] "whatwhatwhat?"

Well, I knew about gsub before at least.

S-PLUS's gsub doesn't recognize the fixed param but otherwise casefold and gsub work the same. (Of course, this is the point of having casefold in R.)

R special cases x^2 and x^0.5

Pulled from src/main/arithmetic.c:

static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall)
{
    int i, n;
    SEXP ans;

    switch (code) {
    case PLUSOP: return s1;
    case MINUSOP:
 ans = duplicate(s1);
 n = LENGTH(s1);
 for (i = 0; i < n; i++)
     REAL(ans)[i] = -REAL(s1)[i];
 return ans;
    default:
 errorcall(lcall, _("invalid unary operator"));
    }
    return s1; /* never used; to keep -Wall happy */
}

Would this be noticeably faster if we redefined the case for MINUSOP as

case MINUSOP:
  {
    const double* y1 = REAL(s1)[i];
    double* yAns;
    ans = duplicate(s1);
    yAns = REAL(ans);
    n = LENGTH(s1);
    for (i = 0; i < n; ++i)
      yAns[i] = -y1[i];
    return ans;
  }

?
We have the following defines:
Rinternals.h: #define REAL(x) ((double *) DATAPTR(x))
Rinternals.h: #define DATAPTR(x) (((SEXPREC_ALIGN *) (x)) + 1)
So, it seems like this could theoretically save 2n additions per negation. Woo! Is gcc smart enough to notice the legality of this move and take advantage of it?

Poking around in arithmetic.c also shows that x^2.0 and x^0.5 are special-cased to use more efficient alternatives than pow while there's another special routine R_pow_di to handle raising a double to an integer power. Which is used by pi^3? Setting breakpoints on R_pow and R_pow_di shows R_pow is used. In fact, pi^as.integer(3) also uses R_pow.

I wonder what the distribution of powers used in production R code looks like. 2 is almost certainly the most popular. Would it be worthwhile special casing other integer powers or even including a generic check for an integer power? Do most implementations of pow do this already?

an assortment of C

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define LINE_SIZE 10000

/**
 * Reads doubles from in until EOF.
 * If le n is not null, *len is set to
 * the number of elements read.
 * The caller is responsible for freeing the 
 * memory of the returned object.
 * Returns null in case of an error.
 */
double* readArray(FILE* in, int* len)
{
  const char* SPACES = " \t\r\n";
  char* line = malloc(LINE_SIZE);
  int capacity = 100;
  double* y = malloc(sizeof(double) * capacity);
  int n = 0;
  while (!ferror(in) && !feof(in)) {
    if (fgets(line, LINE_SIZE, in)) {
      const char* next = line;
      while (*next) {
 next += strspn(next, SPACES);
 if (*next) {
   double v;
   if (sscanf(next, "%lf", &v) == 1) {
     if (capacity == n) {
       capacity *= 2;
       y = realloc(y, capacity * sizeof(double));
     }

     y[n] = v;
     ++n;
   }
   else {
     free(line);
     free(y);
     return 0;
   }
 }
 next += strcspn(next, SPACES);
      }
    }
    else {
      if (!feof(in)) {
 free(line);
 free(y);
 return 0;
      }
    }
  }

  if (len) {
    *len = n;
  }

  return y;
}

/**
 * Reads doubles from file.
 * If len is not null, *len is set to the
 * number of elements read.
 * The caller is responsible for freeing the 
 * memory of the returned object.
 * Returns null in case of an error.
 */
double* readArrayFromFile(const char* file, int* len)
{
  FILE* in = fopen(file, "r");
  double* result = 0;
  if (!in) {
    fprintf(stderr, "Problem opening %s\n", file);
    perror("");
  }
  else {
    result = readArray(in, len);
    fclose(in);
  }
  return result;
}

/**
 * Returns 1 if tail is a suffix of s, 0 otherwise.
 */
int endsWith(const char* s, const char* tail)
{
  const int sLen = strlen(s);
  const int tailLen = strlen(tail);
  if (sLen >= tailLen) {
    return strcmp(s + sLen - tailLen, tail) == 0;
  }
  else {
    return 0;
  }
}

int remove(double* y, const int len,
const double toRemove)
{
  int newLen = 0;
  int i;
  
  for (i = 0; i < len; ++i) {
    const double yi = y[i];
    if (yi != toRemove) {
      y[newLen] = yi;
      ++newLen;
    }
  }
  return newLen;
}

Change the working directory in SAS

Just double click the status bar on the bottom of the main window giving the location of the current working directory!
Change the SAS Working Folder
This was not obvious to me.

These look possibly useful:
proc discrim [PDF]
proc princomp [PDF]

letter counts per paragraph

Count up the number of each letter by paragraph.

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>

int main(int numArgs, char** args)
{
  const int lineSize = 10000;

  if (numArgs != 2) {
    fprintf(stderr, "Usage: par_char_freq <text file>\n");
    exit(1);
  }
  else {
    char* line = malloc(lineSize);
    FILE* in = fopen(args[1], "r");
    int counts[26];
    int paragraphTotal = 0;

    if (!in) {
      fprintf(stderr, "Problem opening %s\n", args[1]);
      perror("");
      exit(1);
    }

    bzero(counts, 26 * sizeof(int));

    while (!feof(in) && !ferror(in)) {
      int lineTotal = 0;
      if (fgets(line, lineSize, in)) {
 int i;
 for (i = 0; line[i]; ++i) {
   const char c = tolower(line[i]);
   if ('a' <= c && c <= 'z') {
     ++counts[c - 'a'];
     ++lineTotal;
   }
 }
 
 if (lineTotal) {
   paragraphTotal += lineTotal;
 }
 else if (paragraphTotal) {
   int i;
   for (i = 0; i < 26; ++i) {
     printf("%d ", counts[i]);
   }
   printf("\n");

   paragraphTotal = 0;
   bzero(counts, 26 * sizeof(int));
 }
      }
    }

    fclose(in);
    free(line);
  }
  return 0;
}

The most permutable word is "trace"?

Which English word has the most permutations which are also English words? (Of course each of the lexible permutations is also an answer. Is "lexible" a sensible term for having the property of being a word?)

Here's some C++ code to look at this:

most_permutable.cpp:
#include <algorithm>
#include <cctype>
#include <fstream>
#include <functional>
#include <iostream>
#include <map>
#include <string>
using namespace std;

template<typename T>
class select2nd {
public:
  typedef T argument_type;
  
  typename T::second_type operator()(const T& t) {
    return t.second;
  }
};

template<typename B, typename U1, typename U2>
class binary_compose {
  public:
  typedef typename B::result_type result_type;
  typedef typename U1::argument_type arg1_type;
  typedef typename U2::argument_type arg2_type;
  binary_compose(const B& b, const U1& u1, const U2& u2) :
    b(b), u1(u1), u2(u2) {}

  result_type operator()(const arg1_type& x1,
    const arg2_type& x2)
  {
    return b(u1(x1), u2(x2));
  }

private:
  B b;
  U1 u1;
  U2 u2;
};

template<typename B, typename U1, typename U2>
binary_compose<B, U1, U2>
compose2(B b, U1 u1, U2 u2)
{
  return binary_compose<B, U1, U2>(b, u1, u2);
}

typedef map<string, int>::iterator iter;

int main(int numArgs, char** args)
{
  map<string, int> permutationCounts;
  ifstream in("/usr/share/dict/words");
  while (in) {
    string word;
    in >> word;
    if (!word.empty() && islower(word[0])) {
      sort(word.begin(), word.end());
      iter p = permutationCounts.find(word);
      if (p != permutationCounts.end()) {
 ++(*p).second;
      }
      else {
 permutationCounts.insert(make_pair(word, 1));
      }
    }
  }

  iter m = max_element(permutationCounts.begin(),
         permutationCounts.end(),
         compose2(less<int>(),
    select2nd<pair<string, int> >(),
    select2nd<pair<string, int> >()));
  cout << (*m).first << " with " << (*m).second << '\n';
}

Drum roll:

$ g++ -Wall most_permutable.cpp
$ ./a.out
acert with 9


What are the 9?

find_perm.cpp:
#include <algorithm>
#include <cstdio>
#include <fstream>
#include <iostream>
#include <string>
using namespace std;

int main(int numArgs, char** args)
{
  if (numArgs != 3) {
    cerr << "Usage: find_perm <file> <letters>\n";
    return 1;
  }

  ifstream in(args[1]);
  if (!in) {
    cerr << "Problem opening " << args[1] << '\n';
    perror("");
    return 1;
  }

  string letters(args[2]);
  sort(letters.begin(), letters.end());
  const string::size_type n = letters.size();

  string word;
  while (in) {
    in >> word;

    if (n == word.size()) {
      string key(word);

      sort(key.begin(), key.end());
      
      if (key == letters) {
 cout << word << '\n';
      }
    }
  }

  in.close();

}

$ g++ -Wall -o find_perm find_perm.cpp
$ ./find_perm oops
Usage: find_perm <file> <letters>
$ ./find_perm oops oops
Problem opening oops
No such file or directory
$ ./find_perm /usr/share/dict/words trace
caret
carte
cater
crate
creat
creta
react
recta
trace


What the heck is "creat"? Hm. Google doesn't know about it. It is really there though:
$ grep -w -n creat /usr/share/dict/words
45105:creat

strings again

getChars <- function(s) {
  n <- nchar(s)
  if (n > 0) substring(s, 1:n, 1:n) else character(0)
}

strip <- function(s, chars) {
  s.chars <- getChars(s)
  paste(s.chars[!(s.chars %in% chars)], collapse="")
}

tr <- function(s, from, to) {
  chars <- getChars(s)
  o <- match(chars, from)
  paste(ifelse(!is.na(o), to[o], chars), collapse="")
}

lower <- function(s) {
  tr(s, from=LETTERS, to=letters)
}

upper <- function(s) {
  tr(s, from=letters, to=LETTERS)
}

Another go: S-PLUS doesn't have strsplit so I use a different (and more efficient?) method for getting at the characters of a string.

> system.time(replicate(10000, strsplit("1234567890", "")[[1]]))
user system elapsed
0.102 0.004 0.105
> system.time(replicate(10000, substring("1234567890", 1:10, 1:10)))
user system elapsed
0.297 0.003 0.299

That's a surprise. Maybe I should try avoiding creating the index list twice? Still strsplit seems so much heavier.
//The source code for strsplit reveals that they make a special case of the pattern "". (See src/main/character.c.)
//Well, this is documented in the help page for strsplit as well.