Renan Cabrera L. PhD

Conservation of Shannon information in classical mechanics

2012-02-12T00:02:00.000-08:00

There are many places where to find the proof of the conservation of the quantum Shannon information under unitary evolution. However, I just could not find the equivalent proof that everybody talks in classical mechanics, so here it goes.

The Lioville equation can be expressed in terms of a self-adjoint operator as observed by Koopman and von Neumann, such that

Given that the Liouvillian follows the Leibniz rule, which is true for conservative systems

it is now easy to prove that

This fact is in complete contrast with respect to the thermodynamical entropy, which can be fundamentally irreversible.

-This entry was powered by http://www.codecogs.com/latex/eqneditor.php

DUST

2012-02-09T13:31:00.000-08:00

DUST stands for Differential Unitary Space Time coding. This coding method is used in the transmission of information from multiple antennas to multiple antennas.

The transmission of information in a block can be modelled with the following equation

where S is the input matrix, H is the channel matrix, \rho is a scalar coefficient that modulates the channel and W is the noise.

DUST encodes the information in a unitary matrix U_n, which is indirectly encoded through a sequence of two S unitary blocks as

If H is full-rank and the channel does not change significantly in the time the two blocks are transmitted. The unitary matrix U can be recovered from the sequence of two received unitary blocks X without! the need to know the channel H as

where the new noise is defined as

-This entry was powered by http://www.codecogs.com/latex/eqneditor.php

Random Unitary Matrices

2012-02-07T11:58:00.000-08:00

The generation of random unitary matrices has many applications in many fields.
The simplest method is based on the QR decomposition of a random complex matrix with independent elements following a Normal distribution.

In Mathematica one can generate a random complex number with Normal distribution for both the real and imaginary parts with the following function

NormalRandomComplex[] := (#[[1]] + #[[2]]*I) &@
RandomReal[MultinormalDistribution[{0, 0}, {{1, 0}, {0, 1}}]]

where {0, 0} specifies the origin and {{1, 0}, {0, 1}} the covariance matrix (standard deviation 1)

The function that generates the random n x n unitary matrix is

RandomUnitaryMatrix[n_] := Module[{z, q, r},
  z = Table[NormalRandomComplex[], {n}, {n}];
  {q, r} = QRDecomposition[z];
  q.DiagonalMatrix[Sign[Tr[r, List]]]

]

In this code, the QR decomposition of the complex random matrix with normal distribution outputs a unitary matrix "q" and a triangular matrix "r". The sought random unitary matrix is a correction of "q", which must be multiplied by a diagonal matrix filled with 1 and -1 according to the sign of the corresponding diagonal element of "r".

New numerical algorithms

2012-01-19T12:42:00.000-08:00

Two recent developments on numerical analysis have called my attention

A new very efficient sparse Fourier transform
http://web.mit.edu/newsoffice/2012/faster-fourier-transforms-0118.html

A significantly faster algorithm for matrix multiplication
http://www.newscientist.com/article/mg21228422.500-mathematical-matrix-multiplier-sees-first-advance-in-24-years.html

De Broglie waves and relativity

2011-12-08T20:31:00.000-08:00

Another extremely interesting connection between relativity and quantum mechanics is that the quantum De Broglie waves can be understood as an effect of the relativistic desynchronization of clocks given that the clocks rotate with a frequency proportional to the total energy of the particle.

Imagine a train with a long chain of synchronized clocks from one extreme to the other with the synchronization established according to an observer moving along with the train. Another person standing still on the train station will conclude that the clocks are actually not synchronized (what he would see is a different thing). If the train is moving to the right, the clocks on the right will be ahead respect to the clocks on the left. The De Broglie wave length will be the length that stretches 12 hours of desynchronization.

[1] W. Baylis, Canadian Journal of Physics, 2007, 85:(12) 1317-1323, 10.1139/p07-121

The ball is round!

2011-12-06T11:48:00.000-08:00

The first thing we learn in a course of relativity is the relativistic length contraction effect for moving objects relative to an observer. However, what one would actually see is a different thing ! because another important effect must be considered; the time required for the light to arrive from the different points of the observed object to our eyes. The combined effect of the Lorentz contraction along with the time delay was investigated reaching to curious results. For example, a relativistic ball remains with a round appearance!, as discussed by Penrose in the following paper

http://adsabs.harvard.edu/abs/1959PCPS...55..137P

A nice visulaization can be found at

http://www.spacetimetravel.org/ueberblick/ueberblick1.html

Thermodynamical entropy and Shannon information

2011-11-07T10:50:00.000-08:00

The moment I read about information theory I thought that it was evident that the thermodynamical entropy was completely equivalent; the same thing, just measured in different units with the Boltzmann constant in the middle. I was shocked when I read that this issue is still under debate today. There are two books that share the idea that the thermodynamical entropy is fundamentally a measure of information and expose the controversy respect with the opposite view:

A Farewell To Entropy: Statistical Thermodynamics based on Information. By Arieh Ben-Naim

E.T. Jaynes: Papers on Probability, Statistics and Statistical Physics. By R.D. Rosenkrantz (Editor)

PyNewtonCUDA

2011-11-01T14:47:00.000-07:00

Denys Bondar and I, just released to the public: PyNewtonCUDA

https://code.google.com/p/py-newton-cuda/

This program is able to propagate a large number of classical particles using GPU technology. The code is written in Python and PyCUDA, which allows to write high-level code for GPU programming while maintaining very high performance.

You can download the code with the following command

svn checkout https://py-newton-cuda.google.com/svn/trunk PyNewtonCUDA

Classical Relativistic Many-Body Dynamics

2011-06-24T23:31:00.000-07:00

I found a book on the interesting and very difficult topic of many-body classical relativistic theory

Classical Relativistic Many-Body Dynamics

Trump, M.A., Schieve, W.C
Springer, Fundamental Theories of Physics, Vol. 103

This relativistic covariant many-body formulation overcomes the no go theorem by Jordan and Sudarshan [1] by giving up the locality of the interactions and the individual mass-shell conditions. These are not simply cosmetic upgrades but have serious consequences on how we should understand classical physics.

"The single particle has no significance; it is the whole system that counts" --Lanczos

[1] FROM CLASSICAL TO
QUANTUM MECHANICS
Giampiero Esposito, Giuseppe Marmo

Geo-coordinates for citations on scientific papers

2011-06-03T05:33:00.000-07:00

Here I draw a map for the places that either published a paper about the Bures measure or were referenced by somebody who made a paper about the Bures measure

The corresponding plot with the citations as links was created by my brother Ruben Cabrera, shown next

The total number of papers in consideration were 662. The extraction of coordinates from heterogeneous references is an expensive computation, so I do not expect to make more of them just for fun.

Mathematica CDF documents

2011-05-19T17:07:00.000-07:00

This is my first Mathematica cdf document on the Internet.
The most remarkable feature of a cdf document is that it can be posted on a web site and evaluated by anybody with a web browser and the freely available cdf player plugin.

Note.- Unfortunately, it only works on Windows and Mac.

Bibliographic Graphs Time-Line

2011-05-13T12:37:00.000-07:00

Different colors for different publication dates. Dark green for the oldest to dark purple to the latest

Bures measue

Haar measure

Bibliographic Graphs

2011-05-10T17:26:00.000-07:00

Continuing with bibliographic graphs.
Here I combine the graphs from the Bures measure in green and the Haar measure in red, where the publications that are common in both are in orange

Bibliographic Graphs

2011-05-08T04:59:00.000-07:00

I just made a Mathematica program to generate the citation graph for scientific academic publications given a specific keyword (and some upper limit on the number of publications to treat). Some examples for three areas of my own interest are shown below

"Bures measure"

The same job done for "Haar measure"

It seems that the research done for the "Haar measure" is made of several unconnected niches. However, I expect that they will find some connection as I let the program run for more time.

One final example for "Quantum control"

Bibliographic research is now easier because I can put the mouse on the top of each node and see what paper is that. Something nice would be to output the whole bibtex file.

Acknowledgements: This program was made possible thanks to the help of my friend Josh Greig.

Extracting data from a pdf file

2011-05-01T13:19:00.000-07:00

There are a few quite expensive commercial programs for extracting data from pdf files, but a few lines of Mathematica code may be able to do the job for most cases. For example, the following file with many pages

results in the following table

Note.- This method works for text-searchable pdf files. Otherwise, this would be an image processing job.

Mathematica as a C code generator

2011-04-21T14:34:00.000-07:00

Mathematica is able to generate C code ready to be compiled and ultimately capable to produce a stand alone executable file. Nice !. There are restrictions on the extent of what can be be processed in this way but it can be really useful in practice as it is right now.

The generation of the executable can be done within Mathematica or outside Mathematica. The Mathematica 8.0 help for my Linux Fedora 64bits requires a minor adjustment in order to work. For Linux systems:

"Libraries" -> "WolframRTL_Static_Minimal.lib"

must be replaced by

"Libraries" -> "WolframRTL_Static_Minimal"

(Thanks to Abdul Dakkak for helping on this)

How do we compile the generated code outside Mathematica? This issue is detailed next

The Mathematica code that generates the C code for a simple example taken from the help is

code =
  CProgram[
   CInclude["stdio.h"],
   CInclude["compute.h"],
   CFunction["int", 
    "main", {{"int", "argc"}, {"char", 
      CDereference[CArray["arg", {}]]}},
    CBlock[{
      CDeclare["double", CAssign["num1", 20.4]],
      CDeclare["double", "num2"],
      CDeclare["WolframLibraryData", 
       CAssign["libData", 
        CCall["WolframLibraryData_new", {"WolframLibraryVersion"}]]],
      CCall["Initialize_compute", {"libData"}],
      CCall["compute", {"libData", "num1", CAddress["num2"]}],
      CCall["printf", {CString["%5.2f\n"], "num2"}],
      CCall["WolframLibraryData_free", {"libData"}],
      CReturn[0]}]]
   ];
cStr = ToCCodeString[code]

mainSource = Export[FileNameJoin[{outDir, "main.c"}], cStr, "Text"]

c = Compile[{{x}}, x^2];
CCodeGenerate[c, "compute", FileNameJoin[{outDir, "compute.c"}], 
  "CodeTarget" -> "WolframRTL"];
CCodeGenerate[c, "compute", FileNameJoin[{outDir, "compute.h"}], 
  "CodeTarget" -> "WolframRTLHeader"];

where

outDir

must be defined according the user's desires

The generated main.c file with a minor manual retouching in the printf function is


#include "compute.h"

#include "WolframRTL.h"


int main(int argc, char *arg[])
{
double num1 = 20.4;

double num2;
WolframLibraryData libData = WolframLibraryData_new(WolframLibraryVersion);
Initialize_compute(libData);

compute(libData, num1, &num2);
printf("%5.2f\n", num2);

WolframLibraryData_free(libData);
return 0;
}

The generated file compute.c is

#include "math.h"

#include "WolframRTL.h"

static WolframCompileLibrary_Functions funStructCompile;

static mbool initialize = 1;

#include "compute.h"

DLLEXPORT int Initialize_compute(WolframLibraryData libData)
{

if( initialize)
{
funStructCompile = libData->compileLibraryFunctions;
initialize = 0;
}

return 0;
}

DLLEXPORT void Uninitialize_compute(WolframLibraryData libData)
{
if( !initialize)
{

initialize = 1;
}
}

DLLEXPORT int compute(WolframLibraryData libData, mreal A1, mreal *Res)
{

mreal R0_0;
mreal R0_1;
R0_0 = A1;
R0_1 = R0_0 * R0_0;
*Res = R0_1;

funStructCompile->WolframLibraryData_cleanUp(libData, 1);
return 0;
}

The generated file compute .h is

#include "WolframLibrary.h"

EXTERN_C DLLEXPORT int Initialize_compute(WolframLibraryData libData);

EXTERN_C DLLEXPORT void Uninitialize_compute(WolframLibraryData libData);

EXTERN_C DLLEXPORT int compute(WolframLibraryData libData, mreal A1, mreal *Res);

Finally, the Makefile that takes care of the compilation is

CXX = /usr/bin/c++

MATHEMATICA_INCLUDE_PATH = /usr/local/Wolfram/Mathematica/8.0/SystemFiles/IncludeFiles/C

MATHEMATICA_LIB_PATH = /usr/local/Wolfram/Mathematica/8.0/SystemFiles/Libraries/Linux-x86-64


a.out: main.c compute.o
 ${CXX} -I${MATHEMATICA_INCLUDE_PATH} -L${MATHEMATICA_LIB_PATH} -lm -lWolframRTL -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -pthread compute.o main.c

compute.o: compute.c
 ${CXX} -c -I${MATHEMATICA_INCLUDE_PATH}  -L${MATHEMATICA_LIB_PATH}  compute.c

Of course, make sure that the libraries are in the right place and eventually you will need to add the path of the Mathematica libraries in ld.so.conf and apply /sbin/ldconfig

A more elaborated example taking a real matrix input and a real matrix output is here

A similar example implemented for complex matrices here

GDP per energy consumption

2011-04-20T23:47:00.000-07:00

Here, the world map for the GDP per unit of energy consumption. The world map has so much variation of this index that I had to make a log plot

In the case of the American continent, the linear scale is better

Energy consumption per capita

2011-04-17T01:46:00.000-07:00

I pulled out some data from Wolfram Alpha to make the following world map of the energy consumption per capita

A close up of America is also interesting

Natural gas consumption

2011-04-08T17:03:00.001-07:00

..............

Energy consumption by country

2011-04-07T15:22:00.000-07:00

This is my first widget from Wolfram Alpha

Wow!! so nice and simple to make

Finances and the rise of machines

2011-03-29T14:24:00.000-07:00

Last Friday March 25th I attended the Princeton University Conference on Quant Trading.

One of the most interesting talks was the analysis of the role of computers in trading
Rise of the Machines: Algorithmic Trading in the Foreign Exchange Market

The conclusion is that computers are becoming more and more prevalent and actually dominate the Forex market already.

Inverse Matrix by the Exchange Method

2011-03-25T16:04:00.000-07:00

After some time off, here I am again.

Here is an interesting implementation of the inverse of a matrix that I found useful when I was developing a CUDA program. This is the most efficient method to implement the inverse of a matrix in terms of memory usage, which is handy when we want to put everything in the fast but limited shared memory.

The Mathematica code is here

InverseExchange.cdf

The cuda code is next
InverseExchange.cu

J. Math. Phys Most Downloaded Articles in August 2010

2011-01-30T23:24:00.000-08:00

My paper about the connection between the Householder decomposition of unitary matrices and the canonical coset representation of unitary groups was the second most downloaded paper in August 2010

http://jmp.aip.org/features/most_downloaded?month=8&year=2010

Ito Calculus with Tensorial

2010-12-06T14:36:00.000-08:00

Here is a short example on how to use the Tensorial package for Ito calculus (stochastic calculus)

ito.pdf

Numerical Linear Algebra and Quantum Control

2010-12-06T14:05:00.000-08:00

Here is a poster about some of my recently published research

Poster pdf

MathLink example 4: Trace of a complex matrix

2010-12-05T21:22:00.000-08:00

In this example, I introduce the ability to treat complex numbers.

-------------------------------------------------------
myTrace.c
-------------------------------------------------------


#include<stdio.h>
#include<string.h>

#include "mathlink.h"
extern void myTrace( void );

void myTrace( void )
{

 int i,j,n;
  
 double *matrix;

 int *dims;
 char **heads;
 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  printf(" MLGetReal64Array error reading data \n");
 };

 n = dims[0];

 double resum=0.;
 double imsum=0.;

 double sum[2];

 int outdims[1];

 char **outheads;
 int outrank = 1; 

 
 for(i=0;i<n;i++){

 resum = resum + matrix[2*i*n + 2*i];

 imsum = imsum + matrix[2*i*n + 2*i + 1];
 }

 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);
 
 MLPutFunction(stdlink,"Complex",2);

 MLPutFloat(stdlink,resum);
 MLPutFloat(stdlink,imsum);


}


