Sample Code in C++ and Ruby

I offer some idiomatic codes often used in scientific world. It is crucial for us to write a code which can be extended to more general purpose. Here I focus on simplicity rather than generality. C++ offers powerful structure but its strict grammar sometimes halt us. Therefore stand-alone sample should accelerate development speed. These codes work without special libraries except smart pointers contained in std::tr1 who deserves native support in the next generation standard C++0x.


PBC in 1D and compact field variables

Generic functions such as iterator on a lattice are not offered. If you need this for large N, then shared_ptr for each lattice point is preferable.

This Turing Machine like example is for Boolean valued state on fixed length ring. But you can easily generalize it for (compact) integer or float valued state models i.e. partial differential equations (PDEs).

#include <iostream>
#include <vector>
#include <bitset>
#include <valarray>
#include <cmath>
using namespace std;

// only works for T = signed int
// signed double
template <typename T>
class universal_mod {
const T m;
public:
universal_mod(const T& base): m(base) {}
T operator()(T x) const {
x -= m*int(x/m);
x += m*(x<0);
return x;
}
};

template <typename T, size_t N>
class Ring {
valarray<T> v;
universal_mod<int> ring;

friend ostream& operator<<(ostream& os, Ring& item) {
for (int i=0; i!=item.size(); ++i)
os << (item[i]? "#": "-") << " ";
return os;
}
public:
Ring(): v(T(),N),ring(N) {}
T& operator[](int i) {
return v[ring(i)];
}
Ring& operator=(const Ring& rhs) { v=rhs.v; return *this; }
int size() const { return N; }
};

// stochastic engine
class Glauber {
const double beta;
double prob(const double& dE) const {
return (1-tanh(beta*dE/2))/2.0;
}
public:
Glauber(double T): beta(1.0/T) {}
bool binary_signal(const double& dE) const {
const double p = prob(dE);
return rand()/(1.0+RAND_MAX) < p;
}
} engine(4*atan(1.0));

// particle tossing machine
class Machine {
bool state; // ACTIVE=true, SLEEP=false
public:
Machine(): state(true) {}
bool is_active() const { return state; }
void excite() { state = true; }
void sleep() { state = false; }
} machine;

int main()
{
srand(time(NULL));

const int N = 40;
Ring<bool,N> f;
for (int i=0; i!=N/2; ++i)
f[i]=true;

for (int t=0; t!=0xFF; ++t) {
cout << f << endl;
// Begin Sweep
int i=N; while (i!=0) {
--i;
if (f[i] && !f[i+1] && machine.is_active() && engine.binary_signal(1.0)) {
f[i] = false;
f[i+1] = true;
machine.sleep();
//cout << "Hopped from " << i << " to " << i+1 << endl;
}
else
machine.excite();
} // End Sweep
}
return 0;
}

lessons we learned


Observer Pattern

Observer Pattern appears very often in scientific OOP. Physical relation between a system and its observer is easy to recognize in this context. However its implementation takes unnatural way. I offer here a minimal translation using a list of weak pointers.

#include<iostream>
#include<list>
#include<tr1/memory> // for Smart Pointers
using namespace std;
using namespace tr1;

class Observer {
public:
Observer() {}
void update() {
cout << "Signal Detected!" << endl;
}
};

class Subject {
typedef weak_ptr<Observer> WPtr;
list<WPtr> obs;
class Eq: unary_function<WPtr,bool> {
WPtr lhs;
public:
Eq(WPtr wp): lhs(wp) {}
bool operator()(WPtr rhs) {
return lhs.lock() == rhs.lock();
}
};
public:
Subject(): obs(0) {}
void attach(WPtr o) { obs.push_back(o); }
void detach(WPtr o) { obs.remove_if(Eq(o)); }
void notify() {
list<WPtr>::iterator oiter = obs.begin();
while (oiter != obs.end()) {
if (shared_ptr<Observer> tmp = oiter->lock())
tmp->update();
else
cout << "No More Exist!" << endl;
++oiter;
}
}
};

