def getAvgAmplitudes(event_array, trace_array, time_range=None):
    """This routine takes an event_array (time x cells) and
    corresponding trace array and returns the average amplitudes of
    events in each cell.

    :param: event_array - 2 or 3d numpy event array (time x cells, or time x cells x trials)
    :param: time_range - optional list of 2 numbers limiting the time range to count events
    :returns: 2d masked numpy array of event average amplitudes. size is cells x largest number of events.
              masked entries are account for variable number of events
    """
    event_array = np.atleast_3d(event_array)
    trace_array= np.atleast_3d(trace_array)

    max_num_events = getCounts(event_array).max()
    time, cells, trials = event_array.shape

    amps = np.zeros((cells, trials, int(max_num_events)))
    amps[:] = np.nan

    for cell in range(cells):
        for trial in range(trials):
            event_ids = np.unique(event_array[:,cell,trial])[1:]
            for i, event_id in enumerate(event_ids):
                amps[cell, trial, i] = trace_array[event_array == event_id].mean()
    amps = np.ma.array(amps, mask=np.isnan(amps))
    amps = np.squeeze(amps)

    return np.ma.masked_array(amps, np.isnan(amps))

def hausdorffnorm(A, B):
    '''
    Finds the hausdorff norm between two matrices A and B.
    INPUTS:
    A: numpy array
    B : numpy array
    OUTPUTS:
    Housdorff norm between matrices A and B
    '''
    # ensure matrices are 3 dimensional, and shaped conformably
    if len(A.shape) == 1:
        A = np.atleast_2d(A)

    if len(B.shape) == 1:
        B = np.atleast_2d(B)

    A = np.atleast_3d(A)
    B = np.atleast_3d(B)

    x, y, z = B.shape
    A = np.reshape(A, (z, x, y))
    B = np.reshape(B, (z, x, y))

    # find hausdorff norm: starting from A to B
    z, x, y = B.shape
    temp1 = np.tile(np.reshape(B.T, (y, z, x)), (max(A.shape), 1))
    temp2 = np.tile(np.reshape(A.T, (y, x, z)), (1, max(B.shape)))
    D1 = np.min(np.sqrt(np.sum((temp1-temp2)**2, 0)), axis=0)

    # starting from B to A
    temp1 = np.tile(np.reshape(A.T, (y, z, x)), (max(B.shape), 1))
    temp2 = np.tile(np.reshape(B.T, (y, x, z)), (1, max(A.shape)))
    D2 = np.min(np.sqrt(np.sum((temp1-temp2)**2, 0)), axis=0)

    return np.max([D1, D2])

    def _viewStatesNGL(self, states, statetype, protein, ligand, mols, numsamples):
        if states is None:
            states = range(self.macronum)
        if isinstance(states, int):
            states = [states]
        if mols is None:
            mols = self.getStates(states, statetype, numsamples=min(numsamples, 15))
        colors = [0, 1, 3, 4, 5, 6, 7, 9]
        if protein is None and ligand is None:
            raise NameError('Please provide either the "protein" or "ligand" parameter for viewStates.')
        if protein:
            mol = Molecule()
        if ligand:
            mol = mols[0].copy()
            mol.remove(ligand, _logger=False)
            mol.coords = np.atleast_3d(mol.coords[:, :, 0])
            mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')
        for i, s in enumerate(states):
            if protein:
                mol.reps.add(sel='segid ST{}'.format(s), style='NewCartoon', color='Index')
            if ligand:
                mol.reps.add(sel='segid ST{}'.format(s), style='Licorice', color=colors[np.mod(i, len(colors))])
                mols[i].filter(ligand, _logger=False)

            mols[i].set('segid', 'ST{}'.format(s))
            tmpcoo = mols[i].coords
            for j in range(mols[i].numFrames):
                mols[i].coords = np.atleast_3d(tmpcoo[:, :, j])
                mol.append(mols[i])

        w = mol.view(viewer='ngl')
        self._nglButtons(w, statetype, states)
        return w

