def _oneEvaluation(self, evaluable):
     """ This method should be called by all optimizers for producing an evaluation. """
     if self._wasUnwrapped:
         self.wrappingEvaluable._setParameters(evaluable)
         res = self.__evaluator(self.wrappingEvaluable)
     elif self._wasWrapped:            
         res = self.__evaluator(evaluable.params)
     else:            
         res = self.__evaluator(evaluable)
         ''' added by JPQ '''
         if self.constrained :
             self.feasible = self.__evaluator.outfeasible
             self.violation = self.__evaluator.outviolation
         # ---
     if isscalar(res):
         # detect numerical instability
         if isnan(res) or isinf(res):
             raise DivergenceError
         # always keep track of the best
         if (self.numEvaluations == 0
             or self.bestEvaluation is None
             or (self.minimize and res <= self.bestEvaluation)
             or (not self.minimize and res >= self.bestEvaluation)):
             self.bestEvaluation = res
             self.bestEvaluable = evaluable.copy()
     
     self.numEvaluations += 1
     
     # if desired, also keep track of all evaluables and/or their fitness.                        
     if self.storeAllEvaluated:
         if self._wasUnwrapped:            
             self._allEvaluated.append(self.wrappingEvaluable.copy())
         elif self._wasWrapped:            
             self._allEvaluated.append(evaluable.params.copy())
         else:            
             self._allEvaluated.append(evaluable.copy())        
     if self.storeAllEvaluations:
         if self._wasOpposed and isscalar(res):
             ''' added by JPQ '''
             if self.constrained :
                 self._allEvaluations.append([-res,self.feasible,self.violation])
             # ---
             else:
                 self._allEvaluations.append(-res)
         else:
             ''' added by JPQ '''
             if self.constrained :
                 self._allEvaluations.append([res,self.feasible,self.violation])
             # ---
             else:
                 self._allEvaluations.append(res)
     ''' added by JPQ '''
     if self.constrained :
         return [res,self.feasible,self.violation]
     else:
     # ---
         return res

def angcomp(ra1,dec1,ra2,dec2,method=3):
    """
    Return the delta_RA and delta_dec (delta = 1-2) components of the angular
    distance between objects.  This is simply an alternate output of the
    angular_distance function above. Distance is returned as degrees, and method chooses a more or less accurate way of determining the distance (with 1 being the fastest/least accurate).
    from astorlib.py

    # UNITS: degrees (output), degrees (input)
    """
    DEGRAD = pi/180.
    import scipy
    if scipy.isscalar(ra1) and scipy.isscalar(ra2):
	from numpy import cos,sin,sqrt
	if ra1-ra2>180:
		ra1 -= 360.
	elif ra2-ra1>180:
		ra2 -= 360.
    else:
	from scipy import cos,sin,sqrt
	if scipy.isscalar(ra1):
		t = ra2.copy()
		ra2 = ra2*0. + ra1
		ra1 = t.copy()
		t = dec2.copy()
		dec2 = dec2*0. + dec1
		dec1 = t.copy()
		del t

	ra1 = scipy.where(ra1-ra2>180,ra1-360.,ra1)
	ra2 = scipy.where(ra2-ra1>180,ra2-360.,ra2)

    ra1 = ra1*DEGRAD
    dec1 = dec1*DEGRAD
    ra2 = ra2*DEGRAD
    dec2 = dec2*DEGRAD

    if method==1:
	deltadec = dec1-dec2
	deltara = (ra1-ra2)*cos(dec2)
    else:
	div = 1.
	if method==3:
		div = sin(dec2)*sin(dec1)+cos(dec2)*cos(dec1)*cos((ra1-ra2))
	deltara = cos(dec2)*sin(ra1-ra2)/div
	deltadec = -(sin(dec2)*cos(dec1)-cos(dec2)*sin(dec1)*cos(ra1-ra2)/div)
	#if sum(div == 0) != 0: #attempt at making array compatable but doesn't
	#work for single integers. Note that could just remove this section of
	#the code since it only here for QA
	if div == 0:
	    import sys
	    print 'tools: div = 0, exiting'
	    sys.exit()

    return deltara/DEGRAD, deltadec/DEGRAD

