Wrapper Scilab

Prerequisite

• Scilab

• OpenTURNS (this tutorial was done with the 0.13.1 version), located at $OPENTURNS_HOME.

• Make sure you have automake, autoconf, libtool, gcc.

Modification of the source files

• Copy the folder $OPENTURNS_HOME/WrapperTemplate/wrapper_calling_shell_command/ in your home folder for example

• Change the wrapper_calling_shell_command folder to ScilabWrapper

• Go to the ScilabWrapper folder

• Delete code_C1.c and code_C1.data

• Adapt the template to your needs:

  ./customize call_scilab
  

  it will substitute all the occurences of wcode by call_scilab in the files and in the file names.

• Save wrapper.c as ScilabWrapper.c:

  /*                                               -*- C -*- */
  /**
   *  @file  ScilabWrapper.c
   *  @brief The wrapper adapts the interface of OpenTURNS and of the wrapped code
   *
   *  (C) Copyright 2005-2007 EDF-EADS-Phimeca
   *
   *  This library is free software; you can redistribute it and/or
   *  modify it under the terms of the GNU Lesser General Public
   *  License as published by the Free Software Foundation; either
   *  version 2.1 of the License.
   *
   *  This library is distributed in the hope that it will be useful
   *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   *  Lesser General Public License for more details.
   *
   *  You should have received a copy of the GNU Lesser General Public
   *  License along with this library; if not, write to the Free Software
   *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
   *
   *  @author: $LastChangedBy: dutka $
   *  @date:   $LastChangedDate: 2008-08-28 17:36:47 +0200 (Thu, 28 Aug 2008) $
   *  Id:      $Id: wrapper.c 916 2008-08-28 15:36:47Z dutka $
   */
  
  #include "WrapperCommon.h"
  #include "WrapperMacros.h"
  
  
  /* WARNING : Please read the following lines
   *
   * In this program, we make the assumption that the end user wishes to
   * call a function (aka NumericalMathFunction) named "call_scilab" (call_scilab stands
   * for wrapped code).
   * In order to individualize the wrapper to the user's needs, we encourage
   * you (as the developer of this wrapper) to renamed any occurence of "call_scilab"
   * (in either case) to the real name of the function.
   * It will also avoid any confusion with other "call_scilab"s written by other entities
   * or developers.
   */
  
  
  /* If you plan to link this wrapper against a FORTRAN library,
   * remember to change the name WCODE for your actual FORTRAN subroutine name.
   * Otherwise you can ignore these lines.
   *
   * Remember that FORTRAN passes its arguments by reference, not by value as C and C++
   * usually do. So you need to pass the pointer to the arguments rather than its value.
   * This is true for single values (integers, reals, etc.) but not for arrays that are
   * already pointers in the C/C++ environment. Those ones can directly be passed as
   * "values" though they are pointers indeed.
   * Be careful that C and C++ arrays start from 0 and FORTRAN from 1!
   * Be also very careful with the size of the value your plan to pass. Integers in C/C++
   * are not INTEGER*8 in many cases. Float or doubles do not necessarily match REAL*4
   * or REAL*8 in FORTRAN.
   *
   * FORTRAN gives no clue for preventing const values from being altered. So you need to
   * protect them by copying them before calling the FORTRAN subroutine if this import
   * to you.
   *
   * Summary: there are only exceptions to the rule and you need to
   * know exactly what you are doing! You may be desappointed at first, but it will keep
   * you away from segmentation violation and other similar fantasies. ;-)
   */
  
  /* If you want to customize this wrapper to your needs, you have to do the following :
   *  - change the current wrapper name 'call_scilab' to any name you chooze. Be careful to respect
   *    the case everywhere you make the change.
   *  - adapt the signatures of the calls. This means that you may want to call :
   *     a) a C function : you have to write it yourself in this file or in another one and
   *        link the wrapper with the corresponding object file.
   *     b) a FORTRAN function : you have to write the function in another file and link
   *        the wrapper with the corresponding object file. But, due to some technical
   *        aspect of C/FORTRAN linking, you have to declare your FORTRAN function to C in
   *        a special manner and use the F77_FUNC macro below. Don't worry, it's easy !
   *        First, change the name 'call_scilab' everywhere in the macro to the name of your FORTRAN
   *        subroutine. Second, declare the arguments passed to the subroutine in a prototype
   *        declaration. Remember that all FORTRAN argumets are passed as reference (aka pointers)
   *        conversely to C. Long and double in C match INTEGER*4 and REAL*8 in FORTRAN
   *        respectively. Avoid using strings.
   *     c) anything else : you are left alone but the preceding rules may help you.
   *  - call your C function or the FORTRAN macro in the execute function of your wrapper.
   */
  
  /* The name of the wrapper's functions is defined in WRAPPERNAME macro */
  #define WRAPPERNAME call_scilab
  
  
  /*
   *  This is the declaration of function named 'call_scilab' into the wrapper.
   */
    
  
  #ifdef __cplusplus
  extern "C" {
  #endif
  
    /*
  *********************************************************************************
  *                                                                               *
  *                             call_scilab function                                    *
  *                                                                               *
  *********************************************************************************
  */
  
    /* The wrapper information informs the NumericalMathFunction object that loads the wrapper of the
     * signatures of the wrapper functions. In particular, it hold the size of the input
     * NumericalPoint (inSize_) and of the output NumericalPoint (outSize_).
     * Those information are also used by the gradient and hessian functions to set the correct size
     * of the returned matrix and tensor.
     */
    
    /* The getInfo function is optional */
    FUNC_INFO( WRAPPERNAME ,
  	     {
  	       GET_EXCHANGED_DATA_FROM( p_state );
  	       SET_INFORMATION_FROM_EXCHANGED_DATA( p_exchangedData );
  	     } )
      
    /* The state creation/deletion functions allow the wrapper to create or delete a memory location
     * that it will manage itself. It can save in this location any information it needs. The OpenTURNS
     * platform only ensures that the wrapper will receive the state (= the memory location) it works
     * with. If many wrappers are working simultaneously or if the same wrapper is called concurrently,
     * this mechanism will avoid any collision or confusion.
     * The consequence is that NO STATIC DATA should be used in the wrapper OR THE WRAPPER WILL BREAKE
     * one day. You may think that you can't do without static data, but in general this is the footprint
     * of a poor design. But if you persist to use static data, do your work correctly and make use
     * of mutex (for instance) to protect your data against concurrent access. But don't complain about
     * poor computational performance! 
     */
      
      
    /* The createState function is optional */
      FUNC_CREATESTATE( WRAPPERNAME ,
  		      {
  			CHECK_WRAPPER_MODE( WRAPPER_FORK );
  			CHECK_WRAPPER_IN(  WRAPPER_FILES );
  			CHECK_WRAPPER_OUT( WRAPPER_FILES );
  		      
  			COPY_EXCHANGED_DATA_TO( p_p_state );
  
  			PRINT( "My message is here" );
  		      } )
      
    /* The deleteState function is optional */
      FUNC_DELETESTATE( WRAPPERNAME ,
  		      {
  			DELETE_EXCHANGED_DATA_FROM( p_state );
  		      } )
      
  
  
  
  
  
    /* Any function declared into the wrapper may declare three actual function prefixed with
     * 'init_', 'exec_' and 'finalize_' folowed by the name of the function, here 'call_scilab'.
     *
     * The 'init_' function is only called once when the NumericalMathFunction object is created.
     * It allows the wrapper to set some internal state, read some external file, prepare the function
     * to run, etc. It takes only one argument, the internal state as created by the 
     *
     * The 'exec_' function is intended to execute what the wrapper is done for: compute an mathematical
     * function or anything else. It takes the internal state pointer as its first argument, the input
     * NumericalPoint pointer as the second and the output NumericalPoint pointer as the third.
     *
     * The 'finalize_' function is only called once when the NumericalMathFunction object is destroyed.
     * It allows the wrapper to flush anything before unloading.
     *
     * Only the 'exec_' function is mandatory.
     */
      
      
    /**
     * Initialization function
     * This function is called once just before the wrapper first called to initialize
     * it, ie create a temparary subdirectory (remember that the wrapper may be called
     * concurrently), store exchanged data in some internal repository, do some
     * pre-computational operation, etc.
     */
      FUNC_INIT( WRAPPERNAME , {} )
        
        
        
  
    /**
     * Execution function
     * This function is called by the platform to do the real work of the wrapper. It may be
     * called concurrently, so be aware of not using shared or global data not protected by
     * a critical section.
     * This function has a mathematical meaning. It operates on one vector (aka point) and
     * returns another vector.
     */
      FUNC_EXEC( WRAPPERNAME,
  	       FUNC_EXEC_BODY_CALLING_COMMAND_IN_TEMP_DIR( "call_scilab" ) )
    
      FUNC_EXEC_SAMPLE_MULTITHREADED( WRAPPERNAME )
  
    /**
     * Finalization function
     * This function is called once just before the wrapper is unloaded. It is the place to flush
     * any output file or free any allocated memory. When this function returns, the wrapper is supposed
     * to have all its work done, so it is not possible to get anymore information from it after that.
     */
      FUNC_FINALIZE( WRAPPERNAME , {} )
  
  
  #ifdef __cplusplus
      } /* end extern "C" */
  #endif
  

