Archive

Archive for February, 2012

New paradigm in science diffusion

February 26, 2012 Leave a comment

Eisenstein October "ten days that shook the world"A recurrent problem among the scientific community is the access to article, this looks to me much more obvious since I arrived within a lab with more physicians and biologist. The situation looks like both a huge loss of time and money for our institutions (where budget tend to decrease) and a net loss for science.

In some research area like mathematics, the problem is even worth because not only access to article is difficult, but the impact factor of mathematical journal remain low (which is annoying in the situation where research is more and more evaluated according to “standard” criteria), as instance the annal of mathematics, according to wikipedia have a Impact factor of 4.1 (rank as the first journal within the mathematics field), which is not very high when you think that it publish works that received the fields medal.  Ways to protest and try to change the business model for science diffusion exist :

Get informed on the problem here, or here (and for those who prefer French link, ici). And eventually support the boycott here

An alternative is to help “open-access” journal, which are more and more renown like PLoS and BMC (to such an extent that they are followed by more conventional editor !), by including them in your must be read journal and or by publishing high quality studies in it…

Use preprint archive systems

Although very popular among the mathematic community, ArXiv.org remain pretty unknown for biologist (although it’s a 20 years old project !). Basically ArXiv is an easy way to publish pre-print (and thus give access ahead of print) for thus who could be interested in your work. A very funny things is that you can find on arXiv works that received prestigious awards before getting published !  A good news, while the biology section of arXiv remain poorly furnished (compared to other domain) Nature publishing group seems to test a similar approach with its Nature Preceding. I must say by now that I am not totally convinced that this will be the exact arXiv for biologist but let wait and see !

So comrade new paradigm are emerging let give them a try and see if world is changing ! 😉

Feeling weird

February 21, 2012 Leave a comment

Sometimes, I just feel weird while I speak to people about statistical/computational/genetics stuff….but you know…as long as you find people that are even weirder…

For the most curious of you, this competition is the MIT integration bee’s, whose purpose is to integrate as fast as possible functions !

Categories: Funny science, Misc

This so small…and so HUGE

February 18, 2012 Leave a comment

I guess this is THE BIG NEWS of the week :

Nanopore announcing their new product

Nanopore DNA sequencing from Oxford Nanopore on Vimeo.

Which of course has already been commented here or here and here, with some extra informations here.

In a nutshell :

  • a small portable disposable device “adaptable” for DNA sequencing
  • a new technology that can reduce several fixed charges in the genotyping process

I am really eager to see how well it will work, to which extent the minION could be use out of the lab, but surely this product will definitely have a huge impact on our day-to-day work !

Using Fortran preprocessor (1)

February 16, 2012 Leave a comment

Library like BLAS or LAPACK are both very efficient and pretty often installed on linux computer.  Anyway, when you have to share a source code, you can’t rely on the fact that those libraries are set up on the host computer. So the safe solution is to find a potentially less efficient alternative, which may sound annoying for people used to BLAS and LAPACK.

Fortunately, for this kind of dilemma there is the PreProcessor, that although not part of the standard, is understood by several fortran compiler (at least gfortran, ifort and xlf).

The basic idea of preprocessing is to allow conditional compilation : depending on certain conditions, one part or another of the code will be compiled.

A first exemple

To illustrate let’s make a  sample program, that should be named “TestPreProcessor.F90”, the extension should be written with capital letters (i.e .F90 and not .f90). That’s the way most compilers automatically know that they should preprocess the code.

program TestPreProcessor
implicit none
#if TEST
 write(*,*)"Compiled with Test Flag"
#ELIF
 write(*,*)"Not Compiled with Test Flag"
#END
end program TestPreProcessor

Now we have two ways  to compile our program.

The classical way

#Compile with no flag
ifort  TestPreProcessor.F90 -o TestPreProcessor
#Test
./TestPreProcessor
#You should see something like this:
Not Compiled with Test Flag

Now the alternative way

#Compile with no flag
ifort  -DTEST TestPreProcessor.F90 -o TestPreProcessor
#Test
./TestPreProcessor
#You should see something like this:
Compiled with Test Flag

The important point is that the flag -D stand for “define”, then during preprocessing, the keyword written after -D will be considered as defined. The keys  #IF…#ELIF…#END will then be examined.”#if TEST” will stand for “is TEST defined ?” if it is true then code below the #if will compiled otherwise code above #elif will be compiled. That ‘s as simple as it seems !

Wrap up and useful example

So we ‘ve explained the basic idea with preprocessing, let apply this to a real problem. Suppose, we want to make a dot product a two vectors, and would like to offer the possibility to use either the dot_product  intrinsic or the equivalent lapack subroutine.

program Dotprod
 implicit none
 integer,parameter::n=1000000
 real(kind=4)::V1(n),V2(n),results
 real,external::sdot

 V1=1.0; V2=1.0
 !//Use BLAS if asked//
#if BLAS
 results=sdot(n,V1,1,V2,1)
#else
 results=dot_product(V1,V2)
