Notes to self

SPSSINC HETCOR

2011-08-20T13:23:00.000+02:00

Computation of heterogeneous correlations are feasible with the R-package polychor, which in version 18 of SPSS has an official interface (Analyze -> Correlate -> heterogeneous correlations). I have experienced numerical instabilities following calls to the SPSSINC HETCOR interface - but a direct call to R and polycor in regular SPSS syntax avoids these problems.

More on convolution

2011-02-20T15:22:00.000+01:00

Often you need to convolve a particular vector with a lot of other vectors of equal length, then the program below can be used. It illustrates the use of Rcpp/inline and FFTW used from within R.

### The program below can be used, when a vector is convolved several times
### with vectors of equal length

### In the first step the FFT of x is stored and FFTW plans are created,
### in subsequent steps the FFT of x is reused and FFTW plans recalled.

library(Rcpp)
require(inline)

### Definition of plugin, makes linkage to FFTW possible
plug <- Rcpp:::Rcpp.plugin.maker(include.before = "#include ",
libs = paste("-L/usr/local/lib/R/site-library/Rcpp/lib -lRcpp",
"-Wl,-rpath,/usr/local/lib/R/site-library/Rcpp/lib",
"-lfftw3 -lm -L/usr/lib"))
registerPlugin("FFTWconv", plug )

### Convolution function used after initiation
convFFTW <- cxxfunction(
signature(x_fftIn = "numeric", yIn = "numeric"),
body = '
Rcpp::NumericMatrix x_fft(x_fftIn);
Rcpp::NumericVector y(yIn);
int ny=y.size();
int n=x_fft.ncol();
Rcpp:NumericVector retPar(n);

double in_1[n];
double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*y_fft;
in_2=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
out_1=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

fftw_plan forward;
fftw_plan backward;

inputfile = fopen("stored.wisdom", "r");
fftw_import_wisdom_from_file(inputfile);
fclose(inputfile);

forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_MEASURE);
backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);

