def visualize2DOF(pred_act1,pred_act2,act,num_bins=10):
    bin_size1 = (numpy.max(pred_act1,axis=0) - numpy.min(pred_act1,axis=0))/num_bins 
    bin_size2 = (numpy.max(pred_act2,axis=0) - numpy.min(pred_act2,axis=0))/num_bins
    	
    of = numpy.zeros((numpy.shape(act)[1],num_bins,num_bins))
    ofn = numpy.zeros((numpy.shape(act)[1],num_bins,num_bins))
	
    for i in xrange(0,numpy.shape(act)[0]):
	idx1 = numpy.round_((pred_act1[i,:]-numpy.min(pred_act1,axis=0)) / bin_size1)    	
	idx2 = numpy.round_((pred_act2[i,:]-numpy.min(pred_act2,axis=0)) / bin_size2)
	
	idx1 = idx1 -(idx1 >= num_bins)
	idx2 = idx2 -(idx2 >= num_bins)
	
	j=0
	for (x,y) in zip(numpy.array(idx1).flatten().tolist(),numpy.array(idx2).flatten().tolist()):
            of[j,x,y] = of[j,x,y] +  act[i,j]
	    ofn[j,x,y] = ofn[j,x,y] + 1 
	    j=j+1
    
    of = of - (ofn <= 0)
    ofn = ofn + (ofn <= 0)
    of = of/ofn
    print of[0]
    print of[1]
    print ofn[0]
    print ofn[1]

    
    
    showRFS(of,joinnormalize=False)

    def function(self, simulation, period):
        period = period.start.offset('first-of', 'month').period('month')
        rfr = simulation.calculate('rfr', period.start.offset('first-of', 'year').period('year').offset(-2))
        age_holder = simulation.compute('age', period)
        scolarite_holder = simulation.compute('scolarite', period)
        P = simulation.legislation_at(period.start).bourses_education.bourse_college

        ages = self.split_by_roles(age_holder, roles = ENFS)
        nb_enfants = zeros(len(rfr))
        for age in ages.itervalues():
            nb_enfants += age >= 0

        plafond_taux_1 = round_(P.plafond_taux_1 + P.plafond_taux_1 * nb_enfants * P.coeff_enfant_supplementaire)
        plafond_taux_2 = round_(P.plafond_taux_2 + P.plafond_taux_2 * nb_enfants * P.coeff_enfant_supplementaire)
        plafond_taux_3 = round_(P.plafond_taux_3 + P.plafond_taux_3 * nb_enfants * P.coeff_enfant_supplementaire)

        eligible_taux_3 = rfr <= plafond_taux_3
        eligible_taux_2 = not_(eligible_taux_3) * (rfr <= plafond_taux_2)
        eligible_taux_1 = not_(or_(eligible_taux_2, eligible_taux_3)) * (rfr <= plafond_taux_1)

        scolarites = self.split_by_roles(scolarite_holder, roles = ENFS)
        nb_enfants_college = zeros(len(rfr))
        for scolarite in scolarites.itervalues():
            nb_enfants_college += scolarite == SCOLARITE_COLLEGE

        montant = nb_enfants_college * (
            eligible_taux_3 * P.montant_taux_3 +
            eligible_taux_2 * P.montant_taux_2 +
            eligible_taux_1 * P.montant_taux_1
            )

        return period, montant / 12

    def function(self, simulation, period):
        period = period.this_month
        rfr = simulation.calculate('rfr', period.n_2)
        age_holder = simulation.compute('age', period)
        scolarite_holder = simulation.compute('scolarite', period)
        P = simulation.legislation_at(period.start).bourses_education.bourse_college

        ages = self.split_by_roles(age_holder, roles = ENFS)
        nb_enfants = sum(
            age >= 0 for age in ages.itervalues()
        )

        scolarites = self.split_by_roles(scolarite_holder, roles = ENFS)

        nb_enfants_college = sum(
            scolarite == SCOLARITE_COLLEGE for scolarite in scolarites.itervalues()
        )

        montant_par_enfant = apply_thresholds(
            rfr,
            thresholds = [
                # plafond_taux_3 est le plus bas
                round_(P.plafond_taux_3 + P.plafond_taux_3 * nb_enfants * P.coeff_enfant_supplementaire),
                round_(P.plafond_taux_2 + P.plafond_taux_2 * nb_enfants * P.coeff_enfant_supplementaire),
                round_(P.plafond_taux_1 + P.plafond_taux_1 * nb_enfants * P.coeff_enfant_supplementaire),
                ],
            choices = [P.montant_taux_3, P.montant_taux_2, P.montant_taux_1]
            )

        montant = nb_enfants_college * montant_par_enfant

        return period, montant / 12

        def loyer_retenu():

            # loyer mensuel réel, multiplié par 2/3 pour les meublés
            L1 = round_(loyer * where(statut_occupation == 5, 2 / 3, 1))

            zone_apl = simulation.calculate('zone_apl_famille', period)

            # Paramètres contenant les plafonds de loyer pour cette zone
            plafonds_by_zone = [[0] + [al.loyers_plafond[ 'zone' + str(zone) ][ 'L' + str(i) ] for zone in range(1, 4)] for i in range(1, 5)]
            L2_personne_seule = take(plafonds_by_zone[0], zone_apl)
            L2_couple = take(plafonds_by_zone[1], zone_apl)
            L2_famille = take(plafonds_by_zone[2], zone_apl) + (al_pac > 1) * (al_pac - 1) * take(plafonds_by_zone[3], zone_apl)

            L2 = select(
                [personne_seule * (al_pac == 0) + chambre, al_pac > 0],
                [L2_personne_seule, L2_famille],
                default = L2_couple
                )

            # taux à appliquer sur le loyer plafond
            coeff_chambre_colloc = select(
                [chambre, coloc],
                [al.loyers_plafond.chambre, al.loyers_plafond.colocation],
                default = 1)

            L2 = round_(L2 * coeff_chambre_colloc, 2)

            # loyer retenu
            L = min_(L1, L2)

            return L

    def indices(self, x, y, clip=False):
        """
        Return the grid pixel indices (i_x, i_y) corresponding to the
        given arrays of grid coordinates. Arrays x and y must have the
        same size. Also return a boolean array of the same length that
        is True where the pixels are within the grid bounds and False
        elsewhere.
		
        If clip is False, a ValueError is raised if any of the pixel
        centers are outside the grid bounds, and array within will be
        all True. If clip is True, then the i_x and i_y values where
        within is False will be nonsense; the safe thing is to use
        only i_x[within] and i_y[within].
        """
        if x.size != y.size:
            raise ValueError("Arrays x and y must have the same length.")
        # This is a workaround for the behavior of int_: when given an
        # array of size 1 it returns an int instead of an array.
        if x.size == 1:
            i_x = np.array([np.int(np.round((x[0] - self.x[0]) / self.dx()))])
            i_y = np.array([np.int(np.round((y[0] - self.y[0]) / self.dy()))])
        else:
            i_x = np.int_(np.round_((x - self.x[0]) / self.dx()))
            i_y = np.int_(np.round_((y - self.y[0]) / self.dy()))
        within = ((0 <= i_x) & (i_x < self.x.size) & (0 <= i_y) & (i_y < self.y.size))
        if not clip and not all(within):
            raise ValueError("Not all points are inside the grid bounds, and clipping is not allowed.")
        return i_x, i_y, within

    def _run_(self):
        '''
        '''
        #print '- Running Mjpeg Decoder...'
        hf = h.HuffCoDec(self.hufftables)
        r, c, chnl = self.R, self.C, self.NCHNL
        Z = self.Z
        
        #hufcd = self.huffcodes#self.fl.readline()[:-1]
        if self.mode == '444':
            for ch in range(chnl):                #hufcd = self.fl.readline()[:-1]            #    print hufcd[0:20]
                nblk, seqrec = hf.invhuff(self.huffcodes[ch], ch)
                for i in range(self.nBlkRows):
                    for j in range(self.nBlkCols):
                        blk = h.zagzig(seqrec[i*self.nBlkCols + j])
                        self.imRaw[r*i:r*i+r, c*j:c*j+c, ch] = np.round_( cv2.idct( blk*Z[:,:,ch] ))
                        
        elif self.mode == '420':
            #import math as m
            if chnl == 1:
                rYmg = self.imRaw
            else:                #Y = self.imRaw[:,:,0]
                Y = np.zeros( (self.M, self.N) )
                dims, CrCb = h.adjImg( downsample(np.zeros( (self.M, self.N, 2) ), self.mode)[1] )
                rYmg = [ Y, CrCb[:,:,0], CrCb[:,:,1] ]
                
            for ch in range(chnl):
                #hufcd = self.fl.readline()[:-1]
                if ch == 0:
                    rBLK = self.nBlkRows
                    cBLK = self.nBlkCols
                else:
                    rBLK, cBLK = int(np.floor(dims[0]/self.R)), int(np.floor(dims[1]/self.C))
            #    print hufcd[0:20]
                nblk, self.seqrec = hf.invhuff(self.huffcodes[ch], ch)
                for i in range(rBLK):
                    for j in range(cBLK):
                        blk = h.zagzig(self.seqrec[i*cBLK + j])
                        #print rYmg[ch][r*i:r*i+r, c*j:c*j+c].shape, ch, i, j
                        rYmg[ch][r*i:r*i+r, c*j:c*j+c] = np.round_( cv2.idct( blk*Z[:,:,ch] ))
            # UPSAMPLE
            if chnl == 1:
                self.imRaw = rYmg #[:self.Mo, : self.No]
            else:
                self.imRaw[:,:,0] = rYmg[0]
                self.imRaw[:,:,1] = upsample(rYmg[1], self.mode)[:self.M, :self.N]
                self.imRaw[:,:,2] = upsample(rYmg[2], self.mode)[:self.M, :self.N]
        
        #self.fl.close()
        
