Pre-computes mask arrays to be used by the 'sharp_round' function
    for roundness computations based on two- and four-fold symmetries.
490009
def precompute_sharp_round(nxk, nyk, xc, yc):
    

    # Create arrays for the two- and four-fold symmetry computations:
    s4m = np.ones((nyk,nxk),dtype=np.int16)
    s4m[yc, xc] = 0

    s2m = np.ones((nyk,nxk),dtype=np.int16)
    s2m[yc, xc] = 0
    s2m[yc:nyk, 0:xc]     = -1;
    s2m[0:yc+1, xc+1:nxk] = -1;

    return s2m, s4m
638433
function getRandomPointOnSphere(r) {
  var u = THREE.Math.randFloat(0, 1);
  var v = THREE.Math.randFloat(0, 1);
  var theta = 2 * Math.PI * u;
  var phi = Math.acos(2 * v - 1);
  var x = r * Math.sin(theta) * Math.sin(phi);
  var y = r * Math.cos(theta) * Math.sin(phi);
  var z = r * Math.cos(phi);
  return {
    x,
    y,
    z,
  };
}
182769
def round_array(array_in):
    
    if isinstance(array_in, ndarray):
        return np.round(array_in).astype(int)
    else:
        return int(np.round(array_in))
1088823
def ZernikeSequence(x,y,maxRadial=5,epsilon=1e-10):
    
    if maxRadial>5:
        raise ValueError('Code for higher radial orders not implemented')
    # Derive radius and exp(i*theta) 
    temp = x + 1j*y
    r1 = np.abs(temp)
    e1 = temp/(r1+epsilon)
     # Generate powers of r recursively
    r2 = r1*r1
    r3 = r2*r1
    r4 = r3*r1
    r5 = r4*r1
    # Generate cos and sin terms recursively from exp(i*theta)
    e2 = e1*e1
    e3 = e2*e1
    e4 = e3*e1
    e5 = e4*e1
    ctheta = e1.real
    stheta = e1.imag
    c2theta = e2.real
    s2theta = e2.imag
    c3theta = e3.real
    s3theta = e3.imag
    c4theta = e4.real
    s4theta = e4.imag
    c5theta = e5.real
    s5theta = e5.imag
    # Generate all the zernikes
    zernike = np.zeros((21,)+x.shape )
    zernike[0] =  1.0
    zernike[1] =  2.0*r1*ctheta
    zernike[2] =  2.0*r1*stheta
    zernike[3] =  sqrt(3.0)*(2.0*r2 - 1.0)
    zernike[4] =  sqrt(6.0)*r2*s2theta
    zernike[5] =  sqrt(6.0)*r2*c2theta
    zernike[6] =  sqrt(8.0)*(3.0*r3 - 2.0*r1)*stheta
    zernike[7] =  sqrt(8.0)*(3.0*r3 - 2.0*r1)*ctheta
    zernike[8] =  sqrt(8.0)*r3*s3theta
    zernike[9] = sqrt(8.0)*r3*c3theta
    zernike[10] = sqrt(5.0)*(6.*r4 - 6.*r2 + 1.)
    zernike[11] = sqrt(10.)*(4.*r4 - 3.*r2)*c2theta
    zernike[12] = sqrt(10.)*(4.*r4 - 3.*r2)*s2theta
    zernike[13] = sqrt(10.)*r4*c4theta
    zernike[14] = sqrt(10.)*r4*s4theta
    zernike[15] = sqrt(12.)*(10*r5-12*r3+3*r1)*ctheta
    zernike[16] = sqrt(12.)*(10*r5-12*r3+3*r1)*stheta
    zernike[17] = sqrt(12.)*(5*r5-4*r3)*c3theta
    zernike[18] = sqrt(12.)*(5*r5-4*r3)*s3theta
    zernike[19] = sqrt(12.)*r5*c5theta
    zernike[20] = sqrt(12.)*r5*s5theta
    # Make zernike zero outside unit circle (useful for dot product)
    zernike = zernike*np.less_equal(r1, 1.0)
    return(zernike[:NumZernike(maxRadial)])