• If you prefer to have a more concise wrapper, make some cleaning. A quite minimal that support multi-threaded evaluations over a sample is:

  #include "WrapperCommon.h"
  #include "WrapperMacros.h"
  
  #define WRAPPERNAME call_scilab
  
  #ifdef __cplusplus
  extern "C" {
  #endif
  
    FUNC_EXEC( WRAPPERNAME,
               FUNC_EXEC_BODY_CALLING_COMMAND_IN_TEMP_DIR( "call_scilab" ) )
  
    FUNC_EXEC_SAMPLE_MULTITHREADED( WRAPPERNAME )
  
  #ifdef __cplusplus
  } /* end extern "C" */
  #endif
  

• Save call_scilab.xml.in as ScilabWrapper.xml.in and adapt it to your needs:

  <?xml version="1.0" encoding="ISO-8859-1"?>
  <!-- 
   @configure_input@
  
   @file  ScilabWrapper.xml.in
   @brief The description file of the wrapper
  
   (C) Copyright 2005-2007 EDF-EADS-Phimeca
  
   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License.
  
   This library is distributed in the hope that it will be useful
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.
  
   You should have received a copy of the GNU Lesser General Public
   License along with this library; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  
   @author: $LastChangedBy: dutka $
   @date:   $LastChangedDate: 2007-10-08 16:24:36 +0200 (lun, 08 oct 2007) $
   Id:      $Id: call_scilab.xml.in 563 2007-10-08 14:24:36Z dutka $
  -->
  <!DOCTYPE wrapper SYSTEM "@OPENTURNS_WRAPPER_DTD_PATH@/wrapper.dtd">
  
  <wrapper>
    <library>
  
      <!-- The path of the shared object -->
      <path>call_scilab.so</path>
  
  
  
      <!-- This section describes all exchanges data between the wrapper and the platform -->
      <description>
  
        <!-- Those variables are substituted in the files above -->
        <!-- The order of variables is the order of the arguments of the function -->
        <variable-list>
  
          <!-- The definition of a variable -->
          <variable id="R" type="in">
            <comment>Yield strength</comment>
            <unit>MPa</unit>
            <regexp>^(\S*)R(\S*)=(\S*)(\R)\S*$</regexp>
            <format>R = %.20e</format>
          </variable>
  
          <!-- The definition of a variable -->
          <variable id="F" type="in">
            <comment>Axial Force</comment>
            <unit>N</unit>
            <regexp>^\S*F\S*=\S*(\R)\S*$</regexp>
            <format>F = %.20e</format>
          </variable>
  
          <!-- The definition of a variable -->
          <variable id="G" type="out">
            <comment>Limit State</comment>
            <unit>MPa</unit>
            <regexp>^\S*G\S*=\S*(\R)\S*$</regexp>
          </variable>
  
        </variable-list>
        
  
  
        <!-- The function that we try to execute through the wrapper -->
        <function provided="yes">call_scilab</function>
  
        <!-- the gradient is  defined  -->
        <gradient provided="no"></gradient>
  
        <!--  the hessian is  defined  -->
        <hessian provided="no"></hessian>
  
      </description>
  
  
    </library>
  
    <external-code>
      <!-- Those data are external to the platform (input files, etc.)-->
      <data>
        <!-- An input file -->
        <file id="data" type="in">
          <name>The input data file</name>
          <path>code_scilab.data</path>
        </file>
  
        <!-- An output file -->
        <file id="result" type="out">
          <name>The output result file</name>
          <path>result.data</path>
  	<subst>G</subst>
        </file>
   
      </data>
  
      <wrap-mode type="fork" state="shared">
        <in-data-transfer mode="files" />
        <out-data-transfer mode="files" />
      </wrap-mode>
  
      <command>scilab -nb -ns -nouserstartup -nw -nwni -nogui -f @prefix@/bin/code_scilab.sce 2&gt;&amp;1 &gt; /dev/null</command>
  
    </external-code>
  
  </wrapper>
  

