def set_numexpr_threads(n=None):
    # if we are using numexpr, set the threads to n
    # otherwise reset
    if _NUMEXPR_INSTALLED and _USE_NUMEXPR:
        if n is None:
            n = ne.detect_number_of_cores()
        ne.set_num_threads(n)

def numexpr_darts(nd=200000, nprocs=1):
    """'nd' is the number of darts. 'nprocs' is number of processors. """
    ne.set_num_threads(nprocs)
    
    # Let us define a numpy version of a throw_dart helper function
    def throw_dart():
        ''' Throw "n" darts at a square of length 1.0, and return
            how many times they landed in a concentric circle of radius
            0.5.
        '''
        x = np.random.uniform(size = nd)
        y = np.random.uniform(size = nd)
        expr1 = ne.evaluate('(x - 0.5)**2 + (y - 0.5)**2')
        radius = np.sqrt(expr1)
        in_or_out = ne.evaluate('radius <= 0.5')
        in_or_out = in_or_out.astype(np.int)
        total_inside = ne.evaluate('sum(in_or_out)')
        return total_inside
    
    number_of_darts_in_circle = 0
    
    # Record the simulation time
    start_time = time()
    number_of_darts_in_circle = throw_dart()
    end_time = time()
    execution_time = end_time - start_time
    
    number_of_darts = nd
    pi_approx = 4 * number_of_darts_in_circle / float(number_of_darts)

    # Return a tuple containing:
    # (0) Number of darts
    # (1) Execution time
    # (2) Darts Thrown per second
    return (number_of_darts, execution_time, number_of_darts / execution_time)

def calc_energy(data):
    """Calculate the energy of the entire system.

    :Parameters: **data** -- the standard python data dictionary
    """
    #name relevant variables
    x = data['x']
    y = data['y']
    nv = data['nv']
    rho = 1000 #density of water
    #initialize EK to zero
    l = nv.shape[0]
    area = np.zeros(l)
    for i in xrange(l):
        #first, calculate the area of the triangle
        xCoords = x[nv[i,:]]
        yCoords = y[nv[i,:]]
        #Compute two vectors for the area calculation.
        v1x = xCoords[1] - xCoords[0]
        v2x = xCoords[2] - xCoords[0]
        v1y = yCoords[1] - yCoords[0]
        v2y = yCoords[2] - yCoords[0]
        #calculate the area as the determinant
        area[i] = abs(v1x*v2y - v2x*v1y)
    #get a vector of speeds.
    sdata = calc_speed(data)
    speed = sdata['speed']

    #calculate EK, use numexpr for speed (saves ~15 seconds)
    ne.set_num_threads(ne.detect_number_of_cores())
    ek = ne.evaluate("sum(rho * area * speed * speed, 1)")
    ek = ek / 4
    data['ek'] = ek
    return data

 def test_changing_nthreads_01_dec(self):
     a = linspace(-1, 1, 1e6)
     b = ((.25*a + .75)*a - 1.5)*a - 2
     for nthreads in range(6, 1, -1):
         numexpr.set_num_threads(nthreads)
         c = evaluate("((.25*a + .75)*a - 1.5)*a - 2")
         assert_array_almost_equal(b, c)

