## Using SDL and OpenGL in R

SDL(Simple DirectMedia Layer) is “A cross-platform multimedia library designed to provide fast access to the graphics framebuffer and audio device”. OpenGL is “The Industry Standard for High Performance Graphics”.

With rdyncall package, which provides "a cross-platform framework for dynamic binding of C libraries using a flexible Foreign Function Interface", we can use SDL and OpenGL in R!

Testing Notes:

1. Install rdyncall package from CRAN: http://cran.r-project.org/web/packages/rdyncall/index.html

Unzip it and copy the DLL files to the system path location (e.g. C:\windows\system32).

3. Test Code (Copied from “rdyncall/demo/SDL.R”)
In R:

# Package: rdyncall
# File: demo/SDL.R
# Description: 3D Rotating Cube Demo using SDL,OpenGL and GLU. (dynport demo)

dynport(SDL)
dynport(GL)
dynport(GLU)

# Globals.

surface <- NULL

# Init.

init <- function()
{
err <- SDL_Init(SDL_INIT_VIDEO)
if (err != 0) error("SDL_Init failed")
surface <<- SDL_SetVideoMode(512,512,32,SDL_DOUBLEBUF+SDL_OPENGL)
}

# GL Display Lists

makeCubeDisplaylist <- function()
{
vertices <- as.double(c(
-1,-1,-1,
1,-1,-1,
-1, 1,-1,
1, 1,-1,
-1,-1, 1,
1,-1, 1,
-1, 1, 1,
1, 1, 1
))

colors <- as.raw( col2rgb( rainbow(8) ) )

triangleIndices <- as.integer(c(
0, 2, 1,
2, 3, 1,
1, 3, 7,
1, 7, 5,
4, 5, 7,
4, 7, 6,
6, 2, 0,
6, 0, 4,
2, 7, 3,
2, 6, 7,
4, 0, 5,
0, 1, 5
))

glEnableClientState(GL_VERTEX_ARRAY)
glVertexPointer(3, GL_DOUBLE, 0, vertices )

glEnableClientState(GL_COLOR_ARRAY)
glColorPointer(3, GL_UNSIGNED_BYTE, 0, colors )

displaylistId <- glGenLists(1)
glNewList( displaylistId, GL_COMPILE )
glPushAttrib(GL_ENABLE_BIT)
glEnable(GL_DEPTH_TEST)
glDrawElements(GL_TRIANGLES, 36, GL_UNSIGNED_INT, triangleIndices)
glPopAttrib()
glEndList()

glDisableClientState(GL_VERTEX_ARRAY)
glDisableClientState(GL_COLOR_ARRAY)

return(displaylistId)
}

# Mainloop.

mainloop <- function()
{
displaylistId <- makeCubeDisplaylist()
evt <- new.struct(SDL_Event)
tbase <- SDL_GetTicks()
quit <- FALSE
while(!quit)
{
tnow <- SDL_GetTicks()
tdemo <- ( tnow - tbase ) / 1000

glClear(GL_COLOR_BUFFER_BIT+GL_DEPTH_BUFFER_BIT)

glMatrixMode(GL_PROJECTION)
aspect <- 512/512
gluPerspective(60, aspect, 3, 1000)

glMatrixMode(GL_MODELVIEW)
gluLookAt(0,0,5,0,0,0,0,1,0)
glRotated(sin(tdemo)*60.0, 0, 1, 0);
glRotated(cos(tdemo)*90.0, 1, 0, 0);

glCallList(displaylistId)

glCallList(displaylistId)

SDL_GL_SwapBuffers()

SDL_WM_SetCaption(paste("time:", tdemo),NULL)
while( SDL_PollEvent(evt) != 0 )
{
if ( evt$type == SDL_QUIT ) quit <- TRUE else if (evt$type == SDL_MOUSEBUTTONDOWN )
{
button <- evt$button cat("button ",button$button," at ",button$x,",",button$y,"\n")
}
}
glerr <- glGetError()
if (glerr != 0)
{
cat("GL Error:", gluErrorString(glerr) )
quit <- 1
}
SDL_Delay(30)
}
glDeleteLists(displaylistId, 1)
}

cleanup <- function()
{
SDL_Quit()
}

run <- function()
{
init()
mainloop()
cleanup()
}

run()


You will see a window showing a rotating cube:

OpenGL: http://www.opengl.org/

Dyncall: http://dyncall.org/index.shtml