def angcomp(ra1,dec1,ra2,dec2,method=3):
    """
    Return the delta_RA and delta_dec (delta = 1-2) components of the angular
    distance between objects.  This is simply an alternate output of the
    angular_distance function above. Distance is returned as degrees, and method chooses a more or less accurate way of determining the distance (with 1 being the fastest/least accurate).
    from astorlib.py
    """
    DEGRAD = pi/180.
    import scipy
    if scipy.isscalar(ra1) and scipy.isscalar(ra2):
	from math import cos,sin,sqrt
	if ra1-ra2>180:
		ra1 -= 360.
	elif ra2-ra1>180:
		ra2 -= 360.
    else:
	from scipy import cos,sin,sqrt
	if scipy.isscalar(ra1):
		t = ra2.copy()
		ra2 = ra2*0. + ra1
		ra1 = t.copy()
		t = dec2.copy()
		dec2 = dec2*0. + dec1
		dec1 = t.copy()
		del t

	ra1 = scipy.where(ra1-ra2>180,ra1-360.,ra1)
	ra2 = scipy.where(ra2-ra1>180,ra2-360.,ra2)

    ra1 = ra1*DEGRAD
    dec1 = dec1*DEGRAD
    ra2 = ra2*DEGRAD
    dec2 = dec2*DEGRAD

    if method==1:
	deltadec = dec1-dec2
	deltara = (ra1-ra2)*cos(dec2)
    else:
	div = 1.
	if method==3:
		div = sin(dec2)*sin(dec1)+cos(dec2)*cos(dec1)*cos((ra1-ra2))
	deltara = cos(dec2)*sin(ra1-ra2)/div
	deltadec = -(sin(dec2)*cos(dec1)-cos(dec2)*sin(dec1)*cos(ra1-ra2)/div)

    if isinstance(div, float) and div == 0:
        raise ValueError("dividing by div = 0")
    elif isinstance(div, numpy.ndarray) and div.any() == 0:
        raise ValueError("dividing by div = 0")


    return deltara/DEGRAD, deltadec/DEGRAD

def generate_tolerances(net, rtol, atol=None):
    if rtol == None:
        rtol = global_rtol
    if scipy.isscalar(rtol):
        rtol = scipy.ones(len(net.dynamicVars)) * rtol

    # We set atol to be a minimum of global_atol to avoid variables with large
    # typical values going negative.
    if (scipy.isscalar(atol) or atol==None):
        typ_vals = [abs(net.get_var_typical_val(id))
                    for id in net.dynamicVars.keys()]
        atol = rtol * scipy.asarray(typ_vals)
        if global_atol:
            atol = scipy.minimum(atol, global_atol)
    return rtol, atol

    def oscillate_hessian(self, params, eps=1e-5,
                          relativeScale=True, stepSizeCutoff=1e-6,
                          verbose=False, f0_zero=False):
        """
        Same as hessian except uses oscillate cost function

        """

	nOv = len(params)
        if scipy.isscalar(eps):
            eps = scipy.ones(len(params), scipy.float_) * eps
	    
        ## compute cost at f(x), set f0 to zero if f0_zero is true
        if f0_zero:
            f0 = 0
        else:
            f0 = self.oscillate_cost(params)

	hess = scipy.zeros((nOv, nOv), scipy.float_)

	    ## compute all (numParams*(numParams + 1))/2 unique hessian elements
        for i in range(nOv):
            for j in range(i, nOv):
                hess[i][j] = self.hessian_elem(self.oscillate_cost, f0,
                                               params, i, j, eps[i], eps[j], 
                                               relativeScale, stepSizeCutoff,
                                               verbose)
                hess[j][i] = hess[i][j]

        return hess

    def periodic_hessian_log_params(self, params, eps=1e-5,
                           relativeScale=False, stepSizeCutoff=1e-6,
                           verbose=False):
        """
        Same as hessian_log_params except uses periodic cost function
        Relative scale is false by default here
        """
	nOv = len(params)
        if scipy.isscalar(eps):
            eps = scipy.ones(len(params), scipy.float_) * eps

	## compute cost at f(x)
	f0 = self.periodic_cost_log_params(scipy.log(params))

	hess = scipy.zeros((nOv, nOv), scipy.float_)

	## compute all (numParams*(numParams + 1))/2 unique hessian elements
        for i in range(nOv):
            for j in range(i, nOv):
                hess[i][j] = self.hessian_elem(self.periodic_cost_log_params, f0,
                                               scipy.log(params), 
                                               i, j, eps[i], eps[j], 
                                               relativeScale, stepSizeCutoff,
                                               verbose)
                hess[j][i] = hess[i][j]

        return hess

    def hessian_log_params(self, params, eps,
                           relativeScale=False, stepSizeCutoff=1e-6,
                           verbose=False):
        """
        Returns the hessian of the model in log parameters.

        eps: Sets the stepsize to try
        relativeScale: If True, step i is of size p[i] * eps, otherwise it is
                       eps
        stepSizeCutoff: The minimum stepsize to take
        vebose: If True, a message will be printed with each hessian element
                calculated
        """
	nOv = len(params)
        if scipy.isscalar(eps):
            eps = scipy.ones(len(params), scipy.float_) * eps

	## compute cost at f(x)
	f0 = self.cost_log_params(scipy.log(params))

	hess = scipy.zeros((nOv, nOv), scipy.float_)

	## compute all (numParams*(numParams + 1))/2 unique hessian elements
        for i in range(nOv):
            for j in range(i, nOv):
                hess[i][j] = self.hessian_elem(self.cost_log_params, f0,
                                               scipy.log(params), 
                                               i, j, eps[i], eps[j], 
                                               relativeScale, stepSizeCutoff,
                                               verbose)
                hess[j][i] = hess[i][j]

        return hess

 def __init__(self,net,prior=S.array([100,1])):
     if SP.isscalar(prior):
         self.fixE1=prior
         CNodeEps.__init__(self,net)
     else:
         self.fixE1 = None
         CNodeEps.__init__(self,net,prior)