def convertRotMatToRisoeU(rMats, U0, symTag='Oh'):
    """
    Makes GrainSpotter gff ouput

    U11 U12 U13 U21 U22 U23 U13 U23 U33

    and takes it into the LLNL/APS frame of reference

    Urows comes from grainspotter's gff output
    U0 comes from XRD.crystallography.latticeVectors.U0
    """
    R = hexrd.XRD.Rotations # formerly import
    
    numU = num.shape(num.atleast_3d(rMats))[0]
    
    Rsamp = num.dot( R.rotMatOfExpMap(piby2*Zl), R.rotMatOfExpMap(piby2*Yl) )
    qin  = R.quatOfRotMat(num.atleast_3d(rMats))
    print "quaternions in (LLNL convention):"
    print qin.T
    qout = num.dot( R.quatProductMatrix( R.quatOfRotMat(Rsamp.T), mult='left' ), \
                    num.dot( R.quatProductMatrix( R.quatOfRotMat(U0), mult='right'),  \
                             qin ).squeeze() ).squeeze()
    if qout.ndim == 1:
        qout = toFundamentalRegion(qout.reshape(4, 1), crysSym=symTag, sampSym=None)
    else:
        qout = toFundamentalRegion(qout, crysSym=symTag, sampSym=None)
    print "quaternions out (Risoe convention, symmetrically reduced)"
    print qout.T
    Uout = R.rotMatOfQuat(qout)
    return Uout

    def depth_image(self):
        self._call_on_changed()

        gl = self.glb
        gl.Clear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        gl.PolygonMode(GL_FRONT_AND_BACK, GL_FILL)
        draw_noncolored_verts(gl, self.camera.v.r, self.f)
        result = np.asarray(deepcopy(gl.getDepth()), np.float64)

        if self.overdraw:
            gl.PolygonMode(GL_FRONT_AND_BACK, GL_LINE)
            draw_noncolored_verts(gl, self.camera.v.r, self.f)
            overdraw = np.asarray(deepcopy(gl.getDepth()), np.float64)
            gl.PolygonMode(GL_FRONT_AND_BACK, GL_FILL)
            boundarybool_image = self.boundarybool_image
            result = overdraw*boundarybool_image + result*(1-boundarybool_image)

        if hasattr(self, 'background_image'):
            if False: # has problems at boundaries, not sure why yet
                bg_px = self.visibility_image == 4294967295
                fg_px = 1 - bg_px
                result = bg_px * self.background_image + fg_px * result
            else:
                tmp = np.concatenate((np.atleast_3d(result), np.atleast_3d(self.background_image)), axis=2)
                result = np.min(tmp, axis=2)

        return result

 def GetPSF(self, vshint = None):
     psfKey = (self.psfType, self.psfFilename, self.lorentzianFWHM, self.beadDiameter)
     
     if not psfKey in self._psfCache.keys():
         if self.psfType == 'file':
             psf, vs = np.load(self.psfFilename)
             psf = np.atleast_3d(psf)
             
             self._psfCache[psfKey] = (psf, vs)        
         elif (self.psfType == 'Laplace'):
             from scipy import stats
             sc = self.lorentzianFWHM/2.0
             X, Y = np.mgrid[-30.:31., -30.:31.]
             R = np.sqrt(X*X + Y*Y)
             
             if not vshint is None:
                 vx = vshint[0]
             else:
                 vx = sc/2.
             
             vs = type('vs', (object,), dict(x=vx/1e3, y=vx/1e3))
             
             psf = np.atleast_3d(stats.cauchy.pdf(vx*R, scale=sc))
                 
             self._psfCache[psfKey] = (psf/psf.sum(), vs)
         elif (self.psfType == 'bead'):
             from PYME.Deconv import beadGen
             psf = beadGen.genBeadImage(self.beadDiameter/2, vshint)
             
             vs = type('vs', (object,), dict(x=vshint[0]/1e3, y=vshint[1]/1e3))
             
             self._psfCache[psfKey] = (psf/psf.sum(), vs)
             
             
     return self._psfCache[psfKey]