def execute(threadCount):
    n = 100000000  # 10 times fewer than C due to speed issues.
    delta = 1.0 / n
    startTime = time()
    set_num_threads(threadCount)
    value = arange(n, dtype=double)
    pi = 4.0 * delta * evaluate('1.0 / (1.0 + ((value - 0.5) * delta) ** 2)').sum()
    elapseTime = time() - startTime
    out(__file__, pi, n, elapseTime, threadCount)

	def __init__(self, endog, exog, **kwargs):
		super(NLRQ, self).__init__(endog, exog, **kwargs)
		ne.set_num_threads(8)
		self._initialize()
		self.PreEstimation = EventHook()
		self.PostVarianceCalculation = EventHook()
		self.PostEstimation = EventHook()
		self.PostInnerStep = EventHook()
		self.PostOuterStep = EventHook()

 def __init__(self, complexity_threshold=2, multicore=True):
     if numexpr_ver is None or numexpr_ver<(2, 0, 0):
         raise ImportError("numexpr version 2.0.0 or better required.")
     self.complexity_threshold = complexity_threshold
     nc = numexpr.detect_number_of_cores()
     if multicore is True:
         multicore = nc
     elif multicore is False:
         multicore = 1
     elif multicore<=0:
         multicore += nc
     numexpr.set_num_threads(multicore)

    def test_multiprocess(self):
        import multiprocessing as mp
        # Check for two threads at least
        numexpr.set_num_threads(2)
        #print "**** Running from main process:"
        _worker()
        #print "**** Running from subprocess:"
        qout = mp.Queue()
        ps = mp.Process(target=_worker, args=(qout,))
        ps.daemon = True
        ps.start()

        result = qout.get()

    def initLayer(self):
        ne.set_num_threads(mp.cpu_count())  # 1 thread per core
        rlayer = (
            self.dlg.comboBox.currentLayer()
        )  # QgsMapLayerRegistry.instance().mapLayersByName(self.dlg.comboBox.currentText())[0]#self.getLayerByName(self.dlg.comboBox.currentText())
        sensorHt = self.dlg.spinBox_sensorHt.value()

        # get list of sun vectors
        vectors = self.skyVectors()
        self.dlg.progressBar.setMaximum(len(vectors))

        scale = rlayer.rasterUnitsPerPixelX()  # assumes square pixels. . .
        bandNum = self.dlg.spinBox_bands.value()
        maxVal = rlayer.dataProvider().bandStatistics(bandNum).maximumValue
        # QgsMessageLog.logMessage("maxVal = %s" % str(maxVal),  "Plugins",  0)
        maxUsrHeight = self.dlg.spinBox_maxHt.value()
        # QgsMessageLog.logMessage("maxUsrHeight = %s" % str(maxUsrHeight),  "Plugins",  0)
        unitZ = maxVal / maxUsrHeight
        # QgsMessageLog.logMessage("unitZ = %s" % str(unitZ),  "Plugins",  0)

        bandCnt = rlayer.bandCount()

        data = self.shaDEM.rasterToArray(rlayer, bandNum)

        # t = time.time()
        a = data["array"].copy()
        adjSensorHt = sensorHt / unitZ
        a = ne.evaluate("a + adjSensorHt")
        # QgsMessageLog.logMessage("Adjusted Sensor Height= %s" % str(adjSensorHt),  "Plugins",  0)
        svfArr = np.zeros(a.shape)
        i = 0

        for vector in vectors:
            # debug - print solar altitude angles
            # QgsMessageLog.logMessage("Vector[%i] solar alt angle: %.2f" % (i+1, math.degrees(math.atan(vector[2]/math.sqrt(vector[0]**2+vector[1]**2)))),  "Profile",  0)

            result = self.shaDEM.ShadowCalc(data, vector, scale, unitZ, maxVal)
            b = result[0]
            dz = result[1]

            svfArr = ne.evaluate("where((b-a) <= 0, svfArr + 1, svfArr)")

            self.dlg.progressBar.setValue(i)
            i += 1

        # t = time.time() - t
        # QgsMessageLog.logMessage("SVF main loop : " + str(t),  "Profile",  0)

        data["array"] = svfArr / self.dlg.spinBox_vectors.value()

        self.saveToFile(data)

def run_experiment(num_iterations):
    previous_threads = set_num_threads(1)

    scratch = zeros(grid_shape)
    grid = zeros(grid_shape)

    block_low = int(grid_shape[0] * .4)
    block_high = int(grid_shape[0] * .5)
    grid[block_low:block_high, block_low:block_high] = 0.005

    start = time.time()
    for i in range(num_iterations):
        evolve(grid, 0.1, scratch)
        grid, scratch = scratch, grid

    set_num_threads(previous_threads)
    return time.time() - start

 def __init__(
         self, expression='in0',
         in_sig=(numpy.complex64,), out_sig=(numpy.complex64,),
         nthreads=None
 ):
     """
     Args:
     expression: either a NumExpr string (in0, in1, ... are the inputs) or
                 a callable (in0, in1 as args) to be used in work()
     in_sig: a list of numpy dtype as input signature
     out_sig: a list of numpy dtype as output signature
     nthreads: how many threads NumExpr should use
     """
     gr.sync_block.__init__(self, "numexpr_evaluate", in_sig, out_sig)
     self._expression = None
     if numexpr and nthreads:
         numexpr.set_num_threads(nthreads)
     self.expression = expression