• Create code_scilab.data and code_scilab.sce:
  code_scilab.data

  R = 300.0
  F = 75000.0
  

  code_scilab.sce implements the same physic example as here:

  // Read the data
  exec("code_scilab.data");
  // Perform the computation
  G = R - F / (100.0 * %pi);
  // Save the result
  f = mopen('result.data', 'wt');
  mfprintf(f, 'G = %.20e', G);
  file('close', f);
  // Exit
  quit;
  

• Save test.py as test_ScilabWrapper.py :

  #                                               -*- Python -*-
  #
  #  @file  test_ScilabWrapper.py
  #  @brief A test file for the wrapper code
  #
  #  (C) Copyright 2005-2007 EDF-EADS-Phimeca
  #
  #  This library is free software; you can redistribute it and/or
  #  modify it under the terms of the GNU Lesser General Public
  #  License as published by the Free Software Foundation; either
  #  version 2.1 of the License.
  #
  #  This library is distributed in the hope that it will be useful
  #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  #  Lesser General Public License for more details.
  #
  #  You should have received a copy of the GNU Lesser General Public
  #  License along with this library; if not, write to the Free Software
  #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  #
  #  @author: $LastChangedBy: dutka $
  #  @date:   $LastChangedDate: 2007-10-05 23:15:55 +0200 (ven, 05 oct 2007) $
  #  Id:      $Id: test_wcode.py 556 2007-10-05 21:15:55Z dutka $
  #
  #  This script imports the OpenTURNS environment in Python, loads
  #  ths wrapper and the wrapped code, and then calls it.
  #
  from openturns import *
  
  f = NumericalMathFunction("ScilabWrapper")
  print "f = ", f
  
  inP = NumericalPoint(f.getInputNumericalPointDimension())
  for i in range(inP.getDimension()):
      inP[i] = 10. * (i + 1)
  print "inP = ", inP
  
  outP = f(inP)
  print "outP = ", outP
  