479286
function getShape(x,y,modules,options){

  //Get the surrounding modules matrix
  var mod_matrix = {};
  mod_matrix.topLeft = (x!=0 && y!=0 && modules[x-1][y-1]);
  mod_matrix.top = (y!=0 && modules[x][y-1]);
  mod_matrix.topRight = (x!=modules.length-1 && y!=0 && modules[x+1][y-1]);
  mod_matrix.left = (x!=0 && modules[x-1][y]);
  mod_matrix.right = (x!=modules.length-1 && modules[x+1][y]);
  mod_matrix.bottomLeft = (x!=0 && y!=modules.length-1 && modules[x-1][y+1]);
  mod_matrix.bottom = (y!=modules.length-1 && modules[x][y+1]);
  mod_matrix.bottomRight = (x!=modules.length-1 && y!=modules.length-1 && modules[x+1][y+1]);

  var surroundingCount = 0;
  if(mod_matrix.top){surroundingCount++;}
  if(mod_matrix.left){surroundingCount++;}
  if(mod_matrix.right){surroundingCount++;}
  if(mod_matrix.bottom){surroundingCount++;}

  //Find out box shape from surroundingCount and the orientation of those surrounding modules
  switch(surroundingCount){
    case 0:
      return getShapeZero(x,y,options,modules);
    case 1:
      return getShapeOne(x,y,options,modules,mod_matrix);
    case 2:
      return getShapeTwo(x,y,options,modules,mod_matrix);
    case 3:
      return getShapeThree(x,y,options,modules,mod_matrix);
    case 4:
      return getShapeFour(x,y,options,modules,mod_matrix);
    default:
      return getShapeZero(x,y,options,modules);
  }
      //This function checks that arrays a1 and a2 are the same:
      /*  a1.length==a2.length && a1.every(function(v,i) { return v === a2[i]}) */

}
751698
def _calc_mask(self):
    """"""
    mask = []
    for row in self._constraints:
      mask.append(tuple(x is None for x in row))
    return tuple(mask)
408007
def circular_mask(size):
    
    r = size/2
    c = (r,r)
    y,x = np.ogrid[-c[0]:size-c[0], -c[1]:size-c[1]]
    mask = ~(x*x + y*y <= r*r)
    return mask
683541
def _create_and_rotate_coordinate_arrays(self, x, y, p):
        

        if p.aspect_ratio == 0 or p.size == 0:
            x = x * 0.0
            y = y * 0.0
        else:
            x = (x*10.0) / (p.size*p.aspect_ratio)
            y = (y*10.0) / p.size

        offset = np.exp(p.size)
        pattern_x = np.add.outer(np.sin(p.orientation)*y, np.cos(p.orientation)*x) + offset
        pattern_y = np.subtract.outer(np.cos(p.orientation)*y, np.sin(p.orientation)*x) + offset

        np.clip(pattern_x, 0, np.Infinity, out=pattern_x)
        np.clip(pattern_y, 0, np.Infinity, out=pattern_y)

        return pattern_x, pattern_y
596760
void constructLinearSystem(Planar<I> input , Se3_F32 g ) {
		int numBands = imageType.getNumBands();

		// first precompute everything that does not depend on pixel values
		inboundsPixels = 0;
		for (int i = 0; i < keypixels.size(); i++) {
			Pixel p = keypixels.data[i];

			// Apply the known warp
			SePointOps_F32.transform(g, p.p3, S);

			if( S.z <= 0 ) {
				p.valid = false;
				continue;
			}

			// Compute projected warped pixel coordinate on image I_1
			p.proj.x = (S.x / S.z) * fx + cx;
			p.proj.y = (S.y / S.z) * fy + cy;

			// make sure it's in the bounds
			if (p.proj.x < 0 || p.proj.x > input.width - 1 || p.proj.y < 0 || p.proj.y > input.height - 1) {
				p.valid = false;
				continue;
			} else {
				p.valid = true;
			}
			inboundsPixels++;

			// pi matrix derivative relative to t at S
			float ZZ = S.z * S.z;

			p.dP11 = fx / S.z;
			p.dP13 = -S.x * fx / ZZ;
			p.dP22 = fy / S.z;
			p.dP23 = -S.y * fy / ZZ;
		}

		// how compute the components which require
		errorOptical = 0;
		int row = 0;
		for (int band = 0; band < numBands; band++) {
			interpDX.setImage(derivX.getBand(band));
			interpDY.setImage(derivY.getBand(band));
			interpI.setImage(input.getBand(band));

			for (int i = 0; i < keypixels.size(); i++) {
				Pixel p = keypixels.data[i];

				if( !p.valid )
					continue;

				// Apply the known warp
				// TODO precompute?
				SePointOps_F32.transform(g, p.p3, S);

				// sample pixel values at warped location in I_1
				// NOTE: This could be highly optimized.  Compute and save interpolation weights once per input
				//       instead of for each band and image (current,dx,dy)
				// TODO create a special bilinear class for this?
				float current = interpI.get( p.proj.x, p.proj.y);
				float dx      = interpDX.get(p.proj.x, p.proj.y);
				float dy      = interpDY.get(p.proj.x, p.proj.y);

				// B = grad^T * dPI/dt = shape(1,3)
				float b1 = dx*p.dP11;
				float b2 = dy*p.dP22;
				float b3 = dx*p.dP13 + dy*p.dP23;

				// C * A(S'(x)) = shape(1,6)
				int indexA = row * 6;
				A.data[indexA++] = -b2*S.z + b3*S.y;
				A.data[indexA++] =  b1*S.z - b3*S.x;
				A.data[indexA++] = -b1*S.y + b2*S.x;
				A.data[indexA++] = b1;
				A.data[indexA++] = b2;
				A.data[indexA]   = b3;

				float error = -(current - p.bands[band]);
				y.data[row] = error;

				errorOptical += Math.abs(error);
				row += 1;
			}
		}
		errorOptical /= row;

		A.numRows = row;
		y.numRows = row;
	}
