def simulatestep(ps, pm, pv, teilchenort, teilchenmobil, number, vel):
    zzv = np.random.random(number)
    zzv2 = zzv < ps
    zzv3 = zzv < pm
        #print "ss"
   # print teilchen, number
    #zzvVel = np.random.randint(-1,2, size=number)
   # zzvVel2 = np.random.random(number)
    
    #zzvVel3 = (zzvVel2 < pv) * zzvVel
    #print zzvVel
    #berechne neuen Zustand für die Teilchen# bin mobil, wenn
    # vorher mobil und es bleibe (zzv3, pm)
    # oder: war nicht mobil und bleibe nicht (invertiert zu oder)
    mobilneu =  np.bitwise_or(np.bitwise_and(teilchenmobil, zzv3),(np.invert(np.bitwise_or(teilchenmobil, zzv2))))
    
    #vel = vel + zzvVel3 * mobilneu
    #vel = vel * mobilneu
    
    
    #print vel
    
    #np.clip(vel, 0, 12, vel)
    
    #print "vel, ss", vel
    
    # wenn mobil, addiere 1 zum Ort
    teilchenortneu = teilchenort + vel*mobilneu
    
    return teilchenortneu, mobilneu, vel

    def _combine_masks(self):
        """Combine multiple mask planes.

        Sets the detected mask bit if any image has a detection,
        and sets other bits only if set in all images.

        Returns
        -------
        np.ndarray
            The combined mask plane.
        """
        mask_arr = (exp.getMaskedImage().getMask() for exp in self.exposures)

        detected_mask = None
        mask_use = None
        for mask in mask_arr:
            mask_vals = mask.getArray()
            if mask_use is None:
                mask_use = mask_vals
            else:
                mask_use = np.bitwise_and(mask_use, mask_vals)

            detected_bit = mask.getPlaneBitMask('DETECTED')
            if detected_mask is None:
                detected_mask = mask_vals & detected_bit
            else:
                detected_mask = np.bitwise_or(detected_mask, (mask_vals & detected_bit))
        mask_vals_return = np.bitwise_or(mask_use, detected_mask)
        return mask_vals_return

 def place_ship(self, position, length, orientation):
     """
     Return None if ship cannot be placed
     """
     ship = None
     if orientation == 'H':
         zeros = np.zeros(self.width * self.height, dtype='int8')
         if (position[0] + length) > self.width:
             return None
         for i in range(length):
             zeros[position[1] * self.width + position[0]+i] = 1
         if np.all(np.bitwise_and(self._layout, zeros) == 0):
             self._layout = np.bitwise_or(self._layout, zeros)
             ship = Ship(position, length, orientation)
     elif orientation == 'V':
         zeros = np.zeros(self.width * self.height, dtype='int8')
         if (position[1] + length) > self.height:
             return None
         for i in range(length):
             zeros[(position[1] + i) * self.width + position[0]] = 1
         if np.all(np.bitwise_and(self._layout, zeros) == 0):
             self._layout = np.bitwise_or(self._layout, zeros)
             ship = Ship(position, length, orientation)
     if ship:
         self._ships.append(ship)
         return ship

def filter(mask, cube, header, clipMethod, threshold, rmsMode, fluxRange, verbose):
	err.message("Running threshold finder.")
	
	# Sanity checks of user input
	err.ensure(
		clipMethod in {"absolute", "relative"},
		"Threshold finder failed. Illegal clip method: '" + str(clipMethod) + "'.")
	err.ensure(
		rmsMode in {"std", "mad", "gauss", "negative"},
		"Threshold finder failed. Illegal RMS mode: '" + str(rmsMode) + "'.")
	err.ensure(
		fluxRange in {"positive", "negative", "all"},
		"Threshold finder failed. Illegal flux range: '" + str(fluxRange) + "'.")
	
	# Scale threshold by RMS if requested
	if clipMethod == "relative":
		threshold *= GetRMS(cube, rmsMode=rmsMode, fluxRange=fluxRange, zoomx=1, zoomy=1, zoomz=1, verbose=verbose)
	
	# Print some information and check sign of threshold
	err.message("  Using threshold of " + str(threshold) + ".")
	err.ensure(threshold >= 0.0, "Threshold finder failed. Threshold value is negative.")
	
	# Run the threshold finder, setting bit 1 of the mask for |cube| >= |threshold|:
	np.bitwise_or(mask, np.greater_equal(np.absolute(cube), threshold), out=mask)
	
	return