Rdyncall Package: http://cran.r-project.org/web/packages/rdyncall/index.html

Posted in programming, r | 1 Comment

## Spaces and Convergences

Metric spaces: a set where a notion of distance (called a metric) between elements of the set is defined.

Vector spaces or Linear space: a mathematical structure formed by a collection of vectors.

Normed vector spaces: a Vector space with a norm defined. The norm is an abstraction of our usual concept of length.

Any normed vector space is a metric space by defining $d(x,y)=||y-x||$.

Metric space may not be a vector space.

Inner product spaces: Vector spaces endowed with inner product.

An inner product naturally induces an associated norm, thus an inner product space is also a normed vector space.

A normed vector spaces is also an inner product spaces if and only if its norm satisfies the parallelogram law.

Banach spaces: complete normed vector spaces.

Hilbert spaces: complete inner product spaces.

Weak convergence(in normed vector space)

Suppose X is a normed vector space, $X^{*}$ is the continuous dual of $X$ , and $x_{1},x_{2},\cdots$ is a sequence in $X$, $x_{o}\in X$. Then we say that {$x_{n}$} converges weakly to $x_{0}$ if $\lim_{n\to\infty}f(x_{n})=f(x_{0})$ for every $f\in X^{*}$. The notation for this is $x_{n}\stackrel{w}{\to}x_{0}$.

Strong convergence(in normed vector space)

$||x_{n}-x_{0}||\to0$ ($n\to\infty$)

The notation for this is $x_{n}\to x_{0}$ or $x_{n}\stackrel{s}{\to}x_{0}$. We also call it convergence in norm.

Weak* convergence

Suppose X is a normed vector space, $\{f_{n}\}\subset X^{*}$, $f_{0}\subset X^{*}$($n=1,2,\cdots$).

It is the same as strong convergence for a series of operators.

For any $x\in X$, $f_{n}(x)\to f_{0}(x)$. The notation for this is $f_{n}\stackrel{w*}{\to}f_{0}$.

Strong convergence

If $\{f_{n}\}$ converse to $f_{0}$ by the norm on $X^{*}$, or $||f_{n}-f_{0}||\to0$($n\to\infty$)

It is the same as uniformly convergence for a series of operators.

References:

Posted in convergence, space | 1 Comment

## Understand the Maximum Likelihood Estimation

Let $X_{1},X_{2}\cdots$, be iid random vector with density or probability function $f(x,\theta)$, where $\theta$ is the unknown parameter and suppose the true value of $\theta$ is $\theta_{0}$.

The likelihood function: $L(x,\theta)=\prod_{i=1}^{n}f(X_{i},\theta)$
The logarithm likelihood function: $l(x,\theta)=\Sigma_{i=1}^{n}logf(X_{i},\theta)$

Note that these two functions are actually joint probability functions of the iid data $X_{1},X_{2},\cdots$; we call them “likelihood” functions instead of “probability” functions because we now consider them as functions of the unknown parameter $\theta$.

The maximum likelihood estimation $\hat{\theta}_{MLE}$ is found by maximize the likelihood function, and we will show the idea behind this procedure.

Jessen inequality states that for a concave function $g$, $Eg(X)\leq g(EX)$ for any random variable $X$. $g(x)=logx$ is concave, so under $\theta_{0}$, $E_{\theta_{0}}logf(X,\theta)-E_{\theta_{0}}logf(X,\theta_{0})=E_{\theta_{0}}log\frac{f(X,\theta)}{f(X,\theta_{0})}\leq log[E_{\theta_{0}}\frac{f(X,\theta)}{f(X,\theta_{0})}]=log[\int\frac{f(X,\theta)}{f(X,\theta_{0})}f(x,\theta_{0})dx]=log1=0.$

That is to say, $max_{\theta}E_{\theta_{0}}logf(X,\theta)=E_{\theta_{0}}logf(X,\theta_{0})$, since $E_{\theta_{0}}logf(X,\theta)\leq E_{\theta_{0}}logf(X,\theta_{0})$.

Now a natural criteria for finding $\theta_{0}$ is to find the parameter which maximizes $h(\theta)\equiv E_{\theta_{0}}logf(X,\theta)$.

However, the function we want to maximize $h(\theta)$ is unknown because $h(\theta)$ has something to do with the unknown $\theta_{0}$.

But note that $\theta_{0}$ only appears in the mean operator $E_{\theta_{0}}$, we can overcome this problem by using sample mean operator $\frac{1}{n}\Sigma$ instead.