def errorScalingFactor(observable, beta):
  """
  Look up the numerical factors to apply to the sky averaged parallax error in order to obtain error
  values for a given astrometric parameter, taking the Ecliptic latitude and the number of transits into
  account.

  Parameters
  ----------

  observable - Name of astrometric observable (one of: alphaStar, delta, parallax, muAlphaStar, muDelta)
  beta       - Values(s) of the Ecliptic latitude.

  Returns
  -------

  Numerical factors to apply to the errors of the given observable.
  """
  if isscalar(beta):
    index=int(floor(abs(sin(beta))*_numStepsSinBeta))
    if index == _numStepsSinBeta:
      return _astrometricErrorFactors[observable][_numStepsSinBeta-1]
    else:
      return _astrometricErrorFactors[observable][index]
  else:
    indices = array(floor(abs(sin(beta))*_numStepsSinBeta), dtype=int)
    indices[(indices==_numStepsSinBeta)] = _numStepsSinBeta-1
    return _astrometricErrorFactors[observable][indices]

def _compute_qth_percentile(sorted, q, axis, out):
    if not isscalar(q):
        p = [_compute_qth_percentile(sorted, qi, axis, None)
             for qi in q]

        if out is not None:
            out.flat = p

        return p

    q = q / 100.0
    if (q < 0) or (q > 1):
        raise ValueError("percentile must be either in the range [0,100]")

    indexer = [slice(None)] * sorted.ndim
    Nx = sorted.shape[axis]
    index = q * (Nx - 1)
    i = int(index)
    if i == index:
        indexer[axis] = slice(i, i + 1)
        weights = array(1)
        sumval = 1.0
    else:
        indexer[axis] = slice(i, i + 2)
        j = i + 1
        weights = array([(j - index), (index - i)], float)
        wshape = [1] * sorted.ndim
        wshape[axis] = 2
        weights.shape = wshape
        sumval = weights.sum()

    # Use add.reduce in both cases to coerce data type as well as
    # check and use out array.
    return add.reduce(sorted[indexer] * weights, axis=axis, out=out) / sumval

    def get_data(self, pos, phase=None):
        """return voltage data for relative position and phase
        
        The relative position vector is matched with neuron_data.horizon and
        
        :Parameters:
            pos : ndarray
                The relative position in the dataset.
            phase : sequence
                The phase within the waveform to retrieve. Either a single
                sample index or a sequence. If None is given, return the
                complete waveform.
                Default=None
        """

        # checking pos
        if vector_norm(pos) > self.horizon:
            raise BeyondHorizonError('norm of relative position [%s] '
                                     'lies beyond the horizon [%s]!'
                                     % (vector_norm(pos), self.horizon))

        # check for phase
        if phase is None:
            phase = xrange(self.intra_v.size)
        if N.isscalar(phase):
            phase = xrange(phase, phase + 1)

        # call get_data implementation
        return self._get_data(pos, phase)

 def callback(x):
     if sp.isscalar(x):
         residuals.append(x)
     else:
         residuals.append(residual_norm(A, x, b))
     if cb is not None:
         cb(x)

 def _bestFound(self):
     """ return the best found evaluable and its associated fitness. """
     bestE = self.bestEvaluable.params.copy() if self._wasWrapped else self.bestEvaluable
     if self._wasOpposed and isscalar(self.bestEvaluation):
         bestF = -self.bestEvaluation
     else:
         bestF = self.bestEvaluation
     return bestE, bestF

 def __init__(self, C, taud, tauf, U):
     if isinstance(C, DelayConnection):
         raise AttributeError, "STP does not handle heterogeneous connections yet."
     NetworkOperation.__init__(self, lambda:None, clock=C.source.clock)
     N = len(C.source)
     P = STPGroup(N, clock=C.source.clock)
     P.x = 1
     P.u = U
     P.ux = U
     if (isscalar(taud) & isscalar(tauf) & isscalar(U)):
         updater = STPUpdater(C.source, P, taud, tauf, U, delay=C.delay * C.source.clock.dt)
     else:
         updater = STPUpdater2(C.source, P, taud, tauf, U, delay=C.delay * C.source.clock.dt)
     self.contained_objects = [updater]
     C.source = P
     C.delay = 0
     C._nstate_mod = 0 # modulation of synaptic weights
     self.vars = P

 def jacobian(self, reference_point):
     jac = self._element.function_gradient(self._mesh_entity.vertex_coords(),
                                           reference_point)
     if sp.isscalar(jac):
         return sp.array([[jac]])
     elif sp.all(sp.array(jac.shape) == 1):
         return jac.reshape((1, 1))
     else:
         return jac