#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

-----------------------------------------------------
myTrace.tm
-----------------------------------------------------

:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Manual

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of complex matrix M"

------------------------------------------------------
Makefile
------------------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@

---------------------------------------------------
myTrace.nb
---------------------------------------------------

link = Install["./myTrace"]

(m = {{1., 2., 1.}, {2., 2., 2.}, {3., 3., 3.}} + I) // MatrixForm

myTrace[m]

MathLink example 3: Trace of a real matrix

2010-12-04T22:45:00.000-08:00

In this example we calculate the trace of a matrix of floating point numbers (double type in C)

---------------------------------------------------
myTrace.c
---------------------------------------------------


#include "mathlink.h"
extern double myTrace( void );

double myTrace( void )
{
 int i,j,n;
  
 double *matrix;
 int *dims;
 char **heads;

 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  return 0.;
 };

 n = dims[0];

 double sum=0;
 for(i=0;i<n;i++) sum = sum + matrix[i*n+i];
 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);

 return sum;
}

#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

----------------------------------------------
myTrace.tm
----------------------------------------------

:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of matrix M"

---------------------------------------------
Makefile
---------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@

Catenary

2010-11-24T23:28:00.000-08:00

The catenary is the curve that is naturally formed when a chain of constant linear density hangs from two fixed points. This is a typical example on variational calculus showing the use of Lagrange multipliers. Next, I develop the problem entirely in Mathematica

Catenary
Catenary(pdf)

You may need MathematicaPlayer if you do not have Mathematica.

This problem can be generalized to find the surface with minimum energy, which is shown next

MathLink Example 2: Sum a list of floating point numbers

2010-11-24T21:23:00.000-08:00

This example considers a sum of of floating point numbers.
----------------------------------------------------------------
addRealList.c
----------------------------------------------------------------

#include "mathlink.h"
extern double addRealList( void );



double addRealList( void )
{
 int i,n;

 double *list; 

 MLGetReal64List(stdlink,&list,&n);

 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
  MLReleaseReal64List(stdlink,list,n);

 return sum;
}



#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

------------------------------------------------------
addRealList.tm
------------------------------------------------------

:Begin:
:Function:       addRealList

:Pattern:        addRealList[ L:{___Real}]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: addRealList::usage = "addRealList[x] gives the sum of the elements on x, given that they are real (double floating numbers)"

------------------------------------------
Makefile
------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

addRealList : addRealListtm.c
 ${CXX} addRealListtm.c addRealList.c -o addRealList -lML64i3 -lm -lpthread -lrt -lstdc++

addRealListtm.c : addRealList.tm
 ${MPREP} addRealList.tm -o $@

clean:

 rm addRealListtm.c
 rm addRealList

--------------------------------------------------------
addRealList.nb
--------------------------------------------------------

link = Install["./addRealList"]

?addRealList

addRealList[{1., 5., 6.}]

MathLink Example 1: Sum a List of Integers

2010-11-23T21:31:00.000-08:00

This example passes a list of integers, performs the sum and returns an integer.
We must note that the length of the list is not type int but type long.

-------------------------------------------------------------------------
addList.c
-------------------------------------------------------------------------

#include "mathlink.h"
extern double addList( double *list, long n);



double addList( double *list, long n)
{

 int i;
 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
 return sum;
}



#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

---------------------------------------------------------------------
addList.tm
---------------------------------------------------------------------

:Begin:
:Function:       addList