def replace_vals(filename, replace, delim=','):
    '''
    Replace the values in filename with specified values in replace_values
    
    Parameters
    ----------
    filename : string 
        Will be read into a rec array

    replace_values : tuple
        First object is value to replace and second object is what to replace
        it with

    
    '''
    data = csv2rec(filename, delimiter=delim, missing=replace[0])
    for nm in data.dtype.names:
        try:
            # Missing float
            isNaN = (np.isnan(data[nm]))
        except:
            isNaN = np.zeros(len(data[nm]), dtype=bool)
        isBlank = np.array([it == '' for it in data[nm]])
        isMinusOne = (data[nm] == -1)# Missing int
        # Missing other
        isNone = np.array([i == None for i in data[nm]])
        ind = np.bitwise_or(isNaN, isBlank)
        ind = np.bitwise_or(ind, isMinusOne)
        ind = np.bitwise_or(ind, isNone)
        data[nm][ind] = replace[1]
    return data

def reset_bad_gain(pdq, gain):
    """
    For pixels in the gain array that are either non-positive or NaN, reset the
    the corresponding pixels in the pixel DQ array to NO_GAIN_VALUE and
    DO_NOT_USE so that they will be ignored.

    Parameters
    ----------
    pdq : int, 2D array
        pixel dq array of input model

    gain : float32, 2D array
        gain array from reference file

    Returns
    -------
    pdq : int, 2D array
        pixleldq array of input model, reset to NO_GAIN_VALUE and DO_NOT_USE
        for pixels in the gain array that are either non-positive or NaN.
    """
    wh_g = np.where( gain <= 0.)
    if len(wh_g[0]) > 0:
        pdq[wh_g] = np.bitwise_or( pdq[wh_g], dqflags.pixel['NO_GAIN_VALUE'] )
        pdq[wh_g] = np.bitwise_or( pdq[wh_g], dqflags.pixel['DO_NOT_USE'] )

    wh_g = np.where( np.isnan( gain ))
    if len(wh_g[0]) > 0:
        pdq[wh_g] = np.bitwise_or( pdq[wh_g], dqflags.pixel['NO_GAIN_VALUE'] )
        pdq[wh_g] = np.bitwise_or( pdq[wh_g], dqflags.pixel['DO_NOT_USE'] )

    return pdq

    def read_ports(self):
        key = psychopy.event.getKeys(keyList=['1','2','3','k'])
        ports = np.asarray([False,False,False])
        if key:
            if not key[0] in self._key_pressed: self._key_pressed.append(key[0])
        if 'k' in self._key_pressed and '1' in self._key_pressed:
            ports = np.bitwise_or(ports,[True,False,False])
            psychopy.event.clearEvents()
            print(self._key_pressed)
            self._key_pressed.remove('k')
            self._key_pressed.remove('1')
        if 'k' in self._key_pressed and '2' in self._key_pressed:
            ports = np.bitwise_or(ports,[False,True,False])
            psychopy.event.clearEvents()
            print(self._key_pressed)
            self._key_pressed.remove('k')
            self._key_pressed.remove('2')
        if 'k' in self._key_pressed and '3' in self._key_pressed:
            ports = np.bitwise_or(ports,[False,False,True])
            psychopy.event.clearEvents()
            print(self._key_pressed)
            self._key_pressed.remove('k')
            self._key_pressed.remove('3')

        return self.get_ports()[ports]

    def __invertibleToRGB(self, rgbVarr, varrIdx, colorLutStruct):
        '''
        Decode an RGB image rendered in INVERTIBLE_LUT mode. The image encodes
        float values as colors, so the RGB value is first decoded into its
        represented float value and then the color table is applied.
        '''
        w0 = np.left_shift(
            rgbVarr[varrIdx[0], varrIdx[1], 0].astype(np.uint32), 16)
        w1 = np.left_shift(
            rgbVarr[varrIdx[0], varrIdx[1], 1].astype(np.uint32), 8)
        w2 = rgbVarr[varrIdx[0], varrIdx[1], 2]

        value = np.bitwise_or(w0, w1)
        value = np.bitwise_or(value, w2)
        value = np.subtract(value.astype(np.int32), 1)
        normalized_val = np.divide(value.astype(float), 0xFFFFFE)

        # Map float value to color lut (use a histogram to support non-uniform
        # colormaps (fetch bin indices))
        colorLut = colorLutStruct.lut
        bins = colorLutStruct.adjustedBins

        idx = np.digitize(normalized_val, bins)
        idx = np.subtract(idx, 1)

        valueImage = np.zeros([rgbVarr.shape[0], rgbVarr.shape[1]],
                              dtype=np.uint32)
        valueImage[varrIdx[0], varrIdx[1]] = idx

        return colorLut[valueImage]