def discrete_data(net, params, pts, interval, vars=None, random=False,
                  uncert_func=typ_val_uncert(0.1, 1e-14)):
    """
    Return a set of data points for the given network generated at the given
    parameters.

    net         Network to generate data for
    params      Parameters for this evaluation of the network
    pts         Number of data points to output
    interval    Integration interval
    vars        Variables to output data for, defaults to all species in net
    random      If False data points are distributed evenly over interval
                If True they are spread randomly and uniformly over each
                variable
    uncert_func Function that takes in a trajectory and a variable id and
                returns what uncertainty should be assumed for that variable,
                either as a scalar or a list the same length as the trajectory.
    """
    # Default for vars
    if vars == None:
        vars = net.species.keys()

    # Assign observed times to each variable
    var_times = {}
    for var in vars:
        if random:
            var_times[var] = scipy.rand(pts) * (interval[1]-interval[0]) + interval[0]
        else:
            var_times[var] = scipy.linspace(interval[0], interval[1], pts)

    # Create a sorted list of the unique times in the var_times dict
    int_times = sets.Set(scipy.ravel(var_times.values()))
    int_times.add(0)
    int_times = list(int_times)
    int_times.sort()

    # Get the trajectory
    traj = Dynamics.integrate(net, int_times, params=params, fill_traj=False)

    # Build up our data dictionary
    data = {}; 
    for var, times in var_times.items():
        var_data = {}
        data[var] = var_data
        
        # Calculate our uncertainties
        var_uncerts = uncert_func(traj, var)
        for time in times:
            val = traj.get_var_val(var, time)
            if scipy.isscalar(var_uncerts):
                uncert = var_uncerts
            else:
                index = traj._get_time_index(time)
                uncert = var_uncerts[index]
            var_data[time] = (val, uncert)
    return data