def calc_speed(data):
    """
    Calculates the speed from ua and va

    :Parameters:
        **data** -- the standard python data dictionary

    .. note:: We use numexpr here because, with multiple threads, it is\
    about 30 times faster than direct calculation.
    """
    #name required variables
    ua = data['ua']
    va = data['va']

    #we can take advantage of multiple cores to do this calculation
    ne.set_num_threads(ne.detect_number_of_cores())
    #calculate the speed at each point.
    data['speed'] = ne.evaluate("sqrt(ua*ua + va*va)")
    return data

def spherical_to_cartesian(ra, dec, threads=1):
    """
    Inputs in degrees.  Outputs x,y,z
    """
    import numexpr as ne
    import math

    ne.set_num_threads(threads)

    pi = math.pi
    rar = ne.evaluate('ra*pi/180.0')
    decr = ne.evaluate('dec*pi/180.0')

    hold1=ne.evaluate('cos(decr)') 

    x = ne.evaluate('cos(rar) * hold1')
    y = ne.evaluate('sin(rar) * hold1')
    z = ne.evaluate('sin(decr)')
 
    return x, y, z

def _great_circle_distance_fast(ra1, dec1, ra2, dec2, threads):
    """
    (Private internal function)
    Returns great circle distance.  Inputs in degrees.
 
    Uses vicenty distance formula - a bit slower than others, but
    numerically stable.

    A faster version than the function above.
    """

    import numexpr as ne
 
    # terminology from the Vicenty formula - lambda and phi and
    # "standpoint" and "forepoint"
    lambs = np.radians(ra1)
    phis = np.radians(dec1)
    lambf = np.radians(ra2)
    phif = np.radians(dec2)
 
    dlamb = lambf - lambs

    #using numexpr
    #nthreads = ne.detect_number_of_cores()
    nthreads = threads
    ne.set_num_threads(nthreads)

    hold1=ne.evaluate('sin(phif)') #calculate these once instead of a few times!
    hold2=ne.evaluate('sin(phis)')
    hold3=ne.evaluate('cos(phif)')
    hold4=ne.evaluate('cos(dlamb)')
    hold5=ne.evaluate('cos(phis)')
    numera = ne.evaluate( 'hold3 * sin(dlamb)')
    numerb = ne.evaluate('hold5 * hold1 - hold2 * hold3 * hold4')
    numer = ne.evaluate('sqrt(numera**2 + numerb**2)')
    denom = ne.evaluate('hold2 * hold1 + hold5 * hold3 * hold4')
    pi=math.pi

    return ne.evaluate('(arctan2(numer, denom))*180.0/pi')

def _spherical_to_cartesian_fast(ra, dec, threads):
    """
    (Private internal function)
    Inputs in degrees.  Outputs x,y,z
    
    A faster version than the function above.
    """
    import numexpr as ne

    #nthreads = ne.detect_number_of_cores()
    nthreads = threads
    ne.set_num_threads(nthreads)

    pi = math.pi
    rar = ne.evaluate('ra*pi/180.0')
    decr = ne.evaluate('dec*pi/180.0')

    hold1=ne.evaluate('cos(decr)') 

    x = ne.evaluate('cos(rar) * hold1')
    y = ne.evaluate('sin(rar) * hold1')
    z = ne.evaluate('sin(decr)')
 
    return x, y, z