for(int i=0;ifor(int i=ny;ifftw_execute(forward);

for(int i=0;i in_2[i][0]=x_fft(0,i)*out_1[i][0]+x_fft(1,i)*out_1[i][1];
in_2[i][1]=x_fft(1,i)*out_1[i][0]-x_fft(0,i)*out_1[i][1];
}
fftw_execute(backward);

for(int i=0;i
fftw_destroy_plan(forward);
fftw_destroy_plan(backward);

return retPar;
', plugin="FFTWconv")

### Convolution function used for initiation / creation of FFTW plans
convFFTWstart <- cxxfunction(
signature(xIn = "numeric", yIn = "numeric"),
body = '
Rcpp::NumericVector x(xIn);
Rcpp::NumericVector y(yIn);
int nx=x.size();
int ny=y.size();
int n=nx+ny-1;
Rcpp:NumericMatrix retPar(3,n);

double in_1[n];
double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*x_fft,*y_fft;
in_2=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
out_1=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

fftw_plan forward;
fftw_plan backward;

forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_MEASURE);
backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);
inputfile = fopen("stored.wisdom", "w");
fftw_export_wisdom_to_file(inputfile);
fclose(inputfile);

for(int i=0;i<(n-nx);i++) in_1[i]=0.0;
for(int i=(n-nx);ifftw_execute(forward);
x_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
for(int i=0;i for(int j=0;j<2;j++) x_fft[i][j]=out_1[i][j];

for(int i=0;ifor(int i=ny;ifftw_execute(forward);

for(int i=0;i in_2[i][0]=x_fft[i][0]*out_1[i][0]+x_fft[i][1]*out_1[i][1];
in_2[i][1]=x_fft[i][1]*out_1[i][0]-x_fft[i][0]*out_1[i][1];
}
fftw_execute(backward);

for(int i=0;iretPar(0,i)=out_2[i]/((double)n);
retPar(1,i)=x_fft[i][0];
retPar(2,i)=x_fft[i][1];
}

fftw_destroy_plan(forward);
fftw_destroy_plan(backward);

return retPar;
', plugin="FFTWconv")

convolve_fftw=function(x,y){
return(convFFTW(x,rev(y)))
}

###R function which handles Rcpp/inline functions above
conv.iteration <<- 0

convolve_fftw=function(x,y){
if(conv.iteration==0){
conv.start <<- convFFTWstart(x,rev(y))
conv.iteration <<- 1
return(conv.start[1,])
}
else return(convFFTW(conv.start[2:3,],rev(y)))
}

An example using Rcpp, inline and FFTW: Convolution

2011-02-20T09:07:00.000+01:00

Following a discussion on the Rcpp mailing list I constructed the example below, which illustrates how to use Rcpp.plugin.maker to assure correct linkage of a library - in this case the FFTW library (Fastest Fourier Transform in the West).
The resulting R function is named convolve_fftw

The example also illustrates the general use of Rcpp and inline

library(Rcpp)
require(inline)

plug <- Rcpp:::Rcpp.plugin.maker(include.before = "#include ",
libs = paste("-L/usr/local/lib/R/site-library/Rcpp/lib -lRcpp",
"-Wl,-rpath,/usr/local/lib/R/site-library/Rcpp/lib",
"-lfftw3 -lm -L/usr/lib"))
registerPlugin("FFTWconv", plug )

convFFTW <- cxxfunction(
signature(xIn = "numeric", yIn = "numeric"),
body = '
Rcpp::NumericVector x(xIn);
Rcpp::NumericVector y(yIn);
int nx=x.size();
int ny=y.size();
int n=nx+ny-1;
Rcpp:NumericVector retPar(n,0.0);

double in_1[n];
double out_2[n];

fftw_complex *in_2,*out_1,*irf_fft,*y_fft;
in_2=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
out_1=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

fftw_plan forward;
fftw_plan backward;

forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_ESTIMATE);
backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_ESTIMATE);

for(int i=0;i<(n-nx);i++) in_1[i]=0.0;
for(int i=(n-nx);ifftw_execute(forward);
irf_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
for(int i=0;i for(int j=0;j<2;j++) irf_fft[i][j]=out_1[i][j];

for(int i=0;ifor(int i=ny;ifftw_execute(forward);
y_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
for(int i=0;i for(int j=0;j<2;j++) y_fft[i][j]=out_1[i][j];

for(int i=0;i in_2[i][0]=irf_fft[i][0]*y_fft[i][0]+irf_fft[i][1]*y_fft[i][1];
in_2[i][1]=irf_fft[i][1]*y_fft[i][0]-irf_fft[i][0]*y_fft[i][1];
}

fftw_execute(backward);

for(int i=0;i
fftw_destroy_plan(forward);
fftw_destroy_plan(backward);

return retPar;
', plugin="FFTWconv")

convolve_fftw=function(x,y){
return(convFFTW(x,rev(y)))
}

Til en kvinde

2011-02-19T17:33:00.000+01:00

Dit lange velduftende krøllede mørke hår
På en krop smukkere end de flestes
Vil på nogle dage trækkes ublidt og besidende op til vores hule på bjergsiden
i nåleskoven langt fra dalen.
Langt oftest kysser jeg dig blidt på din hals til velkomst og afsked
stryger din nakke
- og løfter din krop op mod min.
Hånligt tager du fat i min nakke og fører min tunge
fra dine velpudsede sko op langs dine bløde inderlår til dine bryster
- sådan ligger vi!
vekslende en hel nat
i favntag med mørket og lyset i os.

ModernCV - for black and white printers

2011-02-01T12:01:00.001+01:00

The LaTeX class ModernCV does not provide a black color theme option - and printouts can therefore contain rastered colors if printed on black and white printers.

I am using the classic layout and made the following changes in .../tex/latex/modercv/moderncvthemeclassic.sty to get a black color theme.

First of all I copied the style file and edited it as normal user, then copied the file to its original location as root.

In the code section Package options -> color options I added the lines

% black

\newif\if@colorblack\@colorblackfalse

\DeclareOption{black}{\@colorblacktrue}

immediately after the line

\DeclareOption{grey}{\@colorgreytrue}

Under Style definition -> color I added

\if@colorblack

\definecolor{addresscolor}{rgb}{0,0,0}

\definecolor{sectionrectanglecolor}{rgb}{0,0,0}

\definecolor{sectiontitlecolor}{rgb}{0,0,0}

\definecolor{subsectioncolor}{rgb}{0,0,0}

\definecolor{titlecolor}{rgb}{0,0,0}
\definecolor{quotecolor}{rgb}{0,0,0}
\fi

immediately after the lines

\if@colorgrey

\definecolor{addresscolor}{rgb}{0.35,0.35,0.35}

\definecolor{sectionrectanglecolor}{rgb}{0.75,0.75,0.75}

\definecolor{sectiontitlecolor}{rgb}{0.35,0.35,0.35}

\definecolor{subsectioncolor}{rgb}{0.35,0.35,0.35}

\definecolor{titlecolor}{rgb}{0,0,0}

\definecolor{quotecolor}{rgb}{0,0,0}\fi

Skønhed

2011-01-24T19:34:00.003+01:00

Dine øjne er som midnatsskønne skovsøer i måneskin.

Langs bredden mærker jeg grannåle på de nedfaldne blade.

Vil du elske i nat eller skal jeg drukne ... ?

Du kigger på mig fra bredden,

afventende,

med dine kønne øjne,

skønnere end din smukke krop.

Du finder ro, tryghed har du ikke fundet,

kun vished har du fået til låns.

Din indre usikkerhed tager sine slag af livet,

så utroligt tilstedeværende
og på afstand

Jeg har læst i den danske Shakespeare udgave, Romeo & Julie - og følte en smule inspiration hjulpet godt på vej af drømme om en nærværende kvinde.
Her er Shakespeare:

Elskov er som en røg af hede sukke,
der flammer lutret i en elskers øje,
så rent og klart; men når så modgang rammer,
er den et hav, som tårestrømmen øger.
Hvad er den mer? Et vanvid med forstand,
en bitter drik, som sødt dog lindre kan!

Time-dependent variables and SAS PHREG

2011-01-06T19:02:00.008+01:00

There are several ways to include time-dependent variables in a call to SAS' Proc PHREG: Programming statements within PHREG or data steps combined with PHREG counting process format, i.e. in the model statement we use the syntax

model (t_in,t_out)*censored(0) = X1 ... Xn;

instead of

model t*censored(0) = X1 ... Xn;

Click the title to download the following SAS macro I have written. If you find an error please report it, otherwise you are free to use it. Another macro that can be used for the purpose of splitting a data set according to time-dependent variables is the Lexis macro developed by Bendix Carstensen %LEXIS

/*Macro for splitting data according to 0-1 time dependent variables*/

%macro split_time(
filein, /*Data set file containing td_date among other variables*/
id, /*Patient id - needed since we sort wrt. id and datein - could be removed */
td_var, /*Created 0-1 time-dependent variable*/
td_date, /*Date for update of timedependent variable*/
datein, /*Entry date*/
dateout, /*Exit date*/
cenvar, /*Censored 0 yes 1 no*/
fileout, /*Output file*/
labeltext /*td_var label - plain text - no quotes!*/
);

proc sort data=&filein out=temp; by id &datein; run;

data temp1;
set temp;
if (&dateout GT &td_date) then do;
&td_var=1;
&datein=max(&td_date,&datein);
end;
else delete;
run;

data temp3;
set temp;
if (&dateout LE &td_date) then &td_var=0;
else delete;
run;

data temp2;
set temp;
if (&dateout GT &td_date GT &datein) then do;
&dateout=&td_date;
&cenvar=0;
&td_var=0;
end;
else delete;
run;

data temp;
set temp1 temp2 temp3;
label &td_var="&labeltext";
run;

proc sort data=temp out=&fileout;
by &id &datein;
run;

%mend split_time;

/*Example of its use - three dates included on the input file are used*/
%inc 'Z:\split_time.sas';

proc sort data=lv.litval;
by id dateDep dateCon datePsy;
run;

%split_time(
lv.litval,
id,
conDep,
dateDep,
dateFirstPurchase,
eventTimeHosp,
censoredHosp,
outDataDep,
Concomitant use of antidepressants
);

%split_time(
outDataDep,
id,
conCon,
dateCon,
dateFirstPurchase,
eventTimeHosp,
censoredHosp,
outDataCon,
Concomitant use of anticonvulsants
);

%split_time(
outDataCon,
id,
conPsy,
datePsy,
dateFirstPurchase,
eventTimeHosp,
censoredHosp,
lv.litvalHosp,
Concomitant use of antipsychotics
);

Programmer til dannelse af grundlæggende databaser

2010-11-25T20:40:00.013+01:00

Mange bruger tilsyneladende SAS' proc sql når de gennemfører registerstudier på Danmarks Statistiks servere (DST). Med almindelige data steps kommer man langt, og med en baggrund som statistiker er data steps også mere håndgribelige (i samspil med proc sort, proc tabulate, proc append osv.) Dog skal man tage højde for ikke at benytte autogenererede variable (som f.eks. first.id fra en 'by id' statement) sammen med where options - og tjekke at forståelsen af retain statements er tilstede (ved at kontrollere resulterende datasæt).

Vedlagt et par programmer jeg selv har skrevet og finder nyttige - rettelsesforslag er velkomne og hvis de benyttes bedes man angive kilde.

BDS.zip (Basic Data Steps)

The Pleasure Of Finding Things Out

2009-10-22T22:50:00.002+02:00

Found on YouTube:

Part 1 of 5

Part 2 of 5

Part 3 of 5

Part 4 of 5

Part 5 of 5

Impressive - a tool for pdf-presentations

2009-09-29T23:14:00.005+02:00

Impressive is a presentation tool for pdf files. Click >here< to go to the homepage for the project on SourceForge.net

Note: At time of writing last update of project files was 18/9/2008.

Read the documentation to learn about all the options for the program - including the many advanced options: You can play audio files in the background or view full screen video on entry of specific slides, put a timer on slides, set duration of talk, rotate pages, define transitions for specific pages etc. Many of these features require the use of Info scripts (RTFM)

There is also a useful script gettransitions that will export Beamer transitions to Impressive by creating an appropriate Info file Click here to go to the gettransitions site.

Most often the computer used for the pdf presentation will be a 32 bit Windows box. I therefore suggest to download the file Impressive-0.10.2-win32.zip and put it on the usb-pen with your pdf presentation file (and accompanying audio/video files etc.)

In order to use Impressive simply drop the icon of your pdf file on impressive.exe

Impressive has transitions and background rendering enabled by default - create a shortcut to impressive.exe with the options -t None -b in order to disable transitions and background rendering (Rightclick the impressive.exe icon, select Create Shortcut. Then right-click the newly created shortcut and add -t None -b to the path)

Hyperlinks in your presentation will only work if you use pdftk for rendering: Download the pdftk zip file >here< and extract it in your Impressive folder.

The most basic uses of Impressive are: Pressing the tab key in order to watch a gallery of your slides and use this gallery to jump between slides and highlighting an area on a slide by left-clicking and selecting - view the demo.pdf file (using Impressive!) to learn other basic uses of Impressive.

Windows Remote Desktop on Linux

2009-08-16T15:21:00.012+02:00

In this short guide I will describe how to obtain access to a Windows Remote Desktop session using rdesktop on Linux - assuming you obtain VPN access to your companys network using SSL-explorer.

You are invited to make suggestions and provide alternatives in the comments section - but I will not provide support for SSL-explorer, rdesktop, Linux, java or anything else.

I first tried to connect to a Windows XP box through SSL-Explorer using rdesktop on a Fedora 11 box, but soon realized I needed to change permissions for Sun Java etc. in order to succeed.
Instead I installed Linux Mint and from then on, everything went smooth.

Linux Mint is a distro built on top of Ubuntu - which in itself is built on Debian. Linux Mint provides a smoother experience and easy install options for several non-open-source programs and drivers.

If you are not an open-source fanatic and do not hold the brown color design of Ubuntu dear, but otherwise like the Ubuntu experience, then Linux Mint is your choice.

Here we go:

I login to the SSL-explorer homepage (remember to set a security exception in your browser) - on the initial logon you need to accept some Java-dialogs - the SSL-agent attach itself to the panel and my office computer appears in the applications folder of the SSL-menu.

I right-click the SSL-agent in the panel and choose "Tunnel Monitor", which makes the Tunnel Monitor window appear.

When I click the icon for my office computer, the Tunnel Monitor informs me a connection on local port xxxxx is being attempted.
I then run the following command from a terminal:
rdesktop -g 1600x1150 -a 16 -k da 127.0.0.1:xxxxx

The "-g 1600x1150" option gives me the desired size of the remote desktop window, the "-a 16" option tells rdesktop, we will be using a 16 bit color resolution for the session, and "-k da" makes rdesktop aware, I use a keyboard with Danish layout.

A note on choice of remote desktop window size:
I use the 'xinerama' option provided by the proprietary Nvidia graphical driver. As such I cannot use rdesktop with the "-g workarea" option (I use a setup with two monitors).

On the monitor I use for the remote desktop session, I have a panel on the top of the screen. I right-click the panel and choose "Properties" (I use Linux Mint with the Gnome window manager) The Properties dialog informs me the height of the panel is 24 pixels. Since my monitor has a 1600x1200 resolution, I choose to use the option "-g 1600x1150" (I also need 24 pixels for the remote desktop session window bar).

What Makes Us Happy?

2009-05-17T19:14:00.003+02:00

A recent David Brooks column entitled They had it made - refers to a brilliant essay in The Atlantic on the Grant study - one of the most comprehensive longitudinal studies in history.

Watch the video below and then go read the essay What Makes Us Happy?.

Freeman Dyson

2009-03-28T19:40:00.004+01:00

In our family we do not watch television on a daily basis. Besides The New Yorker, various science journals and the Danish weekly newspaper Weekendavisen I use the Internet as my primary source of (mis)information.
A lot of bad things has been said about The New York Times especially in recent years, but the Op-Eds by David Brooks are brilliant - and I read his columns every week on nytimes.com with great interest.
A couple of days ago The New York Times Magazine had a magnificent profile on Freeman Dyson >Click Here<.
While reading the article I remembered the time before my university studies: I read Feynman's bestselling stories (the profile even has a curious link to a Feynman anecdote) and thought I would find a community of creative non-conformists at the science departments. Well, to a large extend I was wrong - but to read the thoughts of a guy like Freeman Dyson provides comfort - and even HOPE! (and it made me order The Scientist as a Rebel)

Tunnelier and your campus account

2008-12-29T16:38:00.010+01:00

A friend asked me to put these instructions on my blog.

This guide can be used to access the servers at the your campus from your Windows computer via ssh. You can thereby gain access to restricted resources such as your home folder and online journal archives. This method is an alternative to ssh -X and a vpn connection. Using this method reduces the amount of traffic between your computer and the campus server to a minimum. I can also recommend this method, if your ordinary vpn-connection does not allow you to access resources on your home network - quite annoying if your printer is attached to your local network.

Download and install Tunnelier:
http://dl.bitvise.com/Tunnelier-Inst.exe
Run Tunnelier.
In the 'Login' window type 'Host' (some.server), 'Username' (i.e. your campus username)
Select 'password' as Initial method, type your password and check 'Store encrypted password in profile'.
In the 'Options' window uncheck all boxes and choose 'Automatically reconnect if successful connection breaks'.
In the 'Services' window check the box 'Enabled' under 'SOCKS / HTTP Proxy Forwarding and set 'Listen port' to 8080. Furthermore check the box 'Enabled' under 'FTP-to-SFTP Bridge'.
Then click 'Save Profile' to the left and save as 'campusconnect'.
To check that you got the settings right click the 'Login' button.

Then open the 'Computer' folder in the Windows 'Start' menu.
In the 'Computer' folder right-click somewhere in the white area of the window and select 'Add a Network connection'. Click 'Next' a couple of times, and write 'ftp://localhost' in the 'Internet or network address' field and click 'Next'. Uncheck the box 'Log on anonymously', write your campus username and click 'Next'. In the field 'Type a name for this network location' type 'campus folder'. Click 'Next' and then 'Finish'. (We use ftp over ssh to make a secure connection compared to ordinary ftp)
By clicking the campus folder in the Computer folder you can now gain access to your campus home folder.

To keep things seperate I use Firefox to surf restricted internet pages like journal homepages through the secure campus-connection and Internet Explorer for all other purposes.
The details below apply equilly well to Internet Explorer, Opera, etc.
Firefox can be downloaded from this location: http://www.mozilla.com/

Install and launch Firefox. In the 'Tools' menu select 'Options' then select the 'Advanced tab' and 'Settings' next to 'Configure how Firefox connects to the Internet'. Select 'Manual proxy configuration:' Set 'SOCKS Host' to '127.0.0.1' and SOCKS Host 'Port' to '8080' then click 'OK'.

If you automatically want to connect to your campus account on Windows startup (otherwise start Tunnelier and login each time you want access):
Create a shortcut to Tunnelier in the Windows Startup folder (create a shortcut on the Desktop and cut-and-paste it to the Startup folder) On Vista: Click the "Start" button, Click "All Programs" then right click the "Startup" folder. Use Windows search on XP if you have trouble locating the Startup folder.
Edit the Tunnelier shortcut in the Startup folder.
In the 'Target:' field you have to write:
"C:\Program Files (x86)\Bitvise Tunnelier\Tunnelier.exe" -profile=C:\Users\ YOUR_WINDOWS_USER_NAME\Documents\campusconnect.tlp -loginOnStartup
Your path to 'Tunnelier.exe' or 'campusconnect.tlp' might be different, if in doubt use Windows search.

Vista, Visual Studio 2008, Intel Suite and IMSL

2008-03-26T23:40:00.011+01:00

Since the first post on Visual Studio and Intel Fortran I have made some substantial modifications to my setup.

First of all - for a long period of time - I went back to Microsoft XP and installed the Microsoft Visual Studio 2005 Premier Partner Edition, which shipped with my copy of the academic version of Intel Visual Fortran 10.1 Professional Edition with IMSL - no problems.

Microsoft Visual Studio 2008 was released some months ago on MSDNAA, and some weeks ago Intel released updates with the relevant interfaces. Now I have decided to give Vista another shot. I have purchased the Intel® Software Development Suite Student Edition for Windows (because of the included optimization tools), installed Visual Studio 2008 and then the Intel suite (plus my 'old' copy of IMSL).

The interested reader can find screenshots attached, that describe my choice of Project Properties and Tools Options in Visual Studio 2008. http://hellmund.dk/VS08settings.zip

My values in Project Properties -> Linker -> Input -> Additional Dependencies:

imsl.lib imslsuperlu.lib imslhpc_l.lib imsls_err.lib imslmpistub.lib mkl_c.lib mkl_ia32.lib libguide.lib

When I use IMSL procedures I have the following statement in my .F90 files
INCLUDE 'link_fnl_static.h'

Any comments on how to use Intel IPP and MPI would be nice!

Introduction to Bayesian Scientific Computing

2008-03-09T21:45:00.006+01:00

D. Calvetti & E. Somersalo, Springer (2007), is a nice little book, with a remarkable good preface.

Here's some quotes:

"The nature of mathematics is being exact, and its exactness is underlined by the formalism used by mathematicians to write it. This formalism, characterized by theorems and proofs, and syncopated with occasional lemmas, remarks and corollaries, is so deeply ingrained that mathematicians feel uncomfortable when the pattern is broken, to the point of giving the impression that the attitude of mathematicians towards the way mathematics should be written is almost moralistic. There is a definition often quoted, "A mathematician is a person who proves theorems", and a similar, more alchemistic one, credited to Paul Erdös, but more likely going back to Alfréd Rényi, stating that "A mathematician is a machine that transforms coffee into theorems" (Footnote: That said, academic mathematics departments should invest on high quality coffee beans and decent coffee makers, in hope of better theorems. As Paul Turán, a third Hungarian mathematician, remarked, "weak coffee is fit only to produce lemmas".)

Therefore is seems to be the form, not the content, that characterizes mathematics, similarly to what happens in any formal moralistic code wherein form takes precedence over content."

On the reason why they do not use the formal framework of theorems and lemmas:
"when obsessed with a formal language that is void of ambiguities, one loses the capability of expressing emotions exactly, while by liberating the language and making it vague, one creates space for the most exact of all expressions, poetry. Thus the exactness of expression is beyond the language. We feel the same way about mathematics."

"Mathematics is a wonderful tool to express liberally such concepts as qualitative subjective beliefs, but by trying to formalize too strictly how to express them, we may end up creating beautiful mathematics that has a life of its own, in its own academic environment, but which is completely estranged to what we initially set forth."

It cites Peter Lax on the occasion of his receiving the 2005 Abel Prize: "When a mathematician says he has solved the problem he means he knows the solution exits, that it's unique, but very often not much more."

"By following the principle of exclusion and subjective learned opinions, we effectively narrow down the probability distributions of the model parameters so that the model produced plausible results. This process is cumulative: when new information arrives, old information is not rejected, as is often the case in the infamous "model fitting by parameter tweaking", but included as prior information. This mode of building models is not only Bayesian, but also Popperian in the wide sense: data is used to falsify hypotheses thus leading to the removal of impossible events or to assigning them as unlikely, rather than to verify hypotheses, which is in itself a dubious project. (Footnote: See A. Tarantola: Inverse problems, Popper and Bayes, Nature Physics 2, 492-492, (2006)) As the classic philosophic argument goes, producing one white swan, or , for that matter, three, does not prove the theory that all swans are white."

C programming and R

2007-11-26T01:51:00.000+01:00

Check out the slides on this page
http://www.biostat.wisc.edu/~kbroman/teaching/

Karl Broman offers some great advice considering programming in C for use with R amongst other things.
A good and short beginners guide.

How to use Fortran95 code in R

2007-11-03T15:55:00.000+01:00

Aim: Write a small test subroutine in Fortran using basic functionality and some calls to the International Mathematics and Statistics Library (IMSL) - compile a dll and call the subroutine from R taking some memory allocation problems into account.

Prerequisites: A working installation of Microsoft Visual Studio and Intel Fortran Compiler v.10 w. IMSL and R 2.6
This link may also be useful.

Save the following code as test.f90 in C:\Fortrancode

!---begin code---

subroutine test(M,N,S1,V,NU,S2,D,INVX)

!DEC$ ATTRIBUTES DLLEXPORT :: test
!Routines necessary for BSKS,
!used in calculation of
!modified Bessel Function of second kind
use BSKS_INT
use UMACH_INT

!We use linear_operators for elegant
!basic linear algebra syntax
use linear_operators

implicit none

!input:
!V vector with values used as input in Bessel fct.
!M number of numbers in vector V
!XDUM matrix with dimensions (M,N)
!D dimension of dummy square matrix MX
!NU fractional order used in calculation of
!value of Bessel fct. modulo 1. Must be positive.
integer, intent(in) :: M,N,D
double precision, intent(in) :: NU, V(M)
double precision :: INVX((D*D)), MX(D,D), XDUM(M,N)

!output:
!S1 sum of elements in XDUM
!S2 values of Bessel fct.
!INVX is both an input and an output vector,
!converted to a square matrix of dim. D
!the inverse is calculated and returned as INVX
double precision, intent(out) :: S1,S2(M)

!some dummy variables
integer :: i,j,k
double precision :: XX, BS(N)

!Set all M*N elements in XDUM to 0.5
XDUM=0.5

!Sum all elements in XDUM and return the result to S1
S1=sum(XDUM)

!Calculation of Bessel fct. values
!corresponding to order NU+N and the
!M different arguments found in V.
!The M results are returned in S2.
do i = 1, M
XX=V(i)
call BSKS (NU, XX, N, BS)
S2(i)=BS(N)
end do

!Transform the INVX vector into a D*D array, store result in MX
!Note use of forall - which is a fast alternative to a
!ordinary for loop.
forall (i =1:D, j=1:D)
MX(i,j)=INVX(i+(j-1)*D)
end forall

!Calculate inverse of MX
MX = .i. MX

!and return result in INVX
forall (i =1:D, j=1:D)
INVX(i+(j-1)*D)=MX(i,j)
end forall

end subroutine test

!---end code---

Run Microsoft Visual Studio as administrator
Choose File -> New Project
Select Intel(R) Fortran -> Library
and click on Dynamic-link library
Then give the project a name (test) and click OK
In the right hand side frame in the Visual Studio main window right-click on Source Files
Choose Add -> Existing Item...
Browse and select the test.f90 file and click OK

Project-> Properties
Select Fortran->Optimization
Set Optimization to Maximize Speed
Select Fortran->Data
Set SEQUENCE Types Obey Alignment Rules to Yes (/align:sequence)
Select Fortran->External Procedures
Set Name Case Interpretation to Lower Case (/names:lowercase)
and Append Underscore to External Names to Yes (\assume:underscore)
Go to Fortran-> Command Line
Write /heap-arrays in Additional Options:
Click OK

Setting the /heap-arrays option circumvent problems when calling the subroutine from R.
(Dummy arrays are allocated using the heap and not the stack.)

Since we use routines from the IMSL library (and maybe IntelMKL) we need a few more adjustments (using 32 bit versions):
Tools->Options...->Intel Fortran->General
Add C:\Program Files\VNI\imsl\fnl600\IA32\include\STATIC
to Includes: (click ... add the path to the list and click OK) and

Add C:\Program Files\VNI\imsl\fnl600\IA32\lib
to Libraries: and Executables:
In Project->Properties->Linker->Input
in the "Additional Dependencies" line, add the following IMSL libraries and Intel MKL library according to your requirement (I highlighted the ones I use):
Static link for 32 bit application: imsl.lib imslsuperlu.lib imslhpc_l.lib imsls_err.lib imslmpistub.lib mkl_c.lib libguide.lib Dynamic link for 32 bit application: imslmkl_dll.lib mkl_c_dll.lib libguide40.lib Static link Intel 64 bit application: imsl.lib imslsuperlu.lib imslscalar.lib mkl_em64t.lib libguide.lib imsls_err.lib imslmpistub.lib Dynamic link Intel 64 bit application: imslmkl_dll.lib mkl_dll.lib libguide40.lib
In Linker ->Command
write /force:multiple in Additional Options:

In the main window change the option above the file editor frame from Debug to Release

Choose Build -> Build test (and wait)
The file test.dll should now be available in
C:\Users\hellmund\Documents\Visual Studio 2005\Projects\test\test\Release
(if your Windows username is hellmund)

In Intel Fortran Compiler version 10 it is not possible to change the default selection of a multi-threaded runtime library to a single-threaded. We thus have to provide som dlls in order to use our test.dll with dyn.load in R.
Copy the files libifcoremd.dll and libmmd.dll from C:\Program Files\Intel\Compiler\Fortran\10.0.025\IA32\Lib to
C:\Fortrancode

I have used the following R-code to call the subroutine.
When writing R-code, I prefer the editor Tinn-R

#begin code

#Load the dll file

dyn.load("C:/Fortrancode/test.dll");

#Is the dll file loaded?
is.loaded("test");

#Set some parameters
M=100000;
V=sqrt(1:M);
D=3;
INVX=matrix(c(1,2,3,4,1,6,7,8,1),ncol=3,byrow=F);
N=5;
S1=0.0;
NU=0.5;
S2=rep(0,M);

#Call the routine
F=.Fortran("test",as.integer(M),as.integer(N),as.double(S1),as.double(V),as.double(NU),as.double(S2),as.integer(D),as.double(INVX))

#See the output
F[[3]]
F[[6]][1:10]

F[[8]]
#Check the inversion of the matrix INVX
solve(INVX)

#end code