Changes in / [ff431ca:d5014e4] in sasmodels


Ignore:
Files:
3 deleted
7 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/orientation/orientation.rst

    r82592da r5fb0634  
    44================== 
    55 
    6 With two dimensional small angle diffraction data SasView will calculate 
     6With two dimensional small angle diffraction data sasmodels will calculate 
    77scattering from oriented particles, applicable for example to shear flow 
    88or orientation in a magnetic field. 
    99 
    1010In general we first need to define the reference orientation 
    11 of the particles with respect to the incoming neutron or X-ray beam. This 
    12 is done using three angles: $\theta$ and $\phi$ define the orientation of 
    13 the axis of the particle, angle $\Psi$ is defined as the orientation of 
    14 the major axis of the particle cross section with respect to its starting 
    15 position along the beam direction. The figures below are for an elliptical 
    16 cross section cylinder, but may be applied analogously to other shapes of 
    17 particle. 
     11of the particle's $a$-$b$-$c$ axes with respect to the incoming 
     12neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$ 
     13define the orientation of the $c$-axis of the particle, and angle $\Psi$ is 
     14defined as the orientation of the major axis of the particle cross section 
     15with respect to its starting position along the beam direction (or 
     16equivalently, as rotation about the $c$ axis). There is an unavoidable 
     17ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both 
     18serve to rotate the particle about $c$, but this symmetry is destroyed 
     19when $\theta$ is not a multiple of 180. 
     20 
     21The figures below are for an elliptical cross section cylinder, but may 
     22be applied analogously to other shapes of particle. 
    1823 
    1924.. note:: 
    … …  
    2934 
    3035    Definition of angles for oriented elliptical cylinder, where axis_ratio 
    31     b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$ 
     36    b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$ 
    3237    plane, is carried out first, then rotation $\phi$ about the $z$-axis, 
    3338    finally rotation $\Psi$ is around the axis of the cylinder. The neutron 
    34     or X-ray beam is along the $z$ axis. 
     39    or X-ray beam is along the $-z$ axis. 
    3540 
    3641.. figure:: 
    … …  
    4045    with $\Psi$ = 0. 
    4146 
    42 Having established the mean direction of the particle we can then apply 
    43 angular orientation distributions. This is done by a numerical integration 
    44 over a range of angles in a similar way to particle size dispersity. 
    45 In the current version of sasview the orientational dispersity is defined 
    46 with respect to the axes of the particle. 
     47Having established the mean direction of the particle (the view) we can then 
     48apply angular orientation distributions (jitter). This is done by a numerical 
     49integration over a range of angles in a similar way to particle size 
     50dispersity. The orientation dispersity is defined with respect to the 
     51$a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis, 
     52yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 
     53$a$-axis. 
     54 
     55More formally, starting with axes $a$-$b$-$c$ of the particle aligned 
     56with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity 
     57is applied first, using the 
     58`Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_ 
     59$x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$. 
     60The reference orientation then follows, using the 
     61`Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_ 
     62$z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$.  This is implemented 
     63using rotation matrices as 
     64 
     65.. math:: 
     66 
     67    R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\, 
     68        R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi) 
     69 
     70To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the 
     71shape in its canonical orientation, use 
     72 
     73.. math:: 
     74 
     75    [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T 
     76 
     77 
     78The inverse rotation is easily calculated by rotating the opposite directions 
     79in the reverse order, so 
     80 
     81.. math:: 
     82 
     83    R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\, 
     84             R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi) 
     85 
    4786 
    4887The $\theta$ and $\phi$ orientation parameters for the cylinder only appear 
    49 when fitting 2d data. On introducing "Orientational Distribution" in 
    50 the angles, "distribution of theta" and "distribution of phi" parameters will 
     88when fitting 2d data. On introducing "Orientational Distribution" in the 
     89angles, "distribution of theta" and "distribution of phi" parameters will 
    5190appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ 
    52 of the cylinder, the $b$ and $a$ axes of the cylinder cross section. (When 
    53 $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the 
    54 instrument.) The third orientation distribution, in $\Psi$, is about the $c$ 
    55 axis of the particle. Some experimentation may be required to understand the 
    56 2d patterns fully. A number of different shapes of distribution are 
    57 available, as described for polydispersity, see :ref:`polydispersityhelp` . 
     91of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder 
     92cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and 
     93$X$ axes of the instrument.) The third orientation distribution, in $\Psi$, 
     94is about the $c$ axis of the particle. Some experimentation may be required 
     95to understand the 2d patterns fully. A number of different shapes of 
     96distribution are available, as described for size dispersity, see 
     97:ref:`polydispersityhelp`. 
    5898 
    59 Earlier versions of SasView had numerical integration issues in some 
    60 circumstances when distributions passed through 90 degrees. The distributions 
    61 in particle coordinates are more robust, but should still be approached with 
    62 care for large ranges of angle. 
     99Given that the angular dispersion distribution is defined in cartesian space, 
     100over a cube defined by 
     101 
     102.. math:: 
     103 
     104    [-\Delta \theta, \Delta \theta] \times 
     105    [-\Delta \phi, \Delta \phi] \times 
     106    [-\Delta \Psi, \Delta \Psi] 
     107 
     108but the orientation is defined over a sphere, we are left with a 
     109`map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_ 
     110problem, with different tradeoffs depending on how values in $\Delta\theta$ 
     111and $\Delta\phi$ are translated into latitude/longitude on the sphere. 
     112 
     113Sasmodels is using the 
     114`equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_. 
     115In this projection, square patches in angular dispersity become wedge-shaped 
     116patches on the sphere. To correct for the changing point density, there is a 
     117scale factor of $\sin(\Delta\theta)$ that applies to each point in the 
     118integral. This is not enough, though. Consider a shape which is tumbling 
     119freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At 
     120$\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have 
     121a distinct angular preference. If the spin axis is along the beam (which 
     122will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern 
     123should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the 
     124$\sin(\Delta\theta)$ correction. This problem does not appear for a shape 
     125that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in 
     126$[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$ 
     127and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to 
     128symmetry. 
     129 
     130Alternative projections were considered. 
     131The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_ 
     132works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those 
     133points outside $[-180, 180]$. The distortions are a little less for middle 
     134ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$ 
     135and the model is much harder to explain. 
     136The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_ 
     137also improves on the equirectangular projection by extending the range of 
     138reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a 
     139wedge that cuts to the opposite side of the sphere rather than cutting to the 
     140pole. This projection has the nice property that distance from the center are 
     141preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same. 
     142The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ 
     143is like the azimuthal equidistance projection, but it preserves area instead 
     144of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$. 
     145The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_ 
     146has an excellent balance with reasonable distortion in both $\Delta\theta$ 
     147and $\Delta\phi$, as well as preserving small patches. However, it requires 
     148considerably more computational overhead, and we have not yet derived the 
     149formula for the distortion correction, measuring the degree of stretch at 
     150the point $(\Delta\theta, \Delta\phi)$ on the map. 
    63151 
    64152.. note:: 
    65     Note that the form factors for oriented particles are also performing 
    66     numerical integrations over one or more variables, so care should be taken, 
    67     especially with very large particles or more extreme aspect ratios. In such  
    68     cases results may not be accurate, particularly at very high Q, unless the model 
    69     has been specifically coded to use limiting forms of the scattering equations. 
    70      
    71     For best numerical results keep the $\theta$ distribution narrower than the $\phi$  
    72     distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may  
    73     need to reorder the sizes of the three axes to acheive the desired result.  
    74     This is due to the issues of mapping a rectangular distribution onto the  
    75     surface of a sphere. 
     153    Note that the form factors for oriented particles are performing 
     154    numerical integrations over one or more variables, so care should be 
     155    taken, especially with very large particles or more extreme aspect 
     156    ratios. In such cases results may not be accurate, particularly at very 
     157    high Q, unless the model has been specifically coded to use limiting 
     158    forms of the scattering equations. 
    76159 
    77 Users can experiment with the values of *Npts* and *Nsigs*, the number of steps  
    78 used in the integration and the range spanned in number of standard deviations.  
    79 The standard deviation is entered in units of degrees. For a "rectangular"  
    80 distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.  
    81 The new "uniform" distribution avoids this by letting you directly specify the  
     160    For best numerical results keep the $\theta$ distribution narrower than 
     161    the $\phi$ distribution. Thus for asymmetric particles, such as 
     162    elliptical_cylinder, you may need to reorder the sizes of the three axes 
     163    to acheive the desired result. This is due to the issues of mapping a 
     164    rectanglar distribution onto the surface of a sphere. 
     165 
     166Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
     167used in the integration and the range spanned in number of standard deviations. 
     168The standard deviation is entered in units of degrees. For a "rectangular" 
     169distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
     170The new "uniform" distribution avoids this by letting you directly specify the 
    82171half width. 
    83172 
    84 The angular distributions will be truncated outside of the range -180 to +180  
    85 degrees, so beware of using saying a broad Gaussian distribution with large value 
    86 of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would  
    87 give a good integration,as well as becoming rather meaningless. (At some point  
    88 in the future the actual polydispersity arrays may be made available to the user  
    89 for inspection.) 
     173The angular distributions may be truncated outside of the range -180 to +180 
     174degrees, so beware of using saying a broad Gaussian distribution with large 
     175value of *Nsigs*, as the array of *Npts* may be truncated to many fewer 
     176points than would give a good integration,as well as becoming rather 
     177meaningless. (At some point in the future the actual dispersion arrays may be 
     178made available to the user for inspection.) 
    90179 
    91180Some more detailed technical notes are provided in the developer section of 
    92181this manual :ref:`orientation_developer` . 
    93182 
     183This definition of orientation is new to SasView 4.2.  In earlier versions, 
     184the orientation distribution appeared as a distribution of view angles. 
     185This led to strange effects when $c$ was aligned with $z$, where changes 
     186to the $\phi$ angle served only to rotate the shape about $c$, rather than 
     187having a consistent interpretation as the pitch of the shape relative to 
     188the flow field defining the reference orientation.  Prior to SasView 4.1, 
     189the reference orientation was defined using a Tait-Bryan convention, making 
     190it difficult to control.  Now, rotation in $\theta$ modifies the spacings 
     191in the refraction pattern, and rotation in $\phi$ rotates it in the detector 
     192plane. 
     193 
     194 
    94195*Document History* 
    95196 
    96197| 2017-11-06 Richard Heenan 
     198| 2017-12-20 Paul Kienzle 

  • doc/guide/plugin.rst

    rc654160 r7e6bc45e  
    538538If the scattering is dependent on the orientation of the shape, then you 
    539539will need to include *orientation* parameters *theta*, *phi* and *psi* 
    540 at the end of the parameter table.  Shape orientation uses *a*, *b* and *c* 
    541 axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate 
    542 system, with *z* along the beam and *x*-*y* in the detector plane, with *x* 
    543 horizontal and *y* vertical.  The *psi* parameter rotates the shape 
    544 about its *c* axis, the *theta* parameter then rotates the *c* axis toward 
    545 the *x* axis of the detector, then *phi* rotates the shape in the detector 
    546 plane.  (Prior to these rotations, orientation dispersity will be applied 
    547 as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate 
    548 system.)  A particular *qx*, *qy* point on the detector, then corresponds 
    549 to *qa*, *qb*, *qc* with respect to the shape. 
    550  
    551 The oriented C model is called as *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     540at the end of the parameter table.  As described in the section 
     541:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will 
     542be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its 
     543canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the 
     544laboratory frame and beam travelling along $-z$. 
     545 
     546The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
    552547*par1*, etc. are the parameters to the model.  If the shape is rotationally 
    553548symmetric about *c* then *psi* is not needed, and the model is called 

  • explore/jitter.py

    rff10479 r8cfb486  
    165165    # constants in kernel_iq.c 
    166166    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 
     167    'azimuthal_equal_area', 
    167168] 
    168169def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 

  • sasmodels/autoc.py

    r67cc0ff r15be191  
    44from __future__ import print_function 
    55 
    6 import ast 
    76import inspect 
    8 from functools import reduce 
    97 
    108import numpy as np 
    … …  
    9593                constants[name] = obj 
    9694                # Claim all constants are declared on line 1 
    97                 snippets.append('#line 1 "%s"'%escaped_filename) 
    98                 snippets.append(define_constant(name, obj)) 
     95                snippets.append('#line 1 "%s"\n'%escaped_filename) 
     96                snippets.append(py2c.define_constant(name, obj)) 
    9997            elif isinstance(obj, special.Gauss): 
    10098                for var, value in zip(("N", "Z", "W"), (obj.n, obj.z, obj.w)): 
    10199                    var = "GAUSS_"+var 
    102100                    constants[var] = value 
    103                     snippets.append('#line 1 "%s"'%escaped_filename) 
    104                     snippets.append(define_constant(var, value)) 
     101                    snippets.append('#line 1 "%s"\n'%escaped_filename) 
     102                    snippets.append(py2c.define_constant(var, value)) 
    105103                #libs.append('lib/gauss%d.c'%obj.n) 
    106104                source = (source.replace(name+'.n', 'GAUSS_N') 
    … …  
    121119 
    122120    # translate source 
    123     ordered_code = [code[name] for name in ordered_dag(depends) if name in code] 
     121    ordered_code = [code[name] for name in py2c.ordered_dag(depends) if name in code] 
    124122    functions = py2c.translate(ordered_code, constants) 
    125123    snippets.extend(functions) 
    … …  
    127125    # update model info 
    128126    info.source = unique_libs 
    129     info.c_code = "\n".join(snippets) 
     127    info.c_code = "".join(snippets) 
    130128    info.Iq = info.Iqac = info.Iqabc = info.Iqxy = info.form_volume = None 
    131  
    132 def define_constant(name, value): 
    133     if isinstance(value, int): 
    134         parts = ["int ", name, " = ", "%d"%value, ";"] 
    135     elif isinstance(value, float): 
    136         parts = ["double ", name, " = ", "%.15g"%value, ";"] 
    137     else: 
    138         # extend constant arrays to a multiple of 4; not sure if this 
    139         # is necessary, but some OpenCL targets broke if the number 
    140         # of parameters in the parameter table was not a multiple of 4, 
    141         # so do it for all constant arrays to be safe. 
    142         if len(value)%4 != 0: 
    143             value = list(value) + [0.]*(4 - len(value)%4) 
    144         elements = ["%.15g"%v for v in value] 
    145         parts = ["double ", name, "[]", " = ", 
    146                  "{\n   ", ", ".join(elements), "\n};"] 
    147     return "".join(parts) 
    148  
    149  
    150 # Modified from the following: 
    151 # 
    152 #    http://code.activestate.com/recipes/578272-topological-sort/ 
    153 #    Copyright (C) 2012 Sam Denton 
    154 #    License: MIT 
    155 def ordered_dag(dag): 
    156     # type: (Dict[T, Set[T]]) -> Iterator[T] 
    157     dag = dag.copy() 
    158  
    159     # make leaves depend on the empty set 
    160     leaves = reduce(set.union, dag.values()) - set(dag.keys()) 
    161     dag.update({node: set() for node in leaves}) 
    162     while True: 
    163         leaves = set(node for node, links in dag.items() if not links) 
    164         if not leaves: 
    165             break 
    166         for node in leaves: 
    167             yield node 
    168         dag = {node: (links-leaves) 
    169                for node, links in dag.items() if node not in leaves} 
    170     if dag: 
    171         raise ValueError("Cyclic dependes exists amongst these items:\n%s" 
    172                             % ", ".join(str(node) for node in dag.keys())) 

  • sasmodels/modelinfo.py

    r67cc0ff r5ab99b7  
    794794    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    795795    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    796     info.c_code = getattr(kernel_module, 'c_code', None) 
    797796    info.source = getattr(kernel_module, 'source', []) 
    798797    info.c_code = getattr(kernel_module, 'c_code', None) 
    … …  
    810809    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
    811810    # Default single and opencl to True for C models.  Python models have callable Iq. 
    812     info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
    813     info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    814811    info.random = getattr(kernel_module, 'random', None) 
    815812 
    … …  
    827824        except Exception as exc: 
    828825            logger.warn(str(exc) + " while converting %s from C to python"%name) 
     826 
     827    # Needs to come after autoc.convert since the Iq symbol may have been 
     828    # converted from python to C 
     829    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
     830    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    829831 
    830832    if callable(info.Iq) and parameters.has_2d: 

  • sasmodels/models/_cylpy.py

    r67cc0ff rc01ed3e  
    140140py2c = True 
    141141 
     142# TODO: "#define INVALID (expr)" is not supported 
    142143def invalid(v): 
    143144    return v.radius < 0 or v.length < 0 
    … …  
    206207            phi_pd=10, phi_pd_n=5) 
    207208 
    208 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     209qx, qy = 0.2 * cos(2.5), 0.2 * sin(2.5) 
    209210# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
    210211tests = [ 

  • sasmodels/py2c.py

    r0bd0877 r0bd0877  
    1 """ 
    2     py2c 
    3     ~~~~ 
    4  
    5     Convert simple numeric python code into C code. 
    6  
    7     The translate() function works on 
    8  
    9     Variables definition in C 
    10     ------------------------- 
    11     Defining variables within the translate function is a bit of a guess work, 
    12     using following rules: 
    13     *   By default, a variable is a 'double'. 
    14     *   Variable in a for loop is an int. 
    15     *   Variable that is references with brackets is an array of doubles. The 
    16         variable within the brackets is integer. For example, in the 
    17         reference 'var1[var2]', var1 is a double array, and var2 is an integer. 
    18     *   Assignment to an argument makes that argument an array, and the index 
    19         in that assignment is 0. 
    20         For example, the following python code:: 
    21             def func(arg1, arg2): 
    22                 arg2 = 17. 
    23         is translated to the following C code:: 
    24             double func(double arg1) 
    25             { 
    26                 arg2[0] = 17.0; 
    27             } 
    28         For example, the following python code is translated to the 
    29         following C code:: 
    30  
    31             def func(arg1, arg2):          double func(double arg1) { 
    32                 arg2 = 17.                      arg2[0] = 17.0; 
    33                                             } 
    34     *   All functions are defined as double, even if there is no 
    35         return statement. 
    36  
    37 Based on codegen.py: 
    38  
     1r""" 
     2py2c 
     3~~~~ 
     4 
     5Convert simple numeric python code into C code. 
     6 
     7This code is intended to translate direct algorithms for scientific code 
     8(mostly if statements and for loops operating on double precision values) 
     9into C code. Unlike projects like numba, cython, pypy and nuitka, the 
     10:func:`translate` function returns the corresponding C which can then be 
     11compiled with tinycc or sent to the GPU using CUDA or OpenCL. 
     12 
     13There is special handling certain constructs, such as *for i in range* and 
     14small integer powers. 
     15 
     16**TODO: make a nice list of supported constructs*** 
     17 
     18Imports are not supported, but they are at least ignored so that properly 
     19constructed code can be run via python or translated to C without change. 
     20 
     21Most other python constructs are **not** supported: 
     22* classes 
     23* builtin types (dict, set, list) 
     24* exceptions 
     25* with context 
     26* del 
     27* yield 
     28* async 
     29* list slicing 
     30* multiple return values 
     31* "is/is not", "in/not in" conditionals 
     32 
     33There is limited support for list and list comprehensions, so long as they 
     34can be represented by a fixed array whose size is known at compile time, and 
     35they are small enough to be stored on the stack. 
     36 
     37Variables definition in C 
     38------------------------- 
     39Defining variables within the translate function is a bit of a guess work, 
     40using following rules: 
     41*   By default, a variable is a 'double'. 
     42*   Variable in a for loop is an int. 
     43*   Variable that is references with brackets is an array of doubles. The 
     44    variable within the brackets is integer. For example, in the 
     45    reference 'var1[var2]', var1 is a double array, and var2 is an integer. 
     46*   Assignment to an argument makes that argument an array, and the index 
     47    in that assignment is 0. 
     48    For example, the following python code:: 
     49        def func(arg1, arg2): 
     50            arg2 = 17. 
     51    is translated to the following C code:: 
     52        double func(double arg1) 
     53        { 
     54            arg2[0] = 17.0; 
     55        } 
     56    For example, the following python code is translated to the 
     57    following C code:: 
     58 
     59        def func(arg1, arg2):          double func(double arg1) { 
     60            arg2 = 17.                      arg2[0] = 17.0; 
     61                                        } 
     62*   All functions are defined as double, even if there is no 
     63    return statement. 
     64 
     65Debugging 
     66--------- 
     67 
     68*print* is partially supported using a simple regular expression. This 
     69requires a stylized form. Be sure to use print as a function instead of 
     70the print statement. If you are including substition variables, use the 
     71% string substitution style. Include parentheses around the substitution 
     72tuple, even if there is only one item; do not include the final comma even 
     73if it is a single item (yes, it won't be a tuple, but it makes the regexp 
     74much simpler). Keep the item on a single line. Here are three forms that work:: 
     75 
     76    print("x") => printf("x\n"); 
     77    print("x %g"%(a)) => printf("x %g\n", a); 
     78    print("x %g %g %g"%(a, b, c)) => printf("x %g %g %g\n", a, b, c); 
     79 
     80You can generate *main* using the *if __name__ == "__main__":* construct. 
     81This does a simple substitution with "def main():" before translation and 
     82a substitution with "int main(int argc, double *argv[])" after translation. 
     83The result is that the content of the *if* block becomes the content of *main*. 
     84Along with the print statement, you can run and test a translation standalone 
     85using:: 
     86 
     87    python py2c.py source.py 
     88    cc source.c 
     89    ./a.out 
     90 
     91Known issues 
     92------------ 
     93The following constructs may cause problems: 
     94 
     95* implicit arrays: possible namespace collision for variable "vec#" 
     96* swap fails: "x,y = y,x" will set x==y 
     97* top-level statements: code outside a function body causes errors 
     98* line number skew: each statement should be tagged with its own #line 
     99  to avoid skew as comments are skipped and loop bodies are wrapped with 
     100  braces, etc. 
     101 
     102References 
     103---------- 
     104 
     105Based on a variant of codegen.py: 
     106 
     107    https://github.com/andreif/codegen 
    39108    :copyright: Copyright 2008 by Armin Ronacher. 
    40109    :license: BSD. 
    41110""" 
    42 """ 
    43 Update Notes 
    44 ============ 
    45 11/22/2017, O.E.  Each 'visit_*' method is to build a C statement string. It 
    46                     shold insert 4 blanks per indentation level. 
    47                     The 'body' method will combine all the strings, by adding 
    48                     the 'current_statement' to the c_proc string list 
    49    11/2017, OE: variables, argument definition implemented. 
    50    Note: An argument is considered an array if it is the target of an 
    51         assignment. In that case it is translated to <var>[0] 
    52 11/27/2017, OE: 'pow' basicly working 
    53   /12/2017, OE: Multiple assignment: a1,a2,...,an=b1,b2,...bn implemented 
    54   /12/2017, OE: Power function, including special cases of 
    55                 square(x)(pow(x,2)) and cube(x)(pow(x,3)), implemented in 
    56                 translate_power, called from visit_BinOp 
    57 12/07/2017, OE: Translation of integer division, '\\' in python, implemented 
    58                 in translate_integer_divide, called from visit_BinOp 
    59 12/07/2017, OE: C variable definition handled in 'define_c_vars' 
    60               : Python integer division, '//', translated to C in 
    61                 'translate_integer_divide' 
    62 12/15/2017, OE: Precedence maintained by writing opening and closing 
    63                 parenthesesm '(',')', in procedure 'visit_BinOp'. 
    64 12/18/2017, OE: Added call to 'add_current_line()' at the beginning 
    65                 of visit_Return 
    66 01/4/2018, O.E. Added functions 'get_files_names()', , 'print_usage()'. The get_files_names() 
    67                 function retrieves the in/out file names form the command line, and returns 
    68                 true/false if the number of parameters is valid. In case of no input 
    69                 parameters, usage is prompt and the program terminates. 
    70 01/4/2018, O.E. 'translate(functions, constants=None)' returns string, instaed of a list 
    71 01/04/2017 O.E. Fixed bug in 'visit_If': visiting node.orelse in case else exists. 
    72 01/05/2018 O.E. New class C_Vector defined. It holds name, type, size and initial values. 
    73                 Assignment to that class in 'sequence_visit'. 
    74                 Reading and inserting declaration in 'insert_c_vars' 
    75  
    76  
    77 """ 
     111 
     112# Update Notes 
     113# ============ 
     114# 2017-11-22, OE: Each 'visit_*' method is to build a C statement string. It 
     115#                 shold insert 4 blanks per indentation level. The 'body' 
     116#                 method will combine all the strings, by adding the 
     117#                'current_statement' to the c_proc string list 
     118# 2017-11-22, OE: variables, argument definition implemented.  Note: An 
     119#                argument is considered an array if it is the target of an 
     120#                 assignment. In that case it is translated to <var>[0] 
     121# 2017-11-27, OE: 'pow' basicly working 
     122# 2017-12-07, OE: Multiple assignment: a1,a2,...,an=b1,b2,...bn implemented 
     123# 2017-12-07, OE: Power function, including special cases of 
     124#                 square(x)(pow(x,2)) and cube(x)(pow(x,3)), implemented in 
     125#                 translate_power, called from visit_BinOp 
     126# 2017-12-07, OE: Translation of integer division, '\\' in python, implemented 
     127#                 in translate_integer_divide, called from visit_BinOp 
     128# 2017-12-07, OE: C variable definition handled in 'define_c_vars' 
     129#               : Python integer division, '//', translated to C in 
     130#                 'translate_integer_divide' 
     131# 2017-12-15, OE: Precedence maintained by writing opening and closing 
     132#                 parenthesesm '(',')', in procedure 'visit_BinOp'. 
     133# 2017-12-18, OE: Added call to 'add_current_line()' at the beginning 
     134#                 of visit_Return 
     135# 2018-01-03, PK: Update interface for use in sasmodels 
     136# 2018-01-03, PK: support "expr if cond else expr" syntax 
     137# 2018-01-03, PK: x//y => (int)((x)/(y)) and x/y => ((double)(x)/(double)(y)) 
     138# 2018-01-03, PK: True/False => true/false 
     139# 2018-01-03, PK: f(x) was introducing an extra semicolon 
     140# 2018-01-03, PK: simplistic print function, for debugging 
     141# 2018-01-03, PK: while expr: ... => while (expr) { ... } 
     142# 2018-01-04, OE: Fixed bug in 'visit_If': visiting node.orelse in case else exists. 
     143 
     144from __future__ import print_function 
     145 
     146import sys 
    78147import ast 
    79 import sys 
    80148from ast import NodeVisitor 
    81149 
    … …  
    117185 
    118186 
    119 def to_source(tree, constants=None, fname=None, lineno=0): 
    120     """ 
    121     This function can convert a syntax tree into C sourcecode. 
    122     """ 
    123     generator = SourceGenerator(constants=constants, fname=fname, lineno=lineno) 
    124     generator.visit(tree) 
    125 #    c_code = "\n".join(generator.c_proc) 
    126     c_code = "".join(generator.c_proc) 
    127     return (c_code, generator.warnings) 
    128  
     187# TODO: should not allow eval of arbitrary python 
    129188def isevaluable(s): 
    130189    try: 
    … …  
    133192    except Exception: 
    134193        return False 
    135  
    136 class C_Vector(): 
    137     def __init__(self): 
    138         self.size = 0 
    139         self.type = "double" 
    140         self.values = [] 
    141         self.assign_line = -1 
    142         self.name = "" 
    143  
    144     def get_leading_space(self, indent_with="    ", indent_level = 1): 
    145         leading_space = "" 
    146         for n in range(indent_level): 
    147             leading_space += indent_with 
    148         return leading_space 
    149          
    150     def item_declare_string(self, indent_with="    ", indent_level = 1): 
    151         declare_string = self.name + "[" + str(len(self.values)) + "]" 
    152         return declare_string 
    153      
    154     def declare_string(self, indent_with="    ", indent_level = 1): 
    155         declare_string = self.get_leading_space(indent_with, indent_level) 
    156         declare_string += self.type + " " + item_declare_string(indent_with, indent_level) 
    157 #        self.name + "[" + str(len(self.values)) + "];\n" 
    158         return declare_string 
    159      
    160     def get_assign(self, indent_with="    ", indent_level = 1): 
    161         assign_loop = [] 
    162         c_string = self.get_leading_space(indent_with, indent_level - 1) 
    163         c_string += "for (n=0 ; n < " + str(len(self.values)) + " ; n++) {" 
    164         assign_loop.append(c_string) 
    165         for n,value in enumerate(self.values): 
    166             c_string = self.get_leading_space(indent_with, indent_level) 
    167             c_string += self.name + "[" + str(n) + "] = " + str(value) + ";" 
    168             assign_loop.append(c_string) 
    169         c_string = self.get_leading_space(indent_with, indent_level - 1) + "}" 
    170         assign_loop.append(c_string) 
    171         return (assign_loop) 
    172194 
    173195class SourceGenerator(NodeVisitor): 
    … …  
    203225        self.c_vectors = [] 
    204226        self.c_constants = constants if constants is not None else {} 
     227        self.in_expr = False 
    205228        self.in_subref = False 
    206229        self.in_subscript = False 
    … …  
    209232        self.is_sequence = False 
    210233        self.visited_args = False 
    211         self.inside_if = False 
    212         self.C_Vectors = [] 
    213         self.assign_target = "" 
    214         self.assign_c_vector = [] 
    215234 
    216235    def write_python(self, x): 
    … …  
    248267            self.add_unique_var(arg.id) 
    249268 
    250     def write_print(self, node): 
    251         self.write_c ("printf (") 
    252         c_str = self.current_statement 
    253         self.current_statement = "" 
    254         self.visit(node.args[0]) 
    255         print_params = self.current_statement # save print agruments 
    256         self.current_statement = c_str 
    257         print_params = print_params[1:(len(print_params) - 1)] # remove parentheses 
    258         prcnt = print_params.rfind('%') 
    259         if prcnt >= 0: 
    260             format_string = print_params[0:prcnt-1].strip() 
    261             format_string = format_string[1:len(format_string) - 1] 
    262             var_string = print_params[prcnt+1:len(print_params)].strip() 
    263         else: 
    264             format_string = print_params 
    265             var_string = "" 
    266         format_string= format_string[:len(format_string)] 
    267         self.write_c('"' + format_string + '\\n"') 
    268         if (len(var_string) > 0): 
    269             self.write_c(', ' + var_string) 
    270         self.write_c(");") 
    271         self.add_current_line() 
    272  
    273  
    274  
    275269    def newline(self, node=None, extra=0): 
    276270        self.new_lines = max(self.new_lines, 1 + extra) 
    277271        if node is not None and self.add_line_information: 
    278             self.write_c('# line: %s' % node.lineno) 
     272            self.write_c('// line: %s' % node.lineno) 
    279273            self.new_lines = 1 
    280274        if self.current_statement: 
    … …  
    297291        self.body(node.body) 
    298292        if node.orelse: 
     293            self.unsupported(node, "for...else/while...else not supported") 
     294 
    299295            self.newline() 
    300296            self.write_c('else:') 
    … …  
    326322                except AttributeError: 
    327323                    arg_name = arg.id 
    328                 w_str = ("Default Parameters are unknown to C: '%s = %s" \ 
    329                         % (arg_name, str(default.n))) 
     324                w_str = ("C does not support default parameters: %s=%s" 
     325                         % (arg_name, str(default.n))) 
    330326                self.warnings.append(w_str) 
    331327 
    … …  
    366362 
    367363    def add_semi_colon(self): 
    368         semi_pos = self.current_statement.find(';') 
    369         if semi_pos > 0.0: 
    370             self.current_statement = self.current_statement.replace(';', '') 
     364        #semi_pos = self.current_statement.find(';') 
     365        #if semi_pos >= 0: 
     366        #    self.current_statement = self.current_statement.replace(';', '') 
    371367        self.write_c(';') 
    372368 
    373369    def visit_Assign(self, node): 
    374370        self.add_current_line() 
     371        self.in_expr = True 
    375372        for idx, target in enumerate(node.targets): # multi assign, as in 'a = b = c = 7' 
    376373            if idx: 
    … …  
    378375            self.define_c_vars(target) 
    379376            self.visit(target) 
    380             if hasattr(target,'id'): 
    381                 self.assign_target = target.id 
    382377        # Capture assigned tuple names, if any 
    383378        targets = self.tuples[:] 
    … …  
    387382        self.visited_args = False 
    388383        self.visit(node.value) 
    389         if len(self.assign_c_vector): 
    390             for c_statement in self.assign_c_vector: 
    391                 self.add_c_line (c_statement) 
    392 #                self.c_proc.append(c_statement) 
    393             self.assign_c_vector.clear() 
    394         else: 
     384        self.add_semi_colon() 
     385        self.add_current_line() 
     386        # Assign tuples to tuples, if any 
     387        # TODO: doesn't handle swap:  a,b = b,a 
     388        for target, item in zip(targets, self.tuples): 
     389            self.visit(target) 
     390            self.write_c(' = ') 
     391            self.visit(item) 
    395392            self.add_semi_colon() 
    396393            self.add_current_line() 
    397         # Assign tuples to tuples, if any 
    398         # TODO: doesn't handle swap:  a,b = b,a 
    399             for target, item in zip(targets, self.tuples): 
    400                 self.visit(target) 
    401                 self.write_c(' = ') 
    402                 self.visit(item) 
    403                 self.add_semi_colon() 
    404                 self.add_current_line() 
    405             if self.is_sequence and not self.visited_args: 
    406                 for target in node.targets: 
    407                     if hasattr(target, 'id'): 
    408                         if target.id in self.c_vars and target.id not in self.c_dcl_pointers: 
    409                             if target.id not in self.c_dcl_pointers: 
    410                                 self.c_dcl_pointers.append(target.id) 
    411                                 if target.id in self.c_vars: 
    412                                     self.c_vars.remove(target.id) 
     394        if self.is_sequence and not self.visited_args: 
     395            for target in node.targets: 
     396                if hasattr(target, 'id'): 
     397                    if target.id in self.c_vars and target.id not in self.c_dcl_pointers: 
     398                        if target.id not in self.c_dcl_pointers: 
     399                            self.c_dcl_pointers.append(target.id) 
     400                            if target.id in self.c_vars: 
     401                                self.c_vars.remove(target.id) 
    413402        self.current_statement = '' 
    414         self.assign_target = '' 
     403        self.in_expr = False 
    415404 
    416405    def visit_AugAssign(self, node): 
    … …  
    418407            if node.target.id not in self.arguments: 
    419408                self.c_vars.append(node.target.id) 
     409        self.in_expr = True 
    420410        self.visit(node.target) 
    421411        self.write_c(' ' + BINOP_SYMBOLS[type(node.op)] + '= ') 
    422412        self.visit(node.value) 
    423413        self.add_semi_colon() 
     414        self.in_expr = False 
    424415        self.add_current_line() 
    425416 
    … …  
    441432 
    442433    def visit_Expr(self, node): 
     434        #self.in_expr = True 
    443435        self.newline(node) 
    444436        self.generic_visit(node) 
     437        #self.in_expr = False 
    445438 
    446439    def write_c_pointers(self, start_var): 
    … …  
    469462            have_decls = True 
    470463            start_var += 1 
    471         if len(self.C_Vectors) > 0: 
    472             vec_declare = [] 
    473             for c_vec in self.C_Vectors: 
    474                 item_declare = c_vec.item_declare_string() 
    475                 if item_declare not in vec_declare: 
    476                     vec_declare.append (item_declare) 
    477                 if c_vec.name in self.c_vars: 
    478                     self.c_vars.remove(c_vec.name) 
    479             all_declare = "    double " + ", ".join (vec_declare) + ";\n" 
    480             self.c_proc.insert(start_var, all_declare) 
    481             start_var += 1 
    482464 
    483465        if self.c_vars: 
    … …  
    519501        self.current_function = node.name 
    520502 
     503        # remember the location of the next warning that will be inserted 
     504        # so that we can stuff the function name ahead of the warning list 
     505        # if any warnings are generated by the function. 
    521506        warning_index = len(self.warnings) 
     507 
    522508        self.newline(extra=1) 
    523509        self.decorators(node) 
    … …  
    533519        self.insert_signature() 
    534520        self.insert_c_vars(start_vars) 
    535 #        self.c_pointers = [] 
    536         self.c_pointers.clear() 
    537         self.C_Vectors.clear() 
     521        del self.c_pointers[:] 
    538522        self.current_function = "" 
    539         if (warning_index != len(self.warnings)): 
     523        if warning_index != len(self.warnings): 
    540524            self.warnings.insert(warning_index, "Warning in function '" + node.name + "':") 
    541             self.warnings.append("Note: C compilation will fail") 
    542  
    543525 
    544526    def visit_ClassDef(self, node): 
    … …  
    576558 
    577559    def visit_If(self, node): 
    578         self.add_current_line() 
    579         self.inside_if = True 
     560 
    580561        self.write_c('if ') 
     562        self.in_expr = True 
    581563        self.visit(node.test) 
     564        self.in_expr = False 
    582565        self.write_c(' {') 
    583566        self.body(node.body) 
    … …  
    592575                #self.newline() 
    593576                self.write_c('else if ') 
     577                self.in_expr = True 
    594578                self.visit(node.test) 
     579                self.in_expr = False 
    595580                self.write_c(' {') 
    596581                self.body(node.body) 
    … …  
    601586                self.newline() 
    602587                self.write_c('else {') 
    603                 self.body(node.orelse) 
     588                self.body(else_) 
    604589                self.add_c_line('}') 
    605590                break 
    606         self.inside_if = False 
    607591 
    608592    def get_for_range(self, node): 
    … …  
    631615        return start, stop, step 
    632616 
    633     def add_c_int_var (self, iterator): 
    634         if iterator not in self.c_int_vars: 
    635             self.c_int_vars.append(iterator) 
     617    def add_c_int_var(self, name): 
     618        if name not in self.c_int_vars: 
     619            self.c_int_vars.append(name) 
    636620 
    637621    def visit_For(self, node): 
    … …  
    647631                    iterator = self.current_statement 
    648632                    self.current_statement = '' 
    649                     self.add_c_int_var (iterator) 
     633                    self.add_c_int_var(iterator) 
    650634                    start, stop, step = self.get_for_range(node) 
    651                     self.write_c("for(" + iterator + "=" + str(start) + 
     635                    self.write_c("for (" + iterator + "=" + str(start) + 
    652636                                 " ; " + iterator + " < " + str(stop) + 
    653637                                 " ; " + iterator + " += " + str(step) + ") {") 
    … …  
    664648            self.write_c(':') 
    665649            # report the error 
    666             self.unsupported("unsupported " + self.current_statement) 
     650            self.unsupported(node, "unsupported " + self.current_statement) 
    667651 
    668652    def visit_While(self, node): 
    … …  
    670654        self.write_c('while ') 
    671655        self.visit(node.test) 
    672         self.write_c(' { ') 
     656        self.write_c(' {') 
    673657        self.body_or_else(node) 
    674658        self.write_c('}') 
    … …  
    677661    def visit_With(self, node): 
    678662        self.unsupported(node) 
     663 
    679664        self.newline(node) 
    680665        self.write_python('with ') 
    … …  
    691676 
    692677    def visit_Print(self, node): 
    693         # TODO: print support would be nice, though hard to do 
    694678        self.unsupported(node) 
     679 
    695680        # CRUFT: python 2.6 only 
    696681        self.newline(node) 
    … …  
    711696    def visit_Delete(self, node): 
    712697        self.unsupported(node) 
     698 
    713699        self.newline(node) 
    714700        self.write_python('del ') 
    … …  
    720706    def visit_TryExcept(self, node): 
    721707        self.unsupported(node) 
     708 
    722709        self.newline(node) 
    723710        self.write_python('try:') 
    … …  
    728715    def visit_TryFinally(self, node): 
    729716        self.unsupported(node) 
     717 
    730718        self.newline(node) 
    731719        self.write_python('try:') 
    … …  
    737725    def visit_Global(self, node): 
    738726        self.unsupported(node) 
     727 
    739728        self.newline(node) 
    740729        self.write_python('global ' + ', '.join(node.names)) 
    … …  
    746735    def visit_Return(self, node): 
    747736        self.add_current_line() 
     737        self.in_expr = True 
    748738        if node.value is None: 
    749739            self.write_c('return') 
    750740        else: 
    751             self.write_c('return(') 
     741            self.write_c('return ') 
    752742            self.visit(node.value) 
    753         self.write_c(')') 
    754743        self.add_semi_colon() 
     744        self.in_expr = False 
    755745        self.add_c_line(self.current_statement) 
    756746        self.current_statement = '' 
    … …  
    766756    def visit_Raise(self, node): 
    767757        self.unsupported(node) 
     758 
    768759        # CRUFT: Python 2.6 / 3.0 compatibility 
    769760        self.newline(node) 
    … …  
    788779    def visit_Attribute(self, node): 
    789780        self.unsupported(node, "attribute reference a.b not supported") 
     781 
    790782        self.visit(node.value) 
    791783        self.write_python('.' + node.attr) 
    … …  
    807799            elif node.func.id == "SINCOS": 
    808800                self.write_sincos(node) 
    809                 return 
    810             elif node.func.id == "print": 
    811                 self.write_print(node) 
    812801                return 
    813802            else: 
    … …  
    835824                self.visit(node.kwargs) 
    836825        self.write_c(')') 
    837         if (self.inside_if == False): 
    838             self.write_c(';') 
     826        if not self.in_expr: 
     827            self.add_semi_colon() 
     828 
     829    TRANSLATE_CONSTANTS = { 
     830        # python 2 uses normal name references through vist_Name 
     831        'True': 'true', 
     832        'False': 'false', 
     833        'None': 'NULL',  # "None" will probably fail for other reasons 
     834        # python 3 uses NameConstant 
     835        True: 'true', 
     836        False: 'false', 
     837        None: 'NULL',  # "None" will probably fail for other reasons 
     838        } 
    839839 
    840840    def visit_Name(self, node): 
     841        translation = self.TRANSLATE_CONSTANTS.get(node.id, None) 
     842        if translation: 
     843            self.write_c(translation) 
     844            return 
    841845        self.write_c(node.id) 
    842846        if node.id in self.c_pointers and not self.in_subref: 
    … …  
    853857                name not in self.c_constants and not name.isdigit()): 
    854858            if self.in_subscript: 
    855                 self.add_c_int_var (self, node.id) 
     859                self.add_c_int_var(node.id) 
    856860            else: 
    857861                self.c_vars.append(node.id) 
    858862 
     863    def visit_NameConstant(self, node): 
     864        translation = self.TRANSLATE_CONSTANTS.get(node.value, None) 
     865        if translation is not None: 
     866            self.write_c(translation) 
     867        else: 
     868            self.unsupported(node, "don't know how to translate %r"%node.value) 
     869 
    859870    def visit_Str(self, node): 
    860         self.write_c(repr(node.s)) 
     871        s = node.s 
     872        s = s.replace('\\','\\\\').replace('"','\\"').replace('\n','\\n') 
     873        self.write_c('"') 
     874        self.write_c(s) 
     875        self.write_c('"') 
    861876 
    862877    def visit_Bytes(self, node): 
    863         self.write_c(repr(node.s)) 
     878        s = node.s 
     879        s = s.replace('\\','\\\\').replace('"','\\"').replace('\n','\\n') 
     880        self.write_c('"') 
     881        self.write_c(s) 
     882        self.write_c('"') 
    864883 
    865884    def visit_Num(self, node): 
    … …  
    876895        def visit(self, node): 
    877896            self.is_sequence = True 
    878             c_vec = C_Vector() 
    879897            s = "" 
    880898            for idx, item in enumerate(node.elts): 
    … …  
    882900                    s += ', ' 
    883901                if hasattr(item, 'id'): 
    884                     c_vec.values.append(item.id) 
    885 #                    s += item.id 
     902                    s += item.id 
    886903                elif hasattr(item, 'n'): 
    887                     c_vec.values.append(item.n) 
    888 #                    s += str(item.n) 
    889                 else: 
    890                     temp = self.current_statement 
    891                     self.current_statement = "" 
    892                     self.visit(item) 
    893                     c_vec.values.append(self.current_statement) 
    894 #                    s += self.current_statement 
    895                     self.current_statement = temp 
    896             if (self.assign_target): 
    897                 self.add_c_int_var ('n') 
    898                 c_vec.name = self.assign_target 
    899                 self.C_Vectors.append(c_vec) 
    900                 assign_strings = c_vec.get_assign() 
    901                 for assign_str in assign_strings: 
    902                     self.assign_c_vector.append(assign_str) 
    903 #                    self.add_c_line(assign_str) 
    904 #            if s: 
    905 #                self.c_vectors.append(s) 
    906 #                vec_name = "vec"  + str(len(self.c_vectors)) 
    907 #                self.write_c(vec_name) 
     904                    s += str(item.n) 
     905            if s: 
     906                self.c_vectors.append(s) 
     907                vec_name = "vec"  + str(len(self.c_vectors)) 
     908                self.write_c(vec_name) 
    908909        return visit 
    909910 
    … …  
    914915    def visit_Dict(self, node): 
    915916        self.unsupported(node) 
     917 
    916918        self.write_python('{') 
    917919        for idx, (key, value) in enumerate(zip(node.keys, node.values)): 
    … …  
    964966        if is_negative_exp: 
    965967            self.write_c(")") 
    966         self.write_c(" ") 
     968        #self.write_c(" ") 
    967969 
    968970    def translate_integer_divide(self, node): 
    969         self.write_c("(int)(") 
     971        self.write_c("(int)((") 
    970972        self.visit(node.left) 
    971         self.write_c(") /(int)(") 
     973        self.write_c(")/(") 
    972974        self.visit(node.right) 
    973         self.write_c(")") 
     975        self.write_c("))") 
     976 
     977    def translate_float_divide(self, node): 
     978        self.write_c("((double)(") 
     979        self.visit(node.left) 
     980        self.write_c(")/(double)(") 
     981        self.visit(node.right) 
     982        self.write_c("))") 
    974983 
    975984    def visit_BinOp(self, node): 
    … …  
    979988        elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.FloorDiv]: 
    980989            self.translate_integer_divide(node) 
     990        elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.Div]: 
     991            self.translate_float_divide(node) 
    981992        else: 
    982993            self.visit(node.left) 
    … …  
    10431054    def visit_Yield(self, node): 
    10441055        self.unsupported(node) 
     1056 
    10451057        self.write_python('yield ') 
    10461058        self.visit(node.value) 
    … …  
    10481060    def visit_Lambda(self, node): 
    10491061        self.unsupported(node) 
     1062 
    10501063        self.write_python('lambda ') 
    10511064        self.visit(node.args) 
    … …  
    10551068    def visit_Ellipsis(self, node): 
    10561069        self.unsupported(node) 
     1070 
    10571071        self.write_python('Ellipsis') 
    10581072 
    … …  
    10751089    def visit_DictComp(self, node): 
    10761090        self.unsupported(node) 
     1091 
    10771092        self.write_python('{') 
    10781093        self.visit(node.key) 
    … …  
    10831098        self.write_python('}') 
    10841099 
    1085     def visit_NameConstant (self, node): 
    1086         if (hasattr (node, "value")): 
    1087             if (node.value == True): 
    1088                 val = "1" 
    1089             elif node.value == False: 
    1090                 val = "0" 
    1091             else: 
    1092                 val = "" 
    1093             if (len(val) > 0): 
    1094                 self.write_c("(" + val + ")") 
    1095  
    10961100    def visit_IfExp(self, node): 
    1097         self.write_c('(') 
     1101        self.write_c('((') 
    10981102        self.visit(node.test) 
    1099         self.write_c(' ? ') 
     1103        self.write_c(')?(') 
    11001104        self.visit(node.body) 
    1101         self.write_c(' : ') 
     1105        self.write_c('):(') 
    11021106        self.visit(node.orelse) 
    1103         self.write_c(');') 
     1107        self.write_c('))') 
    11041108 
    11051109    def visit_Starred(self, node): 
    … …  
    11171121    def visit_alias(self, node): 
    11181122        self.unsupported(node) 
     1123 
    11191124        self.write_python(node.name) 
    11201125        if node.asname is not None: 
    … …  
    11791184    const = "constant "  # OpenCL needs globals to be constant 
    11801185    if isinstance(value, int): 
    1181         parts = [const + "int ", name, " = ", "%d"%value, ";"] 
     1186        parts = [const + "int ", name, " = ", "%d"%value, ";\n"] 
    11821187    elif isinstance(value, float): 
    1183         parts = [const + "double ", name, " = ", "%.15g"%value, ";"] 
     1188        parts = [const + "double ", name, " = ", "%.15g"%value, ";\n"] 
    11841189    else: 
    11851190        try: 
    … …  
    11951200        elements = ["%.15g"%v for v in value] 
    11961201        parts = [const + "double ", name, "[]", " = ", 
    1197                  "{\n   ", ", ".join(elements), "\n};"] 
     1202                 "{\n   ", ", ".join(elements), "\n};\n"] 
    11981203 
    11991204    return "".join(parts) 
    … …  
    12391244                         % ", ".join(str(node) for node in dag.keys())) 
    12401245 
     1246import re 
     1247PRINT_ARGS = re.compile(r'print[(]"(?P<template>[^"]*)" *% *[(](?P<args>[^\n]*)[)] *[)] *\n') 
     1248SUBST_ARGS = r'printf("\g<template>\\n", \g<args>)\n' 
     1249PRINT_STR = re.compile(r'print[(]"(?P<template>[^"]*)" *[)] *\n') 
     1250SUBST_STR = r'printf("\g<template>\n")' 
    12411251def translate(functions, constants=None): 
    1242     # type: (List[(str, str, int)], Dict[str, any]) -> List[str] 
     1252    # type: (Sequence[(str, str, int)], Dict[str, any]) -> List[str] 
    12431253    """ 
    1244     Convert a set of functions 
     1254    Convert a list of functions to a list of C code strings. 
     1255 
     1256    Returns list of corresponding code snippets (with trailing lines in 
     1257    each block) and a list of warnings generated by the translator. 
     1258 
     1259    A function is given by the tuple (source, filename, line number). 
     1260 
     1261    Global constants are given in a dictionary of {name: value}.  The 
     1262    constants are used for name space resolution and type inferencing. 
     1263    Constants are not translated by this code. Instead, call 
     1264    :func:`define_constant` with name and value, and maybe block_size 
     1265    if arrays need to be padded to the next block boundary. 
     1266 
     1267    Function prototypes are not generated. Use :func:`ordered_dag` 
     1268    to list the functions in reverse order of dependency before calling 
     1269    translate. [Maybe a future revision will return the function prototypes 
     1270    so that a suitable "*.h" file can be generated. 
    12451271    """ 
    12461272    snippets = [] 
    1247     snippets.append("#include <math.h>") 
    1248     snippets.append("") 
    12491273    warnings = [] 
    12501274    for source, fname, lineno in functions: 
    1251         line_directive = '#line %d "%s"'%(lineno, fname.replace('\\', '\\\\')) 
     1275        line_directive = '#line %d "%s"\n'%(lineno, fname.replace('\\', '\\\\')) 
    12521276        snippets.append(line_directive) 
     1277        # Replace simple print function calls with printf statements 
     1278        source = PRINT_ARGS.sub(SUBST_ARGS, source) 
     1279        source = PRINT_STR.sub(SUBST_STR, source) 
    12531280        tree = ast.parse(source) 
    1254         c_code,w = to_source(tree, constants=constants, fname=fname, lineno=lineno) 
    1255         if w: 
    1256             warnings.append ("\n".join(w)) 
     1281        generator = SourceGenerator(constants=constants, fname=fname, lineno=lineno) 
     1282        generator.visit(tree) 
     1283        c_code = "".join(generator.c_proc) 
    12571284        snippets.append(c_code) 
    1258     return "\n".join(snippets),warnings 
    1259 #    return snippets 
    1260  
    1261 def print_usage (): 
    1262         print("""\ 
    1263             Usage: python py2c.py <infile> [<outfile>] 
    1264             if outfile is omitted, output file is '<infile>.c' 
    1265         """) 
    1266  
    1267 def get_files_names (): 
    1268     import os 
    1269     fname_in = fname_out = "" 
    1270     valid_params =  len(sys.argv) > 1 
    1271     if (valid_params == False): 
    1272         print_usage() 
    1273     else: 
    1274         fname_in = sys.argv[1] 
    1275         if len(sys.argv) == 2: 
    1276             fname_base = os.path.splitext(fname_in)[0] 
    1277             fname_out = str(fname_base) + '.c' 
    1278         else: 
    1279             fname_out = sys.argv[2] 
    1280     return valid_params, fname_in,fname_out 
     1285        warnings.extend(generator.warnings) 
     1286    return snippets, warnings 
     1287 
     1288def to_source(tree, constants=None, fname=None, lineno=0): 
     1289    """ 
     1290    This function can convert a syntax tree into C sourcecode. 
     1291    """ 
     1292    c_code = "".join(generator.c_proc) 
     1293    return c_code 
     1294 
    12811295 
    12821296C_HEADER = """ 
    … …  
    12991313} 
    13001314""" 
     1315 
     1316USAGE = """\ 
     1317Usage: python py2c.py <infile> [<outfile>] 
     1318 
     1319if outfile is omitted, output file is '<infile>.c' 
     1320""" 
     1321 
    13011322def main(): 
    1302     print("Parsing...using Python" + sys.version) 
    1303     valid_params, fname_in, fname_out = get_files_names () 
    1304     if (valid_params == False): 
    1305         print("Input parameters error.\nExiting") 
     1323    import os 
     1324    #print("Parsing...using Python" + sys.version) 
     1325    if len(sys.argv) == 1: 
     1326        print(USAGE) 
    13061327        return 
    13071328 
     1329    fname_in = sys.argv[1] 
     1330    if len(sys.argv) == 2: 
     1331        fname_base = os.path.splitext(fname_in)[0] 
     1332        fname_out = str(fname_base) + '.c' 
     1333    else: 
     1334        fname_out = sys.argv[2] 
     1335 
    13081336    with open(fname_in, "r") as python_file: 
    1309             code = python_file.read() 
     1337        code = python_file.read() 
    13101338    name = "gauss" 
    13111339    code = (code 
    … …  
    13151343            .replace('if __name__ == "__main__"', "def main()") 
    13161344           ) 
    1317     translation,warnings = translate([(code, fname_in, 1)]) 
    1318     translation = translation.replace("double main()", "int main(int argc, char *argv[])") 
     1345 
     1346    translation, warnings = translate([(code, fname_in, 1)]) 
     1347    c_code = "".join(translation) 
     1348    c_code = c_code.replace("double main()", "int main(int argc, char *argv[])") 
    13191349 
    13201350    with open(fname_out, "w") as file_out: 
    13211351        file_out.write(C_HEADER) 
    1322         file_out.write(str(translation)) 
     1352        file_out.write(c_code) 
     1353 
    13231354    if warnings: 
    1324         print("".join(str(warnings[0]))) 
    1325     print("...Done") 
     1355        print("\n".join(warnings)) 
     1356    #print("...Done") 
    13261357 
    13271358if __name__ == "__main__": 
    1328     try: 
    1329         main() 
    1330     except Exception as excp: 
    1331         print ("Error:\n" + str(excp.args)) 
     1359    main() 
Note: See TracChangeset for help on using the changeset viewer.