米Yahoo!のサービスに、BOSSというのがある。BOSSはウェブ検索サービスを構築するためのオープンプラットフォームで、APIを通じてサードパーティーが自由にYahoo! Searchのインデックスやランキングアルゴリズムなどの技術を使うことができる。
TechCrunchによると、今年の５月現在でBOSSは1日あたり3000万の検索クエリを扱っており、さらに利用量が伸びている。この規模は、Bingと互角かやや抜いているレベルの規模だ。
BOSSを採用している企業には、セマンティック検索エンジンを手掛けるhaikaやソーシャルサービスのリアルタイム検索エンジンOneRiot、より理解しやすい検索結果の提供を目指すCluuzといったものがある。TechCrunchのpower searchにも使われているようだ。
BOSSを使えば大きな投資資金がなくても検索エンジンを作ることができる。その計算資源を使って、何かビジネスができないだろうか？
Haskellでは、関数の戻り値としてタプルを返すことができます。(defun sendray (pt xr yr zr)first-hitがレイと交差した面とその交点の組を返す関数です。それをsendrayから呼び出して、multiple-value-bindで受け取っています。
(multiple-value-bind (s int) (first-hit pt xr yr zr)
(if s
(* (lambert s int xr yr zr) (surface-color s))
0)))
(defun first-hit (pt xr yr zr)
(let ...略...
(values surface hit)))
sendray :: [Surface] -> Light -> Point ->first_hitでMaybe (Surface, Point)型を返し、それをsendrayでパターンマッチで受け取っています。Haskellでは、タプルを返すのに関数をはさむ必要はなく、そのまま書くことができます。
Double -> Double -> Double -> Color
sendray world light pt xr yr zr =
case first_hit world pt xr yr zr of
Just (s@(Sphere (Color col) _ _ _), int)
-> Color ((lambert s int light) * col)
Nothing
-> Color 0.0
first_hit :: [Surface] -> Point ->
Double -> Double -> Double ->
Maybe (Surface, Point)
first_hit world pt xr yr zr = nearest (map hit world)
where nearest :: [Maybe (Double, Surface, Point)]
-> Maybe (Surface, Point)
nearest hits = case foldr cmp Nothing hits of
Just (d, s, h) -> Just (s, h)
otherwise -> Nothing
...略...
今回は、関数の戻り値に正規の値とそうでない値がある場合に、Haskellではどのように表現するのかを書きます。(defun sphere-intersect (s pt xr yr zr)レイと物体の交点を求めるための方程式の解nがある場合はmake-pointで作った交点を返し、解nがない場合には暗黙的にNilを返しています。
(let* ...略...
(if n
(make-point :x (+ (x pt) (* n xr))
:y (+ (y pt) (* n yr))
:z (+ (z pt) (* n zr))))))
intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Pointminrootで求めた方程式の解nがある場合にはJust Pointを返し、解nがない場合にはNothingを返しています。
intersect (Sphere _ r (Point cx cy cz) _) (Point px py pz) xr yr zr =
let ...略...
in case minroot a b c of
Just n -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
Nothing -> Nothing
(defun intersect (s pr xr yr zr)一方、Haskellではパターンマッチを使うことで、型による分岐を自然に表現することができます。
(funcall (typecase s (sphere #'sphere-intersect))
s pt xr yr zr))
(defun sphere-intersect (s pr xr yr zr)
...)
intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point今回はSphereのみですが、その他の型を追加した場合にはマッチさせるパターンを追記するだけです。
intersect (Sphere _ r (Point cx cy cz)) (Point px py pz) xr yr zr =
let a = (sq xr) + (sq yr) + (sq zr)
b = 2 * ((px-cx)*xr + (py-cy)*yr + (pz-cz)*zr)
c = (sq (px-cx)) + (sq (py-cy)) + (sq (pz-cz)) - (sq r)
n = minroot a b c
in case n of
Just n -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
Nothing -> Nothing
(defun tracer (pathname &optional (res 1))一方、Haskellでは手続き的なループは使わないので、このような書き方はしません。
(with-open-file (p pathname :direction :output)
(format p "P2 ~A ~A 255" (* res 100) (* res 100))
(let ((inc (/ res)))
(do ((y -50 (+ y inc)))
((< (- 50 y) inc))
(do ((x -50 (+ x inc)))
((< (- 50 x) inc))
(print (color-at x y) p)))
tracer :: FilePath -> [Surface] -> Int -> IO()putPGMがpgm画像を出力するアクションで、それにピクセルの色のリストcolsを渡すのです。
tracer pathname world res =
let
cols = [color_at world (pos i) (pos j) | j <- [0..(res*100)-1],
i <- [0..(res*100)-1]]
in
do putPGM pathname (res*100) (res*100) cols
where
putPGM :: FilePath -> Int -> Int -> [Int] -> IO()
putPGM pathname w h cols =
let header = "P2 " ++ show w ++ " " ++ show h ++ " 255\n"
body = foldr (++) "" [(show c) ++ "\n" | c <- cols]
in do writeFile pathname (header ++ body)
pos :: Int -> Double
pos i = (fromIntegral i) / (fromIntegral res) - 50.0
バリ島北部に位置するアグン山。ヒンドゥー教で世界の中心と考えられていることから「世界のへそ」と呼ばれています。とあった。
バリ島では、インドネシアのほかの島々とはちがって、島民の93.2％が、ヒンドゥー教を信仰している。これに、アグン山信仰に代表されるような、むかしから伝わる土着の自然信仰がくわわって、独特のバリ・ヒンドゥー教世界がつくられた。とある。
バリ島では、クディリ朝の支配下に入った11世紀初めごろからヒンドゥー・ジャワ文化の影響が及び始め、その後しばらくジャワの支配を離れるが、マジャパヒト王国がバリを征服した1343年以後、16世紀初めにジャワのイスラム化によって同王国が滅亡するまでにヒンドゥー化が広く浸透した。とある。
Appleは、自分で工場を持ってせっせとiPodを作っているわけではなくて生産を外部に委託している。そういった生産設備を持たない会社が製品の生産を委託する先の会社って何というのだろうと思っていた。
シリコンチップの場合はファブレス(fubless)とファウンダリ(foundary)という単語で調べられるが、電子機器の場合はどうも呼び方が違うようで検索しても見つからない。
では、そういう会社は何と呼ばれているかというと、EMS(Electronics Manufacturing Service)というらしい。
大手では、
といった会社がある。やっぱり台湾が強いのか、２社が台湾にある。あとはアメリカ、カナダ、シンガポール。日本にもあるみたい。
発注元と受託先の役割分担としてはおおまかに２種類あって、
があるらしい。ODMでは設計もやる。
EMSについて疑問に思ったことは、
void output_particles( Particles* ps )
{
static int num = 0;
std::string file_name = \
(boost::format( "result%08d.pov" ) % num++ ).str();
std::cout << "processing " << file_name << " ..." << std::endl;
std::ofstream f( file_name.c_str() );
if ( ! f )
{
std::cerr << "cannot open " << file_name << std::endl;
exit(1);
}
f << "#include ¥"colors.inc¥"¥n"
<< "camera {¥n"
<< " location <10,>¥n"
<< " look_at <10,>¥n"
<< "}¥n"
<< "light_source { <0,> color White }¥n";
FOR_EACH_PARTICLES( ps, p )
{
f << "sphere {¥n"
<< " <" << r(0) << ", " << r(1) << ", " << r(2) << ">, 0.5¥n"
<< " texture {¥n"
<< " pigment { color Yellow }¥n"
<< " }¥n"
<< "}¥n";
}
f << std::endl;
}タイムステップごとに.povファイルを出力しています。
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。void calculate_position( Particles* ps )
{
vec accel(3), g(3), norm(3);
double speed, diff, adj;
g(0) = 0.0; g(1) = -9.8; g(2) = 0.0;
FOR_EACH_PARTICLES( ps, p )
{
accel = p->f * SPH_PMASS;
speed = norm2( accel );
if ( speed > SPH_LIMIT*SPH_LIMIT ) {
accel *= SPH_LIMIT / sqrt(speed);
}
// Z-axis walls
diff = 2.0 * RADIUS - ( p->r(2) - MIN(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = 1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(2) - p->r(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = -1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
// X-axis walls
diff = 2.0 * RADIUS - ( p->r(0) - MIN(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(0) - p->r(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = -1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
// Y-axis walls
diff = 2.0 * RADIUS - ( p->r(1) - MIN(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(1) - p->r(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = -1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
accel += g;
p->v += accel * DT;
p->r += p->v * DT / SPH_SIMSCALE;
// stick on x-y plane
p->r(2) = 0.0;
}
}それぞれの粒子について、粒子にかかる力から加速度を計算したあと、各壁面について境界条件を適用し、重力加速度を加えます。そうやって求まった加速度から、粒子の速度を求め、粒子の位置を更新します。
void calculate_force( Particles* ps )
{
double pterm, vterm, r, c;
vec dr(3), force(3), fcurr(3);
FOR_EACH_PARTICLES( ps, p )
{
force(0) = 0.0; force(1) = 0.0; force(2) = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = ( p->r - pj->r ) * SPH_SIMSCALE;
r = norm( dr );
if ( H > r )
{
c = H - r;
pterm = -0.5 * c * SpikyKern * ( p->p + pj->p ) / r;
vterm = LapKern * SPH_VISC;
fcurr = pterm * dr + vterm * ( pj->v - p->v );
fcurr *= c * p->rho * pj->rho;
force += fcurr;
}
}
p->f = force;
}
}それぞれの粒子について、ナビエ=ストークス方程式の圧力項、粘性項にあたる力を計算しています。
void calculate_density_and_pressure( Particles* ps )
{
double r2, c, sum, H2;
vec dr;
H2 = H*H;
FOR_EACH_PARTICLES( ps, p )
{
sum = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = (p->r - pj->r) * SPH_SIMSCALE;
r2 = norm2( dr );
if ( H2 > r2 )
{
c = H2 - r2;
sum += c * c * c;
}
}
p->rho = sum * SPH_PMASS * Poly6Kern;
p->p = ( p->rho - SPH_RESTDENSITY ) * SPH_INTSTIFF;
p->rho = 1.0 / p->rho; // take inverse for later calculation
}
}
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。typedef boost::numeric::ublas::vector<double> vec;
typedef struct
{
vec r;
vec v;
double rho;
double p;
vec f;
} Particle;
boostライブラリのublasを使い、３次元ベクトルを表現しています。typedef std::vector<particle> Particles;STLのvectorコンテナでParticle構造体のベクタを表現しています。