Finally, we reach that we can find a estimator by maximize $l(x,\theta)=\Sigma_{i=1}^{n}logf(X_{i},\theta)$.

There will be a loss when using sample mean to replace the real one, and the quality of this estimator depends on the law of large numbers.

(Update: 2012/Feb/17) We want to maximize $logP(X;\theta)$, or equivalently, letting

$(logP(X;\theta))'=\frac{P'(X;\theta)}{P(X;\theta)}\to0$.

The above formula can be interpreted as we want

1. $P(X;\theta)\to1$, or as large as possible,

2. $P'(X;\theta)\to0$, or as small as possible.

These two points implies that we hope our estimation $\hat{\theta}$ can make the possibility of the observed data $X$ as large as possible, and at the same time, we want some kind of stability of that possibility, i.e., we hope the rate of change in $P(X;\theta)$ – that is $P'(X;\theta)$, as small as possible.

References:
A Course in Large Sample Theory(Lecture notes), Xianyi Wu
A Course in Large Sample Theory, Thomas S. Ferguson

## Convergence in Probability and Convergence with Probability 1

Definations:

For a random vector $X$ and a sequence of random vectors $X_{1},\cdots X_{n},\cdots$ defined on a probability space $(\Omega,\mathcal{F},P)$.

1. Convergence in Probability or Convergence of probability measures, is kind of convergence in measure, in concept of the measure theory.
$X_{n}\stackrel{P}{\to}X$: $P\{||X_{n}-X||>\epsilon\}\to0$ for any fixed $\epsilon>0$.

2. Convergence with Probability 1, is kind of pointwise convergence in real analysis, it is also called almost surely convergence in the measure theory.
$X_{n}\stackrel{a.s.}{\to}X$: if $P\{X_{n}\to X\}=1$ or $P\{X_{n}\nrightarrow X\}=0$.

Convergence in Probability is weaker than Convergence with Probability 1, as stated below:

Theorem: $X_{n}\stackrel{a.s.}{\to}X\Rightarrow X_{n}\stackrel{P}{\to}X$

Proof: $\{X_{n}\to X\}=${$\omega$: for each $\epsilon>0$ there exists an $N(\omega)>0$ such that $||X_{k}(\omega)-X(\omega)||\leq\epsilon$ for all $k>N(\omega)$}

$=${$\cap_{\epsilon>0}${$\omega:$ there exists an $N(\omega)>0$ such that $||X_{k}(\omega)-X(\omega)||\leq\epsilon$ for all $k>N(\omega)$}}

$=${$\cap_{\epsilon>0}\cup_{n=1}^{\infty}${$\omega:$ $||X_{k}(\omega)-X(\omega)||\leq\epsilon$ for all $k>N(\omega)$}}

$=${$\cap_{\epsilon>0}\cup_{n=1}^{\infty}\cap_{k>n}${$\omega:$ $||X_{k}(\omega)-X(\omega)||\leq\epsilon$}}

So $X_{n}\stackrel{a.s.}{\to}X\Longleftrightarrow P\{X_{n}\to X\}=1$

$\Longleftrightarrow P\{\cap_{\epsilon>0}\cup_{n=1}^{\infty}\cap_{k>n}\{\omega:||X_{k}(\omega)-X(\omega)||\leq\epsilon\}\}=1$

$\Longleftrightarrow P\{\cup_{n=1}^{\infty}\cap_{k>n}\{\omega:||X_{k}(\omega)-X(\omega)||\leq\epsilon\}\}=1$ for all $\epsilon>0$

$\Longleftrightarrow\lim_{n\to\infty}P\{\cap_{k>n}\{\omega:||X_{k}(\omega)-X(\omega)||\leq\epsilon\}\}=1$ for all $\epsilon>0$

$\Longleftrightarrow\lim_{n\to\infty}P\{\cup_{k>n}\{\omega:||X_{k}(\omega)-X(\omega)||>\epsilon\}\}=0$ for all $\epsilon>0$

$\Rightarrow$$\lim_{n\to\infty}P\{||X_{n}-X||>\epsilon\}=0$ for all $\epsilon>0$, or $P\{||X_{n}-X||>\epsilon\}\to0$ for any fixed $\epsilon>0$, that is $X_{n}\stackrel{P}{\to}X$.

