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

error