def convertRotMatToFableU(rMats, U0=num.eye(3), symTag='Oh', display=False):
    """
    Makes GrainSpotter gff ouput

    U11 U12 U13 U21 U22 U23 U13 U23 U33

    and takes it into the hexrd/APS frame of reference

    Urows comes from grainspotter's gff output
    U0 comes from xrd.crystallography.latticeVectors.U0
    """
    numU = num.shape(num.atleast_3d(rMats))[0]

    qin  = quatOfRotMat(num.atleast_3d(rMats))
    qout = num.dot( quatProductMatrix( quatOfRotMat(fableSampCOB.T), mult='left' ), \
                    num.dot( quatProductMatrix( quatOfRotMat(U0), mult='right'),  \
                             qin ).squeeze() ).squeeze()
    if qout.ndim == 1:
        qout = toFundamentalRegion(qout.reshape(4, 1), crysSym=symTag, sampSym=None)
    else:
        qout = toFundamentalRegion(qout, crysSym=symTag, sampSym=None)
    if display:
        print "quaternions in (hexrd convention):"
        print qin.T
        print "quaternions out (Fable convention, symmetrically reduced)"
        print qout.T
        pass
    Uout = rotMatOfQuat(qout)
    return Uout

    def _raw_predict(self, Xnew, full_cov=False, kern=None):
        """
        Make a prediction for the latent function values
        """

        if kern is None: kern = self.kern

        if not isinstance(Xnew, VariationalPosterior):
            Kx = kern.K(self.Z, Xnew)
            mu = np.dot(Kx.T, self.posterior.woodbury_vector)
            if full_cov:
                Kxx = kern.K(Xnew)
                if self.posterior.woodbury_inv.ndim == 2:
                    var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx))
                elif self.posterior.woodbury_inv.ndim == 3:
                    var = Kxx[:,:,None] - np.tensordot(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx).T, Kx, [1,0]).swapaxes(1,2)
                var = var
            else:
                Kxx = kern.Kdiag(Xnew)
                var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T
        else:
            Kx = kern.psi1(self.Z, Xnew)
            mu = np.dot(Kx, self.posterior.woodbury_vector)
            if full_cov:
                raise NotImplementedError, "TODO"
            else:
                Kxx = kern.psi0(self.Z, Xnew)
                psi2 = kern.psi2(self.Z, Xnew)
                var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1)
        return mu, var

    def Project(self, projType):
        import numpy as np
        from PYME.DSView.image import ImageStack
        from PYME.DSView import ViewIm3D
        import os

        if projType == 'mean':        
            filt_ims = [np.atleast_3d(self.image.data[:,:,:,chanNum].mean(2)) for chanNum in range(self.image.data.shape[3])]
        elif projType == 'max':
            filt_ims = [np.atleast_3d(self.image.data[:,:,:,chanNum].max(2)) for chanNum in range(self.image.data.shape[3])]

        fns = os.path.split(self.image.filename)[1]        
        
        im = ImageStack(filt_ims, titleStub = '%s - %s' %(fns, projType))
        im.mdh.copyEntriesFrom(self.image.mdh)
        im.mdh['Parent'] = self.image.filename
        im.mdh['Processing.Projection'] = projType

        if self.dsviewer.mode == 'visGUI':
            mode = 'visGUI'
        else:
            mode = 'lite'

        dv = ViewIm3D(im, mode=mode, glCanvas=self.dsviewer.glCanvas)

        #set scaling to (0,1)
        for i in range(im.data.shape[3]):
            dv.do.Gains[i] = 1.0

    def getXYZ(self):
        """ Get XYZ values in world coordinates for each pixel.

        Usage: XYZ = self.getXYZ()

        Input:
            -NONE-

        Output:
            XYZ - M-by-N-by-3 matrix of [X Y Z] world coordinates for each pixel
        """

        if self.XY is not None:
            return np.c_[np.atleast_3d(self.XY), np.atleast_3d(self)]
        else:
            x = np.arange(0, self.width)
            y = np.arange(0, self.height)
            xx, yy = np.meshgrid(x, y)

            XY = np.zeros((self.height, self.width, 2))

            # From depth map to Point Cloud --> use focal distance
            XY[:, :, 0] = (xx - self.K[0, 2]) / self.K[0, 0]
            XY[:, :, 1] = (yy - self.K[1, 2]) / self.K[1, 1]
            XY = XY * np.atleast_3d(self)
            return np.c_[np.atleast_3d(self.XY), np.atleast_3d(self)]

 def append(self, *args):
     if len(args)<1:
         pass
     else:
         smp=self.mapped_parameters
         print args
         for arg in args:
             #object parameters
             mp=arg.mapped_parameters
             
             if mp.original_filename not in smp.original_files.keys():
                 smp.original_files[mp.original_filename]=arg
                 # add the data to the aggregate array
                 if self.data==None:
                     self.data=np.atleast_3d(arg.data)
                 else:
                     self.data=np.append(self.data,np.atleast_3d(arg.data),axis=2)
                 print "File %s added to aggregate."%mp.original_filename
             else:
                 print "Data from file %s already in this aggregate. \n \
 Delete it first if you want to update it."%mp.original_filename
         # refresh the axes for the new sized data
         self.axes_manager=AxesManager(self._get_undefined_axes_list())
         smp.original_filename="Aggregate Image: %s"%smp.original_files.keys()
         self.summary()

 def GetPSF(self, vshint = None):
     import numpy as np
     from scipy import stats
     
     PSFMode = self.nb2.GetCurrentPage().PSFMode
     #get PSF from file
     if PSFMode == 'File':
         psf, vs = np.load(self.GetPSFFilename())
         psf = np.atleast_3d(psf)
         
         return (self.GetPSFFilename(), psf, vs)        
     elif (PSFMode == 'Laplace'):
         sc = float(self.tLaplaceFWHM.GetValue())/2.0
         X, Y = np.mgrid[-30.:31., -30.:31.]
         R = np.sqrt(X*X + Y*Y)
         
         if not vshint == None:
             vx = vshint*1e3
         else:
             vx = sc/2.
         
         vs = type('vs', (object,), dict(x=vx/1e3, y=vx/1e3))
         
         psf = np.atleast_3d(stats.cauchy.pdf(vx*R, scale=sc))
             
         return 'Generated Laplacian, FWHM=%f' % (2*sc), psf/psf.sum(), vs

    def __init__(self, root, noise, option):
        self.root = root
        self.nFeatures = 4
        self.kernelSize = 3
        self.poolLength = 2
        self.nLambda = 112
        self.batchSize = 64
        self.nClasses = [50] * 12
        self.noise = noise
        self.option = option

        self.labels = ['T0', 'T1', 'T2', 'vmic', 'B0', 'B1', 'v0', 'v1', 'thB0', 'thB1', 'chiB0', 'chiB1']

        self.n_pars = len(self.labels)