856820
def _create_pairwise_gaussian_2d(sx, sy, shape):
    
    feat_size = 2
    feats = np.zeros((feat_size, shape[0], shape[1]), dtype=np.float32)
    for i in range(shape[0]):
        for j in range(shape[1]):
            feats[0, i, j] = i / sx
            feats[1, i, j] = j / sy
    return feats.reshape([feat_size, -1])
508264
def cv_precompute(self, mask, b):
        '''
        

        '''

        # Get current chunk and mask outliers
        m1 = self.get_masked_chunk(b)
        flux = self.fraw[m1]
        K = GetCovariance(self.kernel, self.kernel_params,
                          self.time[m1], self.fraw_err[m1])
        med = np.nanmedian(flux)

        # Now mask the validation set
        M = lambda x, axis = 0: np.delete(x, mask, axis=axis)
        m2 = M(m1)
        mK = M(M(K, axis=0), axis=1)
        f = M(flux) - med

        # Pre-compute the matrices
        A = [None for i in range(self.pld_order)]
        B = [None for i in range(self.pld_order)]
        for n in range(self.pld_order):
            # Only compute up to the current PLD order
            if self.lam_idx >= n:
                X2 = self.X(n, m2)
                X1 = self.X(n, m1)
                A[n] = np.dot(X2, X2.T)
                B[n] = np.dot(X1, X2.T)
                del X1, X2

        if self.transit_model is None:
            C = 0
        else:
            C = np.zeros((len(m2), len(m2)))
            mean_transit_model = med * \
                np.sum([tm.depth * tm(self.time[m2])
                        for tm in self.transit_model], axis=0)
            f -= mean_transit_model
            for tm in self.transit_model:
                X2 = tm(self.time[m2]).reshape(-1, 1)
                C += tm.var_depth * np.dot(X2, X2.T)
                del X2

        return A, B, C, mK, f, m1, m2
329155
def array_map(f, ar):
    ""
    flat_ar = ravel(ar)
    out = zeros(len(flat_ar), flat_ar.typecode())
    for i in range(len(flat_ar)):
        out[i] = f(flat_ar[i])
    out.shape = ar.shape
    return out
242601
def _compute_projection_filters(G, sf, estimated_source):
    
    # epsilon
    eps = np.finfo(np.float).eps

    # shapes
    (nsampl, nchan) = estimated_source.shape
    # handles the case where we are calling this with only one source
    # G should be nsrc X nsrc X nchan X nchan X filters_len X filters_len
    # and sf should be nsrc X nchan X filters_len
    if len(G.shape) == 4:
        G = G[None, None, ...]
        sf = sf[None, ...]
    nsrc = G.shape[0]
    filters_len = G.shape[-1]

    # zero pad estimates and put chan in first dimension
    estimated_source = _zeropad(estimated_source.T, filters_len - 1, axis=1)

    # compute its FFT
    n_fft = int(2**np.ceil(np.log2(nsampl + filters_len - 1.)))
    sef = scipy.fftpack.fft(estimated_source, n=n_fft)

    # compute the cross-correlations between sources and estimates
    D = np.zeros((nsrc, nchan, filters_len, nchan))
    for (j, cj, c) in itertools.product(
        list(range(nsrc)), list(range(nchan)), list(range(nchan))
    ):
        ssef = sf[j, cj] * np.conj(sef[c])
        ssef = np.real(scipy.fftpack.ifft(ssef))
        D[j, cj, :, c] = np.hstack((ssef[0], ssef[-1:-filters_len:-1]))

    # reshape matrices to build the filters
    D = D.reshape(nsrc * nchan * filters_len, nchan)
    G = _reshape_G(G)

    # Distortion filters
    try:
        C = np.linalg.solve(G + eps*np.eye(G.shape[0]), D).reshape(
            nsrc, nchan, filters_len, nchan
        )
    except np.linalg.linalg.LinAlgError:
        C = np.linalg.lstsq(G, D)[0].reshape(
            nsrc, nchan, filters_len, nchan
        )

    # if we asked for one single reference source,
    # return just a nchan X filters_len matrix
    if nsrc == 1:
        C = C[0]
    return C