int main()
{
Subject s;
shared_ptr<Observer> o0(new Observer()); s.attach(o0);
shared_ptr<Observer> o1(new Observer()); s.attach(o1);
shared_ptr<Observer> o2(new Observer()); s.attach(o2);
o0.reset(); // destruct Observer -> No More Exist!
s.detach(o1); // remove from notification list -> Just Unseen from Subject class
s.notify();
return 0;
}

lessons we learned


FFT/FWT

Fast Fourier Transform (FFT)

require 'mathn'

N = 1<<3
f = Array.new(N) {|x| (x<N/2)? -1: 1}

Hadamard = Matrix[[1,1],[1,-1]]/Math.sqrt(2)
def Phase(theta)
Matrix[[1,0],[0,Complex.polar(1,theta)]]
end
def R(k)
lambda {|x,y| (Hadamard*Phase(2*Math::PI*k/N)*Vector[x,y]).to_a}
end

def FFT(arr,s=0,n=arr.size)
abort 'error! not powers of 2' if n&(n-1)!=0
if n==2 then
arr[s],arr[s+1] = R(0).call(arr[s],arr[s+1])
else
arr[s,n] = arr.slice(s,n).partition.with_index {|_,i| i%2==0}.flatten
FFT(arr,s,n/2)
FFT(arr,s+n/2,n/2)
for i in 0...n/2 do
arr[s+i],arr[s+i+n/2] = R(i*N/n).call(arr[s+i],arr[s+i+n/2])
end
end
return
end

puts f
FFT(f)
puts f

Fast Wavelet Transform (FWT)

require 'mathn'

N = 4
f = [9,1,2,0]

Haar = Matrix[[1,1],[1,-1]]/2
R = Proc.new {|x,y| (Haar*Vector[x,y]).to_a}

def FWT(arr,s=0,n=arr.size)
abort 'error! not powers of 2' if n%2!=0
if n==2 then
arr[s],arr[s+1] = R.call(arr[s],arr[s+1])
else
arr[s,n] = arr.slice(s,n).each_slice(2).map(&R).transpose.flatten
FWT(arr,s,n/2)
end
return
end

puts f
FWT(f)
puts f


Matrix Determinant

require 'matrix'
class Matrix
def det_b
abort '!error: not square!' unless square?
n=row_size
a=to_a
(0...n).inject(1) do |det,k|
i=a.transpose[k].slice(k...n).index{|e| e!=0}
return 0 if i.nil? #not full ranked
if i!=0 then
a[k],a[k+i]=a[k+i],a[k]
det*=-1
end
for i in (k+1)...n
for j in (k...n).to_a.reverse
a[i][j]=a[k][k]*a[i][j]-a[i][k]*a[k][j]
a[i][j]/=a[k-1][k-1] if k!=0
end
end
det/=a[k-1][k-1] if k!=0
det*a[k][k]
end
end
end

n=44
require 'mathn'
10.times do
a=Matrix.unit(n).map{|_| rand(10)-5}
#a=Matrix.unit(n).map{|_| 10.0*rand-5.0}
puts '@@@@@@@@@@@@@@@@@@@@@@@@@'
puts 'det: ' << a.det.to_s
puts 'det_e: ' << a.det_e.to_s
puts 'det_b: ' << a.det_b.to_s
end


Inverse Matrix

require 'matrix'
require 'mathn'
A=Matrix[
[1,2,3],
[3,1,2],
[2,3,1]
]
n=A.row_size
abort '!error: not full ranked!' if A.det.zero?
#see (2.6) and Remark 2.1 of [Pan-Reif:1985]
norm=(A*A.transpose).map(&:abs).row_vectors.map{|r| r.to_a.inject(:+)}.max
x=Matrix.scalar(n,1.0/norm)*A.transpose
p x=0xFF.times.inject(x){|x,_| (Matrix.scalar(n,2.0)-x*A)*x}
p Matrix.unit(n)-A*x