• Modify configure.ac:

  #                                               -*- Autoconf -*-
  #
  #  configure.ac
  #
  #  (C) Copyright 2005-2007 EDF-EADS-Phimeca
  #
  #  This library is free software; you can redistribute it and/or
  #  modify it under the terms of the GNU Lesser General Public
  #  License as published by the Free Software Foundation; either
  #  version 2.1 of the License.
  #
  #  This library is distributed in the hope that it will be useful
  #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  #  Lesser General Public License for more details.
  #
  #  You should have received a copy of the GNU Lesser General Public
  #  License along with this library; if not, write to the Free Software
  #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  #
  #  @author: $LastChangedBy: dutka $
  #  @date:   $LastChangedDate: 2008-08-28 17:36:47 +0200 (Thu, 28 Aug 2008) $
  #  Id:      $Id: configure.ac 916 2008-08-28 15:36:47Z dutka $
  #
  #  Process this file with autoconf to produce a configure script.
  #  Set minimal version of Autoconf for configure to work
  #
  
  AC_INIT(ScilabWrapper.c)
  AC_CONFIG_AUX_DIR(.)
  AM_INIT_AUTOMAKE(call_scilab, 0.0)
  AC_DISABLE_STATIC
  AC_PROG_LIBTOOL
  AC_CONFIG_MACRO_DIR([m4])
  
  # We prefer to install the wrapper in the user home directory by default
  # since it is a default search path for Open TURNS library
  AC_PREFIX_DEFAULT([$HOME/openturns])
  
  # You need the C and C++ compilers to build the wrapper
  AC_PROG_CC
  AC_PROG_CXX
  
  OT_DEBUG
  OT_CHECK_THREADS
  OT_CHECK_OPENTURNS
  
  AC_OUTPUT([Makefile ScilabWrapper.xml])
  