# BField, theta, chi, vmac, damping, B0, B1, doppler, kl
        self.lower = np.asarray([-3000.0, -1500.0, -3000.0, 0.0, 0.0, 0.0, -7.0, -7.0, 0.0, 0.0, 0.0, 0.0], dtype='float32')
        self.upper = np.asarray([3000.0, 3000.0, 5000.0, 4.0, 3000.0, 3000.0, 7.0, 7.0, 180.0, 180.0, 180.0, 180.0], dtype='float32')
        
        self.dataFile = "../database/database_sir.h5"

        f = h5py.File(self.dataFile, 'r')
        pars = f.get("parameters")
        stokes = f.get("stokes")
        self.nModels, _ = pars.shape
                
        self.nTraining = int(self.nModels * 0.9)
        self.nValidation = int(self.nModels * 0.1)

# Standardize Stokes parameters
        std_values = np.std(np.abs(stokes[0:self.nTraining,:,:]),axis=0)
        stokes /= std_values[None,:,:]

# Save normalization values
        np.save('{0}_normalization.npy'.format(self.root), std_values)
        
        print("Training set: {0}".format(self.nTraining))

        print("Validation set: {0}".format(self.nValidation))

        self.inTrain = []
        for i in range(4):            
            self.inTrain.append(np.atleast_3d(stokes[0:self.nTraining,i,:]).astype('float32'))

        self.inTest = []
        for i in range(4):            
            self.inTest.append(np.atleast_3d(stokes[self.nTraining:,i,:]).astype('float32'))

        self.outTrain = []
        for i in range(self.n_pars):
            outTrain = np.floor((pars[0:self.nTraining, i] - self.lower[i]) / (self.upper[i] - self.lower[i]) * self.nClasses[i]).astype('int32')            
            self.outTrain.append(np_utils.to_categorical(outTrain, self.nClasses[i]))

        self.outTest = []
        for i in range(self.n_pars):
            outTest = np.floor((pars[self.nTraining:, i] - self.lower[i]) / (self.upper[i] - self.lower[i]) * self.nClasses[i]).astype('int32')
            self.outTest.append(np_utils.to_categorical(outTest, self.nClasses[i]))

        f.close()