1199070
private void getSymmetryAxes(List<Axis> symmAxes, Matrix4d prior, int level, int firstRepeat) {
		if(level >= getNumLevels() ) {
			return;
		}

		Axis elem = axes.get(level);
		Matrix4d elemOp = elem.getOperator();

		// Current axis:
		// elementary maps B -> A
		// prior maps I -> A and J -> B
		// want J -> I = J -> B -> A <- I= inv(prior) * elementary * prior
		Matrix4d currAxisOp = new Matrix4d(prior);
		currAxisOp.invert();
		currAxisOp.mul(elemOp);
		currAxisOp.mul(prior);
		Axis currAxis = new Axis(currAxisOp,elem.getOrder(),elem.getSymmType(),level,firstRepeat);
		symmAxes.add(currAxis);

		//Remember that all degrees are at least 2
		getSymmetryAxes(symmAxes,prior,level+1,firstRepeat);
		//New prior is elementary^d*prior
		Matrix4d newPrior = new Matrix4d(elemOp);
		newPrior.mul(prior);
		int childSize = getNumRepeats(level+1);
		getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize);
		for(int d=2;d<elem.getOrder();d++) {
			newPrior.mul(elemOp,newPrior);
			getSymmetryAxes(symmAxes,newPrior,level+1,firstRepeat+childSize*d);
		}
	}
