public class GenericStatic<T> {
static String o = "";
String oo() {
return o += "o";
}
public static void main(String[] args) {
GenericStatic<Object> a = new GenericStatic<Object>();
GenericStatic<String> b = new GenericStatic<String>();
System.out.println(a.oo());
System.out.println(b.oo());
}
}$ javac GenericStatic.java> (m <- matrix(1:4, nrow=2))This can cause problems when we want to go ahead and continue using the result in matrix operations:
[,1] [,2]
[1,] 1 3
[2,] 2 4
> (v <- m[1, 1:2])
[1] 1 3
> (v <- m[1, 1:2, drop=FALSE])
[,1] [,2]
[1,] 1 3
> v %*% m %*% t(v)
[,1]
[1,] 52
(%i8) display2d: false$
(%i9) %e^t;
(%o9) %e^t
(%i10) display2d: true$
(%i11) %e^t
;
t
(%o11) %e
#!/usr/bin/env pythonHere's a little game to drill you on somewhat similar sounding pairs of English words. I'm sure much could be done to improve it. Is the use of a class here kosher with Python practice? What's contrary to idiom here?
import os
import random
def say(word):
os.system("say " + word)
class Score:
def __init__(self):
self.total = 0
self.correct = 0
if __name__ == '__main__':
pairs = [["bed", "bad"], ["sped", "speed"],
["wary", "very"], ["bull", "ball"],
["mall", "mole"], ["pale", "pull"],
["soup", "stoop"], ["file", "fill"]]
scores = {}
for pair in pairs:
scores["".join(pair)] = Score()
input = ""
print("To quit, enter q")
while input != "q":
pair = random.choice(pairs)
word = random.choice(pair)
print("1) " + pair[0])
print("2) " + pair[1])
say(word)
input = raw_input("What did the computer say?: ")
input = input.strip()
if input != "q":
key = "".join(pair)
score = scores[key]
score.total += 1
if input == word:
score.correct += 1
print("Correct!")
else:
print("Wrong!")
print("You\'ve gotten the right answer for "
+ " vs ".join(pair) + " "
+ str(score.correct) + " out of "
+ str(score.total) + " times")
$ cat neg_dbl_max.c(This is running under Mac OS X 10.4 on an Intel Mac.) It'd be cool if gcc would recognize this constant and print it instead.
#include <float.h>
#include <stdio.h>
int main()
{
printf("%.50g\n", -DBL_MAX);
return 0;
}
$ gcc -Wall neg_dbl_max.c
$ ./a.out
-1.7976931348623157081452742373170435679807056752584e+308
(%i1) e : a^2 + b^2 + c^2;
2 2 2
(%o1) c + b + a
(%i2) ratsubst(r^2, a^2 + b^2 + c^2, e);
2
(%o2) r
(%i3) e : - 3 * a^2 - 3 * b^2 - 3 * c^2;
2 2 2
(%o3) - 3 c - 3 b - 3 a
(%i4) ratsubst(r^2, a^2 + b^2 + c^2, e);
2
(%o4) - 3 r
nan.c:$ gcc nan.c
#include <stdio.h>
#include <stdlib.h>
int main()
{
double x = strtod("NAN", 0);
printf("%f\n", x);
printf("%d\n", (int)(x == x));
return 0;
}
nan1.c:$ gcc nan1.c
#include <math.h>
#include <stdio.h>
int main()
{
printf("%f\n", NAN);
return 0;
}
vincent% cat nan2.c(I'm also working here under Linux. I'm not really sure what string I should drop in for the argument for nan. "hi!" also gives the same results as above. I'll stick with strtod for now.)
#include <math.h>
#include <stdio.h>
int main()
{
double x = nan("");
printf("%f\n", x);
printf("%d\n", (int)(x == x));
return 0;
}
vincent% gcc nan2.c -lm
vincent% ./a.out
0.000000
1
/**Earlier I had implemented this using only the second loop, now only used conditionally, with difference that I first rounded remaining * proportions[i] before truncating it to an int. However, this caused problems with too many points getting allocated, leading to a negative remainder and counts[0] getting an invalid value. Switching to simply truncating to zero, however, led us to sometimes not get exact results even when they were available. Eg, since (int) (97 * 0.5) = 48, the code under the if statement above divides 100 into (26, 49, 25) given proportions (0.25, 0.5, 0.25) rather (25, 50, 25) as one would hope. Hence, I have added the first loop and delegated the "safer" code to be called only in case of emergencies.
* On exit,
* counts[k] = number of points assigned to kth class,
* k = 0:(numClasses - 1) out of numPoints. Ensures that counts[k] > 0
* for every k while also sum(counts) = numPoints.
*/
void proportionsToCounts(const double* proportions,
int numClasses,
int numPoints,
int* counts)
{
int total = 0;
int numEmpty = 0;
int i;
assert(numPoints >= numClasses);
for (i = 0; i < numClasses; ++i) {
counts[i] = (int) round(numPoints * proportions[i]);
total += counts[i];
if (!counts[i]) {
++numEmpty;
}
}
if (numEmpty > 0 || total != numPoints) {
const double remaining = numPoints - numClasses;
int assigned = 0;
for (i = 1; i < numClasses; ++i) {
const int extra = (int) (remaining * proportions[i]);
assigned += extra;
counts[i] = 1 + extra;
}
counts[0] = 1 + numPoints - numClasses - assigned;
}
}
gibbs1d.h:
#ifndef __GIBBS1D_H__
#define __GIBBS1D_H__
#include <vector>
/**
* Produces a sequence of variates with P(x) propto exp(-beta * H(nodes))
* as its stationary distribution, where H(nodes) = the number of consecutive
* pairs with equal state.
*/
class Gibbs1d {
public:
/**
* beta -- larger values imply larger affinities for forming clusters.
* numNodes -- number of points along the line.
* q -- number of states for each point.
*/
Gibbs1d(double beta, int numNodes, int q = 2);
/**
* Update the state of each node by Gibbs sampling.
*/
void iterate();
/**
* Returns a pointer to a numNodes length array giving
* the current state of each node.
*/
const int* current() const;
private:
//Hide the copy ctor and assignment operator
Gibbs1d(const Gibbs1d&);
Gibbs1d& operator=(const Gibbs1d&);
const double beta;
const int numNodes;
const int q;
std::vector<int> states;
const double probNotEq1;
const double probNotEq2;
const double probEq;
const double prob1Out;
};
#endif
gibbs1d.cpp:I think there must be a much more elegant implementation of this. For one thing, the use of "sentinel" values isn't providing much of a payoff here. And the code for computing the probabilities is just very ugly. What's a better alternative to the nested if-statements?
#include <algorithm>
#include <cmath>
#include "gibbs1d.h"
#include "mtrandom.h"
using namespace std;
/**
* beta -- larger values imply larger affinities for forming clusters.
* numNodes -- number of points along the line.
* q -- number of states for each point.
*/
Gibbs1d::Gibbs1d(double beta, int numNodes, int q)
: beta(beta),
numNodes(numNodes),
q(q),
states(numNodes + 2),
probNotEq1(exp(beta) / (2.0 * exp(beta) + q - 2)),
probNotEq2(2.0 * probNotEq1),
probEq(exp(2.0 * beta) / (exp(2.0 * beta) + q - 1)),
prob1Out(exp(beta) / (exp(beta) + q - 1))
{
//Add sentinel values to the ends to avoid special casing
//the boundaries.
states[0] = states[numNodes + 1] = -1;
}
/**
* Returns an integer in the interval [0, q-1]
* != exclude, assigning equal probability to each
* allowed value.
* Assumes exclude is in [0, q-1].
*/
static int uniformNotEq(int q, int exclude)
{
int r = (int) ((q - 1) * genrand_real2());
if (r >= exclude) {
return r + 1;
}
else {
return r;
}
}
/**
* Returns an integer in the interval [0, q-1]
* != exclude1 or exclude2, assigning equal probability to each
* allowed value.
* Assumes exclude1 and exclude2 are in [0, q-1]
* and that exclude1 != exclude2.
*/
static int uniformNotEq(int q, int exclude1, int exclude2)
{
if (exclude1 > exclude2) {
swap(exclude1, exclude2);
}
int r = (int) ((q - 2) * genrand_real2());
if (r >= exclude1) {
r = r + 1;
}
if (r >= exclude2) {
r = r + 1;
}
return r;
}
/**
* Returns the given state with the given probability.
* Otherwise, chooses one of the other q - 1 states
* uniformly.
*/
static int withProbElseUniformNotEq(double prob, int state, int q)
{
const double u = genrand_real2();
if (u <= prob) {
return state;
}
else {
return uniformNotEq(q, state);
}
}
void Gibbs1d::iterate()
{
for (int i = 1; i <= numNodes; ++i) {
const int left = states[i - 1];
const int right = states[i + 1];
if (left >= 0) {
if (right >= 0) {
if (left != right) {
const double u = genrand_real2();
if (u <= probNotEq1) {
states[i] = left;
}
else if (u <= probNotEq2) {
states[i] = right;
}
else {
states[i] = uniformNotEq(q, left, right);
}
}
else { // left == right
states[i] = withProbElseUniformNotEq(probEq, left, q);
}
}
else { //right < 0
states[i] = withProbElseUniformNotEq(prob1Out, left, q);
}
}
else {
if (right >= 0) {
states[i] = withProbElseUniformNotEq(prob1Out, right, q);
}
else {
states[i] = (int) (q * genrand_real2());
}
}
}
}
/**
* Returns a pointer to a numNodes length array giving
* the current state of each node.
*/
const int* Gibbs1d::current() const
{
return &states[1];
}
(%i6) taylor(%e^n, n, 0, 6);What else does Maxima know?
2 3 4 5 6
n n n n n
(%o6)/T/ 1 + n + -- + -- + -- + --- + --- + . . .
2 6 24 120 720
(%i8) sum(k, k, 1, n), simpsum;
2
n + n
(%o8) ------
2
...
(%i10) simplify_sum(sum((-1)^(k+1) / k, k, 1, inf));
- 2 log(2) - %gamma %gamma
(%o10) - ------------------- - ------
2 2
(%i11) ratsimp(%);
(%o11) log(2)
(%i12) simplify_sum(sum(1/k^2, k, 1, inf));
2
%pi
(%o12) ----
6
(%i13) simplify_sum(sum(1/k^3, k, 1, inf));
(%o13) zeta(3)
...
(%i17) simplify_sum(sum(4^(-k) * z^(2*k)/(k!)^2, k, 0, inf));
inf
==== 2 k
\ z
(%o17) > ------
/ k 2
==== 4 k!
k = 0
... No Bessel functions?
(%i19) simplify_sum(sum(1/k^10, k, 1, inf));
10
%pi
(%o19) -----
93555
(%i20) simplify_sum(sum(1/k^50, k, 1, inf));
50
39604576419286371856998202 %pi
(%o20) ---------------------------------------------------
285258771457546764463363635252374414183254365234375
...
(%i22) simplify_sum(sum(1/(k^2 + 1), k, 1, inf));
%i psi (1 - %i) %i psi (%i + 1)
0 0
(%o22) --------------- - ---------------
2 2
... What is this?
#include <algorithm>
#include <tr1/hashtable>
#include <iostream>
#include <iterator>
int main()
{
Internal::X<0> x;
std::copy(x.primes, x.primes + x.n_primes,
std::ostream_iterator<unsigned long>(std::cout, "\n"));
}
Forget indenting. I like to just make my code all one line, with commands separated by ";".
#include <assert.h>I'm not completely happy with this design. I should probably have made "randomPartition" take the counts as an argument (although then it wouldn't've been so "random") since at least for my work, we're almost always interested in the counts (and will just have to recalculate them if they're not available). I wonder what standard implementations of this are available. I suppose I should be checking that the mallocs are succeeding as well. Are there any common platforms where accessing address 0 doesn't lead to a seg fault? Or what platforms don't segfault on accessing 0?
#include <math.h>
#include "partition.h"
#include <Rmath.h>
#include <stdlib.h>
/**
* Returns a array of length numClasses normalized to sum
* to 1.
*/
static double* randomProportions(int numClasses)
{
double* props = malloc(sizeof(double) * numClasses);
int i;
double total = 0.0;
for (i = 0; i < numClasses; ++i) {
props[i] = runif(0.0, 1.0);
total += props[i];
}
for (i = 0; i < numClasses; ++i) {
props[i] /= total;
}
return props;
}
/**
* Ensures at least one element is assigned to each
* class. The caller is responsible for freeing the
* returned object.
*/
int* randomCounts(int len, int numClasses)
{
int* counts = malloc(sizeof(int) * numClasses);
double* props = randomProportions(numClasses);
const double remaining = len - numClasses;
int total = 0;
int i;
for (i = 0; i < numClasses; ++i) {
counts[i] = 1;
}
for (i = 1; i < numClasses; ++i) {
const int extra = (int) round(remaining * props[i]);
total += extra;
counts[i] += extra;
}
counts[0] += remaining - total;
free(props);
return counts;
}
/**
* On exit, classes[i] = the class assigned to the kth class,
* k = 0:(numClasses - 1). Assumes classes has length at least len.
* Assumes that numClasses <= len. Assigns at least one element
* to each class. The classes are assigned in blocks, ie, they
* partition the indeces 0:(len - 1) into numClasses classes,
* ordered from class 0 and up.
*/
void randomPartition(int len, int numClasses, int* classes)
{
int* counts = randomCounts(len, numClasses);
int i;
int currentClass = 0;
assert(numClasses <= len);
for (i = 0; i < len; ++i) {
classes[i] = currentClass;
--counts[currentClass];
if (!counts[currentClass]) {
++currentClass;
}
}
free(counts);
}
case REALSXP:where t and s were previously initialized to 0.0. (Again I wonder if it would make a noticeable difference to cache the result of REAL(x).) So, in exact math, at the end of the loop t = sum(x) - sum(x) = 0. What does this accomplish in finite precision math?
PROTECT(ans = allocVector(REALSXP, 1));
for (i = 0; i < n; i++) s += REAL(x)[i];
s /= n;
if(R_FINITE((double)s)) {
for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
s += t/n;
}
REAL(ans)[0] = s;
break;
SEXP attribute_hidden do_recall(SEXP call, SEXP op, SEXP args, SEXP rho)Renaming a function screws up a recursive function defined in terms of its own name:
{
RCNTXT *cptr;
SEXP s, ans ;
cptr = R_GlobalContext;
/* get the args supplied */
while (cptr != NULL) {
if (cptr->callflag == CTXT_RETURN && cptr->cloenv == rho)
break;
cptr = cptr->nextcontext;
}
args = cptr->promargs;
/* get the env recall was called from */
s = R_GlobalContext->sysparent;
while (cptr != NULL) {
if (cptr->callflag == CTXT_RETURN && cptr->cloenv == s)
break;
cptr = cptr->nextcontext;
}
if (cptr == NULL)
error(_("'Recall' called from outside a closure"));
/* If the function has been recorded in the context, use it
otherwise search for it by name or evaluate the expression
originally used to get it.
*/
if (cptr->callfun != R_NilValue)
PROTECT(s = cptr->callfun);
else if( TYPEOF(CAR(cptr->call)) == SYMSXP)
PROTECT(s = findFun(CAR(cptr->call), cptr->sysparent));
else
PROTECT(s = eval(CAR(cptr->call), cptr->sysparent));
ans = applyClosure(cptr->call, s, args, cptr->sysparent, R_BaseEnv);
UNPROTECT(1);
return ans;
}
> embed(1:5, dim=3)> as.integer(intToBits(as.integer(3)))
[,1] [,2] [,3]
[1,] 3 2 1
[2,] 4 3 2
[3,] 5 4 3
> embed(1:5, dim=2)
[,1] [,2]
[1,] 2 1
[2,] 3 2
[3,] 4 3
[4,] 5 4
intToBitString <- function(x, bits=32) {
substring(paste(rev(as.integer(intToBits(as.integer(x)))),
collapse=""),
32 - bits + 1, 32)
}
doubleToBitString <- function(x) {
stopifnot(is.finite(x))
e <- floor(log2(x))
m <- (abs(x * 2 ^ -e) - 1)
shift <- 8
s.bits <- ""
remaining <- .Machine$double.digits
while (remaining > 0) {
u <- min(remaining, shift)
m <- m * 2 ^ u
k <- floor(m)
s.bits <- paste(s.bits, intToBitString(k, u), sep="")
m <- m - k
remaining <- remaining - u
}
mexp <- .Machine$double.exponent
paste(if (x > 0) 0 else 1,
intToBitString(e + 2^(mexp - 1) - 1, mexp),
s.bits)
}(We could easily avoid depending on the base being 2, but is it worth it?) The output for pi seems to match that given for pi on Wikipedia page.create.so <- function() {
writeLines(c("void doubleToBitString(const double* x,"
" int* length, char** result)",
"{",
" const char* bytes = (const void*)x;",
" int numBytes = sizeof(double);",
" static char buffer[1024];",
" int i;",
" for (i = 0; i != numBytes; ++i) {",
" int j;",
" for (j = 0; j != 8; ++j) {",
" int mask = 1 << j;",
" buffer[8*i + j] = (mask & bytes[i]) ? \'1\' : \'0\';",
" }",
" }",
" buffer[8*numBytes] = \'\\0\';",
" *result = buffer;",
" *length = numBytes;",
"}"),
con = "bitstring.c")
system("R CMD SHLIB bitstring.c", ignore.stderr=T)
}
doubleToBitString <- function(x) {
if (!file.exists("bitstring.so")) {
create.so()
}
dyn.load("bitstring.so")
out <- .C("doubleToBitString",
x=as.double(x),
len=as.integer(1),
bits=as.character('junk'))
return(out$bits)
}Whoa. Does this write into memory it shouldn't be touching?
Call:
survreg(formula = Surv(y, z) ~ x + frailty.gaussian(v), data = d)
coef se(coef) se2 Chisq DF p
(Intercept) 1.20 0.254 0.156 22.3 1.0 2.4e-06
x 2.96 0.166 0.162 319.3 1.0 0.0e+00
frailty.gaussian(v) 176.6 63.2 1.2e-12
Scale= 4.7
Iterations: 7 outer, 23 Newton-Raphson
Variance of random effect= 4
Degrees of freedom for terms= 0.4 1.0 63.2 1.0
Likelihood ratio test=463 on 63.5 df, p=0 n= 1000
#include <stdlib.h>
/**
* 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)
{
int capacity = 100;
double* y = malloc(sizeof(double) * capacity);
int n = 0;
while (!ferror(in) && !feof(in)) {
double v;
int status = fscanf(in, "%lf", &v);
if (status == 1) {
if (capacity == n) {
capacity *= 2;
y = realloc(y, capacity * sizeof(double));
}
y[n] = v;
++n;
}
else if (status == 0) {
fprintf(stderr,
"# Invalid input? (read %d values so far)\n",
n);
free(y);
return 0;
}
else if (status == EOF && !feof(in)) {
fprintf(stderr,
"# Got error (read %d values so far)\n", n);
perror("");
free(y);
return 0;
}
}
if (len) {
*len = n;
}
return y;
}
mset.data <- function(width=200, iterations=100) {
mset <- as.vector(outer(seq(-2, 2, len=width),
seq(-2, 2, len=width) * 1i, FUN="+"))
cmset <- NULL
z <- numeric(length(mset))
for (i in 1:iterations) {
z <- z^2 + mset
outside <- Mod(z) > 2
cmset <- append(cmset, mset[outside])
inside <- !outside
mset <- mset[inside]
z <- z[inside]
}
all <- c(cmset, mset)
all.next <- all + all^2
inside <- rep(c(0, 1), c(length(cmset), length(mset)))
x <- Re(all)
y <- Im(all)
data.frame(inside=inside, x1=Re(all), y1=Im(all),
x2=Re(all.next), y2=Im(all.next))
}
mset.prob <- function(width=200) {
d <- mset.data(width=width)
f <- glm(inside ~ I(x1^2) + I(y1^2) +
x1*y1 + I(x2^2) + I(y2^2) + x2*y2,
data=d, family=binomial())
z1 <- as.vector(outer(seq(-2, 2, len=width),
seq(-2, 2, len=width) * 1i, FUN="+"))
z2 <- z1 + z1^2
p <- predict(f, type="response",
newdata=data.frame(x1=Re(z1), y1=Im(z1),
x2=Re(z2), y2=Im(z2)))
matrix(p, nrow=width, ncol=width)
}$ cat int_div_by_zero.cKaboom! Nothing. This result immediately raises the question, "Is integer math division as slow as or slower than floating point division since it seems to be implemented in terms of it?"
#include <stdio.h>
int main()
{
printf("%d\n", 3 / 0);
return 0;
}
$ gcc int_div_by_zero.c
int_div_by_zero.c: In function 'main':
int_div_by_zero.c:5: warning: division by zero
$ ./a.out
Floating point exception
##Generate n variates from the standard (mu=0, sigma=1)> x <- rnorm(10000)
##Smallest Extreme Value distribution.
rsev <- function(n) log(-log(runif(n)))
##g -- RNG for the underlying standard distribution
##for the log-location-scale family.
##Defaults to rsev to produce Weibull variates.
gen.data <- function(x, beta0=3, beta1=9, sigma=4, g=rsev) {
n <- length(x)
exp(beta0 + beta1 * x + sigma * g(n))
}
Could be quite a bit more efficient. It should work though.
/**
* Compute a cumulative sum of y, ie,
* result[j] = sum(i=0:j) y[i].
* The caller is responsible for freeing the returned
* object.
*/
double* normedCumSum(const double* y, int len)
{
double* result = malloc(len * sizeof(double));
if (len > 0) {
double norm;
int i;
result[0] = y[0];
for (i = 1; i < len; ++i) {
result[i] = result[i - 1] + y[i];
}
norm = 1.0 / result[len - 1];
for (i = 0; i < len; ++i) {
result[i] *= norm;
}
}
return result;
}
/**
* Sample reps elements from y (with replacement) with the given
* weights w. The caller is responsible for freeing the returned
* object.
*/
double* sample(const double* y, const double* w, int len, int reps)
{
double* result = malloc(reps * sizeof(double));
double* probs = normedCumSum(w, len);
int i;
for (i = 0; i < reps; ++i) {
int j;
const double u = runif(0.0, 1.0);
for (j = 0; j < len; ++j) {
if (u <= probs[j]) {
result[i] = y[j];
break;
}
}
}
free(probs);
return result;
}
int main()
{
double y[2] = {0, 1};
double w[2] = {0.6, 0.4};
double* x = sample(y, w, 2, 1000);
double sum = 0.0;
int i;
for (i = 0; i < 1000; ++i) {
sum += x[i];
}
printf("mean: %g\n", sum / 1000);
return 0;
}
RANDSAMPLE Random sample, with or without replacement.I'm still learning how to use the help system. It's pretty silly, but I found this via Google search rather than via the internal search for the program. ("sample", "random", and "random sample" all returned relatively huge lists of results, with the one I was interested settled down in the middle. Ah, putting quote marks around "random sample" brings randsample to near the top.)
Y = RANDSAMPLE(N,K) returns Y as a 1-by-K vector of values sampled
uniformly at random, without replacement, from the integers 1:N.
Y = RANDSAMPLE(POPULATION,K) returns K values sampled uniformly at
random, without replacement, from the values in the vector POPULATION.
Y = RANDSAMPLE(...,REPLACE) returns a sample taken with replacement if
REPLACE is true, or without replacement if REPLACE is false (the default).
Y = RANDSAMPLE(...,true,W) returns a weighted sample, using positive
weights W, taken with replacement. W is often a vector of probabilities.
This function does not support weighted sampling without replacement.
Example: Generate a random sequence of the characters ACGT, with
replacement, according to specified probabilities.
R = randsample('ACGT',48,true,[0.15 0.35 0.35 0.15])
See also rand, randperm.
row_major.c:
#include <stdio.h>
int main() {
char a[4][3][2];
int i;
for (i = 0; i < 4; ++i) {
int j;
for (j = 0; j < 3; ++j) {
int k;
for (k = 0; k < 2; ++k) {
printf("(%d, %d, %d) ==> %d\n", i, j, k,
(int)(&a[i][j][k] - &a[0][0][0]));
}
}
}
return 0;
}
charClass2.c:$ gcc -Wall charClass2.c
#include <stdio.h>
#include <string.h>
int main(int numArgs, char** args)
{
int i;
for (i = 0; i < numArgs; ++i) {
char x[3];
char y[1000];
if (sscanf(args[i], "%2[a]%999[a]", x, y) == 2) {
printf("%zu %zu\n", strlen(x), strlen(y));
}
else {
printf(":-(\n");
}
}
return 0;
}
$ ghciThis looks like fun.
___ ___ _
/ _ \ /\ /\/ __(_)
/ /_\// /_/ / / | | GHC Interactive, version 6.6, for Haskell 98.
/ /_\\/ __ / /___| | http://www.haskell.org/ghc/
\____/\/ /_/\____/|_| Type :? for help.
Loading package base ... linking ... done.
Prelude> :l stripHtml.hs
[1 of 1] Compiling Main ( stripHtml.hs, interpreted )
Ok, modules loaded: Main.
*Main> stripHtml "<b><blink>What's</blink></b><h1> a better</h1> way to d<tt>o this</tt>?"
"What's a better way to do this?"
#include <stdio.h>
#include <string.h>
/**
* If str is of the form "<m>:<n>",
* set *min = m
* and *max = n. Otherwise,
* if str if of the form "<n>",
* set *min = defaultMin and
* *max = n.
* Return the number of numbers
* successfully parsed (1 or 2 indicating success,
* 0 failure).
*/
int charToIntRange(const char* str,
int defaultMin,
int* min,
int* max)
{
if (strchr(str, ':')) {
if (sscanf(str, "%d:%d", min, max) == 2) {
return 2;
}
else {
return 0;
}
}
else {
*min = defaultMin;
return sscanf(str, "%d", max);
}
}
int main(int numArgs, char** args)
{
int i;
for (i = 1; i < numArgs; ++i) {
int m;
int n;
if (charToIntRange(args[i], -1666, &m, &n)) {
printf("%d:%d\n", m, n);
}
else {
printf("failed to parse <<%s>>\n", args[i]);
}
}
return 0;
}
##Rotate a square matrix 90 degrees clockwise.
rot90 <- function(a) {
n <- dim(a)[1]
stopifnot(n==dim(a)[2])
t(a[n:1, ])
}
f = fopen('test.txt', 'w');
fprintf(f, '%g\n', 1:50);
fclose(f);(I'm sure there's a better way to do this.)
> d <- as.matrix(dist(1:4, diag=T, upper=T))Is it surprising to find the sample variance such a relatively poor approximation to the true variance even with 10000 observations? (I guess the answer is "no!")
> s <- 0.3^d
> a <- chol(s)
> m <- matrix(nrow=10000,ncol=4)
> for (i in 1:10000) { x <- rnorm(4); m[i,] <- a %*% x}
> var(m)
[,1] [,2] [,3] [,4]
[1,] 1.09839578 0.30672920 0.1062927 0.03032918
[2,] 0.30672920 1.01055057 0.3031668 0.09055224
[3,] 0.10629268 0.30316682 1.0044863 0.27050059
[4,] 0.03032918 0.09055224 0.2705006 0.88942696
> s
1 2 3 4
1 1.000 0.30 0.09 0.027
2 0.300 1.00 0.30 0.090
3 0.090 0.30 1.00 0.300
4 0.027 0.09 0.30 1.000
##Produce elements that are marginally Rician distributed but for
##which the noise between elements is correlated with the
##structure cor(X_i, X_j) = cor(Y_i, Y_j) = rho^abs(i-j)
##but {X_i} independent of {Y_i}, i = 1:length(mu) and
##result[i] = sqrt(X_i^2 + Y_i^2).
##The value is a matrix of draws from the given distribution,
##each row a single draw so that the ith column corresponds
##to mu[i]
corRice <- function(mu, n, rho, sigma) {
k <- length(mu)
d <- as.matrix(dist(1:k, diag=TRUE, upper=TRUE))
r <- (rho ^ d)
a <- sigma * chol(r)
result <- matrix(nrow=n, ncol=k)
for (i in 1:n) {
x <- (a %*% rnorm(k)) + mu
y <- a %*% rnorm(k)
result[i, ] <- sqrt(x^2 + y^2)
}
result
}
static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall)Would this be noticeably faster if we redefined the case for MINUSOP as
{
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 */
}
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;
}
#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;
}
#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;
}
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';
}
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();
}
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.> library(SASmixed)
> names(Semiconductor)
> str(Semiconductor)
'data.frame': 48 obs. of 5 variables:
$ resistance: num 5.22 5.61 6.11 6.33 6.13 6.14 5.6 5.91 5.49 4.6 ...
$ ET : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$ Wafer : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
$ position : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
$ Grp : Factor w/ 12 levels "1/1","1/2","1/3",..: 1 1 1 1 2 2 2 2 3 3 ...
...
> options(contrasts=c("contr.treatment", "contr.poly"))
> lm1 <- lm(resistance ~ ET * position, data=Semiconductor)
> summary(lm1)
Call:
lm(formula = resistance ~ ET * position, data = Semiconductor)
Residuals:
Min 1Q Median 3Q Max
-1.01333 -0.25750 0.04333 0.28333 0.74667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.61333 0.26891 20.874 <2e-16 ***
ET2 0.38000 0.38030 0.999 0.325
ET3 0.52333 0.38030 1.376 0.178
ET4 0.72667 0.38030 1.911 0.065 .
position2 -0.16333 0.38030 -0.429 0.670
position3 -0.06000 0.38030 -0.158 0.876
position4 0.27333 0.38030 0.719 0.478
ET2:position2 0.35667 0.53782 0.663 0.512
ET3:position2 0.37333 0.53782 0.694 0.493
ET4:position2 0.37667 0.53782 0.700 0.489
ET2:position3 -0.16667 0.53782 -0.310 0.759
ET3:position3 -0.30333 0.53782 -0.564 0.577
ET4:position3 -0.38333 0.53782 -0.713 0.481
ET2:position4 -0.35000 0.53782 -0.651 0.520
ET3:position4 -0.31667 0.53782 -0.589 0.560
ET4:position4 -0.07333 0.53782 -0.136 0.892
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4658 on 32 degrees of freedom
Multiple R-squared: 0.4211, Adjusted R-squared: 0.1498
F-statistic: 1.552 on 15 and 32 DF, p-value: 0.1449
> options(contrasts=c("contr.helmert", "contr.poly"))
> lm2 <- lm(resistance ~ ET * position, data=Semiconductor)
> summary(lm2)
Call:
lm(formula = resistance ~ ET * position, data = Semiconductor)
Residuals:
Min 1Q Median 3Q Max
-1.01333 -0.25750 0.04333 0.28333 0.74667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.00292 0.06723 89.292 < 2e-16 ***
ET1 0.17000 0.09507 1.788 0.08324 .
ET2 0.09722 0.05489 1.771 0.08606 .
ET3 0.10986 0.03881 2.830 0.00797 **
position1 0.05667 0.09507 0.596 0.55535
position2 -0.11000 0.05489 -2.004 0.05360 .
position3 0.03542 0.03881 0.912 0.36834
ET1:position1 0.08917 0.13446 0.663 0.51197
ET2:position1 0.03250 0.07763 0.419 0.67826
ET3:position1 0.01667 0.05489 0.304 0.76337
ET1:position2 -0.05750 0.07763 -0.741 0.46427
ET2:position2 -0.03528 0.04482 -0.787 0.43700
ET3:position2 -0.02444 0.03169 -0.771 0.44617
ET1:position3 -0.05167 0.05489 -0.941 0.35363
ET2:position3 -0.01111 0.03169 -0.351 0.72818
ET3:position3 0.01125 0.02241 0.502 0.61909
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4658 on 32 degrees of freedom
Multiple R-squared: 0.4211, Adjusted R-squared: 0.1498
F-statistic: 1.552 on 15 and 32 DF, p-value: 0.1449
> options(contrasts=c("contr.treatment", "contr.poly"))
> lmr1 <- lm(resistance ~ ET + position, data=Semiconductor)
> anova(lm1, lmr1)
Analysis of Variance Table
Model 1: resistance ~ ET * position
Model 2: resistance ~ ET + position
Res.Df RSS Df Sum of Sq F Pr(>F)
1 32 6.9421
2 41 7.7515 -9 -0.8095 0.4146 0.9178
> options(contrasts=c("contr.helmert", "contr.poly"))
> lmr2 <- lm(resistance ~ ET + position, data=Semiconductor)
> anova(lm2, lmr2)
Analysis of Variance Table
Model 1: resistance ~ ET * position
Model 2: resistance ~ ET + position
Res.Df RSS Df Sum of Sq F Pr(>F)
1 32 6.9421
2 41 7.7515 -9 -0.8095 0.4146 0.9178> a <- as.factor(rep(1:3, each=4))
> b <- as.factor(rep(1:2, each=6))
> y <- rnorm(12)
> options(contrasts=c("contr.treatment", "contr.poly"))
> model.matrix(y ~ a*b)
(Intercept) a2 a3 b2 a2:b2 a3:b2
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 0 0
4 1 0 0 0 0 0
5 1 1 0 0 0 0
6 1 1 0 0 0 0
7 1 1 0 1 1 0
8 1 1 0 1 1 0
9 1 0 1 1 0 1
10 1 0 1 1 0 1
11 1 0 1 1 0 1
12 1 0 1 1 0 1
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
> options(contrasts=c("contr.helmert", "contr.poly"))
> model.matrix(y ~ a*b)
(Intercept) a1 a2 b1 a1:b1 a2:b1
1 1 -1 -1 -1 1 1
2 1 -1 -1 -1 1 1
3 1 -1 -1 -1 1 1
4 1 -1 -1 -1 1 1
5 1 1 -1 -1 -1 1
6 1 1 -1 -1 -1 1
7 1 1 -1 1 1 -1
8 1 1 -1 1 1 -1
9 1 0 2 1 0 2
10 1 0 2 1 0 2
11 1 0 2 1 0 2
12 1 0 2 1 0 2
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.helmert"
attr(,"contrasts")$b
[1] "contr.helmert"
> b <- as.ordered(b)
> model.matrix(y ~ a*b)
(Intercept) a1 a2 b.L a1:b.L a2:b.L
1 1 -1 -1 -0.7071068 0.7071068 0.7071068
2 1 -1 -1 -0.7071068 0.7071068 0.7071068
3 1 -1 -1 -0.7071068 0.7071068 0.7071068
4 1 -1 -1 -0.7071068 0.7071068 0.7071068
5 1 1 -1 -0.7071068 -0.7071068 0.7071068
6 1 1 -1 -0.7071068 -0.7071068 0.7071068
7 1 1 -1 0.7071068 0.7071068 -0.7071068
8 1 1 -1 0.7071068 0.7071068 -0.7071068
9 1 0 2 0.7071068 0.0000000 1.4142136
10 1 0 2 0.7071068 0.0000000 1.4142136
11 1 0 2 0.7071068 0.0000000 1.4142136
12 1 0 2 0.7071068 0.0000000 1.4142136
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.helmert"
attr(,"contrasts")$b
[1] "contr.poly"
##Find the correlation between each pair of columns
##in d != base conditional on base being in the
##corresponding subset.
##
##d -- 3 col data.frame or matrix
##subset -- list of vectors
##base -- column to use in conditioning
cor.on.subsets <- function(d, subsets, base=2) {
stopifnot(dim(d)[2] == 3)
stopifnot(base %in% 1:3)
x <- setdiff(1:3, base)
n <- length(subsets)
r <- numeric(n)
for (i in 1:n) {
d1 <- d[d[, base] %in% subsets[[i]], ]
r[i] <- cor(d1[,x[1]], d1[, x[2]])
}
r
}
##Plot a qq-plot of the sq Mahalanobis distanceI learned about a couple of R functions in writing this: mahalanobis and colMeans. I wonder what else I'm missing. (On a hunch, I tried help.search("mahalanobis").)
##vs. the corresponding chi-sq quantiles.
##Return a p-value for the sq. corr. efficient
##between the two based on simulation.
chisq.qqplot <- function(a.matrix, nsamples=10000) {
d.sq <- mahalanobis(a.matrix, colMeans(a.matrix), cov(a.matrix))
d.sq <- sort(d.sq)
rows <- dim(a.matrix)[1]
cols <- dim(a.matrix)[2]
q <- qchisq(((1:rows) - 0.5) / rows, df=cols)
plot(q, d.sq, xlab="theoretical", ylab="sq. distance")
qq.cor <- cor(d.sq, q)
sim.cors <- numeric(nsamples)
for (i in 1:nsamples) {
x <- matrix(rnorm(cols * rows), ncol=cols)
d2 <- mahalanobis(x, colMeans(x), cov(x))
d2 <- sort(d2)
sim.cors[i] <- cor(q, d2)
}
list(d.sq = d.sq,
qq.cor=qq.cor,
qq.cor.sq=qq.cor^2,
sim.cors=sim.cors,
p.value=sum(sim.cors^2 < qq.cor^2) / nsamples)
}
printMatrix <- function(m, r=2) {
cat("\\left(\\begin{array}{",
rep("c", nrow(m)), "}\n", sep="")
cat(apply(m, 1,
function(row) paste(round(row, r),
collapse=" & ")),
sep=" \\\\\n");
cat("\\end{array}\\right)\n");
}
foo <- function(s3, r=2) {
s <- matrix(c(s3[1], s3[2], s3[2], s3[3]), ncol=2);
cat("$\n")
printMatrix(s, r=r)
cat("$\n\n");
cat("cor: $", cor(s)[1,2], "$\n\n", sep="");
cat("generalized variance: $", det(s), "$\n\n", sep="");
cat("total variance: $", sum(diag(s)), "$\n", sep="");
}
printMatrix(a) := block([as, body],There's just got to be a function in Maxima somewhere that already does this. (There's also just got to be a function that creates a list based on repeating an element. makelist(x, i, 1, n) does the job though.)
as : matrixmap(lambda([x],
printf(false, "~a", x)), a),
body: simplode(maplist(lambda([row],
simplode(row, " & ")),
as),
sconc(" ", "\\\\", newline)),
sconc("\\left(\\begin{array}{",
simplode(makelist("r", i, 1,
length(a[1]))), "}",
newline,
body,
newline,
"\\end{array}\\right)"))$
#include <iostream>$ g++ array_init_every_time.cpp
using namespace std;
int main()
{
for (int i = 0; i < 3; ++i) {
int a[1] = {13};
cout << a[0] << endl;
}
return 0;
}
R:Okay. That's not what I was expecting.
> a <- "Yes"
> f <- function() { return(a) }
> a <- "No"
> f()
[1] "No"
Scheme:Hm. But (in R):
> (define a "Yes")
> (define f (lambda () a))
> (define a "No")
> (f)
"No"
printMatrix(a) := block([as],Naturally, again I have to wonder very strongly if there isn't already a better way to do this.
as : matrixmap(lambda([x],
printf(false, "~a", x)), a),
simplode(maplist(lambda([row],
simplode(row, " & ")),
as),
sconc(" ", "\\\\", newline)))$
diag(d) := block([n], n : length(d),
genmatrix(lambda([i, j],
if i = j then d[i] else 0),
n, n))$
(%i90) diag_matrix(1,2,3);
[ 1 0 0 ]
[ ]
(%o90) [ 0 2 0 ]
[ ]
[ 0 0 3 ]