def sky2pix(header,ra,dec):
	hdr_info = parse_header(header)
	if scipy.isscalar(ra):
		ra /= raddeg
		dec /= raddeg
	else:
		ra = ra.astype(scipy.float64)/raddeg
		dec = dec.astype(scipy.float64)/raddeg
	if hdr_info[3]=="DEC":
		ra0 = hdr_info[0][1]/raddeg
		dec0 = hdr_info[0][0]/raddeg
	else:
		ra0 = hdr_info[0][0]/raddeg
		dec0 = hdr_info[0][1]/raddeg

	phi = pi + arctan2(-1*cos(dec)*sin(ra-ra0),sin(dec)*cos(dec0)-cos(dec)*sin(dec0)*cos(ra-ra0))
	argtheta = sin(dec)*sin(dec0)+cos(dec)*cos(dec0)*cos(ra-ra0)
	if scipy.isscalar(argtheta):
		theta = arcsin(argtheta)
	else:
		argtheta[argtheta>1.] = 1.
		theta = arcsin(argtheta)

	if hdr_info[5]=="TAN":
		r_theta = raddeg/tan(theta)
		x = r_theta*sin(phi)
		y = -1.*r_theta*cos(phi)
	elif hdr_info[5]=="SIN":
		r_theta = raddeg*cos(theta)
		x = r_theta*sin(phi)
		y = -1.*r_theta*cos(phi)
	if hdr_info[3]=="DEC":
		a = x.copy()
		x = y.copy()
		y = a.copy()
	inv = linalg.inv(hdr_info[2])
	x0 = inv[0,0]*x + inv[0,1]*y
	y0 = inv[1,0]*x + inv[1,1]*y

	x = x0+hdr_info[1][0]-1
	y = y0+hdr_info[1][1]-1

	return x,y

 def inverse_jacobian(self, reference_point):
     jac = self.jacobian(reference_point)
     if sp.isscalar(jac):
         return sp.array([[1 / jac]])
     elif sp.all(sp.array(jac.shape) == 1):
         return 1 / jac.reshape((1, 1))
     elif len(jac.shape) == 1 or jac.shape[0] == 1 or jac.shape[1] == 1:
         raise NotImplementedError("Maybe a projection approach is needed here.")
     elif jac.shape[0] == jac.shape[1]:
         print(jac.shape, jac)
         return spl.inv(jac)

    def transformSkyCoordinateErrors(self, phi, theta, sigPhiStar, sigTheta, rhoPhiTheta=0):
        """
        Converts the sky coordinate errors from one reference system to another, including the covariance
        term. Equations (1.5.4) and (1.5.20) from section 1.5 in the Hipparcos Explanatory Volume 1 are used.

        Parameters
        ----------

        phi         - The longitude-like angle of the position of the source (radians).
        theta       - The latitude-like angle of the position of the source (radians).
        sigPhiStar  - Standard error in the longitude-like angle of the position of the source (radians or
                      sexagesimal units, including cos(latitude) term)
        sigTheta    - Standard error in the latitude-like angle of the position of the source (radians or
                      sexagesimal units)

        Keywords (optional)
        -------------------

        rhoPhiTheta - Correlation coefficient of the position errors. Set to zero if this keyword is not
                      provided.

        Retuns
        ------
        
        sigPhiRotStar  - The transformed standard error in the longitude-like angle (including
                         cos(latitude) factor)
        sigThetaRot    - The transformed standard error in the latitude-like angle.
        rhoPhiThetaRot - The transformed correlation coefficient.
        """
        if isscalar(rhoPhiTheta) and not isscalar(sigTheta):
            rhoPhiTheta=zeros_like(sigTheta)+rhoPhiTheta
        c, s = self._getJacobian(phi,theta)
        cSqr = c*c
        sSqr = s*s
        covar = sigPhiStar*sigTheta*rhoPhiTheta
        varPhiStar = sigPhiStar*sigPhiStar
        varTheta = sigTheta*sigTheta
        varPhiStarRot = cSqr*varPhiStar+sSqr*varTheta+2.0*covar*c*s
        varThetaRot = sSqr*varPhiStar+cSqr*varTheta-2.0*covar*c*s
        covarRot = (cSqr-sSqr)*covar+c*s*(varTheta-varPhiStar)
        return sqrt(varPhiStarRot), sqrt(varThetaRot), covarRot/sqrt(varPhiStarRot*varThetaRot)

 def BlockInd2SubWithoutBand(self,ind):
     if sp.isscalar(ind):
         if ind < 0 or ind > self.numTotalBlocks:
             raise exceptions.IndexError('BlockInd2SubWithouBand')
     else:
         if ind.any() < 0 or ind.any() >self.numTotalBlocks:
             raise exceptions.IndexError('BlockInd2SubWithouBand')
     sub = []
     for i in range(0,self.Dim):
         sub.append( ind%self.M )
         ind //= self.M
     return a(sub,dtype='int')