import telescope
import receiver
from .. import constants
from ..constants import DTOR, RTOD
from .. import float_type
from ..util import tools, math, fits, hdf, time
from ..util.tools import struct
from ..data import c_interface
from copy import deepcopy
from mpl_toolkits.axes_grid1 import make_axes_locatable
"""

DTOR = np.pi/180.
RTOD = 180./np.pi

ne.set_num_threads(int(os.environ.get('OMP_NUM_THREADS',4))) # Don't use all the processors!

tqu=['T','Q','U']
teb=['T','E','B']
#pol_index = {'T':0, 'Q':1, 'U':2, 'I':0, 'V':3, 'T':0, 'E':1, 'B':2, 'I':0} # For translating letters into array indices

proj_name_to_index = {'Sanson-Flamsteed':0, 'CAR':1, 'SIN':2, 'Healpix':3, 'Sterographic':4, 'Lambert':5, 'CAR00':7, 'rot00':8, 'BICEP':9}
__proj_name_to_index = dict( [(key.lower(), value) for key, value in proj_name_to_index.iteritems()]) # Lower-case version for comparisons.

# Projections tested in testAllProj:
proj_index_in_idl = [0,1,2,4,5]
proj_to_test = [0,1,2,4,5,7]

##########################################################################################
def read(filename, file_type=None, verbose=True):
    """

    from bcolz.py2help import check_output

    # make an absolute path if required, for example when running in a clone
    if not path.isabs(path_):
        path_ = path.join(os.getcwd(), path_)
    # look up the commit using subprocess and git describe
    try:
        # redirect stderr to stdout to make sure the git error message in case
        # we are not in a git repo doesn't appear on the screen and confuse the
        # user.
        label = check_output(["git", "describe"], cwd=path_,
                             stderr=subprocess.STDOUT).strip()
        return label
    except OSError:  # in case git wasn't found
        pass
    except subprocess.CalledProcessError:  # not in git repo
        pass

git_description = _get_git_descrtiption(__path__[0])

# Initialization code for the Blosc and numexpr libraries
_blosc_init()
ncores = detect_number_of_cores()
blosc_set_nthreads(ncores)
# Benchmarks show that using several threads can be an advantage in bcolz
blosc_set_nthreads(ncores)
if numexpr_here:
    numexpr.set_num_threads(ncores)
import atexit
atexit.register(_blosc_destroy)