if $X_{n}\stackrel{a.s.}{\to}X$, let the set of unconvergence points as $U_{a.s.}=\{\omega:X_{n}(\omega)\nrightarrow X(\omega)\}$ then $P\{U_{a.s.}\}=0$.
if $X_{n}\stackrel{P}{\to}X$, let the set of unconvergence points as $U_{P}=\{\omega:X_{n}(\omega)\nrightarrow X(\omega)\}$ then we donot have $P\{U_{P}\}=0$.

Convergence in Probability do not require $P\{U_{P}\}=0$ but release the condition to $\lim_{n\to\infty}P\{U_{P}(n)\}=0$, where $U_{P}(n)$ is a set which will go smaller as $n\to\infty$.

From the proof above we can see that $U_{P}=\{\cup_{\epsilon}\cap_{n=1}^{\infty}\cup_{k}\{\omega:||X_{k}(\omega)-X(\omega)||>\epsilon\}\}$ and $U_{P}(n)=\{\omega:||X_{n}-X||>\epsilon\}$, that is to say $U_{P}$ is bigger than $U_{a.s.}$

References:
A Course in Large Sample Theory(Lecture notes), Xianyi Wu
A Course in Large Sample Theory, Thomas S. Ferguson

Posted in analysis, convergence, probability | 24 Comments

## Stochastic o and O symbols

Stochastic $o$ and $O$ symbols are the basic symbols for Asymptotic Statistics or Large Sample Theory.

(i) $A_{n}=o_{p}(B_{n})$: if $|\frac{A_{n}}{b_{n}}|\stackrel{P}{\to}0$.

sequence of random variables $A_{n}$ is of smaller order in probability than a sequence $B_{n}$.

In particular, $A_{n}=o_{p}(1)$, “small oh-P-one”, if and only if $A_{n}\stackrel{P}{\to}0$; so $A_{n}=o_{p}(B_{n})$ means $A_{n}=Y_{n}B_{n}$ and $Y_{n}\stackrel{P}{\to}0$.

Example: $X=o_{p}(1)$ means $X\stackrel{P}{\to}0$, and $X=o_{p}(n^{-1/2})$ means $n^{1/2}X\stackrel{P}{\to}0$, or $X$ goes to $0$ faster than $\frac{1}{n^{1/2}}$ in probability(such as $X=\frac{1}{n}$).

(ii) $A_{n}=O_{p}(B_{n})$ : if given $\epsilon>0$, there exists a constant $M=M(\epsilon)$ and an integer $n_{0}=n_{0}(\epsilon)$ such that $P(|A_{n}|\leq M|B_{n}|)\geq1-\epsilon$ for all $n>n_{0}$.

sequence $A_{n}$ is be of order less than or equal to that of $B_{n}$ in probability.

In particular, $A_{n}=O{}_{p}(1)$, “big oh-P-one”, if for any$\epsilon>0$, there exists a constant $M$ and an integer $n_{0}$ such that $P(|A_{n}|\leq M)\geq1-\epsilon$ for all $n>n_{0}$, $A_{n}$ is said to be bounded in probability(or tight); so $A_{n}=O_{p}(B_{n})$ means $A_{n}=Y_{n}B_{n}$ and $Y_{n}=O_{p}(1)$.

It’s easy to see from the definition that $O_{p}(1)=O_{p}(C)$ for any constant $0.