def jaccard_distance(u, v):
    """return jaccard distance"""
    u = numpy.asarray(u)
    v = numpy.asarray(v)
    return (numpy.double(numpy.bitwise_and((u != v),
            numpy.bitwise_or(u != 0, v != 0)).sum())
            /  numpy.double(numpy.bitwise_or(u != 0, v != 0).sum()))

    def pattern_distance_jaccard(a, b):
        """
        Computes a distance measure for two binary patterns based on their
        Jaccard-Needham distance, defined as

        .. math::

            d_J(a,b) = 1 - J(a,b) = \\frac{|a \\cup b| - |a \\cap b|}{|a \\cup b|}.

        The similarity measure takes values on the closed interval [0, 1],
        where a value of 1 is attained for disjoint, i.e. maximally dissimilar
        patterns a and b and a value of 0 for the case of :math:`a=b`.

        Parameters
        ----------
        a : list or array, int or bool
            Input pattern
        b : list or array, int or bool
            Input pattern

        Returns
        -------
        dist : double
            Jaccard distance between `a` and `b`.
        """
        # Note: code taken from scipy. Duplicated as only numpy references wanted for base functionality
        a = np.atleast_1d(a).astype(bool)
        b = np.atleast_1d(b).astype(bool)
        dist = (np.double(np.bitwise_and((a != b), np.bitwise_or(a != 0, b != 0)).sum())
                / np.double(np.bitwise_or(a != 0, b != 0).sum()))
        return dist

def partBy2(x):
    x = np.bitwise_and(x, 0x1f00000000ffff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 32)), 0x1f00000000ffff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 16)), 0x1f0000ff0000ff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 8)), 0x100f00f00f00f00f)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 4)), 0x10c30c30c30c30c3)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 2)), 0x1249249249249249)
    return x

def mas_combining4(row):
    rows = ['DER_mass_MMC','DER_mass_vis','DER_mass_transverse_met_lep']
    frame = row[rows[0]] + row[rows[1]] + row[rows[2]]
    mask = np.bitwise_or(row[rows[0]]==-999, row[rows[1]]==-999.)
    mask = np.bitwise_or(mask, row[rows[1]]==-999.)
    frame[mask] = -999.
    frame.name = 'MNew_mc4'
    return frame