def regression_plot(Z,X,band_names=None,visible_only=True,figsize=(12,7)):
    """
    Produce a figure with a plot for each image band that displays the
    relationship between depth and radiance and gives a visual representation
    of the regression carried out in the `slopes` and `regressions` methods.

    Notes
    -----
    This method doesn't come directly from Lyzenga 1978 but the author of this
    code found it helpful.

    Parameters
    ----------
    Z : np.ma.MaskedArray
        Array of depth values repeated for each band so that Z.shape==X.shape.
        The mask needs to be the same too so that Z.mask==X.mask for all the
        bands.
    X : np.ma.MaskedArray
        The array of log transformed radiance values from equation B1 of
        Lyzenga 1978.

    Returns
    -------
    figure
        A matplotlib figure.
    """
    if band_names is None:
        band_names = ['Band'+str(i+1) for i in range(X.shape[-1])]
    nbands = X.shape[-1]
    if np.atleast_3d(Z).shape[-1] == 1:
        Z = np.repeat(np.atleast_3d(Z), nbands, 2)
    if visible_only:
        fig, axs = plt.subplots( 2, 3, figsize=figsize)
    else:
        fig, axs = plt.subplots( 2, 4, figsize=figsize )
    regs = regressions(Z,X)
    for i, ax in enumerate(axs.flatten()):
        if i > nbands-1:
            continue
        slp, incpt, rval = regs[:,i]
        # print X.shape, Z.shape
        x, y = equalize_array_masks(Z[...,i], X[...,i])
        if x.count() < 2:
            continue
        x, y = x.compressed(), y.compressed()
        # print "i = {}, x.shape = {}, y.shape = {}".format(i, x.shape, y.shape)
        ax.scatter( x, y, alpha=0.1, edgecolor='none', c='gold' )
        smth = lowess(y,x,frac=0.2)
        # ax.plot(smth.T[0],smth.T[1],c='black',alpha=0.5)
        ax.plot(smth.T[0],smth.T[1],c='black',alpha=0.5,linestyle='--')
        reglabel = "m=%.2f, r=%.2f" % (slp,rval)
        f = lambda x: incpt + slp * x
        ax.plot( x, f(x), c='brown', label=reglabel, alpha=1.0 )
        ax.set_title( band_names[i] )
        ax.set_xlabel( r'Depth (m)' )
        ax.set_ylabel( r'$X_i$' )
        ax.legend(fancybox=True, framealpha=0.5)
    plt.tight_layout()
    return fig

 def load_texture(self, filename, gray=False, blur=False):
     print "Loading texture from " + filename
     self.pixels = np.atleast_3d(scipy.misc.imread(filename, flatten=gray))
     if blur:
         self.pixels = \
             np.atleast_3d(scipy.misc.imfilter(self.pixels.squeeze(), 
                                                'blur'))
     print "Done loading texture"