(iii) $A_{n}\asymp_{p}B_{n}$: if given $\epsilon>0$, there exist constants $0 and an integer $n_{0}$ such that $P[m<|\frac{A_{n}}{B_{n}}| for all $n>n_{0}$.

sequence $A_{n}$ is said to be of the same order as $B_{n}$ in probability.

Some facts:

$o_{p}(1)+o_{p}(1)=o_{p}(1)$: If $X_{n}\stackrel{P}{\to}0$ and $Y_{n}\stackrel{P}{\to}0$, then $Z_{n}=X_{n}+Y_{n}\stackrel{P}{\to}0$. (example of continuous-mapping theorem)

$o_{p}(1)+O_{p}(1)=O_{p}(1)$

$O_{p}(1)o_{p}(1)=o_{p}(1)$: If the sequence $\{Y_{n},n=1,2,\cdots\}$ is bounded in probability and if $\{C_{n}\}$ is a sequence of random variables tending to 0 in probability, then $C_{n}Y_{n}\stackrel{P}{\to}0$.

$(1+o_{p}(1))^{-1}=O_{p}(1)$

$o_{p}(R_{n})=R_{n}o_{p}(1)$

$O_{p}(R_{n})=R_{n}O_{p}(1)$

$o_{p}(O_{p}(1))=o_{p}(1)$

Lemma: Let $R$ be a function defined on domain in $\mathcal{R}^{k}$ such that $R(0)=0$. Let $X_{n}$ be a sequence of random vectors with values in the domain of $R$ that converges in probability to zero. Then, for every $p>0$,

(i) if $R(h)=o(||h||^{p})$ as $h\to0$, then $R(X_{n})=o_{p}(||X_{n}||^{p})$;

(ii) if $R(h)=O(||h||^{p})$ as $h\to0$, then $R(X_{n})=O_{p}(||X_{n}||^{p})$;

Result: For a random variable $S$, $S=ES+O_{p}(\sqrt{Var(S}))$.

Proof:

We only needs to prove that $(S-ES)/\sqrt{Var(S)}=O_{p}(1)$ or equally, for any$\epsilon>0$, there exists a constant $M$ and an integer $n_{0}$ such that $P(|(S-ES)/\sqrt{Var(S)}|\leq M)\geq1-\epsilon$ for all $n>n_{0}$.

Let $NS=(S-ES)/\sqrt{Var(S)}$,

According to Markov inequality, $P(|NS|\leq M)\geq ENS^{2}/M^{2}=[E(S-ES)^{2}/Var(S)]/M^{2}=1/M^{2}\to0$ as $M\to\infty$.

From the proof above we know that for any normalized random variable $NS=(S-ES)/\sqrt{Var(S)}$, we have $NS=O_{p}(1)$, or $NS$ is bounded in probability- the reason is natural, if any random variable is not bounded, either its mean is too large($E(S_{n})\to\infty$) or it varies too much($Var(S_{n})\to\infty$), and normalization will eliminate those two possibilities. On the other hand, for a specified random variable $S$, if $ES<\infty$ and $Var(S)<\infty$, then $S=ES+O_{p}(1)$, especially, when $ES=0$, $S=O_{p}(1)$ .

Example: from center limit theorem we know that $\sqrt{n}(\bar{X}-EX)\to N(0,DX)$, then we have

//$\sqrt{n}(\bar{X}-EX)=N(0,DX)+o_{p}(1)=\sqrt{DX}O_{p}(1)+o_{p}(1)=O_{p}(1)$//,

//$\bar{X}=EX+O_{p}(1)\times n^{-1/2}=EX+O_{p}(n^{-1/2})$//.

$P(\frac{\sqrt{n}(\bar{X}-EX)}{DX}>M )\to P(Z>M),\ Z\sim N(0,1)$. $P(Z>M)$ can be smaller than $\forall\epsilon$ as long as $M$ is large enough, so $\frac{\sqrt{n}(\bar{X}-EX)}{DX}=O_p(1)$, or $\bar{X}=EX+O_{p}(n^{-1/2})$.

The weak law of large numbers states that $\bar{X}\stackrel{P}{\to}EX$, so we have

$\bar{X}-EX=o_{p}(1)$.

(Update: 2012/Feb/17) Similarly, let $X_{n}$ be a sequence of random vectors, using Markov inequality $P(|X_{n}|>M)\leq\frac{E|X_{n}|^{k}}{M^{k}}$, we have

1. If there is a number $k>0$ such that $E|X_{n}|^{k}$ is bounded, then $X_{n}=O_{p}(1)$;
similarly, if $E|X_{n}|^{k}\leq ca_{n}$, where $c$ is a constant and $a_{n}$ is a sequence of positive numbers,
then $X_{n}=O_{p}(a_{n}^{1/k})$.

2. If there is a number $k>0$ such that $E|X_{n}|^{k}\to0$ (So $M$ can be $\epsilon$), then $X_{n}=o_{p}(1)$;
similarly, if $E|X_{n}|^{k}\leq ca_{n}$, where $c$ is a constant and $a_{n}$ is a sequence of positive numbers,
then $X_{n}=o_{p}(b_{n})$ for any sequence $b_{n}>0$ such that $b_{n}^{-1}a_{n}^{1/k}\to0$.

3. If there are sequences of vectors $\{\mu_{n}\}$ and singularization matrices $\{A_{n}\}$ such that $A_{n}(X_{n}-\mu_{n})$ converges in distribution, then $X_{n}=\mu_{n}+O_{p}(||A_{n}^{-1}||)$.

References:

## Using Lib FANN in R via Rcpp

Test Environment: OS: Windows XP, R: R 2.12.2(G:\Program Files\R\R-2.12.2\), Rcpp: 0.9.4, Rtools: 212 (G:\Rtools), FANN: 2.1.0
Lib FANN(Fast Artificial Neural Network Library) is a free open source neural network library written in C. In this tutorial I will show you how to use the Rcpp package to wrap FANN to R.

2. Unzip the source code to some folder(“G:\libFANN\fann-2.1.0\“)

3. Make sure that you installed Rtools and setted the path for MinGW properly. You can check this post for that. The path for my MinGW(Rtools) is “G:\Rtools\MinGW

4. Install MSYS: Download and install MSYS 1.0.11 from http://www.mingw.org/wiki/MSYS. I installed MSYS into “G:\msys“. Make sure to tell msys where your installed MinGW(Rtools) is. This can be done in the post-install process of MSYS, or manually create a file named “fstab” in folder “G:\msys\1.0\etc“, with the following line:
G:/Rtools/MinGW /mingw #tell msys the path of MinGW

5. Run msys(there should be a shortcut on your desktop: “G:\msys\1.0\msys.bat -norxvt“), do the config for libFANN:
cd G:\libFANN\fann-2.1.0
./configure
Now you should find the Makefile created in “G:\libFANN\fann-2.1.0\src, let’s compile the source code of FANN:
cd src
make
After the compilation, you will find some files such as “fixedfann.o, floatfann.o, doublefann.o” created in the src folder.

6. Let’s first test the compiled lib: Copy the files “xor.data, xor_train.c, xor_test.c” from “examples” folder to “src” folder.
ar rc libFANNLib.a doublefann.o fixedfann.o floatfann.o
ranlib libFANNLib.a
gcc -o xor_train xor_train.c -Iinclude libFANNLib.a
gcc -o xor_test xor_test.c -Iinclude libFANNLib.a
Now will can find “xor_train.exe, xor_test.exe” in the “src” folder, click one after the other to run the test(to view the results, you should run them both from cmd).

7. It’s time to wrap FANN to R, we will write a simple package “RcppFANN” using Rcpp Module: go to R

library(Rcpp);
Rcpp.package.skeleton("RcppFANN");


8. Go to the folder “C:\Documents and Settings\Administrator\My Documents\RcppFANN\src” and edit the two files “rcpp_hello_world.h, rcpp_hello_world.cpp“, we simply wrap the example code “xor_train.c, xor_test.c” into “rcpp_hello_world.cpp” using Rcpp Module:

//rcpp_hello_world.h
#ifndef _RcppFANN_RCPP_HELLO_WORLD_H
#define _RcppFANN_RCPP_HELLO_WORLD_H

#include

/*
* note : RcppExport is an alias to extern "C" defined by Rcpp.
*
* It gives C calling convention to the rcpp_hello_world function so that
* it can be called from .Call in R. Otherwise, the C++ compiler mangles the
* name of the function and .Call can't find it.
*
* It is only useful to use RcppExport when the function is intended to be called
* on Rcpp-devel for a misuse of RcppExport
*/

#endif


//rcpp_hello_world.cpp
#include "rcpp_hello_world.h"

using namespace Rcpp ;

#include

#include "fann.h"

int xor_test()
{
fann_type *calc_out;
unsigned int i;
int ret = 0;

struct fann *ann;
struct fann_train_data *data;

printf("Creating network.\n");

#ifdef FIXEDFANN
ann = fann_create_from_file("xor_fixed.net");
#else
ann = fann_create_from_file("xor_float.net");
#endif

if(!ann)
{
printf("Error creating ann --- ABORTING.\n");
return 0;
}

fann_print_connections(ann);
fann_print_parameters(ann);

printf("Testing network.\n");

#ifdef FIXEDFANN
#else
#endif

for(i = 0; i < fann_length_train_data(data); i++) 	{ 		fann_reset_MSE(ann); 		calc_out = fann_test(ann, data->input[i], data->output[i]);
#ifdef FIXEDFANN
printf("XOR test (%d, %d) -> %d, should be %d, difference=%f\n",
data->input[i][0], data->input[i][1], calc_out[0], data->output[i][0],
(float) fann_abs(calc_out[0] - data->output[i][0]) / fann_get_multiplier(ann));

if((float) fann_abs(calc_out[0] - data->output[i][0]) / fann_get_multiplier(ann) > 0.2)
{
printf("Test failed\n");
ret = -1;
}
#else
printf("XOR test (%f, %f) -> %f, should be %f, difference=%f\n",
data->input[i][0], data->input[i][1], calc_out[0], data->output[i][0],
(float) fann_abs(calc_out[0] - data->output[i][0]));
#endif
}

printf("Cleaning up.\n");
fann_destroy_train(data);
fann_destroy(ann);

return ret;
}

////

int FANN_API test_callback(struct fann *ann, struct fann_train_data *train,
unsigned int max_epochs, unsigned int epochs_between_reports,
float desired_error, unsigned int epochs)
{
printf("Epochs     %8d. MSE: %.5f. Desired-MSE: %.5f\n", epochs, fann_get_MSE(ann), desired_error);
return 0;
}

int xor_train()
{
fann_type *calc_out;
const unsigned int num_input = 2;
const unsigned int num_output = 1;
const unsigned int num_layers = 3;
const unsigned int num_neurons_hidden = 3;
const float desired_error = (const float) 0;
const unsigned int max_epochs = 1000;
const unsigned int epochs_between_reports = 10;
struct fann *ann;
struct fann_train_data *data;

unsigned int i = 0;
unsigned int decimal_point;

printf("Creating network.\n");
ann = fann_create_standard(num_layers, num_input, num_neurons_hidden, num_output);

fann_set_activation_steepness_hidden(ann, 1);
fann_set_activation_steepness_output(ann, 1);

fann_set_activation_function_hidden(ann, FANN_SIGMOID_SYMMETRIC);
fann_set_activation_function_output(ann, FANN_SIGMOID_SYMMETRIC);

fann_set_train_stop_function(ann, FANN_STOPFUNC_BIT);
fann_set_bit_fail_limit(ann, 0.01f);

fann_init_weights(ann, data);

printf("Training network.\n");
fann_train_on_data(ann, data, max_epochs, epochs_between_reports, desired_error);

printf("Testing network. %f\n", fann_test_data(ann, data));

for(i = 0; i < fann_length_train_data(data); i++) 	{ 		calc_out = fann_run(ann, data->input[i]);
printf("XOR test (%f,%f) -> %f, should be %f, difference=%f\n",
data->input[i][0], data->input[i][1], calc_out[0], data->output[i][0],
fann_abs(calc_out[0] - data->output[i][0]));
}

printf("Saving network.\n");

fann_save(ann, "xor_float.net");

decimal_point = fann_save_to_fixed(ann, "xor_fixed.net");
fann_save_train_to_fixed(data, "xor_fixed.data", decimal_point);

printf("Cleaning up.\n");
fann_destroy_train(data);
fann_destroy(ann);

return 0;
}

////

RCPP_MODULE(mod){
function("xor_test", &xor_test);
function("xor_train", &xor_train);
}


9. Edit the “Makevars.win” to:
## Use the R_HOME indirection to support installations of multiple R version
PKG_LIBS = $(shell “${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe” -e “Rcpp:::LdFlags()”) G:/libFANN/fann-2.1.0/src/libFANNLib.a PKG_CPPFLAGS = -I”G:\libFANN\fann-2.1.0\src\include” 10. Now we can compile and install our test package: Windows Start -> run -> cmd -> cd C:\Documents and Settings\Administrator\My Documents C:\Documents and Settings\Administrator\My Documents> “G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” check RcppFANN “G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” INSTALL RcppFANN 11. Finally, let test our package: copy the data files “xor.data“(located in “G:\libFANN\fann-2.1.0\examples“) to folder “C:\Documents and Settings\Administrator\My Documents\“(or your working directory). In R: library("Rcpp"); library("RcppFANN"); mod mod$xor_train();



Result:

[1] 0


If everything goes well, you can find three files “xor_fixed.data, xor_fixed.net, xor_float.net“(the trained network) created in “C:\Documents and Settings\Administrator\My Documents\“.

mod$xor_test();  Result: [1] 0  Well, FANN is successfully running in R, but we can’t see anything useful in R because the results of the original test function are printed using the C “printf” function. To make the lib really useful, you can write some further wrapper functions to retrieve those results and supply data to FANN in R. Source code of the test package “RcppFANN”: https://bzstat.googlecode.com/svn/trunk/blogpost/RcppFANN Posted in programming, r | Leave a comment ## The Rcpp package test notes Test Environment: OS: Windows XP, R: R 2.12.2(G:\Program Files\R\R-2.12.2\), Rcpp: 0.9.4 1. Download and install the Rtools(for windows): http://www.murdoch-sutherland.com/Rtools/ I installed the Rtools into “G:\Rtools“. 2. Set Enviroment Path for windows: Desktop -> right-click My Computer -> click Properties -> click the Advanced tab -> click Environment Variables button -> highlight the “Path” variable in the Systems Variable section -> click the Edit button -> addG:\Rtools\bin;G:\Rtools\MinGW\bin;” there (“C:\MiKTeX\miktex\bin;” may needed too.) and in the Systems Variable section -> click the New button -> add a new variable named “CYGWIN” with value “nodosfilewarning” (Turn off MS-DOS style path warnings with cygwin utilities) 3. Install the R package “Rcpp” and “inline”. 4. Quick test using inline library(inline); library("Rcpp"); inc <- ' double norm(double x, double y){ return sqrt(x*x+y*y); } RCPP_MODULE(mod){ function("norm", &norm); } ' fx <- cxxfunction( signature(), "" , include = inc, plugin = "Rcpp" ); mod <- Module( "mod", getDynLib(fx) ); mod$norm(3,4);


Result:

[1] 5


5. Quick test of writing a simple package

Rcpp.package.skeleton("RcppTestPackage");

Creating directories ...
Creating DESCRIPTION ...
Creating NAMESPACE ...
Saving functions and data ...
Making help files ...
Done.
Further steps are described in './RcppTestPackage/Read-and-delete-me'.

>> added useDynLib directive to NAMESPACE
>> added Makevars file with Rcpp settings
>> added Makevars.win file with Rcpp settings
>> added example src file using Rcpp classes
>> added example R file calling the C++ example
>> added Rd file for rcpp_hello_world


Now you can find the created file in “C:\Documents and Settings\Administrator\My Documents\RcppTestPackage” (or your working directory)

Now we can compile and install the test package:
Windows Start -> run -> cmd -> cd C:\Documents and Settings\Administrator\My Documents
“G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” check RcppTestPackage
“G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” INSTALL RcppTestPackage
(“G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” build RcppTestPackage (This is the file you can upload to CRAN.))

Finally, go to R and test your package:

library(RcppTestPackage);
rcpp_hello_world();


Result:

[[1]]
[1] "foo" "bar"

[[2]]
[1] 0 1


6. Create a simple package using Rcpp Module:

Rcpp.package.skeleton("RcppTestModule");


Find the created file in “C:\Documents and Settings\Administrator\My Documents\RcppTestPacakge” (or your working directory)
Go to the “src” folder and make changes to the files “rcpp_hello_world.h” and “rcpp_hello_world.cpp”,
Let’s write something meaningful:

//rcpp_hello_world.h
#ifndef _RcppTestPacakge_RCPP_HELLO_WORLD_H
#define _RcppTestPacakge_RCPP_HELLO_WORLD_H

#include <Rcpp.h>

/*
* note : RcppExport is an alias to extern "C" defined by Rcpp.
*
* It gives C calling convention to the rcpp_hello_world function so that
* it can be called from .Call in R. Otherwise, the C++ compiler mangles the
* name of the function and .Call can't find it.
*
* It is only useful to use RcppExport when the function is intended to be called
* on Rcpp-devel for a misuse of RcppExport
*/

#endif

//rcpp_hello_world.cpp
#include "rcpp_hello_world.h"

using namespace Rcpp;

double norm(double x, double y){
return sqrt(x*x+y*y);
}

RCPP_MODULE(mod){
function("norm",
&norm,
List::create( _["x"] = 0.0, _["y"] = 0.0),
"Provides a simple vector norm");
}


Similarly, compile and install the test package:
Windows Start -> run -> cmd -> cd C:\Documents and Settings\Administrator\My Documents
“G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” check RcppTestModule
“G:\Program Files\R\R-2.12.2\bin\i386\Rcmd” INSTALL RcppTestModule

Finally, test our package in R:

library("RcppTestModule");
mod <- Module( "mod", PACKAGE = "RcppTestModule");
mod$norm(3,4);  Result: [1] 5  and mod$norm();


Result:

[1] 0


also

show(mod\$norm);


Result:

internal C++ function <0241A0A8>
docstring : Provides a simple vector norm
signature : double norm(double, double)


Source code of package “RcppTestModule”: https://bzstat.googlecode.com/svn/trunk/blogpost/RcppTestModule