def get_keyv(iarr, level):
    i1, i2, i3 = (v.astype("int64") for v in iarr)
    i1 = spread_bitsv(i1, level)
    i2 = spread_bitsv(i2, level) << 1
    i3 = spread_bitsv(i3, level) << 2
    np.bitwise_or(i1, i2, i1)
    np.bitwise_or(i1, i3, i1)
    return i1

 def __call__(self):
     if hasattr(self.ctx.params, 'cleaner') and self.ctx.params.cleaner != "None":
         #load module
         mod = importlib.import_module(self.ctx.params.cleaner)        
 
         rfi_mask_vx, rfi_mask_vy = mod.rm_rfi(self.ctx)
         
         self.ctx.tod_vx.mask = np.bitwise_or(self.ctx.tod_vx.mask, rfi_mask_vx)
         self.ctx.tod_vy.mask = np.bitwise_or(self.ctx.tod_vy.mask, rfi_mask_vy)

 def colRowIsOnSciencePixelList(self, col, row, padding=DEFAULT_PADDING):
     """similar to colRowIsOnSciencePixelList() but takes lists as input"""
     out = np.ones(len(col), dtype=bool)
     col_arr = np.array(col)
     row_arr = np.array(row)
     mask = np.bitwise_or(col_arr < 12. - padding, col_arr > 1111 + padding)
     out[mask] = False
     mask = np.bitwise_or(row_arr < 20. - padding, row_arr > 1043 + padding)
     out[mask] = False
     return out

    def _numpy(self, data, weights, shape):
        q = self.quantity(data)
        self._checkNPQuantity(q, shape)
        self._checkNPWeights(weights, shape)
        weights = self._makeNPWeights(weights, shape)
        newentries = weights.sum()

        import numpy

        selection = numpy.isnan(q)
        numpy.bitwise_not(selection, selection)
        subweights = weights.copy()
        subweights[selection] = 0.0
        self.nanflow._numpy(data, subweights, shape)

        # avoid nan warning in calculations by flinging the nans elsewhere
        numpy.bitwise_not(selection, selection)
        q = numpy.array(q, dtype=numpy.float64)
        q[selection] = 0.0
        weights = weights.copy()
        weights[selection] = 0.0

        if all(isinstance(v, Count) and v.transform is identity for c, v in self.bins) and numpy.all(numpy.isfinite(q)) and numpy.all(numpy.isfinite(weights)):

            h, _ = numpy.histogram(q, [float("-inf")] + [(c1 + c2)/2.0 for (c1, v1), (c2, v2) in zip(self.bins[:-1], self.bins[1:])] + [float("inf")], weights=weights)

            for hi, (c, v) in zip(h, self.bins):
                v.fill(None, float(hi))

        else:
            selection = numpy.empty(q.shape, dtype=numpy.bool)
            selection2 = numpy.empty(q.shape, dtype=numpy.bool)

            for index in xrange(len(self.bins)):
                if index == 0:
                    high = (self.bins[index][0] + self.bins[index + 1][0])/2.0
                    numpy.greater_equal(q, high, selection)

                elif index == len(self.bins) - 1:
                    low = (self.bins[index - 1][0] + self.bins[index][0])/2.0
                    numpy.less(q, low, selection)

                else:
                    low = (self.bins[index - 1][0] + self.bins[index][0])/2.0
                    high = (self.bins[index][0] + self.bins[index + 1][0])/2.0
                    numpy.less(q, low, selection)
                    numpy.greater_equal(q, high, selection2)
                    numpy.bitwise_or(selection, selection2, selection)

                subweights[:] = weights
                subweights[selection] = 0.0
                self.bins[index][1]._numpy(data, subweights, shape)

        # no possibility of exception from here on out (for rollback)
        self.entries += float(newentries)