def update(fig):
    """Fit new pointing model and update plots."""
    # Perform early redraw to improve interactivity of clicks (which typically change state of target dots)
    # Target state: 0 = flagged, 1 = unflagged, 2 = highlighted
    target_state = keep * ((target_index == fig.highlighted_target) + 1)
    # Specify colours of flagged, unflagged and highlighted dots, respectively, as RGBA tuples
    dot_colors = np.choose(target_state, np.atleast_3d(np.vstack([(1,1,1,1), (0,0,1,1), (1,0,0,1)]))).T
    for ax in fig.axes[:7]:
        ax.dots.set_facecolors(dot_colors)
    fig.canvas.draw()

    # Fit new pointing model and update results
    params, sigma_params = new_model.fit(az[keep], el[keep], measured_delta_az[keep], measured_delta_el[keep],
                                         std_delta_az[keep], std_delta_el[keep], enabled_params)
    new.update(new_model)

    # Update rest of figure
    fig.texts[3].set_text("$\chi^2$ = %.1f" % new.chi2)
    fig.texts[4].set_text("all sky rms = %.3f' (robust %.3f')" % (new.sky_rms, new.robust_sky_rms))
    new.metrics(target_index == fig.highlighted_target)
    fig.texts[5].set_text("target sky rms = %.3f' (robust %.3f')" % (new.sky_rms, new.robust_sky_rms))
    new.metrics(keep)
    fig.texts[-1].set_text(unique_targets[fig.highlighted_target])
    # Update model parameter strings
    for p, param in enumerate(display_params):
        fig.texts[2*p + 6].set_text(param_to_str(new_model, param) if enabled_params[param] else '')
        # HACK to convert sigmas to arcminutes, but not for P9 and P12 (which are scale factors)
        # This functionality should really reside inside the PointingModel class
        std_param = rad2deg(sigma_params[param]) * 60. if param not in [8, 11] else sigma_params[param]
        std_param_str = ("%.2f'" % std_param) if param not in [8, 11] else ("%.0e" % std_param)
        fig.texts[2*p + 7].set_text(std_param_str if enabled_params[param] and opts.use_stats else '')
        # Turn parameter string bold if it changed significantly from old value
        if np.abs(params[param] - old_model.values()[param]) > 3.0 * sigma_params[param]:
            fig.texts[2*p + 6].set_weight('bold')
            fig.texts[2*p + 7].set_weight('bold')
        else:
            fig.texts[2*p + 6].set_weight('normal')
            fig.texts[2*p + 7].set_weight('normal')
    daz_az, del_az, daz_el, del_el, quiver, before, after = fig.axes[:7]
    # Update quiver plot
    quiver_scale = 0.1 * fig.quiver_scale_slider.val * np.pi / 6 / deg2rad(old.robust_sky_rms / 60.)
    quiver.quiv.set_segments(quiver_segments(new.residual_az, new.residual_el, quiver_scale))
    quiver.quiv.set_color(np.choose(keep, np.atleast_3d(np.vstack([(0.3,0.3,0.3,0.2), (0.3,0.3,0.3,1)]))).T)
    # Update residual plots
    daz_az.dots.set_offsets(np.c_[rad2deg(az), rad2deg(new.residual_xel) * 60.])
    del_az.dots.set_offsets(np.c_[rad2deg(az), rad2deg(new.residual_el) * 60.])
    daz_el.dots.set_offsets(np.c_[rad2deg(el), rad2deg(new.residual_xel) * 60.])
    del_el.dots.set_offsets(np.c_[rad2deg(el), rad2deg(new.residual_el) * 60.])
    after.dots.set_offsets(np.c_[np.arctan2(new.residual_el, new.residual_xel), new.abs_sky_error])
    resid_lim = 1.2 * max(new.abs_sky_error.max(), old.abs_sky_error.max())
    daz_az.set_ylim(-resid_lim, resid_lim)
    del_az.set_ylim(-resid_lim, resid_lim)
    daz_el.set_ylim(-resid_lim, resid_lim)
    del_el.set_ylim(-resid_lim, resid_lim)
    before.set_ylim(0, resid_lim)
    after.set_ylim(0, resid_lim)
    # Redraw the figure
    fig.canvas.draw()

    def __init__(self, root, noise, option):
        self.root = root
        self.nFeatures = 100
        self.kernelSize = 3
        self.poolLength = 2
        self.nLambda = 50
        self.batchSize = 256
        self.nClasses = [50, 50, 50, 50, 10, 20, 20, 20, 20]
        self.noise = noise
        self.option = option

                                # BField, theta, chi, vmac, damping, B0, B1, doppler, kl
        self.lower = np.asarray([0.0,      0.0,   0.0, -7.0, 0.0,  0.15, 0.15, 0.20,  1.0], dtype='float32')
        self.upper = np.asarray([3000.0, 180.0, 180.0,  7.0, 0.5,   1.2,  1.2, 0.80,  5.0], dtype='float32')
        
        self.dataFile = "/net/duna/scratch1/aasensio/deepLearning/milne/database/database_6301_hinode_1component.h5"

        f = h5py.File(self.dataFile, 'r')
        pars = f.get("parameters")
        stokes = f.get("stokes")
        self.nModels, _ = pars.shape

        std_values = np.std(np.abs(stokes),axis=0)
        stokes /= std_values[None,:,:]

	self.sigma_noise = 1e-3 / np.mean(std_values, axis=0)