• Modify Makefile.am :

  #                                               -*- Makefile -*-
  #
  #  Makefile.am
  #
  #  (C) Copyright 2005-2007 EDF-EADS-Phimeca
  #
  #  This library is free software; you can redistribute it and/or
  #  modify it under the terms of the GNU Lesser General Public
  #  License as published by the Free Software Foundation; either
  #  version 2.1 of the License.
  #
  #  This library is distributed in the hope that it will be useful
  #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  #  Lesser General Public License for more details.
  #
  #  You should have received a copy of the GNU Lesser General Public
  #  License along with this library; if not, write to the Free Software
  #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  #
  #  @author: $LastChangedBy: dutka $
  #  @date:   $LastChangedDate: 2007-10-05 23:15:55 +0200 (ven, 05 oct 2007) $
  #  Id:      $Id: Makefile.am 556 2007-10-05 21:15:55Z dutka $
  #
  ACLOCAL_AMFLAGS = -I m4
  
  wrapperdir          = $(prefix)/wrappers
  
  wrapper_LTLIBRARIES = call_scilab.la
  call_scilab_la_SOURCES   = ScilabWrapper.c
  call_scilab_la_CPPFLAGS  = $(OPENTURNS_CPPFLAGS)
  call_scilab_la_LDFLAGS   = -module -no-undefined -version-info 0:0:0
  call_scilab_la_LDFLAGS  += $(OPENTURNS_LDFLAGS)
  call_scilab_la_LIBADD    = $(OPENTURNS_LIBS)
  
  bin_SCRIPTS = code_scilab.sce
  
  XMLWRAPPERFILE = ScilabWrapper.xml
  wrapper_DATA   = $(XMLWRAPPERFILE)
  EXTRA_DIST     = $(XMLWRAPPERFILE).in test.py code_scilab.sce code_scilab.data
  

Building the wrapper

• Run the bootstrap script:

  ./bootstrap
  

• Configure the wrapper:

  ./configure --prefix=/home/%your home folder%/openturns --with-openturns=$OPENTURNS_HOME
  

• Build the wrapper :

  make
  

• Install your wrapper in /home/%your home folder%/openturns/wrappers :

  make install
  

• The wrapper will be seen by OpenTURNS automatically, as it was configured to be placed in a standard folder, where OpenTURNS checks by default the presence of wrappers.

Checking the wrapper

The test_ScilabWrapper.py script should give:

python test_ScilabWrapper.py
Welcome to OpenTURNS version 0.13.1
USR - (func_createState_call_scilab) My message is here
f =  class=NumericalMathFunction name=ScilabWrapper implementation=class=NumericalMathFunctionImplementation name=ScilabWrapper description=[R,F,G] evaluationImplementation=class=ComputedNumericalMathEvaluationImplementation name=ScilabWrapper gradientImplementation=class=CenteredFiniteDifferenceGradient name=Unnamed epsilon=class=NumericalPoint name=Unnamed dimension=2 implementation=class=NumericalPointImplementation name=Unnamed dimension=2 values=[0.00316228,0.00316228] evaluation=class=ComputedNumericalMathEvaluationImplementation name=ScilabWrapper hessianImplementation=class=CenteredFiniteDifferenceHessian name=Unnamed epsilon=class=NumericalPoint name=Unnamed dimension=2 implementation=class=NumericalPointImplementation name=Unnamed dimension=2 values=[0.01,0.01] evaluation=class=ComputedNumericalMathEvaluationImplementation name=ScilabWrapper
inP =  [10,20]
outP =  [9.93634]

Running a stochastic analysis using the wrapper

The analysis is done using the TUI. Note the use of a BlockSize parameter of 2, in order to benefit from a dual-core processor. It improves the overall simulation time by a factor of 40% on the machine used for the test:

from openturns import *
from math import * 

stochasticDimension = 2
RandomGenerator().SetSeed(0)

# Definition of the function to use : Scilab Coupling
LimitState = NumericalMathFunction("ScilabWrapper")
dim = LimitState.getInputNumericalPointDimension()

# Mean
mean = NumericalPoint(dim, 0.0)
mean[0] = 300.0
mean[1] = 75000.0

# Standard deviation
sigma = NumericalPoint(dim, 0.0)
sigma[0] = 30.0
sigma[1] = 5000.0

component = Description(1)
BorneInf = 0.0

# Initialization of the distribution collection:
aCollection = DistributionCollection()

# Initialization of the distribution collection:
aCollection = DistributionCollection()

# Create a first marginal : distribution 1D
marginal = LogNormal(mean[0], sigma[0], BorneInf, LogNormal.MUSIGMA)
marginal.setName("Yield strength")
component[0] = "R"
marginal.setDescription(component)
# Graphical output of the PDF
pdfgraph=marginal.drawPDF()
pdfgraph.draw("/tmp", "pdf_R", 640, 480)
# Fill the first marginal of aCollection
aCollection.add( Distribution(marginal, "Yield strength") )

# Create a second marginal : distribution 1D
marginal = Normal(mean[1], sigma[1])
marginal.setName("Traction_load")
component[0] = "F"
marginal.setDescription(component)
# Graphical output of the PDF
pdfgraph=marginal.drawPDF()
pdfgraph.draw("/tmp", "pdf_F", 640, 480)
# Fill the second marginal of aCollection
aCollection.add( Distribution(marginal, "Traction_load") )

# Create a copula : IndependentCopula (no correlation)
aCopula = IndependentCopula(aCollection.getSize())
aCopula.setName("Independent copula")

# Instanciate one distribution object
myDistribution = ComposedDistribution(aCollection, Copula(aCopula))
myDistribution.setName("myDist")

# We create a 'usual' RandomVector from the Distribution
vect = RandomVector(Distribution(myDistribution))

# We create a composite random vector
G = RandomVector(LimitState, vect)

# We create an Event from this RandomVector
myEvent = Event(G, ComparisonOperator(Less()), 0.0, "Event 1")

cv = 0.10
NbSim = 50

myAlgo = MonteCarlo(myEvent)
myAlgo.setMaximumOuterSampling(NbSim)
myAlgo.setBlockSize(2)
myAlgo.setMaximumCoefficientOfVariation(cv)

Log.Show(Log.INFO)
myAlgo.run()

result = myAlgo.getResult()
probability = result.getProbabilityEstimate()
print "MonteCarlo result=" , result
print "Number of executed iterations =" , result.getOuterSampling()
print "Pf = " , probability
print "CV =" , sqrt(result.getVarianceEstimate())