def color(img_on, img_off):
    img_off = cv2.imread(img_off)
    img_on = cv2.imread(img_on)

    img = np.array((np.array(img_on, dtype=np.int16)-np.array(img_off, dtype=np.int16)).clip(0,255), dtype=np.uint8)
    img = cv2.medianBlur(img,3)
    cv2.imshow("soustracted",cv2.resize(img, (img.shape[1]/2,img.shape[0]/2)))
    cv2.waitKey(0)

    lower = np.array([0, 0, 20], dtype=np.uint8)
    upper = np.array([5, 5, 255], dtype=np.uint8)
    mask0 = cv2.inRange(img, lower, upper)
    cv2.imshow("color_threshold_red",cv2.resize(mask0, (img.shape[1]/2,img.shape[0]/2)))
    cv2.waitKey(0)

    img_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

    lower = np.array([0, 170, 20], dtype=np.uint8)
    upper = np.array([10, 255, 255], dtype=np.uint8)
    mask1 = cv2.inRange(img_hsv, lower, upper)

    lower = np.array([170, 170, 20], dtype=np.uint8)
    upper = np.array([180, 255, 255], dtype=np.uint8)
    mask2 = cv2.inRange(img_hsv, lower, upper)

    mask = np.bitwise_or(mask1, mask2)
    cv2.imshow("hsv_threshold_low_brightness",cv2.resize(mask, (mask.shape[1]/2,mask.shape[0]/2)))
    cv2.waitKey(0)

    lower = np.array([80, 0, 200], dtype=np.uint8)
    upper = np.array([100, 255, 255], dtype=np.uint8)
    mask3 = cv2.inRange(img_hsv, lower, upper)
    cv2.imshow("hsv_threshold_high_brightness",cv2.resize(mask3, (mask3.shape[1]/2,mask3.shape[0]/2)))
    cv2.waitKey(0)

    mask = np.bitwise_or(mask, mask3)
    mask = np.bitwise_or(mask0, mask)
    mask = cv2.GaussianBlur(mask,(3,3),0)
    mask = cv2.inRange(mask, np.array([250]), np.array([255]))

    res = cv2.bitwise_and(img_on, img_on, mask=mask)
    tmp = np.zeros(res.shape)

    for line in range(res.shape[0]):
        moments = cv2.moments(res[line,:,2])
        if(moments['m00'] != 0):
            tmp[line][round(moments['m01']/moments['m00'])] = [0,255,0]


    cv2.imshow("hsv_or_color",cv2.resize(mask, (mask.shape[1]/2,mask.shape[0]/2)))
    cv2.waitKey(0)
    cv2.imshow("result",cv2.resize(tmp, (res.shape[1]/2,res.shape[0]/2)))
    cv2.waitKey(0)

def gameoflife(W,H,ITER,DIST,random_state):
    random.setstate(random_state)
    LIVING_LOW= 2
    LIVING_HIGH = 3
    ALIVE = 3

    full      = np.zeros((W+2,H+2), dtype=np.long, dist=DIST)
    dead      = np.zeros((W,H),     dtype=np.long, dist=DIST)
    live      = np.zeros((W,H),     dtype=np.long, dist=DIST)
    live2     = np.zeros((W,H),     dtype=np.long, dist=DIST)
    neighbors = np.zeros((W,H),     dtype=np.long, dist=DIST)

    cells = full[1:W+1,1:H+1]
    ul = full[0:W, 0:H]
    um = full[0:W, 1:H+1]
    ur = full[0:W, 2:H+2]
    ml = full[1:W+1, 0:H]
    mr = full[1:W+1, 2:H+2]
    ll = full[2:W+2, 0:H]
    lm = full[2:W+2, 1:H+1]
    lr = full[2:W+2, 2:H+2]

    for i in xrange(W):
      for j in range(H):
          if random.random() > .8:
              cells[i][j] = 1
    for i in xrange(ITER):
        # zero neighbors
        np.bitwise_and(neighbors,0,neighbors)
        # count neighbors
        neighbors += ul
        neighbors += um
        neighbors += ur
        neighbors += ml
        neighbors += mr
        neighbors += ll
        neighbors += lm
        neighbors += lr
        # extract live cells neighbors
        np.multiply(neighbors, cells, live)
        # find all living cells among the already living
        np.equal(live, LIVING_LOW, live2)
        np.equal(live, LIVING_HIGH, live)
        # merge living cells into 'live'
        np.bitwise_or(live, live2, live)
        # extract dead cell neighbors
        np.equal(cells, 0, dead)
        dead *= neighbors
        np.equal(dead,ALIVE,dead)
        # make sure all threads have read their values
        np.bitwise_or(live, dead, cells)
    return full

