| Andrew Cooke | Contents | Latest | RSS | Twitter | Previous | Next

C[omp]ute

Welcome to my blog, which was once a mailing list of the same name and is still generated by mail. Please reply via the "comment" links.

Always interested in offers/projects/new ideas. Eclectic experience in fields like: numerical computing; Python web; Java enterprise; functional languages; GPGPU; SQL databases; etc. Based in Santiago, Chile; telecommute worldwide. CV; email.

Personal Projects

Lepl parser for Python.

Colorless Green.

Photography around Santiago.

SVG experiment.

Professional Portfolio

Calibration of seismometers.

Data access via web services.

Cache rewrite.

Extending OpenSSH.

C-ORM: docs, API.

Last 100 entries

Crypto AG DID work with NSA / GCHQ; UNUMS (Universal Number Format); MOOCs (Massive Open Online Courses); Interesting Looking Game; Euler's Theorem for Polynomials; Weeks 3-6; Reddit Comment; Differential Cryptanalysis For Dummies; Japanese Graphic Design; Books To Be Re-Read; And Today I Learned Bugs Need Clear Examples; Factoring a 67 bit prime in your head; Islamic Geometric Art; Useful Julia Backtraces from Tasks; Nothing, however, is lost with less discomfort than that which, when lost, cannot be missed; Article on Didion; Cost of Living by City; British Slavery; Derrida on Metaphor; African SciFi; Traits in Julia; Alternative Japanese Lit; Pulic Key as Address (Snow); Why Information Grows; The Blindness Of The Chilean Elite; Some Victoriagate Links; This Is Why I Left StackOverflow; New TLS Implementation; Maths for Physicists; How I Am 8; 1000 Word Philosophy; Cyberpunk Reading List; Detailed Discussion of Message Dispatch in ParserCombinator Library for Julia; FizzBuzz in Julia w Dependent Types; kokko - Design Shop in Osaka; Summary of Greece, Currently; LLVM and GPUs; See Also; Schoolgirl Groyps (Maths); Japanese Lit; Another Example - Modular Arithmetic; Music from United; Python 2 and 3 compatible alternative.; Read Agatha Christie for the Plot; A Constructive Look at TempleOS; Music Thread w Many Recommendations; Fixed Version; A Useful Julia Macro To Define Equality And Hash; k3b cdrom access, OpenSuse 13.1; Week 2; From outside, the UK looks less than stellar; Huge Fonts in VirtualBox; Keen - Complex Emergencies; The Fallen of World War II; Some Spanish Fiction; Calling C From Fortran 95; Bjork DJ Set; Z3 Example With Python; Week 1; Useful Guide To Starting With IJulia; UK Election + Media; Review: Reinventing Organizations; Inline Assembly With Julia / LLVM; Against the definition of types; Dumb Crypto Paper; The Search For Quasi-Periodicity...; Is There An Alternative To Processing?; CARDIAC (CARDboard Illustrative Aid to Computation); The Bolivian Case Against Chile At The Hague; Clear, Cogent Economic Arguments For Immigration; A Program To Say If I Am Working; Decent Cards For Ill People; New Photo; Luksic And Barrick Gold; President Bachelet's Speech; Baltimore Primer; libxml2 Parsing Stream; configure.ac Recipe For Library Path; The Davalos Affair For Idiots; Not The Onion: Google Fireside Chat w Kissinger; Bicycle Wheels, Inertia, and Energy; Another Tax Fraud; Google's Borg; A Verion That Redirects To Local HTTP Server; Spanish Accents For Idiots; Aluminium Cans; Advice on Spray Painting; Female View of Online Chat From a Male; UX Reading List; S4 Subgroups - Geometric Interpretation; Fucking Email; The SQM Affair For Idiots; Using Kolmogorov Complexity; Oblique Strategies in bash; Curses Tools; Markov Chain Monte Carlo Without all the Bullshit; Email Para Matias Godoy Mercado; The Penta Affair For Idiots; Example Code To Create numpy Array in C; Good Article on Bias in Graphic Design (NYTimes); Do You Backup github?

© 2006-2015 Andrew Cooke (site) / post authors (content).

Uniformly Random, Correlated Numbers in Matlab/Octave

From: andrew cooke <andrew@...>

Date: Sat, 23 Jul 2011 11:27:19 -0400

Someone just asked me how to implement the work at
http://acooke.org/random.pdf in Matlab.  Here's the code for Octave, which is
very similar:


n = 1000;  # the number of points
t = 0.1; # the 't' parameter from the paper

sq2 = sqrt(2);

# first, generate points in the u,v coordinates as in figure 4.
# u goes from 0 to sq2, v from -t*sq2 to t*sq2
uv = rand(n,2) * [sq2,0;0,2*t*sq2] - repmat([0,t*sq2],n,1);
# now rotate that by 45 degrees to xy
xy = uv * [1/sq2,1/sq2;-1/sq2,1/sq2];
# finally, "fold in" points outside the unit square
for i = 1:n
  #printf("%d x %f y %f\n", i, xy(i,1), xy(i,2));
  if (xy(i,1) < 0) xy(i,1) = -xy(i,1); endif
  if (xy(i,2) < 0) xy(i,2) = -xy(i,2); endif
  if (xy(i,1) > 1) xy(i,1) = 2-xy(i,1); endif
  if (xy(i,2) > 1) xy(i,2) = 2-xy(i,2); endif
endfor;

#xy
plot(xy(:,1),xy(:,2),'*')


If anyone with more experience than me know of a better way of folding the
points (ie without the explicit "for" loop), please let me know!

Andrew

Avoiding For Loops

From: Jonathan Dursi <ljdursi@...>

Date: Sat, 23 Jul 2011 12:26:57 -0400

I think this is pretty easily done with the vector-friendly min & max 
operations:

function xy = randminmax(n)
     t = 0.1; # the 't' parameter from the paper

     sq2 = sqrt(2);

     # first, generate points in the u,v coordinates as in figure 4.
     # u goes from 0 to sq2, v from -t*sq2 to t*sq2
     uv = rand(n,2) * [sq2,0;0,2*t*sq2] - repmat([0,t*sq2],n,1);
     # now rotate that by 45 degrees to xy
     xy = uv * [1/sq2,1/sq2;-1/sq2,1/sq2];

     # finally, "fold in" points outside the unit square
     xy = max(xy,-xy);
     xy = min(xy,2-xy);
endfunction


t1=now(); xy=randfor(100000); t2=now(); t2-t1
ans =  1.8261e-05
t1=now(); xy2=randminmax(100000); t2=now(); t2-t1
ans =  1.4284e-07

-- 
Jonathan Dursi | SciNet, Compute/Calcul Canada

Re: Avoiding For Loops

From: andrew cooke <andrew@...>

Date: Sat, 23 Jul 2011 13:50:57 -0400

Thanks!  I was looking for something like a map - never ocurred to me to look
for specific functions.  I'll forward that to the person who asked.

(And even when I read your suggestion, before I looked at the code, I couldn't
see how the upper bound would work.  I'm clearly not getting the matrix-based
approach at all... guess I better practice!)

Andrew

Comment on this post