Changes in / [d5014e4:ff431ca] in sasmodels


Ignore:
Files:
3 added
7 edited

Legend:

Unmodified
Added
Removed
  • doc/guide/orientation/orientation.rst

    r5fb0634 r82592da  
    44================== 
    55 
    6 With two dimensional small angle diffraction data sasmodels will calculate 
     6With two dimensional small angle diffraction data SasView 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 particle's $a$-$b$-$c$ axes with respect to the incoming 
    12 neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$ 
    13 define the orientation of the $c$-axis of the particle, and angle $\Psi$ is 
    14 defined as the orientation of the major axis of the particle cross section 
    15 with respect to its starting position along the beam direction (or 
    16 equivalently, as rotation about the $c$ axis). There is an unavoidable 
    17 ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both 
    18 serve to rotate the particle about $c$, but this symmetry is destroyed 
    19 when $\theta$ is not a multiple of 180. 
    20  
    21 The figures below are for an elliptical cross section cylinder, but may 
    22 be applied analogously to other shapes of particle. 
     11of the particles with respect to the incoming neutron or X-ray beam. This 
     12is done using three angles: $\theta$ and $\phi$ define the orientation of 
     13the axis of the particle, angle $\Psi$ is defined as the orientation of 
     14the major axis of the particle cross section with respect to its starting 
     15position along the beam direction. The figures below are for an elliptical 
     16cross section cylinder, but may be applied analogously to other shapes of 
     17particle. 
    2318 
    2419.. note:: 
     
    3429 
    3530    Definition of angles for oriented elliptical cylinder, where axis_ratio 
    36     b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$ 
     31    b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$ 
    3732    plane, is carried out first, then rotation $\phi$ about the $z$-axis, 
    3833    finally rotation $\Psi$ is around the axis of the cylinder. The neutron 
    39     or X-ray beam is along the $-z$ axis. 
     34    or X-ray beam is along the $z$ axis. 
    4035 
    4136.. figure:: 
     
    4540    with $\Psi$ = 0. 
    4641 
    47 Having established the mean direction of the particle (the view) we can then 
    48 apply angular orientation distributions (jitter). This is done by a numerical 
    49 integration over a range of angles in a similar way to particle size 
    50 dispersity. 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, 
    52 yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 
    53 $a$-axis. 
    54  
    55 More formally, starting with axes $a$-$b$-$c$ of the particle aligned 
    56 with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity 
    57 is 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$. 
    60 The 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 
    63 using 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  
    70 To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the 
    71 shape 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  
    78 The inverse rotation is easily calculated by rotating the opposite directions 
    79 in 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  
     42Having established the mean direction of the particle we can then apply 
     43angular orientation distributions. This is done by a numerical integration 
     44over a range of angles in a similar way to particle size dispersity. 
     45In the current version of sasview the orientational dispersity is defined 
     46with respect to the axes of the particle. 
    8647 
    8748The $\theta$ and $\phi$ orientation parameters for the cylinder only appear 
    88 when fitting 2d data. On introducing "Orientational Distribution" in the 
    89 angles, "distribution of theta" and "distribution of phi" parameters will 
     49when fitting 2d data. On introducing "Orientational Distribution" in 
     50the angles, "distribution of theta" and "distribution of phi" parameters will 
    9051appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ 
    91 of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder 
    92 cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and 
    93 $X$ axes of the instrument.) The third orientation distribution, in $\Psi$, 
    94 is about the $c$ axis of the particle. Some experimentation may be required 
    95 to understand the 2d patterns fully. A number of different shapes of 
    96 distribution are available, as described for size dispersity, see 
    97 :ref:`polydispersityhelp`. 
     52of 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 
     54instrument.) The third orientation distribution, in $\Psi$, is about the $c$ 
     55axis of the particle. Some experimentation may be required to understand the 
     562d patterns fully. A number of different shapes of distribution are 
     57available, as described for polydispersity, see :ref:`polydispersityhelp` . 
    9858 
    99 Given that the angular dispersion distribution is defined in cartesian space, 
    100 over 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  
    108 but the orientation is defined over a sphere, we are left with a 
    109 `map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_ 
    110 problem, with different tradeoffs depending on how values in $\Delta\theta$ 
    111 and $\Delta\phi$ are translated into latitude/longitude on the sphere. 
    112  
    113 Sasmodels is using the 
    114 `equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_. 
    115 In this projection, square patches in angular dispersity become wedge-shaped 
    116 patches on the sphere. To correct for the changing point density, there is a 
    117 scale factor of $\sin(\Delta\theta)$ that applies to each point in the 
    118 integral. This is not enough, though. Consider a shape which is tumbling 
    119 freely 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 
    121 a distinct angular preference. If the spin axis is along the beam (which 
    122 will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern 
    123 should 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 
    125 that 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$ 
    127 and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to 
    128 symmetry. 
    129  
    130 Alternative projections were considered. 
    131 The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_ 
    132 works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those 
    133 points outside $[-180, 180]$. The distortions are a little less for middle 
    134 ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$ 
    135 and the model is much harder to explain. 
    136 The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_ 
    137 also improves on the equirectangular projection by extending the range of 
    138 reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a 
    139 wedge that cuts to the opposite side of the sphere rather than cutting to the 
    140 pole. This projection has the nice property that distance from the center are 
    141 preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same. 
    142 The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ 
    143 is like the azimuthal equidistance projection, but it preserves area instead 
    144 of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$. 
    145 The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_ 
    146 has an excellent balance with reasonable distortion in both $\Delta\theta$ 
    147 and $\Delta\phi$, as well as preserving small patches. However, it requires 
    148 considerably more computational overhead, and we have not yet derived the 
    149 formula for the distortion correction, measuring the degree of stretch at 
    150 the point $(\Delta\theta, \Delta\phi)$ on the map. 
     59Earlier versions of SasView had numerical integration issues in some 
     60circumstances when distributions passed through 90 degrees. The distributions 
     61in particle coordinates are more robust, but should still be approached with 
     62care for large ranges of angle. 
    15163 
    15264.. note:: 
    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. 
     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. 
    15976 
    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  
    166 Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
    167 used in the integration and the range spanned in number of standard deviations. 
    168 The standard deviation is entered in units of degrees. For a "rectangular" 
    169 distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
    170 The new "uniform" distribution avoids this by letting you directly specify the 
     77Users can experiment with the values of *Npts* and *Nsigs*, the number of steps  
     78used in the integration and the range spanned in number of standard deviations.  
     79The standard deviation is entered in units of degrees. For a "rectangular"  
     80distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.  
     81The new "uniform" distribution avoids this by letting you directly specify the  
    17182half width. 
    17283 
    173 The angular distributions may be truncated outside of the range -180 to +180 
    174 degrees, so beware of using saying a broad Gaussian distribution with large 
    175 value of *Nsigs*, as the array of *Npts* may be truncated to many fewer 
    176 points than would give a good integration,as well as becoming rather 
    177 meaningless. (At some point in the future the actual dispersion arrays may be 
    178 made available to the user for inspection.) 
     84The angular distributions will be truncated outside of the range -180 to +180  
     85degrees, so beware of using saying a broad Gaussian distribution with large value 
     86of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would  
     87give a good integration,as well as becoming rather meaningless. (At some point  
     88in the future the actual polydispersity arrays may be made available to the user  
     89for inspection.) 
    17990 
    18091Some more detailed technical notes are provided in the developer section of 
    18192this manual :ref:`orientation_developer` . 
    18293 
    183 This definition of orientation is new to SasView 4.2.  In earlier versions, 
    184 the orientation distribution appeared as a distribution of view angles. 
    185 This led to strange effects when $c$ was aligned with $z$, where changes 
    186 to the $\phi$ angle served only to rotate the shape about $c$, rather than 
    187 having a consistent interpretation as the pitch of the shape relative to 
    188 the flow field defining the reference orientation.  Prior to SasView 4.1, 
    189 the reference orientation was defined using a Tait-Bryan convention, making 
    190 it difficult to control.  Now, rotation in $\theta$ modifies the spacings 
    191 in the refraction pattern, and rotation in $\phi$ rotates it in the detector 
    192 plane. 
    193  
    194  
    19594*Document History* 
    19695 
    19796| 2017-11-06 Richard Heenan 
    198 | 2017-12-20 Paul Kienzle 
  • doc/guide/plugin.rst

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

    r8cfb486 rff10479  
    165165    # constants in kernel_iq.c 
    166166    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 
    167     'azimuthal_equal_area', 
    168167] 
    169168def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
  • sasmodels/autoc.py

    r15be191 r67cc0ff  
    44from __future__ import print_function 
    55 
     6import ast 
    67import inspect 
     8from functools import reduce 
    79 
    810import numpy as np 
     
    9395                constants[name] = obj 
    9496                # Claim all constants are declared on line 1 
    95                 snippets.append('#line 1 "%s"\n'%escaped_filename) 
    96                 snippets.append(py2c.define_constant(name, obj)) 
     97                snippets.append('#line 1 "%s"'%escaped_filename) 
     98                snippets.append(define_constant(name, obj)) 
    9799            elif isinstance(obj, special.Gauss): 
    98100                for var, value in zip(("N", "Z", "W"), (obj.n, obj.z, obj.w)): 
    99101                    var = "GAUSS_"+var 
    100102                    constants[var] = value 
    101                     snippets.append('#line 1 "%s"\n'%escaped_filename) 
    102                     snippets.append(py2c.define_constant(var, value)) 
     103                    snippets.append('#line 1 "%s"'%escaped_filename) 
     104                    snippets.append(define_constant(var, value)) 
    103105                #libs.append('lib/gauss%d.c'%obj.n) 
    104106                source = (source.replace(name+'.n', 'GAUSS_N') 
     
    119121 
    120122    # translate source 
    121     ordered_code = [code[name] for name in py2c.ordered_dag(depends) if name in code] 
     123    ordered_code = [code[name] for name in ordered_dag(depends) if name in code] 
    122124    functions = py2c.translate(ordered_code, constants) 
    123125    snippets.extend(functions) 
     
    125127    # update model info 
    126128    info.source = unique_libs 
    127     info.c_code = "".join(snippets) 
     129    info.c_code = "\n".join(snippets) 
    128130    info.Iq = info.Iqac = info.Iqabc = info.Iqxy = info.form_volume = None 
     131 
     132def 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 
     155def 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

    r5ab99b7 r67cc0ff  
    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) 
    796797    info.source = getattr(kernel_module, 'source', []) 
    797798    info.c_code = getattr(kernel_module, 'c_code', None) 
     
    809810    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
    810811    # 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)) 
    811814    info.random = getattr(kernel_module, 'random', None) 
    812815 
     
    824827        except Exception as exc: 
    825828            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)) 
    831829 
    832830    if callable(info.Iq) and parameters.has_2d: 
  • sasmodels/models/_cylpy.py

    rc01ed3e r67cc0ff  
    140140py2c = True 
    141141 
    142 # TODO: "#define INVALID (expr)" is not supported 
    143142def invalid(v): 
    144143    return v.radius < 0 or v.length < 0 
     
    207206            phi_pd=10, phi_pd_n=5) 
    208207 
    209 qx, qy = 0.2 * cos(2.5), 0.2 * sin(2.5) 
     208qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    210209# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
    211210tests = [ 
  • sasmodels/py2c.py

    r0bd0877 r0bd0877  
    1 r""" 
    2 py2c 
    3 ~~~~ 
    4  
    5 Convert simple numeric python code into C code. 
    6  
    7 This code is intended to translate direct algorithms for scientific code 
    8 (mostly if statements and for loops operating on double precision values) 
    9 into C code. Unlike projects like numba, cython, pypy and nuitka, the 
    10 :func:`translate` function returns the corresponding C which can then be 
    11 compiled with tinycc or sent to the GPU using CUDA or OpenCL. 
    12  
    13 There is special handling certain constructs, such as *for i in range* and 
    14 small integer powers. 
    15  
    16 **TODO: make a nice list of supported constructs*** 
    17  
    18 Imports are not supported, but they are at least ignored so that properly 
    19 constructed code can be run via python or translated to C without change. 
    20  
    21 Most 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  
    33 There is limited support for list and list comprehensions, so long as they 
    34 can be represented by a fixed array whose size is known at compile time, and 
    35 they are small enough to be stored on the stack. 
    36  
    37 Variables definition in C 
    38 ------------------------- 
    39 Defining variables within the translate function is a bit of a guess work, 
    40 using 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  
    65 Debugging 
    66 --------- 
    67  
    68 *print* is partially supported using a simple regular expression. This 
    69 requires a stylized form. Be sure to use print as a function instead of 
    70 the print statement. If you are including substition variables, use the 
    71 % string substitution style. Include parentheses around the substitution 
    72 tuple, even if there is only one item; do not include the final comma even 
    73 if it is a single item (yes, it won't be a tuple, but it makes the regexp 
    74 much 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  
    80 You can generate *main* using the *if __name__ == "__main__":* construct. 
    81 This does a simple substitution with "def main():" before translation and 
    82 a substitution with "int main(int argc, double *argv[])" after translation. 
    83 The result is that the content of the *if* block becomes the content of *main*. 
    84 Along with the print statement, you can run and test a translation standalone 
    85 using:: 
    86  
    87     python py2c.py source.py 
    88     cc source.c 
    89     ./a.out 
    90  
    91 Known issues 
    92 ------------ 
    93 The 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  
    102 References 
    103 ---------- 
    104  
    105 Based on a variant of codegen.py: 
    106  
    107     https://github.com/andreif/codegen 
     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 
     37Based on codegen.py: 
     38 
    10839    :copyright: Copyright 2008 by Armin Ronacher. 
    10940    :license: BSD. 
    11041""" 
    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  
    144 from __future__ import print_function 
    145  
     42""" 
     43Update Notes 
     44============ 
     4511/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] 
     5211/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 
     5712/07/2017, OE: Translation of integer division, '\\' in python, implemented 
     58                in translate_integer_divide, called from visit_BinOp 
     5912/07/2017, OE: C variable definition handled in 'define_c_vars' 
     60              : Python integer division, '//', translated to C in 
     61                'translate_integer_divide' 
     6212/15/2017, OE: Precedence maintained by writing opening and closing 
     63                parenthesesm '(',')', in procedure 'visit_BinOp'. 
     6412/18/2017, OE: Added call to 'add_current_line()' at the beginning 
     65                of visit_Return 
     6601/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. 
     7001/4/2018, O.E. 'translate(functions, constants=None)' returns string, instaed of a list 
     7101/04/2017 O.E. Fixed bug in 'visit_If': visiting node.orelse in case else exists. 
     7201/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""" 
     78import ast 
    14679import sys 
    147 import ast 
    14880from ast import NodeVisitor 
    14981 
     
    185117 
    186118 
    187 # TODO: should not allow eval of arbitrary python 
     119def 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 
    188129def isevaluable(s): 
    189130    try: 
     
    192133    except Exception: 
    193134        return False 
     135 
     136class 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) 
    194172 
    195173class SourceGenerator(NodeVisitor): 
     
    225203        self.c_vectors = [] 
    226204        self.c_constants = constants if constants is not None else {} 
    227         self.in_expr = False 
    228205        self.in_subref = False 
    229206        self.in_subscript = False 
     
    232209        self.is_sequence = False 
    233210        self.visited_args = False 
     211        self.inside_if = False 
     212        self.C_Vectors = [] 
     213        self.assign_target = "" 
     214        self.assign_c_vector = [] 
    234215 
    235216    def write_python(self, x): 
     
    267248            self.add_unique_var(arg.id) 
    268249 
     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 
    269275    def newline(self, node=None, extra=0): 
    270276        self.new_lines = max(self.new_lines, 1 + extra) 
    271277        if node is not None and self.add_line_information: 
    272             self.write_c('// line: %s' % node.lineno) 
     278            self.write_c('# line: %s' % node.lineno) 
    273279            self.new_lines = 1 
    274280        if self.current_statement: 
     
    291297        self.body(node.body) 
    292298        if node.orelse: 
    293             self.unsupported(node, "for...else/while...else not supported") 
    294  
    295299            self.newline() 
    296300            self.write_c('else:') 
     
    322326                except AttributeError: 
    323327                    arg_name = arg.id 
    324                 w_str = ("C does not support default parameters: %s=%s" 
    325                          % (arg_name, str(default.n))) 
     328                w_str = ("Default Parameters are unknown to C: '%s = %s" \ 
     329                        % (arg_name, str(default.n))) 
    326330                self.warnings.append(w_str) 
    327331 
     
    362366 
    363367    def add_semi_colon(self): 
    364         #semi_pos = self.current_statement.find(';') 
    365         #if semi_pos >= 0: 
    366         #    self.current_statement = self.current_statement.replace(';', '') 
     368        semi_pos = self.current_statement.find(';') 
     369        if semi_pos > 0.0: 
     370            self.current_statement = self.current_statement.replace(';', '') 
    367371        self.write_c(';') 
    368372 
    369373    def visit_Assign(self, node): 
    370374        self.add_current_line() 
    371         self.in_expr = True 
    372375        for idx, target in enumerate(node.targets): # multi assign, as in 'a = b = c = 7' 
    373376            if idx: 
     
    375378            self.define_c_vars(target) 
    376379            self.visit(target) 
     380            if hasattr(target,'id'): 
     381                self.assign_target = target.id 
    377382        # Capture assigned tuple names, if any 
    378383        targets = self.tuples[:] 
     
    382387        self.visited_args = False 
    383388        self.visit(node.value) 
    384         self.add_semi_colon() 
    385         self.add_current_line() 
     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: 
     395            self.add_semi_colon() 
     396            self.add_current_line() 
    386397        # Assign tuples to tuples, if any 
    387398        # 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) 
    392             self.add_semi_colon() 
    393             self.add_current_line() 
    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) 
     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) 
    402413        self.current_statement = '' 
    403         self.in_expr = False 
     414        self.assign_target = '' 
    404415 
    405416    def visit_AugAssign(self, node): 
     
    407418            if node.target.id not in self.arguments: 
    408419                self.c_vars.append(node.target.id) 
    409         self.in_expr = True 
    410420        self.visit(node.target) 
    411421        self.write_c(' ' + BINOP_SYMBOLS[type(node.op)] + '= ') 
    412422        self.visit(node.value) 
    413423        self.add_semi_colon() 
    414         self.in_expr = False 
    415424        self.add_current_line() 
    416425 
     
    432441 
    433442    def visit_Expr(self, node): 
    434         #self.in_expr = True 
    435443        self.newline(node) 
    436444        self.generic_visit(node) 
    437         #self.in_expr = False 
    438445 
    439446    def write_c_pointers(self, start_var): 
     
    462469            have_decls = True 
    463470            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 
    464482 
    465483        if self.c_vars: 
     
    501519        self.current_function = node.name 
    502520 
    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. 
    506521        warning_index = len(self.warnings) 
    507  
    508522        self.newline(extra=1) 
    509523        self.decorators(node) 
     
    519533        self.insert_signature() 
    520534        self.insert_c_vars(start_vars) 
    521         del self.c_pointers[:] 
     535#        self.c_pointers = [] 
     536        self.c_pointers.clear() 
     537        self.C_Vectors.clear() 
    522538        self.current_function = "" 
    523         if warning_index != len(self.warnings): 
     539        if (warning_index != len(self.warnings)): 
    524540            self.warnings.insert(warning_index, "Warning in function '" + node.name + "':") 
     541            self.warnings.append("Note: C compilation will fail") 
     542 
    525543 
    526544    def visit_ClassDef(self, node): 
     
    558576 
    559577    def visit_If(self, node): 
    560  
     578        self.add_current_line() 
     579        self.inside_if = True 
    561580        self.write_c('if ') 
    562         self.in_expr = True 
    563581        self.visit(node.test) 
    564         self.in_expr = False 
    565582        self.write_c(' {') 
    566583        self.body(node.body) 
     
    575592                #self.newline() 
    576593                self.write_c('else if ') 
    577                 self.in_expr = True 
    578594                self.visit(node.test) 
    579                 self.in_expr = False 
    580595                self.write_c(' {') 
    581596                self.body(node.body) 
     
    586601                self.newline() 
    587602                self.write_c('else {') 
    588                 self.body(else_) 
     603                self.body(node.orelse) 
    589604                self.add_c_line('}') 
    590605                break 
     606        self.inside_if = False 
    591607 
    592608    def get_for_range(self, node): 
     
    615631        return start, stop, step 
    616632 
    617     def add_c_int_var(self, name): 
    618         if name not in self.c_int_vars: 
    619             self.c_int_vars.append(name) 
     633    def add_c_int_var (self, iterator): 
     634        if iterator not in self.c_int_vars: 
     635            self.c_int_vars.append(iterator) 
    620636 
    621637    def visit_For(self, node): 
     
    631647                    iterator = self.current_statement 
    632648                    self.current_statement = '' 
    633                     self.add_c_int_var(iterator) 
     649                    self.add_c_int_var (iterator) 
    634650                    start, stop, step = self.get_for_range(node) 
    635                     self.write_c("for (" + iterator + "=" + str(start) + 
     651                    self.write_c("for(" + iterator + "=" + str(start) + 
    636652                                 " ; " + iterator + " < " + str(stop) + 
    637653                                 " ; " + iterator + " += " + str(step) + ") {") 
     
    648664            self.write_c(':') 
    649665            # report the error 
    650             self.unsupported(node, "unsupported " + self.current_statement) 
     666            self.unsupported("unsupported " + self.current_statement) 
    651667 
    652668    def visit_While(self, node): 
     
    654670        self.write_c('while ') 
    655671        self.visit(node.test) 
    656         self.write_c(' {') 
     672        self.write_c(' { ') 
    657673        self.body_or_else(node) 
    658674        self.write_c('}') 
     
    661677    def visit_With(self, node): 
    662678        self.unsupported(node) 
    663  
    664679        self.newline(node) 
    665680        self.write_python('with ') 
     
    676691 
    677692    def visit_Print(self, node): 
     693        # TODO: print support would be nice, though hard to do 
    678694        self.unsupported(node) 
    679  
    680695        # CRUFT: python 2.6 only 
    681696        self.newline(node) 
     
    696711    def visit_Delete(self, node): 
    697712        self.unsupported(node) 
    698  
    699713        self.newline(node) 
    700714        self.write_python('del ') 
     
    706720    def visit_TryExcept(self, node): 
    707721        self.unsupported(node) 
    708  
    709722        self.newline(node) 
    710723        self.write_python('try:') 
     
    715728    def visit_TryFinally(self, node): 
    716729        self.unsupported(node) 
    717  
    718730        self.newline(node) 
    719731        self.write_python('try:') 
     
    725737    def visit_Global(self, node): 
    726738        self.unsupported(node) 
    727  
    728739        self.newline(node) 
    729740        self.write_python('global ' + ', '.join(node.names)) 
     
    735746    def visit_Return(self, node): 
    736747        self.add_current_line() 
    737         self.in_expr = True 
    738748        if node.value is None: 
    739749            self.write_c('return') 
    740750        else: 
    741             self.write_c('return ') 
     751            self.write_c('return(') 
    742752            self.visit(node.value) 
     753        self.write_c(')') 
    743754        self.add_semi_colon() 
    744         self.in_expr = False 
    745755        self.add_c_line(self.current_statement) 
    746756        self.current_statement = '' 
     
    756766    def visit_Raise(self, node): 
    757767        self.unsupported(node) 
    758  
    759768        # CRUFT: Python 2.6 / 3.0 compatibility 
    760769        self.newline(node) 
     
    779788    def visit_Attribute(self, node): 
    780789        self.unsupported(node, "attribute reference a.b not supported") 
    781  
    782790        self.visit(node.value) 
    783791        self.write_python('.' + node.attr) 
     
    799807            elif node.func.id == "SINCOS": 
    800808                self.write_sincos(node) 
     809                return 
     810            elif node.func.id == "print": 
     811                self.write_print(node) 
    801812                return 
    802813            else: 
     
    824835                self.visit(node.kwargs) 
    825836        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         } 
     837        if (self.inside_if == False): 
     838            self.write_c(';') 
    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 
    845841        self.write_c(node.id) 
    846842        if node.id in self.c_pointers and not self.in_subref: 
     
    857853                name not in self.c_constants and not name.isdigit()): 
    858854            if self.in_subscript: 
    859                 self.add_c_int_var(node.id) 
     855                self.add_c_int_var (self, node.id) 
    860856            else: 
    861857                self.c_vars.append(node.id) 
    862858 
    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  
    870859    def visit_Str(self, node): 
    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('"') 
     860        self.write_c(repr(node.s)) 
    876861 
    877862    def visit_Bytes(self, node): 
    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('"') 
     863        self.write_c(repr(node.s)) 
    883864 
    884865    def visit_Num(self, node): 
     
    895876        def visit(self, node): 
    896877            self.is_sequence = True 
     878            c_vec = C_Vector() 
    897879            s = "" 
    898880            for idx, item in enumerate(node.elts): 
     
    900882                    s += ', ' 
    901883                if hasattr(item, 'id'): 
    902                     s += item.id 
     884                    c_vec.values.append(item.id) 
     885#                    s += item.id 
    903886                elif hasattr(item, 'n'): 
    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) 
     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) 
    909908        return visit 
    910909 
     
    915914    def visit_Dict(self, node): 
    916915        self.unsupported(node) 
    917  
    918916        self.write_python('{') 
    919917        for idx, (key, value) in enumerate(zip(node.keys, node.values)): 
     
    966964        if is_negative_exp: 
    967965            self.write_c(")") 
    968         #self.write_c(" ") 
     966        self.write_c(" ") 
    969967 
    970968    def translate_integer_divide(self, node): 
    971         self.write_c("(int)((") 
     969        self.write_c("(int)(") 
    972970        self.visit(node.left) 
    973         self.write_c(")/(") 
     971        self.write_c(") /(int)(") 
    974972        self.visit(node.right) 
    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("))") 
     973        self.write_c(")") 
    983974 
    984975    def visit_BinOp(self, node): 
     
    988979        elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.FloorDiv]: 
    989980            self.translate_integer_divide(node) 
    990         elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.Div]: 
    991             self.translate_float_divide(node) 
    992981        else: 
    993982            self.visit(node.left) 
     
    10541043    def visit_Yield(self, node): 
    10551044        self.unsupported(node) 
    1056  
    10571045        self.write_python('yield ') 
    10581046        self.visit(node.value) 
     
    10601048    def visit_Lambda(self, node): 
    10611049        self.unsupported(node) 
    1062  
    10631050        self.write_python('lambda ') 
    10641051        self.visit(node.args) 
     
    10681055    def visit_Ellipsis(self, node): 
    10691056        self.unsupported(node) 
    1070  
    10711057        self.write_python('Ellipsis') 
    10721058 
     
    10891075    def visit_DictComp(self, node): 
    10901076        self.unsupported(node) 
    1091  
    10921077        self.write_python('{') 
    10931078        self.visit(node.key) 
     
    10981083        self.write_python('}') 
    10991084 
     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 
    11001096    def visit_IfExp(self, node): 
    1101         self.write_c('((') 
     1097        self.write_c('(') 
    11021098        self.visit(node.test) 
    1103         self.write_c(')?(') 
     1099        self.write_c(' ? ') 
    11041100        self.visit(node.body) 
    1105         self.write_c('):(') 
     1101        self.write_c(' : ') 
    11061102        self.visit(node.orelse) 
    1107         self.write_c('))') 
     1103        self.write_c(');') 
    11081104 
    11091105    def visit_Starred(self, node): 
     
    11211117    def visit_alias(self, node): 
    11221118        self.unsupported(node) 
    1123  
    11241119        self.write_python(node.name) 
    11251120        if node.asname is not None: 
     
    11841179    const = "constant "  # OpenCL needs globals to be constant 
    11851180    if isinstance(value, int): 
    1186         parts = [const + "int ", name, " = ", "%d"%value, ";\n"] 
     1181        parts = [const + "int ", name, " = ", "%d"%value, ";"] 
    11871182    elif isinstance(value, float): 
    1188         parts = [const + "double ", name, " = ", "%.15g"%value, ";\n"] 
     1183        parts = [const + "double ", name, " = ", "%.15g"%value, ";"] 
    11891184    else: 
    11901185        try: 
     
    12001195        elements = ["%.15g"%v for v in value] 
    12011196        parts = [const + "double ", name, "[]", " = ", 
    1202                  "{\n   ", ", ".join(elements), "\n};\n"] 
     1197                 "{\n   ", ", ".join(elements), "\n};"] 
    12031198 
    12041199    return "".join(parts) 
     
    12441239                         % ", ".join(str(node) for node in dag.keys())) 
    12451240 
    1246 import re 
    1247 PRINT_ARGS = re.compile(r'print[(]"(?P<template>[^"]*)" *% *[(](?P<args>[^\n]*)[)] *[)] *\n') 
    1248 SUBST_ARGS = r'printf("\g<template>\\n", \g<args>)\n' 
    1249 PRINT_STR = re.compile(r'print[(]"(?P<template>[^"]*)" *[)] *\n') 
    1250 SUBST_STR = r'printf("\g<template>\n")' 
    12511241def translate(functions, constants=None): 
    1252     # type: (Sequence[(str, str, int)], Dict[str, any]) -> List[str] 
     1242    # type: (List[(str, str, int)], Dict[str, any]) -> List[str] 
    12531243    """ 
    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. 
     1244    Convert a set of functions 
    12711245    """ 
    12721246    snippets = [] 
     1247    snippets.append("#include <math.h>") 
     1248    snippets.append("") 
    12731249    warnings = [] 
    12741250    for source, fname, lineno in functions: 
    1275         line_directive = '#line %d "%s"\n'%(lineno, fname.replace('\\', '\\\\')) 
     1251        line_directive = '#line %d "%s"'%(lineno, fname.replace('\\', '\\\\')) 
    12761252        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) 
    12801253        tree = ast.parse(source) 
    1281         generator = SourceGenerator(constants=constants, fname=fname, lineno=lineno) 
    1282         generator.visit(tree) 
    1283         c_code = "".join(generator.c_proc) 
     1254        c_code,w = to_source(tree, constants=constants, fname=fname, lineno=lineno) 
     1255        if w: 
     1256            warnings.append ("\n".join(w)) 
    12841257        snippets.append(c_code) 
    1285         warnings.extend(generator.warnings) 
    1286     return snippets, warnings 
    1287  
    1288 def 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  
     1258    return "\n".join(snippets),warnings 
     1259#    return snippets 
     1260 
     1261def print_usage (): 
     1262        print("""\ 
     1263            Usage: python py2c.py <infile> [<outfile>] 
     1264            if outfile is omitted, output file is '<infile>.c' 
     1265        """) 
     1266 
     1267def 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 
    12951281 
    12961282C_HEADER = """ 
     
    13131299} 
    13141300""" 
    1315  
    1316 USAGE = """\ 
    1317 Usage: python py2c.py <infile> [<outfile>] 
    1318  
    1319 if outfile is omitted, output file is '<infile>.c' 
    1320 """ 
    1321  
    13221301def main(): 
    1323     import os 
    1324     #print("Parsing...using Python" + sys.version) 
    1325     if len(sys.argv) == 1: 
    1326         print(USAGE) 
     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") 
    13271306        return 
    13281307 
    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  
    13361308    with open(fname_in, "r") as python_file: 
    1337         code = python_file.read() 
     1309            code = python_file.read() 
    13381310    name = "gauss" 
    13391311    code = (code 
     
    13431315            .replace('if __name__ == "__main__"', "def main()") 
    13441316           ) 
    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[])") 
     1317    translation,warnings = translate([(code, fname_in, 1)]) 
     1318    translation = translation.replace("double main()", "int main(int argc, char *argv[])") 
    13491319 
    13501320    with open(fname_out, "w") as file_out: 
    13511321        file_out.write(C_HEADER) 
    1352         file_out.write(c_code) 
    1353  
     1322        file_out.write(str(translation)) 
    13541323    if warnings: 
    1355         print("\n".join(warnings)) 
    1356     #print("...Done") 
     1324        print("".join(str(warnings[0]))) 
     1325    print("...Done") 
    13571326 
    13581327if __name__ == "__main__": 
    1359     main() 
     1328    try: 
     1329        main() 
     1330    except Exception as excp: 
     1331        print ("Error:\n" + str(excp.args)) 
Note: See TracChangeset for help on using the changeset viewer.