def simulatestep(p_in, teilchenort, teilchenmobil, number, vel):
    psvektor = np.array([ps] * number)
    zzv = np.random.random(number)
    #bleibe stationär
    zzv2 = zzv < psvektor
    #bleibe mobil
    zzv3 = zzv < p_in
    #print "pin", p_in
    
    '''print "vektoren"
    print p_in
    print zzv
    print zzv2
    print zzv3'''
        
    #berechne neuen Zustand für die Teilchen
    # bin mobil, wenn
    # vorher mobil und es bleibe (zzv3, bestimmt durch p_in, variable WKeiten)
    # oder: nicht mobil und nicht bleibe (zzv2, bestimmt durch gammavektor, fest TODO: geht das geschickter?)
    mobilneu =  np.bitwise_or(np.bitwise_and(teilchenmobil, zzv3),(np.invert(np.bitwise_or(teilchenmobil, zzv2)))) 
    # wenn mobil, addiere vel zum Ort
    teilchenortneu = teilchenort + vel*mobilneu
    
    #print "mobneu", mobilneu
    
    vel = np.clip(vel, velmin, velmax)
    # berechne neue Geschwindigkeiten; falls mobil addiere was dazu, clip davor und danach zum Schutz
    velneu = (vel + ((velmax-vel)/veldivisor)) * mobilneu # wenn nicht mobil, schneller, wenn mobil
    #print "velneu", velneu
    #vel ist mindestens 1, für den neuen Ort wird mobil ja eh noch mal befragt
    velneu = np.clip(velneu, 0, velmax)
    #print max(velneu)
    
    #print "p_in", p_in
    #time.sleep(1)
    
    # wo ich nicht mobil bin: gamma0 -> p1: [0, g0, 0,0,0, g0, g0, 0...]
    p1 = np.invert(mobilneu) * gamma0
    
    #wo ich mobil bin, erhöht sich p TODO p2: [p, 0, p,p,p,0,0,p,...]
    p2 = mobilneu * (gamma0 + vel/gammadivisor)
    
    #print "veldings", 1/vel
    #print gamma0/gammadivisor
    #print "p2", p2
    p_out = p1 + p2
    #print "pout", p_out
    p_out = np.clip(p_out, gamma0, 0.9999)
    #print velneu,'\n',  p_out

    return teilchenortneu, mobilneu, velneu, p_out

def apply_flat_field (science, flat ):
    """
    Short Summary
    -------------
    Flat fields the data and error arrays, and updates data quality array
    based on bad pixels in flat field arrays. Applies portion of flat field
    corresponding to science image subarray.

    Parameters
    ----------
    science: JWST data model
        input science data model

    flat: JWST data model
        flat field data model

    Returns
    -------
    None
    """

    # If the input science data model is a subarray, extract the same
    # subarray from the flatfield model
    if ref_matches_sci (flat, science):
        flat_data = flat.data
        flat_dq   = flat.dq
    else: 
        log.info("Extracting matching subarray from flat")
        flat_data = get_subarray(flat.data, science)
        flat_dq   = get_subarray(flat.dq, science)

    # For pixels whose flat is either NaN or NO_FLAT_FIELD, update their DQ to
    # indicate that no flat is applied to those pixels
    flat_dq[np.isnan(flat_data)] = np.bitwise_or(flat_dq[np.isnan(flat_data)],
                                                 dqflags.pixel['NO_FLAT_FIELD'])

    # Replace NaN's in flat with 1's
    flat_data[np.isnan(flat_data)] = 1.0 

    # Reset flat values of pixels having DQ values containing NO_FLAT_FIELD
    # to 1.0, so that no flat fielding correction is made 
    wh_dq = np.bitwise_and( flat_dq, dqflags.pixel['NO_FLAT_FIELD'])
    flat_data[ wh_dq == dqflags.pixel['NO_FLAT_FIELD'] ] = 1.0  

    # Flatten data and error arrays
    science.data /= flat_data
    science.err  /= flat_data

    # Combine the science and flat DQ arrays
    science.dq = np.bitwise_or(science.dq, flat_dq)