:Pattern:        addList[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { RealList  }
:ReturnType:     Real

:End:

:Evaluate: addList::usage = "addList[x] gives the sum of the elements of the list x"

---------------------------------------------------------------------
Makefile
---------------------------------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

BINARIES = addList

all : $(BINARIES)

addList : addListtm.c addList.c
 ${CXX}  addList.c addListtm.c  -lML64i3 -lm -lpthread -lrt -lstdc++ -o $@

addListtm.c : addList.tm
 ${MPREP} addList.tm -o addListtm.c

-------------------------------------------
addList.nb
-------------------------------------------

link = Install["./addList"]

addList[{1, 5, 6}]

MathLink Example 0

2010-11-22T23:00:00.000-08:00

MathLink is a protocol on how to communicate Mathematica with external programs. For example, one would like to run compiled code in C within the Mathematica environment. I am going to post a sequence of MathLink programs for Linux with an increasing level of complexity, starting from a simple example already found in the Mathematica help (included for completeness) finishing with an example where I send and receive a matrix of complex numbers processed by the GPU using CUDA

Most examples are made of 4 files:

*.c : where the computation occurs.

*.tm : MathLink template where the input and output are specified and the documentation of the function is placed.
Makefile : where the compilation sequence is specified
*.nb : Mathematica notebook that executes the whole program

These files have to be placed in the same folder and compiled with: $make , which may need editing in order to provide the path of the executable file mprep that comes as part of Mathematica.

As seen in Makefile, mprep generates a proper *.c code from the *.tm file, which is compiled along with the *.c file that contains the computational process.

The first example adds two integers

---------------------------------------------------------------------------------
addtwo.c
---------------------------------------------------------------------------------

#include "mathlink.h"
extern int addtwo( int i, int j);



int addtwo( int i, int j)
{

 return i+j;
}



#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

-----------------------------------------------------------------
addtwo.tm
-----------------------------------------------------------------

int addtwo P(( int, int));

:Begin:
:Function:       addtwo

:Pattern:        AddTwo[i_Integer, j_Integer]
:Arguments:      { i, j }
:ArgumentTypes:  { Integer, Integer }
:ReturnType:     Integer

:End:

:Evaluate: AddTwo::usage = "AddTwo[x, y] gives the sum of two machine integers x and y."

------------------------------------------------------------
Makefile
------------------------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

BINARIES = addtwo

all : $(BINARIES)

addtwo : addtwotm.c addtwo.c
 ${CXX}  addtwotm.c addtwo.c  -lML64i3 -lm -lpthread -lrt -lstdc++ -o $@

addtwotm.c : addtwo.tm
 ${MPREP} addtwo.tm -o addtwotm.c

------------------------------------------------
addtwo.nb
------------------------------------------------

SetDirectory["/home/rcabrera/Documents/source/mathematica/MathLinkExamples/mytest2/alternative"]

link = Install["./addtwo"]

?AddTwo

AddTwo[2, 3]

Subversion

2010-10-14T14:03:00.000-07:00

Subversion is a version control software very useful for software developers
Here I found a nice video tutorial
http://showmedo.com/videotutorials/video?name=950000&fromSeriesID=95

Efimov states

2010-10-13T16:37:00.000-07:00

I've heard about Efimov states some time ago and now I decided to put them in my shopping list

http://www.sciencedaily.com/releases/2009/12/091211131526.htm

The Challenge to make quantum computers

2010-10-03T02:28:00.000-07:00

The challenge to implement a completely scalable quantum computer is tied with the understanding of the quantum/classical transition where quantum mechanics dominates the explanation of the microscopic world (up to the molecular scale more or less), but classical mechanics explains very well the ordinary macroscopic world. A completely scalable quantum computer would eventually give us a macroscopic taste of the imaginable strange quantum effects that so far have been seen only in the microscopic world.

"Because there are no known fundamental obstacles to such scalability (practical quantum computer with large number of qubits), it has been suggested that failure to achieve it would reveal new physics" -Emanuel Knill

I feel that the most recent papers are becoming more conservative about their predictions on the feasibility of quantum computing despite fact that there is work stating that fault tolerant quantum computation is possible with two ingredients:

Maximum Error/Gate about 10^-4 to 10^-5
Effective error correction codes with ancillary qubits.

The last ingredient seems to be more or less accomplished and the former one does not seem to be fundamentally unattainable despite the fact that we are currently very far. A paper that summarizes these facts with a high degree of scepticism is

Is Fault-Tolerant Quantum Computation Really Possible?

which I also like for its entertaining and straightforward writing style.

Humble opinion about God and science

2010-09-21T01:08:00.000-07:00

For some reason many physicists in the past and present such as Einstein and Hawking used the term God when they wanted to describe something essential and fundamental such as the laws of physics or the Big Bang theory. For most people this was an indication of their religiosity. However, anybody who knows the context of those words can tell that truth is different. Here is an article explaining this issue.

My personal point of view is that science does not require God by definition and not because science can prove that God does not exist. Defining science as the rational understanding of nature based on fundamental laws (and axioms), the acknowledgement of the existence of God would contradict this premise. Using God to explain something in science would be equivalent to admit an irrational element in the theory. In other words, admitting something that is not based on fundamental laws or axioms.

However, I think there is still the inconsistent possibility to put God outside science as the something beyond the laws of physics and axioms. Some would place God as the creator of the laws of physics and axioms. I do not see any consistent and rational way to put God as the explanation of absolutely anything. Having said that, I do not consider myself as completely consistent and rational being and I do not think anybody was, is or will ever be, so, here it lies a window for the endless debate for the existence of God.

Will we ever find a rational and consistent explanation for our consciousness and existence? Something tells me that the most likely answer is NOT.

Inverse Tomography

2010-08-25T05:43:00.000-07:00

A few years ago I had the pleasure to attend Pedro Goldman's presentation (Argentinian scientist) about his research in inverse tomography and applications to the calculation of the optimal dosage of radiation for radiotherapy.

Fast Optimization or the Radiation Therapy of Tumors -the Impossible Possible

Sitting down on my chair I would not expect that this was going to be one of the most shocking and interesting presentations I ever attended.

This research field is mathematically interesting and extremely important for medical physics.

Computational Law

2010-08-20T12:22:00.000-07:00

I did not know about this this research field until my last trip to Beijing China where I met Set Chandler, who develops his programs in Mathematica

http://www.wolfram.com/products/mathematica/portraits/sethchandler/

The general idea is to develop models as response of the implementation of the laws. The idea is to minimize the subjectivity in today's approach.

A much more ambitious project is to analyse the self consistency of the laws and their level of complexity. The complexity does not only depend on the document's length but also in how intricate is the relation with itself and other laws.

Videos created in Mathematica

2010-08-18T15:31:00.000-07:00

I just learned about a video channel for videos created in Mathematica.

http://vimeo.com/channels/mathematica

Medical Imaging

2010-08-10T17:32:00.000-07:00

Medical Imaging is today a very active multidisciplinary research field where differential geometry, high performance computing and diverse mathematics meet together.

In my trip to China I had the opportunity to see some of the latests developments at professor Bart M. ter Haar Romeny's group at the University of Eindhoven in Holland.

Some of the lastest thesis in his group that can be donwloaded are

E.M. Franken, Enhancement of Crossing Elongated Structures in Images, PhD. Thesis, 2008, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romenij, R. Duits

M.A.van Almsick, Context Models of Lines and Contours, PhD. Thesis, 2007, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romeny, L.M.J. Florack

E. Balmashnova, Scale-Euclidean Invariant Object Retrieval, PhD. Thesis, 2007, Technische Universiteit Eindhoven.Advisors: B.M. ter Haar Romenij,co-advisor: L.M.J. Florack

F.M.W. Kanters, Towards Object Based Image Editing, PhD. Thesis, 2007, TU/e. Advisors: B.M. ter Haar Romenij, L.M.J. Florack

B. Platel, Exploring the Deep Structure of Images, PhD. Thesis, 2007, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romenij, Co-advisor: L.M.J. Florack

Additionally, professor Romeny has a book entitled

Front-End Vision and Multi-Scale Image Analysis: Multi-scale Computer Vision Theory and Applications, written in Mathematica

Online version

The complete package can be found at
http://extras.springer.com/2003/978-1-4020-1507-6

If you are not afraid of modestly advanced mathematics and want something really applied, try this.

The canonical coset decomposition of unitary matrices through Householder transformations

2010-08-03T14:35:00.000-07:00

My paper describing the connection between the Householder decomposition of and the canonical coset parametrization of unitary operators was published today in the Journal of Mathematical Physics

J. Math. Phys. 51, 082101 (2010)

http://arxiv.org/abs/1008.2477

This parametrization is very useful for calculating the parameters that uniquely characterises a given unitary operator U(N). I show that instead of solving a complicated systems of non-linear equations one can perform very stable and efficient Householder transformations. Right now I am working on some practical applications based on that result.

Independent Component Analysis

2010-07-30T08:24:00.000-07:00

I was already familiar with the Principal component analysis, which is a method for finding a linear combination of variables such that the possible correlation between them disappears i.e. elimination of the diagonal elements of the covariance matrix. Statistical independence implies the elimination of correlation but the opposite is not necessarily true because there is possible association between the variables on higher orders. The Independent Component Analysis seems to be the generalization needed to go beyond the simple correlation.

Bartlett, M.S. and Movellan, J.R. and Sejnowski, T.J., Face recognition by independent component analysis, IEEE Transactions on neural networks,13, pp 1450--1464, 2002

Hyvarinen, A. and Oja, E. Independent component analysis: a tutorial, Neural Networks, 13, pp 411--430 2000

10th International Mathematica Symposium

2010-07-30T07:20:00.000-07:00

I just returned from the 10th International Mathematica Symposium held in Beijing China, where I gave a talk. It was a great opportunity to see the work of other people around the world and the new features of Mathematica 8.

Some of the topics included medical imaging, stochastic simulations, video editing, real-time control, automatic C code generation, CUDA and more.

Quantum entanglement and classical chaos

2010-07-07T16:31:00.000-07:00

I recently attended a talk about the connection between quantum entanglement and classical chaos

http://www.nature.com/nature/journal/v461/n7265/full/nature08396.html

This tells me that the underlying reason of the sensitivity to initial conditions in classical mechanics may be the rise of entanglement in a sequence with the collapse of the wave function.

Venn diagrams in R

2010-07-05T17:05:00.000-07:00

This is my first attempt to use R, the well known program for statistics.

I am using the package LIMMA and the instructions found here in order to draw Venn diagrams. My code is

-----------------------------------------
library(limma)

data1<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h1.csv", sep=',', header=T)
data2<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h2.csv", sep=',', header=T)
data3<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h3.csv", sep=',', header=T)

set1<-data1$x
set2<-data2$x
set3<-data3$x

set123<-union(set1,union(set2,set3))

bool1 = logical(length(set123))
bool2 = logical(length(set123))
bool3 = logical(length(set123))

for(i in 1:length(set123)) bool1[i]<-is.element(set123[i],set1)
for(i in 1:length(set123)) bool2[i]<-is.element(set123[i],set2)
for(i in 1:length(set123)) bool3[i]<-is.element(set123[i],set3)

c3<-cbind(bool1, bool2, bool3)

a <- vennCounts(c3)

vennDiagram(a)
-----------------------------------------------

Note.- The data is stored in a column with label "x"

Mersenne Twister in CUDA

2010-07-02T21:05:00.000-07:00

The Mersenne Twister algorithm (what a cool name) for generating random numbers is implemented as an example in the SDK folder provided by NVIDIA. However, as it happens with all the examples, it depends on the CUTIL library, which is not part of the standard CUDA API. There is another problem, the compilation of each example is carried out by a general set of makefiles customized for compiling all the examples. These are big problems if we want to generate a stand alone version of the examples. This is desirable because in this way we could adapt the code for our own purposes and eventually link with extra code.

In the rest of this post I will give the makefiles that accomplish the stand alone compilation independent of CUTIL for the Mersenne Twister example.

You will have to follow the following steps :

- Make a copy of the Mersenne Twister folder from the SDK (I renamed as MersenneTwisterCUDA),

- Remove all the references for the CUTIL library

- Add the following lines on MersenneTwister.cu

char raw_path[] = "/YOUR_PATH/MersenneTwisterCUDA/data/MersenneTwister.raw";
char dat_path[] = "/YOUR_PATH/MersenneTwisterCUDA/data/MersenneTwister.dat";

while removing the corresponding lines that originally define raw_path and dat_path.

- Finally, paste the make files I will provide and run $make on the command line
(Define YOUR_PATH according to the structure of your directories)

file: common
----------------------------------------------------------------

CUDA_PATH = /usr/local/cuda

C := gcc
CC := c++
NVCC := $(CUDA_PATH)/bin/nvcc

CUDA_INCLUDE_PATH = /usr/local/cuda/include
CUDA_LIB_PATH := $(CUDA_PATH)/lib64
CUDA_FLAGS = -lcuda -lcudart

===================================================

file: Makefile:
----------------------------------------------------------------------------
include ./common.mk

MersenneTwister.out: genmtrand.c.o MersenneTwister_gold.cpp.o MersenneTwister.cu.o
$(NVCC) -L$(CUDA_LIB_PATH) $(CUDA_FLAGS) genmtrand.c.o MersenneTwister_gold.cpp.o MersenneTwister.cu.o -o MersenneTwister.out

MersenneTwister.cu.o: MersenneTwister.cu dci.h MersenneTwister.h MersenneTwister_kernel.cu
$(NVCC) -c -I$(CUDA_INCLUDE_PATH) MersenneTwister.cu -o MersenneTwister.cu.o

MersenneTwister_gold.cpp.o: MersenneTwister_gold.cpp
$(CC) -c MersenneTwister_gold.cpp -o MersenneTwister_gold.cpp.o

genmtrand.c.o: genmtrand.c
$(C) -c genmtrand.c -o genmtrand.c.o

Classical Origin of Fermion Spin

2010-05-30T21:52:00.000-07:00

The spin was postulated by Pauli from experimental evidence, but it was only with the arrival of the Dirac equation that the spin appears naturally. This leaded to many people to consider the spin as fundamentally "quantum". In the following paper we argue that the spin appears naturally from classical relativistic mechanics alone

W. E. Baylis, R Cabrera, J.D. Keselica, Quantum/Classical Interface: Classical Geometric Origin of Fermion Spin, Advances in Applied Clifford Algebras, 2010

Maxwell's demon

2010-05-12T22:23:00.000-07:00

Mark G Raizen recently gave a talk about techniques for trapping and cooling atoms. In this way, he actually implemented Maxwell's demon in practice!

http://www.sciencemag.org/cgi/content/abstract/324/5933/1403

Yes, this demon exists and does not violate the second law of thermodynamics because information is a thermodynamic variable with an essential role.

Attention span

2010-05-05T13:58:00.000-07:00

The amount of information we can absorb not only depends on the time we actually engage our attention. It also depends on how we distribute this time in subintervals. It seems that we can only pay our highest degree of attention for a few seconds [wikipedia article], so what do we do in the rest?

Related to this topic is the strategy followed by Hollywood, who are interested in maintaining their movies as appealing as possible. In this next article there is an interesting analysis of the duration and distribution of the shoot lengths.

James E. Cutting, Attention and the Evolution of Hollywood Film Psychological Science, 2010

One of the main conclusions is that the distribution of power in the frequency domain must obey 1/f , where f is the frequency. I am sure this is particularly important for teachers and students.

More comments about this article can be found here.

Eigenvalues: C/Lapack

2010-05-02T20:22:00.000-07:00

Here is another example on the use of Lapack. This time the objective is to calculate the eigenvectors and eigenvalues of a complex matrix

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//.......................................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

  for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//......................................................................................................
//  MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
//  The eigenvectors are stored in columns
//.....................................................................................................
void MatrixComplexEigensystem( double complex *eigenvectorsVR, double complex *eigenvaluesW, double complex *A, int N)
{

 int i;
  
  double complex *AT = (double complex*) malloc( N*N*sizeof(double complex) );

  zgeTranspose( AT, A , N);
  
  char JOBVL ='N';   // Compute Right eigenvectors

  char JOBVR ='V';   // Do not compute Left eigenvectors

  double complex VL[1];

  int LDVL = 1; 
  int LDVR = N;

 int LWORK = 4*N; 

 double complex *WORK =  (double complex*)malloc( LWORK*sizeof(double complex));

 double complex *RWORK = (double complex*)malloc( 2*N*sizeof(double complex));

 
 int INFO;

 zgeev_( &JOBVL, &JOBVR, &N, AT ,  &N , eigenvaluesW ,

   VL, &LDVL,
   eigenvectorsVR, &LDVR, 
   WORK, 
   &LWORK, RWORK, &INFO );

 zgeTranspose( AT, eigenvectorsVR , N);

 for(i=0;i<N*N;i++) eigenvectorsVR[i]=AT[i];

 
  free(WORK);
  free(RWORK);
  free(AT);
}


int main()
{
  int i,j;
  const int N = 3;

  
  double complex A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

  double complex eigenVectors[N*N];
  double complex eigenValues[N];

  
  MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);
  
  printf("\nEigenvectors\n");

  for(i=0;i<N;i++){
    for(j=0;j<N;j++) printf(" (%f,%f) \t", eigenVectors[i*N + j]);

    printf("\n");
  }
  
  printf("\nEigenvalues \n");
  for(i=0;i<N;i++) printf("\n (%f, %f) \t",  eigenValues[i] );

  
printf("\n------------------------------------------------------------\n"); 
return 0;

}