#.........这里部分代码省略.........
        :math:`f(x,z) = 
        k_m^2 \\left(\\left(\\frac{n(x,z)}{n_m}\\right)^2 -1\\right)`


    See Also
    --------
    odt_to_ri : conversion of the object function :math:`f(\mathbf{r})` 
        to refractive index :math:`n(\mathbf{r})`.


    Notes
    -----
    This implementation can deal with projection angles that are not
    distributed along a circle  about the rotational axis. If there are
    slight deviations from this circle, simply pass the 3D rotational
    positions instead of the 1D angles to the `angles` argument. In
    principle, this should improve the reconstruction. The general
    problem here is that the backpropagation algorithm requires a
    ramp filter in Fourier space that is oriented perpendicular to the
    rotational axis. If the sample does not rotate about a single axis,
    then a 1D parametric representation of this rotation must be found
    to correctly determine the filter in Fourier space. Such a
    parametric representation could e.g. be a spiral between the poles
    of the unit sphere (but this kind of rotation is probably difficult
    to implement experimentally).
    
    Do not use the parameter `lD` in combination with the Rytov
    approximation - the propagation is not correctly described.
    Instead, numerically refocus the sinogram prior to converting
    it to Rytov data (using e.g. :func:`odtbrain.sinogram_as_rytov`)
    with a numerical focusing algorithm (available in the Python
    package :py:mod:`nrefocus`).
    """
    ne.set_num_threads(num_cores)

    if copy:
        uSin = uSin.copy()
        angles = angles.copy()

    # `tilted_axis` is required for several things:
    # 1. the filter |kDx*v + kDy*u| with (u,v,w)==tilted_axis
    # 2. the alignment of the rotational axis with the y-axis
    # 3. the determination of the point coordinates if only
    #    angles in radians are given.
    
    # For (1) we need the exact axis that corresponds to our input data.
    # For (2) and (3) we need `tilted_axis_yz` (see below) which is the
    # axis `tilted_axis` rotated in the detector plane such that its
    # projection onto the detector coincides with the y-axis.
    
    # Normalize input axis
    tilted_axis=norm_vec(tilted_axis)

    # `tilted_axis_yz` is computed by performing the inverse rotation in
    # the x-y plane with `angz`. We will again use `angz` in the transform
    # within the for-loop to rotate each projection according to its 
    # acquisition angle.
    angz = np.arctan2(tilted_axis[0], tilted_axis[1])
    rotmat = np.array([ 
                       [np.cos(angz), -np.sin(angz), 0],
                       [np.sin(angz),  np.cos(angz), 0],
                       [0           ,             0, 1],
                       ])
    # rotate `tilted_axis` onto the y-z plane.
    tilted_axis_yz = norm_vec(np.dot(rotmat, tilted_axis))
    

import os
import numpy as np
import numexpr as ne
ne.set_num_threads(ne.ncores) # inclusive HyperThreading cores
import sys

sys.path.append(os.path.expanduser('~/devel/mapalign/mapalign'))
sys.path.append(os.path.expanduser('~/devel/hcp_corr'))

import embed
import hcp_util

def fisher_r2z(R):
    return ne.evaluate('arctanh(R)')

def fisher_z2r(Z):
    X = ne.evaluate('exp(2*Z)')
    return ne.evaluate('(X - 1) / (X + 1)')

# here we go ...

## parse command line arguments
# first arg is output prefix, e.g. /ptmp/sbayrak/fisher/fisher_
cliarg_out_prfx = sys.argv[1]
# the rest args are the subject path(s), e.g. /ptmp/sbayrak/hcp/*
cliarg_rest = sys.argv[2:]

# list of all subjects as numpy array
subject_list = np.array(cliarg_rest) # e.g. /ptmp/sbayrak/hcp/*

cnt_files = 4

#.........这里部分代码省略.........
        
    jmc, jmm : instance of :func:`multiprocessing.Value` or ``None``
        The progress of this function can be monitored with the 
        :mod:`jobmanager` package. The current step ``jmc.value`` is
        incremented ``jmm.value`` times. ``jmm.value`` is set at the 
        beginning.
    verbose : int
        Increment to increase verbosity.


    Returns
    -------
    f : ndarray of shape (Nx, Ny, Nx), complex if ``onlyreal==False``
        Reconstructed object function :math:`f(\mathbf{r})` as defined
        by the Helmholtz equation.
        :math:`f(x,z) = 
        k_m^2 \\left(\\left(\\frac{n(x,z)}{n_m}\\right)^2 -1\\right)`


    See Also
    --------
    odt_to_ri : conversion of the object function :math:`f(\mathbf{r})` 
        to refractive index :math:`n(\mathbf{r})`.

    Notes
    -----
    Do not use the parameter `lD` in combination with the Rytov
    approximation - the propagation is not correctly described.
    Instead, numerically refocus the sinogram prior to converting
    it to Rytov data (using e.g. :func:`odtbrain.sinogram_as_rytov`)
    with a numerical focusing algorithm (available in the Python
    package :py:mod:`nrefocus`).
    """
    ne.set_num_threads(num_cores)

    if copy:
        uSin = uSin.copy()

    A = angles.shape[0]
    # jobmanager
    if jmm is not None:
        jmm.value = A + 2
    
    # check for dtype
    dtype = np.dtype(dtype)
    assert dtype.name in ["float32", "float64"], "dtype must be float32 or float64!"

    assert num_cores <= _ncores, "`num_cores` must not exceed number " +\
                                 "of physical cores: {}".format(_ncores)

    assert uSin.dtype == np.complex128, "uSin dtype must be complex128."

    dtype_complex = np.dtype("complex{}".format(
        2 * np.int(dtype.name.strip("float"))))

    # set ctype
    ct_dt_map = {np.dtype(np.float32): ctypes.c_float,
                 np.dtype(np.float64): ctypes.c_double
                 }

    assert len(uSin.shape) == 3, "Input data `uSin` must have shape (A,Ny,Nx)."
    assert len(uSin) == A, "`len(angles)` must be  equal to `len(uSin)`."
    assert len(list(padding)) == 2, "Parameter `padding` must be boolean tuple of length 2!"
    assert np.array(padding).dtype is np.dtype(bool), "Parameter `padding` must be boolean tuple."
    assert coords is None, "Setting coordinates is not yet supported."