#endif
 !//only a test//
 write(*,*)results
end program Dotprod

And the corresponding makefile could look like :

PROG=Dotprod
COMP=ifort
#--------------------------------------------------------------------
#Compilations options + Library path
#--------------------------------------------------------------------
BLAS= -DBLAS -lblas
FF= -O3 $(BLAS) # Comment or uncomment $(BLAS)
#--------------------------------------------------------------------
#Main rules
#--------------------------------------------------------------------
$(PROG): $(PROG).F90
$(COMP) $(FF) -o $(PROG) $(PROG).F90

These should  allow to use either a BLAS subroutine or a standard internal subroutine throughout your code, and ask for very few changes to the end-user who would like to compile your code.

There are several other very useful preprocessor reserved key…but we will see them in an other post 🙂

Categories: Fortran, Linux Tags: ,

Fortran command line arguments

February 9, 2012 2 comments

helloOne little detail that can make your program looks professional is
the ability to pass arguments on the command line.

Surprisingly a lot of Fortran programmers don’t know that this can be done easily, this post will give some tips to use commands arguments.

First, until recently Fortran standard didn’t include any instructions to handle command line argument. Nevertheless, several compilers (in fact all the compilers,  I’ve worked with) like ifort, gfortran, g95 were accepting the use of C function iargc  and the subroutine getarg. And I must confess that until recently this was the tools I used.

Anyways since more and more fortran compilers are compliant with the F03 standard, the official (and recommended) way to deal with command line arguments is the function :
command_argument_count()
which retrieve the number of argument passed to the command line

and the subroutine
get_command_argument(i,aname)
which put character string present at the “i”th argument into variable “aname”.

First step

To illustrate use of command arguments let’s make a  sample program

program TestArg
 integer::narg,cptArg !#of arg & counter of arg
 character(len=20)::name !Arg name

!Check if any arguments are found
 narg=command_argument_count()
!Loop over the arguments
 if(narg>0)then
!loop across options
 do cptArg=1,narg
  call get_command_argument(cptArg,name)
   select case(adjustl(name))
    case("--help","-h")
     write(*,*)"This is program TestArg : Version 0.1"
    case default
     write(*,*)"Option ",adjustl(name),"unknown"
   end select
 end do
 end if
end program TestArg

Compile and test.

ifort TestArg.f90 -o TestArg ; ./TestArg --help

This is for the basic case.  A good hint is trying to use “classical” flags, in fact many software use somehow the same flag for classical purpose “verbosity”, “help” etc… I do believe that the arguments are really helpful because they allow you to offer a variety of options that won’t bother the user by default. In fact, without optional arguments, we tend to ask the user a lot of information that are for common users either obvious or useless.

A convenient use of the command arguments, is to  create a debug flag that will help to correct your program by enabling the printing of more information, or version flag that can tell you which version of the software are you using…these flags can reduce tremendously your debugging time !

Adding some difficulties….

I generally use arg to activate an option, and later in the program I execute some pieces of code to be executed. As instance, if I want to write down a certain file to be read, I add a flag and later in the program execution, I will ask for the name of this file. Anyway, you probably know some software (plink, beagle  as instance) which looks for all the information they need in the command arguments : With plink, you would as instance write :

plink --ped mydata.ped --map autosomal.map

A program such as the previous one, won’t handle such situation…but be can program something similar (although I believe this is not so nice).

program TestComplexArg
 integer::narg,cptArg
 character(len=20)::name,pedFile,mapFile
 logical::lookForPed=.FALSE.
 logical::lookForMap=.FALSE.
 logical::fileExist

!Check if arguments are found
 narg=command_argument_count()

 if(narg>0)then
!loop across options
 do cptArg=1,narg
  call get_command_argument(cptArg,name)
   select case(adjustl(name))
!First known args
    case("--ped")
     lookForPed=.TRUE. !change logical value
    case("--map")
     lookForMap=.TRUE.
    case default
!Treat the second arg of a serie
     if(LookForPed)then
      pedFile=adjustl(name) !assign a value to pedfile
      inquire(file=pedFile,exist=fileExist)!check if it exist
      if(.not.fileExist)then
       write(*,*)'file ',pedFile,' not found'
       stop
      endif
      LookForPed=.FALSE. !put the logical variable to its initial value
     elseif(LookForMap)then
      mapFile=adjustl(name)
      inquire(file=mapFile,exist=fileExist)
      if(.not.fileExist)then
       write(*,*)'file ',mapFile,' not found'
       stop
      endif
      LookForMap=.FALSE.
     else
      write(*,*)"Option ",adjustl(name),"unknown"
     endif
    end select
   end do
  endif
end program TestComplexArg

Compile and test.

#Compile
ifort TestArg.f90 -o TestArg
#Create empty files
touch mydata.ped  autosomal.map
#Test
./TestArg --ped mydata.ped --map autosomal.map

I hope these examples can be useful, and allow you to code “better” (either less verbose or more user-friendly) software.

Categories: Fortran, Linux, Shell Tags: ,