#        imrec = cv2.cvtColor((self.imRaw[:self.Mo, :self.No]+128), cv2.COLOR_YCR_CB2BGR)
#        imrec = self.imRaw[:self.Mo, :self.No]+128
        imrec = self.imRaw+128.0
#        imrec[imrec>255.0]=255.0
#        imrec[imrec<0.0]=0.0
								
        #print 'Mjpeg Decoder Complete...'
        
        return imrec

def rotz(ang):
    """Generate a homogenous trransform for ang radians around the z axis"""
    s = N.round_(sin(ang), decimals=14); c = N.round_(cos(ang), decimals=14)
    return N.array([
        [c,-s, 0, 0],
        [s, c, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ])

def is_equidistant(x):
    '''
    >>> is_equidistant((0,1,2))
    True
    >>> is_equidistant((0,1,2.5))
    False
    '''
    d = np.diff(x)
    return (np.round_(d,8)==np.round_(d[0],8)).all()

def process_histogram(PabsFlip, N1, uCut, lCut, angleInc, radStep):
    """
    Create orientation Histogram
    Sum pixel intensity along different angles
    :param PabsFlip:
    :param N1:
    :param uCut: upper-cut parameter from the settings.SettingsWindow
    :param lCut: lower-cut parameter form the settings.SettingsWindow
    :param angleInc: angle-increment from the
    :param radStep: radial-step
    :return:
    """
    n1 = np.round(N1 / 2) - 1
    freq = np.arange(-n1, n1 + 1, 1)
    x, y = freq, freq

    # Variables for settings
    CO_lower = lCut
    CO_upper = uCut
    angleInc = angleInc
    radStep = radStep

    #  Set up polar coordinates prior to summing the spectrum
    theta1Rad = np.linspace(0.0, 2 * math.pi, num=360/angleInc)
    f1 = np.round_(N1 / (2 * CO_lower))
    f2 = np.round_(N1 / (2 * CO_upper))

    rho1 = np.linspace(f1, f2, num=(f2 - f1)/radStep)  # frequency band
    PowerX = np.zeros((theta1Rad.size, theta1Rad.size))
    PowerY = np.zeros((theta1Rad.size))

    # Interpolate using a Spine
    PowerSpline = scipy.interpolate.RectBivariateSpline(y=y, x=x, z=PabsFlip)
    n_dx = 0.001

    for p in range(0, theta1Rad.size):
        # converting theta1Rad and rho1 to cartesian coordinates
        xfinal = rho1 * math.cos(theta1Rad[p])
        yfinal = rho1 * math.sin(theta1Rad[p])

        # Evaluate spin on path
        px = PowerSpline.ev(yfinal, xfinal)
        PowerY[p] = np.sum(px)

    # Only use the data in the first two quadrants (Spectrum is symmetric)
    num = len(theta1Rad)
    PowerYFinal = PowerY[0:num // 2]
    theta1RadFinal = theta1Rad[0:num // 2]

    power_area = np.trapz(PowerYFinal, theta1RadFinal)
    normPower = PowerYFinal / power_area

    # TODO: Ask Rici what those are
    return normPower, theta1RadFinal

 def _binary_preds(self, model_preds, mean_preds, stack_preds):
     """
     """
     stack_preds_bin = []
     model_preds_bin = np.round_(model_preds, decimals=0)
     mean_preds_bin = np.round_(mean_preds, decimals=0)
     stack_preds_bin = np.round_(stack_preds, decimals=0) \
                     if self.stack else 0
             
     return model_preds_bin, mean_preds_bin, stack_preds_bin
     

def restarize_events(events, durations, dt, t_max):
    """ build a binary sequence of events. Each event start is approximated
    to the nearest time point on the time grid defined by dt and t_max.
    """
    smpl_events = np.array(np.round_(np.divide(events, dt)), dtype=int)
    smpl_durations = np.array(np.round_(np.divide(durations, dt)), dtype=int)
    smpl_events = extend_sampled_events(smpl_events, smpl_durations)
    if np.allclose(t_max % dt, 0):
        bin_seq = np.zeros(int(t_max / dt) + 1)
    else:
        bin_seq = np.zeros(int(np.round((t_max + dt) / dt)))
    bin_seq[smpl_events] = 1

    return bin_seq

def create_curve_states_2D():
    R90CW,FH,FV,R90CW_FH,R90CW_FV,R180CW = get_standard_matrixes()
    state0m = np.matrix([[0, 0, 0, 1],
        [0, 1, 0, 1],
        [1, 1, 0, 1],
        [1, 0, 0, 1]])


    state0 = np.round_(state0m.getT()[:2:].getT().astype(int)).astype(int)
    state1 = np.round_(R90CW_FH.dot(state0m.getT())[:2:].getT()).astype(int)
    state2 = np.round_(R90CW_FV.dot(state0m.getT())[:2:].getT()).astype(int)
    state3 = np.round_(R180CW.dot(state0m.getT())[:2:].getT()).astype(int)

    return np.array([state0,state1, \
    state2, state3])

def general_axis_rotation(axis,  angle):
    """Generates a rotation matrix around <axis> by <angle>, using the right-hand
    rule.
    Arguments: 
        axis - a 3-component 1D array representing a unit vector
        angle - rotation counterclockwise in radians around the axis when the axis 
            points to the viewer.
    Returns: A 3x3 array representing the matrix of rotation.
    Reference: [1] p.47
    """
    N.round_(sin(ang), decimals=14); c = N.round_(cos(ang), decimals=14); v = 1 - c
    add = N.array([[0,          -axis[2], axis[1]],  
                            [axis[2],  0,          -axis[0]], 
                            [-axis[1], axis[0],  0        ] ])
    return N.multiply.outer(axis,  axis)*v + N.eye(3)*c + add*s

    def test_round_(self):
        self.assertQuantityEqual(
            np.round_([.5, 1.5, 2.5, 3.5, 4.5] * pq.J),
            [0., 2., 2., 4., 4.] * pq.J
            )

        self.assertQuantityEqual(
            np.round_([1,2,3,11] * pq.J, decimals=1),
            [1, 2, 3, 11] * pq.J
            )

        self.assertQuantityEqual(
            np.round_([1,2,3,11] * pq.J, decimals=-1),
            [0, 0, 0, 10] * pq.J
            )

def log_der_13(z, nstop):
    '''
    Calculate logarithmic derivatives of Riccati-Bessel functions psi
    and xi for complex arguments.  Riccati-Bessel conventions follow
    Bohren & Huffman.

    See Mackowski et al., Applied Optics 29, 1555 (1990).

    Parameters
    ----------
    z: complex number
    nstop: maximum order of computation
    '''
    z = np.complex128(z) # convert to double precision

    # Calculate Dn_1 (based on \psi(z)) using downward recursion.
    # See Mackowski eqn. 62
    nmx = np.maximum(nstop, int(np.round_(np.absolute(z)))) + 15
    dn1 = log_der_1(z, nmx, nstop)

    # Calculate Dn_3 (based on \xi) by up recurrence
    # initialize
    dn3 = zeros(nstop+1, dtype = 'complex128')
    psixi = zeros(nstop+1, dtype = 'complex128')
    dn3[0] = 1.j
    psixi[0] = -1j*exp(1.j*z)*sin(z)
    for dindex in arange(1, nstop+1):
        # Mackowski eqn 63
        psixi[dindex] = psixi[dindex-1] * ( (dindex/z) - dn1[dindex-1]) * (
            (dindex/z) - dn3[dindex-1])
        # Mackowski eqn 64
        dn3[dindex] = dn1[dindex] + 1j/psixi[dindex]

    return dn1, dn3

def compute_allegement_fillon(simulation, period):
    """
        Exonération Fillon
        http://www.securite-sociale.fr/comprendre/dossiers/exocotisations/exoenvigueur/fillon.htm
    """
    assiette = simulation.calculate_add('assiette_allegement', period)
    smic_proratise = simulation.calculate_add('smic_proratise', period)
    taille_entreprise = simulation.calculate('taille_entreprise', period)
    majoration = (taille_entreprise <= 2)  # majoration éventuelle pour les petites entreprises
    # Calcul du taux
    # Le montant maximum de l’allègement dépend de l’effectif de l’entreprise.
    # Le montant est calculé chaque année civile, pour chaque salarié ;
    # il est égal au produit de la totalité de la rémunération annuelle telle
    # que visée à l’article L. 242-1 du code de la Sécurité sociale par un
    # coefficient.
    # Ce montant est majoré de 10 % pour les entreprises de travail temporaire
    # au titre des salariés temporaires pour lesquels elle est tenue à
    # l’obligation d’indemnisation compensatrice de congés payés.
    Pf = simulation.legislation_at(period.start).cotsoc.exo_bas_sal.fillon
    seuil = Pf.seuil
    tx_max = (Pf.tx_max * not_(majoration) + Pf.tx_max2 * majoration)
    if seuil <= 1:
        return 0
    ratio_smic_salaire = smic_proratise / (assiette + 1e-16)
    # règle d'arrondi: 4 décimales au dix-millième le plus proche
    taux_fillon = round_(tx_max * min_(1, max_(seuil * ratio_smic_salaire - 1, 0) / (seuil - 1)), 4)

    # Montant de l'allegment
    return taux_fillon * assiette

def compute_taux_exoneration(assiette_allegement, smic_proratise, taux_max, seuil_max, seuil_min = 1):
    ratio_smic_salaire = smic_proratise / (assiette_allegement + 1e-16)
    # règle d'arrondi: 4 décimales au dix-millième le plus proche ( # TODO: reprise de l'allègement Fillon unchecked)
    return round_(
        taux_max * min_(1, max_(seuil_max * seuil_min * ratio_smic_salaire - seuil_min, 0) / (seuil_max - seuil_min)),
        4,
        )

    def findNonLinear2D(self, t=0, npxls=af.MIN_PXLS_YX, phaseContrast=True):
        """
        calculate local cross correlation of projection images

        set self.mapyx and self.regions
        """
        #if self.refyz is None or self.refyx is None:
        self.setRefImg(removeEdge=False)

        # preparing the initial mapyx
        # mapyx is not inherited to avoid too much distortion
        self.mapyx = N.zeros((self.img.nt, self.img.nw, 2, self.img.ny, self.img.nx), N.float32)

        if N.all((N.array(self.mapyx.shape[-2:]) - self.img.shape[-2:]) >= 0):
            slcs = imgGeo.centerSlice(self.mapyx.shape[-2:], win=self.img.shape[-2:], center=None)
            self.mapyx = self.mapyx[slcs] # Ellipsis already added
        else:
            self.mapyx = imgFilters.paddingValue(self.mapyx, shape=self.img.shape[-2:], value=0)

        self.last_win_sizes = N.zeros((self.nt, self.nw), N.uint16)
        # calculation
        for w in range(self.img.nw):
            if w == self.refwave:
                self.mapyx[t,w] = 0
                continue

            self.echo('Projection local alignment -- W: %i' % w)
            if self.img.nz > 1:
                img = self.img.get3DArr(w=w, t=t)
                #img = af.fixSaturation(img, self.getSaturation(w=w, t=t))

                zs = N.round_(self.refzs-self.alignParms[t,w,0]).astype(N.int)
                if zs.max() >= self.img.nz:
                    zsbool = (zs < self.img.nz)
                    zsinds = N.nonzero(zsbool)[0]
                    zs = zs[zsinds]

                imgyx = af.prep2D(img, zs=zs, removeEdge=False)
                del img
            else:
                imgyx = N.squeeze(self.img.get3DArr(w=w, t=t))
                #imgyx = af.fixSaturation(imgyx, self.getSaturation(w=w, t=t))

            affine = self.alignParms[t,w]

            imgyx = imgyx.astype(N.float32)
            self.refyx = self.refyx.astype(N.float32)

            yxs, regions, arr2, win = af.iterWindowNonLinear(imgyx, self.refyx, npxls, affine=affine, initGuess=self.mapyx[t,w], phaseContrast=self.phaseContrast, maxErr=self.maxErrYX, cthre=self.cthre, echofunc=self.echofunc)

            self.mapyx[t,w] = yxs
            if self.regions is None or self.regions.shape[-2:] != regions.shape:
                self.regions = N.zeros((self.img.nt, self.img.nw)+regions.shape, N.uint16)
            self.regions[t,w] = regions

            self.last_win_sizes[t,w] = win
            
            self.progress()

        self.echo('Projection local alignment done')

def compute_subset_score(indices, pred_set, y):
    subset = [vect for i, vect in enumerate(pred_set) if i in indices]
    mean_preds = sp.mean(subset, axis=0)
    mean_preds = np.round_(mean_preds, decimals=0)
    mean_score = compute_score(y, mean_preds)

    return mean_score, indices

    def __get_figure_bar(self, dataframe_main, column_name, xlabel, ylabel):
        if not is_all_zeros(dataframe_main[column_name]):
            mean = np.round_(dataframe_main[column_name].mean(), 4)

            # Save to file
            if ENABLE_SEABORN:
                output_file_name = self.__report_file.split('.')[0]
                output_file_name += '_' + xlabel + '_' + column_name + '.eps'

                fig, ax = plt.subplots(1, 1)
                ax.get_xaxis().set_visible(False)

                title = 'Avg: ' + column_name + '=' + str(mean)
                bar_plot = dataframe_main[column_name].plot(ax=ax,
                                                            kind='bar',
                                                            table=True,
                                                            title=title)
                bar_plot.figure.savefig(output_file_name)
                bar_plot.figure.clf()

            # Bokeh draw
            bar_ratio = Bar(dataframe_main,
                            values=column_name,
                            label='ModuleId',
                            xlabel=xlabel,
                            ylabel=ylabel,
                            title='Avg: ' + column_name + '=' + str(mean))
            return bar_ratio