Complex matrix inverse: C++/Lapack

2010-05-02T01:58:00.000-07:00

So, here is an example on how to call lapack from c++

#include<iostream>
#include<math.h>

#include<complex>
#include <stdlib.h>

using namespace std;

extern "C" void zgetrf_( int*, int* , complex<double>* , int*, int* , int* );

extern "C" void zgetri_( int*, complex<double>* , int*, int* , complex<double>*, int* , int* );

//........................................................................................
void zgeTranspose( complex<double> *Transposed, complex<double> *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

 for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//.........................................................................................
void MatrixComplexInverse(complex<double> *invA, complex<double> *A, int n)
{

 int LWORK=10*n;

 int *permutations;

 complex<double> *WORK, *tempA;

 tempA = new complex<double>[n*n];

 permutations =  new int[2*n];
 WORK = new complex<double>[n*n];


 int INFO;

 zgeTranspose(tempA,A,n);


 zgetrf_( &n, &n, tempA , &n, permutations , &INFO );

 if (INFO != 0) {
   cout<<"ComplexMatrixInverse: Error at zgetrf  \n"; exit(0);
   }



 zgetri_( &n, tempA , &n, permutations , WORK, &LWORK, &INFO );

 if (INFO != 0) {
   cout<<"ComplexMatrixInverse: Error at zgetri  \n"; exit(0);
   }

 zgeTranspose(invA,tempA,n);

 delete [] WORK;

 delete [] tempA;
 delete [] permutations;

}

/////////////////////////////////////////////////////////////////////////////////////////

int main()
{
 int i,j;
 const int N = 3;

 complex<double> I(0.,1.);

 complex<double> A[] = { 1. + I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

 complex<double> invA[N*N];

 MatrixComplexInverse(invA,A,N);


 for(i=0;i<N;i++){
   for(j=0;j<N;j++) cout << invA[i*N + j]<<"\t";

   cout<<"\n";
 }

cout<<"---------------------------\n";

return 0;

}

////////////////////=====================////////////////////

and the make file is

#
CC = c++

#edit LAPACK_PATH if necessary
LAPACK_PATH = /usr/lib64/atlas

a,out: main.cpp
 $(CC)  main.cpp  -L$(LAPACK_PATH) -llapack -lblas  -lgfortran -lm

My basic rules of programming

2010-04-28T00:51:00.001-07:00

Here I am listing some general rules of programming that I learned to follow.
Some of them are enforced by certain languages or programming environments but some are very permissive and require to be more careful. These rules are already well known but here they are again

Define very well the input from the output arguments of a function. The input arguments should remain unmodified in the routine.
Never write a routine that cannot fit in one page of computer screen in the worst case.
Develop and validate a program in reasonable small steps. It may take longer to write the code but it will take less time to debug.
Use a maximum of two nested loops in a certain routine.
Use meaningful names for the variables even if they become long.
Do not repeat the same task in different places.

Complex matrix inverse: C/Lapack

2010-04-26T23:09:00.000-07:00

Here a simple example on how to use Lapack within C to invert matrices using the LU decomposition.

The C++ code is different and I will post that another day

[The html code for this post was generated by http://www.bedaux.net/cpp2html/]

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//........................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//.........................................................................................
void MatrixComplexInverse(double complex *invA, double complex *A, int n)
{

int LWORK=10*n;

int *permutations;

double complex *WORK, *tempA;

tempA = (double complex*) malloc( n*n*sizeof(double complex) );

permutations = (int*) malloc( 2*n*sizeof(int) );

WORK = (double complex *)malloc(LWORK*sizeof(double complex));


int INFO;

zgeTranspose(tempA,A,n);


zgetrf_( &n, &n, tempA , &n, permutations , &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetrf  \n"); exit(0);
 }



zgetri_( &n, tempA , &n, permutations , WORK, &LWORK, &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetri  \n"); exit(0);
 }

zgeTranspose(invA,tempA,n);

free(WORK);

free(tempA);
free(permutations);

}


/////////////////////////////////////////////////////////////////////////////////////////


int main()
{
int i,j;
const int N = 3;


double complex A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

double complex invA[N*N];

MatrixComplexInverse(invA,A,N);


for(i=0;i<N;i++){
 for(j=0;j<N;j++) printf(" (%f,%f) \t", invA[i*N + j]);

 printf("\n");
}

printf("------------------------------------------------------------\n");
return 0;

}


The make file is

#
CC = gcc

#edit LAPACK_PATH if necessary
LAPACK_PATH = /usr/lib64/atlas

a,out: main.c
$(CC)  main.c  -L$(LAPACK_PATH) -llapack -lblas  -lgfortran -lm

Complex matrix inverse III: CUDA

2010-04-25T15:31:00.000-07:00

Here is cuinverse_kernel.cu

__global__ void cgeMatrixInverse_kernel(cuFloatComplex *invA , cuFloatComplex *A , int N , cuFloatComplex *Work)
{

int i;
int idx = blockDim.x * blockIdx.x + threadIdx.x;

cuFloatComplex * ThreadWorkSpace;

ThreadWorkSpace = Work + idx*cgeMatrixInverse_WorkSpace()*N*N;

for(i=0;i<N*N;i++) A[ i + idx*N*N ] =  A[i];

A[ idx*N*N ] = make_cuFloatComplex( (float) idx , 1./sqrtf(2));

cgeMatrixInverse(invA + idx*N*N , A + idx*N*N , N , ThreadWorkSpace);


}

_______________________

and my makefile is

#  edit the path of your working directory
BASE_PATH = /home/rcabrera/Documents/source/c/inverseCUDA

#

CUDA_PATH = /usr/local/cuda

CC = gcc
NVCC    := $(CUDA_PATH)/bin/nvcc

CUDA_INCLUDE_PATH = /usr/local/cuda/include
CUDA_LIB_PATH := $(CUDA_PATH)/lib64
CUDA_LIBS = -lcuda -lcudart

a.out: main.cu cuinverse_kernel.cu cumatrixtools.h
$(NVCC)  -I$(CUDA_INCLUDE_PATH) -L$(CUDA_LIB_PATH) $(CUDA_LIBS) -lm  main.cu -o a.out

# $(NVCC) -deviceemu  -I$(CUDA_INCLUDE_PATH) -L$(CUDA_LIB_PATH) $(CUDA_LIBS) -lm  main.cu -o a.out

clean :
rm -f a.out

Complex matrix inverse II: CUDA

2010-04-25T15:28:00.000-07:00

So, here is cumatrixtools.h

#include<cuComplex.h>
#ifndef _CUMATRIXTOOLS_H_
#define _CUMATRIXTOOLS_H_


////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void
cgeSquareMatrixProduct(cuFloatComplex *out, cuFloatComplex *A, cuFloatComplex *B, int N  )
{

int i,j,k;
cuFloatComplex sum;

for(i=0;i<N;i++)

for(j=0;j<N;j++){
  //sum = 0. + 0.*I;
    sum  = make_cuFloatComplex(0.,0.);

  for(k=0;k<N;k++){
    sum =  cuCaddf(sum  , cuCmulf(A[i*N+k],B[k*N+j]) );
  }

  out[i*N+j] = sum;
}


}
//////////////////////////////////////////////////////////////////////

__device__ __host__ static __inline__ void cgeTranspose(cuFloatComplex  *At, cuFloatComplex  *A, int N)
{

int i,j;
for( i = 0; i<N; i++)

for( j = 0; j<N; j++)
  At[i+j*N] = A[i*N+j] ;
}

//////////////////////////////////////////////////////////////////////////////////////
//                LU decomposition of complex matrices
/////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeDoolittle_LU_Decomposition(cuFloatComplex *LU, cuFloatComplex *A, int n)
{

 int i, j, k, p;
 cuFloatComplex *p_k, *p_row, *p_col;

for(k=0; k<n*n; k++) LU[k]=A[k];


for (k = 0, p_k = LU; k < n; p_k += n, k++) {

    for (j = k; j < n; j++) {

       for (p = 0, p_col = LU; p < k; p_col += n,  p++)

//           *(p_k + j) -= *(p_k + p) * *(p_col + j);
       *(p_k + j) = cuCsubf( *(p_k + j) , cuCmulf( *(p_k + p) , *(p_col + j) ));
    }

    if ( cuCabsf(*(p_k + k)) != 0.0 )  //return -1;


    for (i = k+1, p_row = p_k + n; i < n; p_row += n, i++) {

       for (p = 0, p_col = LU; p < k; p_col += n, p++)

         // *(p_row + k) -= *(p_row + p) * *(p_col + k);
    *(p_row + k) = cuCsubf( *(p_row + k) , cuCmulf(  *(p_row + p) , *(p_col + k) ));
         //*(p_row + k) /= *(p_k + k);


  *(p_row + k) = cuCdivf( *(p_row + k) , *(p_k + k)  );
    }
 }

}

//////////////////////////////////////////////////////////////////////////////////////////////////
// Back substitution for lower triangular matrices assuming that the diagonal is filled with 1's
// Given  T x = b ,
// where T is a lower triangula matrix NxN with 1's in the diagonal
// b is a known vector
// x is an unknown vector
// the routine otputs x = xout
//////////////////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeBackSubstitutionLowerDiagOne(cuFloatComplex *xout, cuFloatComplex *T , cuFloatComplex *b ,int N  )
{

cuFloatComplex bTemp;
int i,j;

for(i=0;i<N;i++){

  bTemp = b[i];
  for(j=0;j<i;j++) bTemp = cuCsubf( bTemp , cuCmulf(T[ i*N + j], xout[j] ) );

  xout[i] = bTemp ;    
}
}

/////////////////////////////////////////////////////////////////////////////////////
// Back substitution for upper triangular matrice
////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeBackSubstitutionUpper(cuFloatComplex *xout, cuFloatComplex *T , cuFloatComplex *b ,int N  )
{

cuFloatComplex bTemp;
int i,j;

for(i=N-1;i>=0;i--){

  bTemp = b[i];
  for(j=i+1;j<N;j++) bTemp = cuCsubf( bTemp , cuCmulf(T[ i*N + j], xout[j] ) );

  xout[i] = cuCdivf( bTemp , T[ i*N + i] );    
}
}

///////////////////////
__device__ __host__ static __inline__ void cgeIdentity(cuFloatComplex *Id, int N)
{

int i,j;
for(i=0;i<N;i++)

for(j=0;j<N;j++){
Id[i*N+j] = make_cuFloatComplex(0.f,0.f);
}


for(i=0;i<N;i++){
  Id[i*N+i] = make_cuFloatComplex(1.f,0.f);
}
}

//////////////////////////////////////////////////////////////////////////////////////
//  Inverse of a matrix using the triangularization of matrices
//  Warning:
//          It does not destroys the original matrix A
//          A work space for 3 complex matrices must be supplied in W
//////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ int cgeMatrixInverse_WorkSpace(){ return 3;}

__device__ __host__ static __inline__ void cgeMatrixInverse(cuFloatComplex *inv , cuFloatComplex *A , int N , cuFloatComplex *W)
{

int i;
cuFloatComplex *Id;
cuFloatComplex *invL, *invU;


Id = W;      //Double purpose work space
  invL = W + N*N;

invU = W + 2*N*N;


cgeIdentity( Id, N);

 
 cgeDoolittle_LU_Decomposition(inv,A,N);


for(i=0;i<N;i++)   cgeBackSubstitutionLowerDiagOne( invL + i*N , inv , Id + i*N ,  N  );


for(i=0;i<N;i++) cgeBackSubstitutionUpper( invU + i*N , inv , Id + i*N ,  N  );


cgeSquareMatrixProduct(Id , invL , invU , N  );


cgeTranspose(inv , Id , N);

}


#endif

Complex matrix Inverse I: CUDA

2010-04-17T21:58:00.000-07:00

Here is a short program in CUDA to invert complex matrices through the LU decomposition (without permutations) with back substitution in __device__ mode. This means that the whole inversion is designed to run in each thread.

This code was adapted from a C code in real double precision found here

The code is standalone in the sense that it only requires the Standard tools provided by CUDA (without CuBlas) and tested in the NVIDIA GeForce 240M
on a Linux Fedora 11.

The program can be organized in three files

cumatrixtools.h: which contains the method itself.
cuinverse_kernel.cu: which contains the kernel function, which maps the method on the threreads.
main.cu: which administers the memory and calls kernel function.

The purpose of this program is to invert several variations of a matrix A. The matrix A is passed to the kernel and the inverse is evaluated after modifying the first element of A. Each branch executes a different matrix. The kernel is called with <<<2,32>>>, which means that the threads are organized in two blocks, each bloch with 32 threads. The more physical processors one has, the more threads per block one could assign in order to gain performance

The whole program is distributed in three posts (I,II and III)

main.cu looks like

#include <stdio.h>
#include <complex.h>

#include<cuComplex.h>

#include"cumatrixtools.h"
#include"cuinverse_kernel.cu"

void cuPrintMatrix( cuFloatComplex  *C ,int N, int M )
{

int i,j;
for(i=0;i<N;i++){

for(j=0;j<M;j++) printf(" (%f,%f)\t ", cuCrealf(C[i*N + j]) , cuCimagf(C[i*N + j]) );

printf(" \n ");
}
}

///////////////////////////////////////////////////////////////////////////////
const int N=3;     //Matrix dimension

const int NBranches = 64;
///////////////////////////////////////////////////////////////////////////////
// Main program
///////////////////////////////////////////////////////////////////////////////
int main(int argc, char **argv){

int i;
float  complex A[9] =  //  Base matrix
    {

0.f+I/sqrtf(2.f), 0.0f+ I/sqrt(2.0f),  0.+ 0.*I,

0. -I/2.     , 0.+I/2       ,   0.+ I/sqrt(2.),
-0.5+0.*I ,  0.5+0.*I  ,  -1./sqrt(2.)+.0*I

};

cuFloatComplex  *h_A, *h_invA;
cuFloatComplex  *d_A, *d_invA,   *d_WorkSpace;


printf("...allocating CPU memory.\n");
h_A =          (cuFloatComplex *)   malloc(                              N*N*sizeof(cuFloatComplex ));
h_invA =       (cuFloatComplex *)   malloc(                    NBranches*N*N*sizeof(cuFloatComplex ));

printf("...allocating GPU memory.\n");

cudaMalloc((void **)&d_A,                                      NBranches*N*N*sizeof(cuFloatComplex ));

cudaMalloc((void **)&d_invA,                                   NBranches*N*N*sizeof(cuFloatComplex ));

cudaMalloc((void **)&d_WorkSpace, NBranches*cgeMatrixInverse_WorkSpace()*N*N*sizeof(cuFloatComplex ));


printf("...Copying memory.\n ");
for(i=0;i<N*N;i++ ) h_A[i] = make_cuFloatComplex( crealf(A[i]) , cimagf(A[i])  );

cudaMemcpy(d_A, h_A, N*N*sizeof(cuFloatComplex) , cudaMemcpyHostToDevice);


printf("...The base matrix is:\n");
cuPrintMatrix( h_A , N, N );


printf("\n...Calling the kernel.\n");

cudaThreadSynchronize();
cgeMatrixInverse_kernel<<<2,32>>>(d_invA, d_A , N ,d_WorkSpace);   // Divinding the 64 branches in 2 blocks of 32 threads

    cudaThreadSynchronize();

cudaMemcpy(h_invA, d_invA, NBranches*N*N*sizeof(float), cudaMemcpyDeviceToHost);


printf("\n The inverse of the first branch is \n");
cuPrintMatrix( h_invA , N, N );


printf("\n The inverse of the second branch is \n");
cuPrintMatrix( h_invA + N*N , N, N );


printf("\n and so on ..\n");

free(h_A);
free(h_invA);


cudaFree(d_A);
cudaFree(d_invA);

cudaThreadExit();
printf("\n-------------------------------------------------------\n");
}

Multidimensional Arrays

2010-02-01T21:17:00.000-08:00

Linear algebra is very well developed with a vast literature in numerical methods. This not true for multi-dimensional arrays higher than 2, which correspond to the familiar matrices. In general it is difficult to generalize matrix operations to higher dimensions. One paper that gives an introduction to the topic and treats the problem of diagonalizing higher dimensional arrays is
Tensor Diagonalization

The literature usually refers multi-dimensional arrays as as Tensors, but I prefer to use the term array, because in physics a tensor is more than a multidimensional array.

The following figures correspond to an 9x9x9 array after and before a diagonalization attempt, with the application of orthogonal transformations.

Benford's law

2009-12-19T06:33:00.000-08:00

This is the probabilistic law that rules the distribution of de digits that appear in for example, our bills, bank statements, national debts, ..., where the digit 1 is seen as the most probable.

The reason of this counter-intuitive law is that the natural underlying distribution is logarithmic. Benford's law

When probabilities are involved, our intuition may catastrophically lead us to wrong conclusions.

The Black-Scholes equation

2009-12-13T20:29:00.000-08:00

I do not remember the first time I heard about the Black-Scholes equation, but I finally fulfilled my curiosity. The Black-Scholes equation is a differential equation for the value of a contract that establishes the optional right to buy (call options) or sell stocks (put options). The typical call option is a contract that gives the option to buy the stock at a certain point in the future T for a given strike price k. The value of the stock changes in time so that if I buy a call option written with a strike price k, I am betting that the price of the stock at time T will be higher than k (plus the cost of the contract itself), so that I will be able to execute the transaction recovering the difference. Otherwise, I could let the contract expire, losing all the money I paid for the call option.

The value of the contract is well know in the last day T, when the option can be executed because it is exactly the difference of price of the stock minus the strike price of the contract if the stock price is higher than k, otherwise, the contract's value is ZERO. The big problem is to estimate the fair value of the contract as a function of the stock price S and time t. Under certain assumptions, the value of the contract obeys the Black-Scholes partial differential equation.

The figure below shows a typical solution of the Black-Scholes partial differential equation. The curve in red corresponds to the value of the contract at time T as a function of the stock price, acting as boundary condition. The remaining curves correspond to times farther and farther in the past. The strike price is k=120, the maturity time is T=1, the interest rate is r=0.05 and the volatility is sigma=0.5

This is an example adapted from the book:
Computational Financial Mathematics using Mathematica, by Srdjan Stojanovic.

Now, having a reasonable idea of the value of the contract, I can even buy and sell contracts!.

Functional Programming: MapReduce

2009-12-08T17:05:00.000-08:00

Functional programming is a programming paradigm that was introduced by LISP. Many modern languages computational systems were inspired by LISP and one them is Mathematica.

It seems that the next biggest impact of functional programming is going to come from the application to parallel/distributed programming. One example of this emerging technology is the MapReduce framework developed by Google. In similar way, Yahoo is developing Haddop for the same purposes.
Another example is of course the implementation of WolframAlpha mostly developed in Mathematica.

There are many tutorials available including youtube videos such as
MapReduce Cluster Computing

One of the implementations that called my attention is MARS, which is developed on the top of CUDA.

Computational Linear Algebra

2009-12-04T20:11:00.000-08:00

How do we implement programs that require linear algebra?
For most of the cases there are no reasons to spent effort in developing our own libraries for linear algebra because we can use very good and freely available libraries.

The first layer is the Fortran BLAS library that implements matrix products. There are many variants of BLAS developed by certain companies and research institutions.
The second layer is the Fortran Lapack library, which implements more advanced routines such as matrix decompositions and computation of eigenvalues on the top of BLAS. However, there are no high level routines such as matrix inverse, determinants, etc. The reason of the absence of these high level routines is that there are many ways to implement them in terms of the Lapack routines but one has to choose a particular method according to the personal requirements for maximum efficiency. The implementation of programs in Lapack has high potential for efficiency optimization because one can specifies the type of matrix in the operation ie. real, complex, symmetric, etc and the decomposition that fits better for the final purpose.
The third layer is the Blas/Lapack wrapper, which implements the operations that are taught in a first curse of Linear algebra. The routines in Lapack are usually taught in a second course of linear algebra. The wrappers are easier to use but some of the flexibility is lost.
For those who work with C, GSL must be the best option. For C++, one has the choice to use the Boost library, which includes a lot more than linear algebra. Personally, the library I really feel as the most friendly is the Armadillo C++ library for linear algebra.

Magnetic sensors

2009-11-19T01:02:00.000-08:00

One of the most important technologies with the most important applications is the ability to detect small magnetic fields. This allows us to make more compact hard disks and more recently more sensitive detectors, including bio-detectors.

Gigant Magneto Resistance (GMR) is one of the important effects responsible for these technologies that exploit the difference of electrical resistivity when the current traverses two different magnetic materials (separated by a very thin non-magnetic material) and the spins of the materials are parallel or anti-parallel. Giant Magnetoresistance: The Really Big Idea Behind a Very Tiny Tool

The Tunnel Magnetoresistance effect (TM) , which relies in the quantum tunnel effect, is another important similar effect that can be even more sensitive. Tunnel Magnetoresistance effect

There transmission of light in a gas such as Rubidium can be also very sensitive to very small external magnetic fields. Laser magnetometer

A link with some of these technologies applied to sensors is
GMR for sensors

Geometric Optics with Lie Groups

2009-10-28T23:52:00.001-07:00

Donald Barnharth, from Optical Software introduced me with
another very interesting application of Lie groups in geometric optics.

This should not be a surprise because geometric optics as well as classical mechanics can be expressed in terms of a variational principle. In classical mechanics we have the minimal action and in geometric optics we have Fermat's principle, which states that the trajectory of light rays minimize the optical path length defined as

where n is the refraction index, ds is the differential arc length and C is the path. The whole machinery of classical mechanics can be translated to geometric optics with the corresponding Hamiltonian formulation as the phase space formulation of geometric optics.

The key idea of Hamiltonian mechanics is that the trajectory of the particles can be seen as active continuous symplectic transformations. For example, one has the following time-evolution operator

where H is the Hamiltonian (independent of time), with the Lie operator defined in terms of the Poisson brackets as

and

The same formalism can be applied to geometric optics with some particular adjustments. The time-evolution operator is replaced by the transformation that propagates the optical phase-state in the space. For practical applications this transformation is factorized in a perturbative-like product expansion that reminds the Fer expansion as

where, the first factor represents the para-axial approximation while the rest represent the corresponding higher order corrections.

References

LaTeX was powered by MathTran

Relativistic Many-body dynamics

2009-10-25T21:08:00.000-07:00

An intuitive generalization of the explicitly covariant single-body relativistic dynamical equation to multiple multiple interacting particles is difficult. The no-go theorem of Currie Jordan and Sudarshan is an example of the difficulties to device such explicitly covariant version of relativistic interacting particles. I do not know if a satisfactory elegant solution was found.

A book that explains this problem is

FROM CLASSICAL TO
QUANTUM MECHANICS
Giampiero Esposito, Giuseppe Marmo

Also see
Form of relativistic dynamics with world lines

The Lorentz-Dirac Force

2009-10-25T20:18:00.000-07:00

There is a naive general idea that we already understand classical mechanics very well, but I do not think that is the case. There are many fundamental questions in classical mechanics without a satisfactory answer. One of them is the description of the trajectory of a charged particle in the presence of an electromagnetic field. Yes, we have the Lorentz force, but it does not take into account the effect of the radiation of the charged particle in the acceleration process. The Lorentz-Dirac force is an attempt to include the effect of the radiation but unfortunately, the solutions of this equation are pathological. Most books only mention this situation but Baylis's book in electrodynamics devotes a complete chapter this subject.
Electrodynamics: a modern geometric approach
By William Eric Baylis (chapter 12)

Returning from IMUC

2009-10-25T06:09:00.000-07:00

I am finally writing again after my return from the magnificent International Mathematica User Conference 2009, where I learned about the features of the future Mathematica 8 and met many people. I am going to devote a series of entries about IMUC 2009, but some general observations and things that I learned are

Mathematica is growing and improving at an increasingly faster rate.
The new Mathematica notebook will have nearly the same capacity found in LaTeX for creating static documents. However, is the dynamic capability what makes it revolutionary and much better than LaTeX.
The Mathematica kernel is eventually going to support efficient tensor computation. These routines are going to be implemented in separated modules, so that Tensorial (the tensor package I develop with David Park and Jean-Francois Gouyet) will benefit from these new features. Moreover, now I have more ideas on how to independently improve Tensorial.
The rendering of 3D images with CUDA is truly amazing and I am waiting to play with it when it gets implemented in the kernel. Now I also know that my next laptop has to have an NVidia graphics card.
I am also waiting for the Mathematica plugin for web browsers, so that I can directly publish mathematica notebooks on the web without the need to export any html.

Hyperdeterminants

2009-10-13T18:22:00.000-07:00

The concept of hyperdeterminant is as new to me as a few days and is introducing me into a fascinating new branch of algebra that I did not even suspect. I am very familiar with tensors and how they can be seen as multidimensional arrays that generalize matrices, but the concept of hypermatrices goes further.

The hyperdeterminant was invented (discovered) by the famous mathematician Cayle, who gave us many things including the Cayle transform, useful to approximate the exponential of anti-Hermitian matrices, and the amazing Cayle-Hamilton theorem of linear algebra.

There is active research today and I even found a blog
hyperdeterminant.wordpress.com

More recently I found that Thomas Wolf and Sergey Tsarev are implemening serious computer work in order to find hyperdeterminants of higher order
Hyperdeterminants as integrable discrete systems.

What I had in mind was the hypertederminant of at least 5x5x5 cubic matrices, but now I known this is an exceedingly difficult problem. Now my question is if there is one way to find a good approximation.

MMM,.... I just, learned about other alternative definitions of hyperdeterminants, so the possibilities are still open.

Inverted Retina

2009-10-13T11:50:00.000-07:00

If you did not know, we have our retinas inverted in the sense that the light that enters our eyes has to pass through many layers of nerves, blood vessels and all the wiring before reaching the photo detectors themselves. Even more troublesome, is the fact that it seems that the images are completely distorted and blurred when they arrive the photo-detectors. However we surely know we can see very well. Why did we evolve this feature? and how do we really see? Some of the answers can be found in this absolutely amazing paper

V. D. Svet and A. M. Khazen, About the formation of an image in the inverted retina of the eye, Biophysics, Volume 54, Number 2 / April, 2009, pp 193-203

The key point of this article is that we process the images in blocks. This means that we collect a sequence of images, which we process together in order to improve the signal/noise ratio.

There is a recent independent article here, suggesting that we actually see images in discrete sequences, which is completely consistent with Svet's theory.

I met Dr Svet in Windsor, Canada, where he gave a few seminars about topics concerning sonar imaging and detection. I expect to write more entries in my blog about other topics of his research.

What really matters is to challenge our brains

2009-10-11T13:01:00.000-07:00

Our brains thrive when we learn new challenging things. It is not enough to practice what we already know well
Effect of challenging our brains in our brains

Rodolfo Sanchez Ph.D.

2009-10-10T21:52:00.000-07:00

I was very glad to hear again from my old Bolivian friend Rodolfo Sanchez, who earned his Ph.D. degree in physics in Germany and now is working as a postdoctoral researcher at the GSI institute

Rodolfo Sanchez PhD

This world in small because it happens that he is now collaborating with Dr Gordon Drake, who is professor at Windsor, Canada.

Tree Fractal

2009-10-10T01:21:00.000-07:00

I was a little bit surprised to find one of my Mathematica fractals as part of a collection of other fractals, as you can see in the October 17, 2004 entry

http://nylander.wordpress.co/2004/10/

Calculation of the unitary part of the Bures measure for N-level quantum systems

2009-10-08T11:28:00.000-07:00

My paper about the Bures measure was published

Calculation of the unitary part of the Bures measure for N-level quantum systems

The Bures measure can be written as the product of two factors. One corresponding to the populations and the other corresponding to the unitary transformation. Many people arrived to formulas for the Bures volume and the volume of the unitary part of the Bures measure, but in this paper we give an explicit expression of the measure itself in terms of even balls.

Much more about measures is going to come the public very soon!

Fisher Information In Natural Selection

2009-10-07T22:52:00.000-07:00

I introduced myself with the concept of Fisher information in my readings about quantum estimation and the Cramer-Rao bound. It was when I read the book

B. Roy Frieden, Physics from Fisher Information: A Unification,

which showed me that the Fisher information may even play an important role in the foundations of physics itself. Unexpectedly, I run into the Fisher information again when I was reading about Jeffreys' prior for Bayesian estimation. Now, I came across with another interesting paper about the role of Fisher information in Natural selection

Frank, S. A. 2009. Natural selection maximizes Fisher
information. Journal of Evolutionary Biology 22:231–244

Can we use this in order to refine genetic/evolutionary optimization algorithms?

Quantum estimation

2009-10-07T14:37:00.000-07:00

Qantum estimation/metrology is a research field closely tied with quantum information theory. Everything was born with the Heisenberg uncertainty principle, when we realized that in general measurements can be non-commutative in contrast with the classical world where everything is commutative i.e. the order of the measurements is irrelevant. The consequence of the non-commutativity is that two given observables may be incompatible for simultaneous arbitrarily precise measurements.

The next step was done by Helstrom in the 70's who introduced the Symmetric Logarithmic Derivative (SLD) operator to obtain the Quantum Fisher information. It was later in the 80's when Wootters came with the refined and elegant concept of quantum statistical distinguishability, which was tied with the work of Helstrom by Braunstein and Caves (1994) in their famous paper Statistical distance and the geometry of quantum states. This paper is remarkable because it also established the connection with the completely independent work of Uhlmann about the Bures metric, which was born in the generalization of the Berry phase for mixed states.

More recently, we witnessed many important developments in the quantum single parameter estimation with one unexpected twist. Now we know that quantum mechanics and in particular quantum entanglement have the potential to defeat their classical estimation methods in the quest for higher precision. For example, you can see my older post about NOON states and quantum interferometry.

However, this new knowledge cannot be applied in the real world until we learn how to generate a considerable number of entangled states in a reliable way.

You can also see my wikipedia article about the Bures metric

Learning from being wrong

2009-10-06T23:14:00.000-07:00

Being wrong seems to be something to avoid, but it may be the price to pay if we want to improve or optimize something. This is the general philosophy behind Genetic/Evolutionary algorithms and even statistical sampling.

In Genetic/Evolutionar algorithms, the key ingredient is the introduction of mutation. In most cases a mutation will be destructive, but there will be occasions when the mutation will be beneficial. We improve by discarding the destructive mutations, while keeping the beneficial ones.

It was Thomas Bäck from Leiden University who told us how Toyota promotes mutation in the production system in order to make more efficient cars.

In statistical sampling, the objective is to explore the space of a probability distribution (explore the opportunities). Here we have the Metropolis (and Metropolis-Hastings) algorithm, where one is forced to be wrong sometimes in order to explore the probability distribution.

Ben Schumacher, the author of the quantum noiseless coding theorem has something to tell about this
Ben Schumacher on "Being wrong"

Disclaimer: The only way to profit from being wrong (mutate) is to be ready to correct ourselves as fast as possible.

A scientific approach to science education

2009-10-04T17:12:00.000-07:00

In 2006, when I was in Windsor, I had the opportunity to meet Carl Wieman, the winner of the physics Nobel price in 2001 for his experimental work in the production of the Bose-Einstein condensate. However, something that called my attention even more was his research and discoveries about science education. He showed us with experimental data how our current approach is far for being the most efficient way to teach science. Even worse, sometimes teaching science as we do it today, leads to students with even more misconceptions!

More about his research in education can be found HERE

Quantum/Classical transition

2009-10-03T14:18:00.000-07:00

The nature of the quantum/classical transition is according to my opinion, the most important fundamental question of physics. There are many approaches and many insights that have consequences on the continuing debate on the interpretation of quantum mechanics and its mysterious features such as the collapse of the wave function. Some of the approaches and related topics, not necessarily independent from each other are

The Ehrenfest theorem (My favorite).
The Koopman- von Neumann equation for the classical wave function.
The Wigner function for a formulation in the phase space.
Very realed with the previous approach is the Moyal bracket and Weyl quantization.
Feynman Path integrals and decoherence.
The classical spinor formalism. Electrodynamics: a modern geometric approach By William Eric Baylis
The limit to the Thomas-Fermi model by Elliott H. Lieb
The analogy of classical statistical mechanics with quantum mechanics inspired on the Wigner function by Blokhintsev. Another researcher, A. O. Bolivar pursues this idea further but it seems that he was unaware of the previous work by Blokhintsev.
Related with the previous approach is also the work by Amir Caldeira and Anthony J. Leggett, who proposed a quantum dissipation model.
The very intriguing derivation of the Schrodinger equation from classical mechanics with complexified Brownian motion by Edward Nelson.
The p-mechanics formalism, which exploits the representations of the Heisemberg. This method is also closely related with the Weyl quantization.
Geometric quantization, is another method inspired on the Weyl quantization trying to maintain a coordinate-free procedure based on differential geometry.

The creation of the universe in 6 days

2009-10-02T23:57:00.000-07:00

And God created the universe in 6 days, but in order to give us a nice sky full of beautiful shiny stars, he decided to simulate the light as though as it was coming from thousands, millions and billions of years ago in the past. This was necessary because the speed of light is so fast that in 6 days there would be no chances for us to see anything, not even the light from our closest star, excepting our own sun. We would not be able to see the Milky Way as we see it today, not even in 6 thousand years.

However, God did his job so well that there is no way for anybody to find any hint of crack in the simulation. Unfortunately, this implies that humans will never be able to prove that the universe was created in 6 days, no matter how hard they could try.

Shannons's noisy-chanel coding theorem

2009-10-02T21:31:00.000-07:00

This is an amazing theorem that completely defies intuition and lies at the core of classical information theory. One way to defeat the noise is to encode the original information with added redundancy. The price to pay is a reduction on the efficiency of transmission, also called transmission rate. Intuition says that as the communication error goes to zero, through redundancy coding, the transmission rate should go to zero as well. NOT TRUE!. Shannon actually found a finite bound for the transmission rate, even as the communication error goes to zero.
http://en.wikipedia.org/wiki/Noisy-channel_coding_theorem

Princeton Physics Plasma Laboratory

2009-10-02T14:11:00.000-07:00

The director of the Princeton Physics Plasma Laboratory, Stewart Prager, gave today an overview presentation about the research carried out in the institution. They design and construct tokamaks for nuclear fusion and collaborate with other institutions and projects such as ITER.. After seeing the steady progress in this research field, I am much more optimistic about the feasibility to have one day a commercial energy generator plant based in nuclear fusion. Something that also called me the attention was the importance of Lithium as coating material for the inner surfaces of the plasma chamber.

Quantum Interferometer

2009-09-30T12:13:00.000-07:00

Jonathan Dowling from Louisiana State University gave a presentation about the use of entangled photons for high precision interferometry [1] showing that they can lead to much better measurements than it is possible with classical light. These entangled photon states are of the form

Some problems remain such as the difficulty to produce NOON states with high power, but they are very promising anyway.

So far they devised how to make ingenious conditional measurements in order to generate NOON states.

Another important topic was the emulation of non-linear effects by conditional measurements within linear optics [2].

Refererences

J Dowling arXiv:0904.0163
G. G. Lapaire1, et al, Conditional linear-optical measurement schemes generate effective photon nonlinearities DOI: 10.1103/PhysRevA.68.042314
Jonathan Dowling ppt presentation

LaTeX was possible in this entry with MathTran

Human chromosome 2

2009-09-26T09:11:00.000-07:00

Humans have 23 chromosomes while all the other hominids have 24 chromosomes. The question is, how can this be possible if we are suppose to have a common ancestor. Now we know that human chromosome 2 is the result of the fusion of two chromosomes found in our relative hominids.
YouTube:Chromosome 2

Quantum Feedback Control

2009-09-25T14:18:00.000-07:00

Kurt Jacobs gave a presentation about Quantum Feedback Control. He showed how to introduce the action of weak continuous measurements into a master equation that. This equation resembles the Lindblad equation with an extra stochastic term, where the feedback can be introduced in the Hamiltonian. One application was to create a cat state of a harmonic oscillator.

Adiabatic quantum complexity

2009-09-25T00:11:00.000-07:00

The adiabatic quantum computer was proposed in arXiv:quant-ph/0001106v1 in order to solve certain type of problems.

For low dimensions it was shown that the adiabatic quantum computer is polynomial, for a certain problem with exponential complexity in a classical computer. Peter Young, who gave a presentation at Princeton yesterday, is using Quantum Montecarlo Simuations, as used in statistical mechanics for the partition function, in order to engage this problem for higher dimensions where direct simulations are not feasible.
The Complexity Of The Quantum Adiabatic Algorithm

Intelligence gene

2009-09-23T09:09:00.000-07:00

The following article seems to indicate that there is a gene that is associated with higher intelligence, but those who have it are in disadvantage in high pressure exams
memory gene

Epigenetics

2009-09-22T18:18:00.000-07:00

Epigenetics studies the effect of the environment on the expression of genes. A saw the following NOVA tv program
http://www.pbs.org/wgbh/nova/genes/issa.html
and I was very surprised.

LaTeX in my blog

2009-09-22T15:34:00.000-07:00

I want to write equations in my blog, so I found mimetex, but I could not run it from the cgi Princeton server. I do not know why. I may need to set up complicated permissions or there is a compatibility issue between mimetex and the cgi server.

Here I am testing the cgi provided by mimetex

It works, but I can only use it as a test.

Another method, less elegant though, is using MathTran, which I used to make the following equation

Quantum mechanics in plants??

2009-09-21T21:37:00.000-07:00

Profesor Gregory D. Scholes gave a presentation in Princeton about the possible role of quantum mechanical effects in the absorption of light and energy transfer in photosynthetic processes in plants.

www.nature.com/nature/journal/v431/n7006/full/431256a.html

This seems to be very controversial because they claim high coherence at room temperature. Personally I find it difficult to believe it but I'll wait and see more developments.

PkPd

2009-09-21T20:40:00.000-07:00

PKPD stands for Pharmacokinetic/Pharmacodynamic and it seems to be one of the most fruitful research areas in the boundaries between mathematics, statistics, chemistry and BIOLOGY.

http://en.wikipedia.org/wiki/Pharmacodynamics

Particularly, I am mostly interested in the modeling with differential equations and the use of control theory.

Alpha

2009-09-21T20:21:00.000-07:00

Alpha is what they call knowledge data base

http://www.wolframalpha.com/

The difference from regular search engines is is designed to deliver well formatted data.

For example, using "Dow Jones" it gives a very nice page with graphs, which can be downloaded along with the raw data from a Mathematica notebook.

The Dow Jones Industrial Average is one of the most important indexes in Wall street.

In the same way one can download the information for specific companies such as Toyota or Microsoft.

International Mathematica User Conference 2009

2009-09-21T16:09:00.000-07:00

I am going to give an oral presentation at the International Mathematica User Conference 2009 in Champaign Illinois. October 22-24.

List of presentations

My presentation is going to be about my package developed for performing Magnus expansions. This expansion can be used to find approximate solutions for systems of linear equations and linear operators in general

Magnus expansion in Wikipedia

Entanglement measure

2009-09-20T10:48:00.000-07:00

I developed a very nice method to measure the entanglement of pure states

Renan Cabrera, Herschel Rabitz, The landscape of quantum transitions driven by single-qubit unitary transformations with implications for entanglement, J. Phys. A: Math. Theor. 42 (2009) 275303.

doi: 10.1088/1751-8113/42/27/275303

This method is based in the measurement of the Bures distance between a given state and the closest separable state, which allows one to calculate the generalized Schmidt state. The entanglement can be measured from the coefficients of this Schmidt state.

Bures measure

2009-09-20T10:41:00.000-07:00

My paper about the calculation of the unitary part of the Bures measure was recently accepted for publication

Renan Cabrera, Herschel Rabitz, Calculation of the Unitary part of the Bures Measure for N-level Quantum Systems , Accepted at J. Phy A: Math. Theor.

This paper shows how the Bures measure can be expressed as the product of the measure of even euclidean balls. This means that we now have a simple and easy formula for sampling. The ultimate application will be in Bayesian quantum estimation.

My wikipedia article that explains the basics of the Bures metric is

http://en.wikipedia.org/wiki/Bures_metric

Group theory for wireless communications

2009-09-20T10:24:00.000-07:00

The technology of transmission of information between multiple antennas began to develop 10 years ago. The purpose is to increase the fidelity and transmission rate between multiple sources and multiple receptors.

The particular technology that I am interested mostly is in the Unitary Space-Time Codes (UST). In this technology, the information is encoded in blocks spanning space and time. In the simplest case, we could think of N antennas with a sequence of N pulses, so that the information block can be represented as an NxN complex matrix. It is mathematically convenient to use unitary matrices, so the name UST is justified in this case.

Inside the UST field I am paying attention to the Cayley encoding [1,2,3] because it seems to be elegant, simple and easy to understand.

References

Cayley Differential Unitary Space-Time Codes by B. Hochwald, B. Hassibi
Unitary space-time modulation via Cayley transform by Y Jing, B Hassibi
Yindi Jing Thesis
Representation Theory for High-Rate Multiple-Antenna Code Design by B. Hassibi, B. Hochwald, A. Shokrollahi, W. Sweldens
Differential Unitary Space-Time Modulation by B. Hochwald, W. Sweldens
Random Matrices for Wireless Communications
A. Tulino, S. Verdu (at Princeton)
Circuits for Wireless Communications: Selected Readingsby Banlue Srisuchinwong (Editor), Wanlop Surakampontorn (Editor), Sawasd Tantaratana (Editor)
Space-time coding for broadband wireless communications by By Georgios B. Giannakis, Zhiqiang Liu
Random Matrices For Wireless Communications I
Random Matrices For Wireless Communications II

More information can be found at Bell labs

http://mars.bell-labs.com/

More applications of group theory in engineering can be found at

http://www.usna.edu/Users/math/wdj/repn_thry_appl.htm

Greetings

2009-09-20T06:29:00.000-07:00

Hi, My name is Renan Cabrera L.

I am a physicist interested in quantum mechanics, group theory, and science in general.

My web page is

Renan Cabrera's Web Page