868671
function createUniformWrapper(gl, program, uniforms, locations) {

  function makeGetter(index) {
    var proc = new Function('gl', 'prog', 'locations', 
      'return function(){return gl.getUniform(prog,locations[' + index + '])}') 
    return proc(gl, program, locations)
  }

  function makePropSetter(path, index, type) {
    switch(type) {
      case 'bool':
      case 'int':
      case 'sampler2D':
      case 'samplerCube':
        return 'gl.uniform1i(locations[' + index + '],obj' + path + ')'
      case 'float':
        return 'gl.uniform1f(locations[' + index + '],obj' + path + ')'
      default:
        var vidx = type.indexOf('vec')
        if(0 <= vidx && vidx <= 1 && type.length === 4 + vidx) {
          var d = type.charCodeAt(type.length-1) - 48
          if(d < 2 || d > 4) {
            throw new Error('gl-shader: Invalid data type')
          }
          switch(type.charAt(0)) {
            case 'b':
            case 'i':
              return 'gl.uniform' + d + 'iv(locations[' + index + '],obj' + path + ')'
            case 'v':
              return 'gl.uniform' + d + 'fv(locations[' + index + '],obj' + path + ')'
            default:
              throw new Error('gl-shader: Unrecognized data type for vector ' + name + ': ' + type)
          }
        } else if(type.indexOf('mat') === 0 && type.length === 4) {
          var d = type.charCodeAt(type.length-1) - 48
          if(d < 2 || d > 4) {
            throw new Error('gl-shader: Invalid uniform dimension type for matrix ' + name + ': ' + type)
          }
          return 'gl.uniformMatrix' + d + 'fv(locations[' + index + '],false,obj' + path + ')'
        } else {
          throw new Error('gl-shader: Unknown uniform data type for ' + name + ': ' + type)
        }
      break
    }
  }

  function enumerateIndices(prefix, type) {
    if(typeof type !== 'object') {
      return [ [prefix, type] ]
    }
    var indices = []
    for(var id in type) {
      var prop = type[id]
      var tprefix = prefix
      if(parseInt(id) + '' === id) {
        tprefix += '[' + id + ']'
      } else {
        tprefix += '.' + id
      }
      if(typeof prop === 'object') {
        indices.push.apply(indices, enumerateIndices(tprefix, prop))
      } else {
        indices.push([tprefix, prop])
      }
    }
    return indices
  }

  function makeSetter(type) {
    var code = [ 'return function updateProperty(obj){' ]
    var indices = enumerateIndices('', type)
    for(var i=0; i<indices.length; ++i) {
      var item = indices[i]
      var path = item[0]
      var idx  = item[1]
      if(locations[idx]) {
        code.push(makePropSetter(path, idx, uniforms[idx].type))
      }
    }
    code.push('return obj}')
    var proc = new Function('gl', 'prog', 'locations', code.join('\n'))
    return proc(gl, program, locations)
  }

  function defaultValue(type) {
    switch(type) {
      case 'bool':
        return false
      case 'int':
      case 'sampler2D':
      case 'samplerCube':
        return 0
      case 'float':
        return 0.0
      default:
        var vidx = type.indexOf('vec')
        if(0 <= vidx && vidx <= 1 && type.length === 4 + vidx) {
          var d = type.charCodeAt(type.length-1) - 48
          if(d < 2 || d > 4) {
            throw new Error('gl-shader: Invalid data type')
          }
          if(type.charAt(0) === 'b') {
            return dup(d, false)
          }
          return dup(d)
        } else if(type.indexOf('mat') === 0 && type.length === 4) {
          var d = type.charCodeAt(type.length-1) - 48
          if(d < 2 || d > 4) {
            throw new Error('gl-shader: Invalid uniform dimension type for matrix ' + name + ': ' + type)
          }
          return dup([d,d])
        } else {
          throw new Error('gl-shader: Unknown uniform data type for ' + name + ': ' + type)
        }
      break
    }
  }

  function storeProperty(obj, prop, type) {
    if(typeof type === 'object') {
      var child = processObject(type)
      Object.defineProperty(obj, prop, {
        get: identity(child),
        set: makeSetter(type),
        enumerable: true,
        configurable: false
      })
    } else {
      if(locations[type]) {
        Object.defineProperty(obj, prop, {
          get: makeGetter(type),
          set: makeSetter(type),
          enumerable: true,
          configurable: false
        })
      } else {
        obj[prop] = defaultValue(uniforms[type].type)
      }
    }
  }

  function processObject(obj) {
    var result
    if(Array.isArray(obj)) {
      result = new Array(obj.length)
      for(var i=0; i<obj.length; ++i) {
        storeProperty(result, i, obj[i])
      }
    } else {
      result = {}
      for(var id in obj) {
        storeProperty(result, id, obj[id])
      }
    }
    return result
  }

  //Return data
  var coallesced = coallesceUniforms(uniforms, true)
  return {
    get: identity(processObject(coallesced)),
    set: makeSetter(coallesced),
    enumerable: true,
    configurable: true
  }
}
392228
public void roundRectangle(float x, float y, float w, float h, float r) {
        if (w < 0) {
            x += w;
            w = -w;
        }
        if (h < 0) {
            y += h;
            h = -h;
        }
        if (r < 0)
            r = -r;
        float b = 0.4477f;
        moveTo(x + r, y);
        lineTo(x + w - r, y);
        curveTo(x + w - r * b, y, x + w, y + r * b, x + w, y + r);
        lineTo(x + w, y + h - r);
        curveTo(x + w, y + h - r * b, x + w - r * b, y + h, x + w - r, y + h);
        lineTo(x + r, y + h);
        curveTo(x + r * b, y + h, x, y + h - r * b, x, y + h - r);
        lineTo(x, y + r);
        curveTo(x, y + r * b, x + r * b, y, x + r, y);
    }
164087
def _reduction_output_shape(x, output_shape, reduced_dim):
  """"""
  if output_shape is None:
    if reduced_dim is None:
      return Shape([])
    else:
      if reduced_dim not in x.shape.dims:
        raise ValueError(
            "reduced_dim=%s not in x.shape.dims=%s" % (reduced_dim, x.shape))
      return x.shape - reduced_dim
  if reduced_dim is not None:
    if [reduced_dim] != [d for d in x.shape.dims if d not in output_shape.dims]:
      raise ValueError(
          "reduced_dim contradicts output_shape:"
          "x=%s output_shape=%s reduced_dim=%s" %
          (x, output_shape, reduced_dim))
  return output_shape