# Save normalization values        
        np.save('{0}_normalization.npy'.format(self.root), std_values)
                
        self.nTraining = int(self.nModels * 0.9)
        self.nValidation = int(self.nModels * 0.1)
        
        print("Training set: {0}".format(self.nTraining))

        print("Validation set: {0}".format(self.nValidation))

        self.inTrain = []
        for i in range(4):            
            self.inTrain.append(np.atleast_3d(stokes[0:self.nTraining,:,i]).astype('float32'))

        self.inTest = []
        for i in range(4):            
            self.inTest.append(np.atleast_3d(stokes[self.nTraining:,:,i]).astype('float32'))

        self.outTrain = []
        for i in range(9):
            outTrain = np.floor((pars[0:self.nTraining, i] - self.lower[i]) / (self.upper[i] - self.lower[i]) * self.nClasses[i]).astype('int32')            
            self.outTrain.append(np_utils.to_categorical(outTrain, self.nClasses[i]))

        self.outTest = []
        for i in range(9):
            outTest = np.floor((pars[self.nTraining:, i] - self.lower[i]) / (self.upper[i] - self.lower[i]) * self.nClasses[i]).astype('int32')
            self.outTest.append(np_utils.to_categorical(outTest, self.nClasses[i]))

        f.close()

def merge3D(A, B, position):
    A = np.atleast_3d(A)
    B = np.atleast_3d(B)

    mat_temp = np.nan * np.ones([max(A.shape[0] + position[0], B.shape[0]), max(position[1] + A.shape[1], B.shape[1]),
                                 max(position[2] + A.shape[2], B.shape[2])])
    mat_temp[0:B.shape[0], 0:B.shape[1], 0:B.shape[2]] = B
    mat_temp[position[0]:position[0] + A.shape[0], position[1]:position[1] + A.shape[1],
    position[2]:position[2] + A.shape[2]] = A

    return mat_temp

 def covariance(self):
     """
     Posterior covariance
     $$
     K_{xx} - K_{xx}W_{xx}^{-1}K_{xx}
     W_{xx} := \texttt{Woodbury inv}
     $$
     """
     if self._covariance is None:
         #LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
         self._covariance = (np.atleast_3d(self._K) - np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K, [1,0]).T).squeeze()
         #self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K)
     return self._covariance

 def get_transparent_item_heights_and_mask(self, low_limit, high_limit):
     low_limit_3d = numpy.atleast_3d(low_limit)
     high_limit_3d = numpy.atleast_3d(high_limit)
     max_height = self.blocks.shape[2]
     shape = self.blocks.shape
     trimmed_shape = (shape[0], shape[1], shape[2]-1)
     cell_depth = numpy.indices(trimmed_shape)[2]
     cell_is_selected = numpy.logical_and(cell_depth>=low_limit_3d, cell_depth<high_limit_3d)
     selectable_substance = numpy.logical_and(tileid_is_transparent[self.blocks[:,:,:-1]], self.blocks[:,:,:-1] != 0)
     potential_blocks = numpy.logical_and(selectable_substance, cell_is_selected)
     floor_heights = (max_height-2)-numpy.argmax(potential_blocks[:,:,::-1], axis=2)
     mask = get_cells_using_heightmap(potential_blocks, floor_heights)
     return numpy.clip(floor_heights, low_limit, high_limit), mask