dataesr/scientific-paragraphs-categorization huggingface dataset:
Lastly we highlight as always for the Europe dataset that the average original error is lower than the relative one. As in the repeated RED case, we remark that ST-HOSVD is convenient when the first two components of multilinear are strictly lower than 500, otherwise T-HOSVD provides better results. Then we  present the final statistics for original and relative errors of Rao estimates for Earth dataset. As always in vector epO and in epR are respectively stored the original and relative errors per pixel. In this case comparing Table  9  with Table  11 , we notice that starting tensors with repeated NIR band provides on average better results. Indeed for both the decomposition techniques at the fifth approximation the average original and relative error per pixel is around 17. 5% in the first case and 16.
In Figure  11  there are approximated Rao estimates at different multilinear rank, for this NIR repeated case. As before, our eyes hardly perceive differences between Rao computed over NASA or self-made NDVI and over the last three approximated NDVIs. Example 4. 24. Firstly we remark that May 2014 do not realise the minimum original and relative error starting from tensors with repeated NIR band. Therefore these are not the worst approximations in the Earth dataset, in this second case. In Figure  12  there are approximated Rao estimates at different multilinear rank, for this NIR repeated case. As in Figure  10  the problematic areas are the North pole and the Pacific ocean. Indeed we notice that there is a biodiversity overestimation in the Arctic regions, which decreases for increasing multilinear rank components.
Conclusion We have shown that the presented approach is extremely convenient for Rényi index estimation to save storage memory. As reported in Tables 4,6, 5 and 7 the average error per pixel is around 5. 5% for Earth dataset and is close to 7. 6% for the Europe's one when using respectively 16. 4% and 14. 8% of the total tensor information. Moreover we have seen that starting from tensor with repeated NIR bands is more convenient than with RED repeated, in the Rényi case. For Rao case starting with twice NIR rasters is more convenient only for Earth dataset. Lastly we have remarked In the Rao case we face greater on average original and relative error per pixel. Indeed in Tables  8  and 9  we have an average error per pixel close to 19. 5% and 17% for Europe and for Earth dataset respectively.
If tw_values is an empty vector, we return missing value. We do a first Parallel call of IndexOut inside IndexOP and a second one of IndexOP outside it. Even in this case the presented code will work only when the moving window is entirely contained in the raster. To speed up the computations, we develop a special version of this code which works when the moving window is not fully contained. It is available here  [Ian19]. 6. 2. Approximation codes for T-HOSVD. In order to implement the approximation Algorithm 1 of T-HOSVD we need a python library expressly developed to manage tensors. Moreover we also want that this library to interact properly with NumPy the python most used library for numerical computation, cf. [Oli06]. Therefore we choose TensorLy, cf. [KPAP19]. This recently developed library, firstly presented in 2016, wants to make tensor study and manipulation easy and accessible.
Besides their creators projected TensorLy to perfectly match with other famous python libraries, as NumPy, SciPy. They developed most of the main tensor operations and related functions. To compute the singular value decomposition of tensors flattening, the central step in T-HOSVD and ST-HOSVD algorithm, we need a python function able to manage huge arrays. This is not a trivial task, since development of python functions is left to singular initiative. We decide to use svds function from SciPy sparse linear algebra function, cf. [JOP ]. This implementation of SVD takes advantage of matrix sparsity in performing the matrix-vector multiplication. The input arguments of ST-HOSVD and T-HOSVD implementations are the same. The first difference with T_hosvd is the initial declaration of core tensor, set equal to T.
, 2011). Furthermore, hippocampal activity has been documented during rest periods of early learning in humans  (Albouy et al. , 2015; Jacobacci et al. , 2020). Thus, it is possible that waking replay, through offline recapitulation of prior practice  (Carr et al. , 2011; Diba and Buzsaki, 2007; Foster and Wilson, 2006; Singer and Frank, 2009) , promotes wakeful consolidation of skill, an issue that has not been investigated in humans or animal models  (Carr et al. , 2011; Ramanathan et al. , 2015). In this study, we posited a contribution of neural replay in wakeful memory consolidation. We tested this by applying a multivariate analytical approach to magnetoencephalography (MEG) brain activity data acquired during skill learning to decode spontaneous, waking replay during the periods of rest in which sequences of simple discrete actions learned during interspersed practice are bound into a consolidated skill.
Small asterisks appeared above a sequence item when a keypress was recorded, providing feedback to the participant about their current location within the sequence. After completion of a full iteration of the sequence, the asterisks were removed. The asterisk display did not provide error feedback since they appeared for both correct and incorrect keypresses. During the 10 s interleaved rest periods, the sequence was replaced with a string of five ''X'' symbols, which participants were instructed to fixate on. Visual stimuli and task instructions were presented, and keypress responses recorded using a custom script running in E-Prime 2 (Psychology Software Tools, Inc. ).
Magnetoencephalography (MEG) MEG acquisition Continuous MEG was recorded at a sampling frequency of 600Hz with a CTF 275 MEG system (CTF Systems, Inc. , Canada) while participants were seated in an upright position. The system is composed of a whole head array of 275 radial 1 st -order gradiometer/ SQUID channels housed in a magnetically shielded room (Vacuumschmelze, Germany). Two of the gradiometers were malfunctioning and were not used. A third channel was removed (MLT16-1609) after visual inspection of the data indicated the presence of high noise levels for all recordings, resulting in 272 total channels of MEG data. Synthetic 3rd order gradient balancing was used to remove background noise online. A TTL trigger sent from the task computer to the MEG data acquisition computer was used to temporally align the behavioral and MEG data.
Determination of the early learning trial cutoff Here, early learning corresponds to the set of trials between trial 1 and a cutoff trial, T, on which 95% of the performance gains are first achieved at the group level. We determined T as follows. First, we fit the trial-by-trial group average correct sequence speed data with an exponential model of the form: LðtÞ = C 1 + C 2 $ À 1 À e Àlt Á ; (1) where L(t) represents the group average learning state on practice trial, t, C 1 and C 2 control the pre-training performance and asymptote (i. e. , -extent of learning), and l controls the learning rate. Parameters C 1 (boundary constraints = [0,5]), C 2  ([0,15] ) and l ([0,2]) were estimated using a constrained nonlinear least-squares method (MATLAB's lsqcurvefit, trust-region-reflective algorithm).
Magnetoencephalography (MEG) MEG data preparation and pre-processing Each participant's MEG data were organized into the following samples: pre-training rest (single segment of 5 min), inter-trial practice rest (36 segments of 10 s), practice (36 segments of 10 s) and post-training rest (single segment of 5 min). All MEG samples were band-pass filtered between 1-100Hz, notch-filtered at 60Hz and subjected to Independent Component Analysis (ICA). ICA was fit using MNE Python's  (Gramfort et al. , 2013)  implementation of FastICA  (Fan et al. , 2016)  after pre-whitening with Principal Component Analysis (set to retain the top 15 components). Spatial topographies and time series of all components were visually inspected and those determined to be strongly associated with known electrocardiogram, eye movement and blinks, and head movement artifacts were removed from the MEG sample.
MEG source reconstruction and parcellation Source reconstruction of the MEG data was performed using the standard pipeline in MNE Python  (Gramfort et al. , 2014). A single forward solution and inverse solution was calculated for each MEG sample. A boundary element model (BEM) required for generation of the forward solution was constructed for individual subjects from inner-skull and pial layer surfaces obtained by segmentation of the subject's T1-MRI volume in FreeSurfer  (Fischl, 2012). Hippocampal surface labels were obtained by parcellating the participant's MRI according to the coordinates in the Brainnetome atlas  (Fan et al. , 2016)  (https://atlas. brainnetome. org/). Neocortical and hippocampal source spaces were obtained by sampling the corresponding surfaces on a recursive octahedral grid with 4.
9mm between dipoles, and spatially aligned to the MEG sensor data with the co-registration data obtained from BrainSight. The BEM was then used to generate the forward solution at each source dipole location. The Linearly Constrained Minimum-Variance (LCMV) beamformer operator was used to compute the inverse solution from the recorded MEG sensor space data. The entire duration of each MEG sample was used to calculate the inverse solution data covariance matrix. All inverse solutions used the same sample noise covariance matrix, which was calculated over the first 20 s of the pre-training rest MEG sample for individual subjects. MNI-space transformations obtained from the FreeSurfer segmentation pipeline were then used to spatially register the resulting source reconstruction time series to the Brainnetome atlas, enabling quantitative estimation of millisecond-resolved brain activity for 212 neocortical and 4 hippocampal parcels (i.
e. , -216 total parcels). Parcel activity was estimated by averaging the time series of all source dipoles falling within a given parcel boundary. A process called mean-flipping was used to prevent within-parcel source cancellation when averaging. That is, the sign was flipped for all sources with a sign that differed from that of the average source. Since the sensor space data had been bandpass filtered with a cut-off frequency of 100Hz, the source data were downsampled from 600Hz to 200Hz to increase computation speed in subsequent processing steps. All remaining analyses were performed directly on the parcellated source-space MEG time series. Training and evaluation of keypress decoders Individualized keypress decoders were trained for each participant using Python's Scikit-learn library  (Pedregosa et al.
, 2011). The training data of a subject's decoder was assembled as follows. We time-averaged the parcellated source for each correct keypress in all 36 practice MEG samples between t ± Dt ms, where t is the keypress timestamp and Dt is treated as a cross-validated hyperparameter described in detail below. The timestamp for each keypress was obtained from behavioral data files output by E-prime. This produced the predictor matrix X [keypresses x parcels] and the target array y [keypresses x 1] where keypresses varied between subjects (some typed more sequences than others) and parcels = 216 as earlier described. We used X and y to train four one-versus-all Radial-basis Support Vector Machines (Python Scikit-learn's SVC) with a pipeline parameterized with a feature centering/scaling policy, a time window for averaging the MEG source around a keypress (Dt), a regularization term (C) and a kernel coefficient (g).
Grid search under a particular criterion worked as follows. First, the training set was divided into 5 mutually exclusive stratified partitions so that the proportions of all keypresses were preserved as in the original data (since the training sequence was 41324, each folder would have approximately twice as many keys = 4). Next, we used the random oversampler in Python's Imbalanced-learn toolbox  (Lemaitre et al. , 2017)  (with a ''not majority'' sampling policy) to independently oversamples keypresses 1, 2 and 3 within each training partition. Note that test partitions were never oversampled, and training partitions did not share common oversampled patterns as this would lead to double-dipping and ultimately overfitting. Finally, we performed the same 5-fold cross-validation procedure on each of the 420 configurations in the grid.
Note that the expected key at each position is always the keypress with the top detection probability. All participants had at least one decoder with a non-zero detection rate. The selected decoder was the one used to interrogate waking rest MEG data samples for replay events. Sequence replay detection and replay rate calculation Given a rest period MEG sample, X rest [T x parcels] , of T time-points x parcels recorded from a given subject, we estimated the probability that a given sequence, seq = ABCDE, was replayed at a time 0 % t < T within that sample to be the product of the probabilities of each of the sequence's constituent keys A, B, C, D and E ˛[1, 2, 3, 4] being observed both in the right order and ending at time t, that is: p seq = ABCDE t ½ = p key = A X resampled t À 4; : ½ À Á $p key = B X resampled t À 3; : ½ À Á /p key = E X resampled t; : ½ À Á ; (2) Where X resampled = resample X rest t À 5$tscl + 1 : t; : ½ À Á is an MEG sample segment extracted from t À 5$tscl + 1 to t at all 216 parcels and compressed down to 5 time-points.
A downside to SK is a dependence on large degrees of replication in order to support the MoM method. Binois, Gramacy, and Ludkovski (2018)  provide similar results but relaxing the necessity of having a large degree of replication, and place inference fully within a likelihood framework. That methodology is the primary focus of the implementation described here. Other approaches for achieving a degree of input-dependent variance include quantile kriging  (Plumlee and Tuo 2014) ; the use of pseudoinputs  (Snelson and Ghahramani 2005) , also sometimes called a predictive process  (Banerjee, Gelfand, Finley, and Sang 2008) ; and (non-GP-based) tree methods  (Pratola, Chipman, George, and McCulloch 2017a). Although the hetGP package has many aspects in common, and in some cases is directly inspired by these approaches, none of these important methodological contributions are, to our knowledge, coupled with an open source R implementation.
In the computer experiments and machine learning communities, sequential design of experiments/active learning (e. g. ,  Seo, Wallat, Graepel, and Obermayer 2000; Gramacy and Lee 2009)  and Bayesian optimization (e. g. ,  Snoek, Larochelle, and Adams 2012)  are a primary goal of modeling. The hetGP package provides hooks for learning and design to reduce predictive variance, organically determining an input-dependent degree of replication required to efficiently separate signal from noise  (Binois, Huang, Gramacy, and Ludkovski 2019) ; for Bayesian optimization; for contour finding  (Lyu, Binois, and Ludkovski 2018) , and for leptokurtic responses  (Shah, Wilson, and Ghahramani 2014; Chung, Binois, Gramacy, Bardsley, Moquin, Smith, and Smith 2019). A description and detailed examples are provided herein.
The DiceKriging  (Roustant, Ginsbourger, and Deville 2012)  package is often used for SK experiments; however, users must preprocess the data and calculate their own moment-based variance estimates. The mlegp package  (Dancik and Dorman 2008)  offers a more hands-off SK experience, facilitating replicate detection, but without coupled modeling and sequential design features. Other R packages include GPfit (MacDonald,  Ranjan, and Chipman 2015) , which focuses expressly on deterministic modeling, whereas those like kergp  (Deville, Ginsbourger, and Roustant 2018)  target flexible covariance kernels, including additive or qualitative versions, in the homoskedastic setting. The remainder of this paper is organized as follows. Section 2 reviews hetGP's approach to ordinary (homoskedastic) GP regression, and linear algebra identities that enable substantial speedups when the design has a high degree of replication.
Gaussian process modeling under replication Here we address computer experiments involving a function f (x) : x ∈ R d → R requiring expensive simulation to approximate a real-world (physical) relationship. An emulator fN is a regression or response surface fit to input-output pairs (x 1 , y 1 ),. , (x N , y N ), where y i = f (x i ) + ǫ i , with a centered (Gaussian) random noise ǫ i , e. g. , ǫ i ∼ N (0, r(x i )). Predictive equations from fN , notated as fN (x ′ ) for a new location x ′ , may serve as a cheap surrogate for f (x ′ ) at new inputs x ′ for visualization, sensitivity analysis, optimization, and so forth. The hetGP package targets thrifty GP surrogates for stochastic computer model simulations whose variance and mean are both believed to vary nonlinearly. Although our emphasis and vocabulary are tilted toward computer surrogate modeling throughout, many of the techniques we describe are equally well suited to applications in machine learning and geostatistics or anywhere else GP regression may be applied.
The data-generating function coded below was introduced by  (Gramacy and Lee 2009)  to illustrate challenges in GP regression when the response surface is nonstationary. (Heteroskedasticity is a form of nonstationarity; and although here the nonstationarity is in the mean,  Binois et al. (2018, see Appendix C)  showed that nevertheless hetGP methods can offer an appropriate quantification of predictive uncertainty on these data. ) In order to help separate signal from noise, degree 2 replication is used. Replication is an important focus of this manuscript, with further discussion in Section 2. 2. Initial designs are Latin hypercube samples, generated by using the lhs package  (Carnell 2018).
In 1D applications the error bars derived from predictive mean and variance surfaces take on a sausage shape, being fat away from the training data sites and thin, or "pinched," at the training locations. Many libraries, such as those cited in our introduction, offer automations of these procedures. In this manuscript we highlight features of the hetGP package, which offers similar calculations as a special case. The code below provides an illustration. R> fit  y, rep(Lwr, 2), rep(Upr, 2), known = list(beta0 = 0), + init = c(list(theta = rep(0. 1, 2) , g = 0. 1 * var(y)))) R> c(fit$theta, fit$g) [1] 0. 631256079 2. 125039287 0. 007229522                                                                               These estimated parameters are nearly identical to the ones above, obtained "by hand".
Arguments lower and upper (3rd and 4th arguments) as well as init of mleHomGP are used to define the same optimization problem as above. By default mleHomGP estimates a constant trend β 0 , while it is fixed with known here in order to facilitate an apples-to-apples comparison. x 1 x 2 y mean surface -2 -1 0 1 2 3 4 -2 -1 0 1 2 3 4 variance x1 x2 There are subtle differences between the objective and optimization being performed, which accounts for the slight differences in higher significant digits. The hetGP package offers an S3 predict method that can facilitate prediction exactly in the manner illustrated above; but rather than duplicate Figure  1  here, we shall delay an illustration until later. So to conclude this warm-up, GPs are relatively simple to specify, and inference involves a few dozen lines of code to define an objective (log-likelihood) and its gradient, and a library routine for optimization.
Besides the SK feature offered by the mleGP package, we are not aware of any other software package, for R or otherwise, that automatically preprocesses the data (e. g. , via find_reps) into a form that can exploit these Woodbury identities to speed calculations in the face of heavy replication, or with comparable computational advantage. Replication is a common feature in the sampling of noisy processes, and there is a substantial computational benefit to handling this form of data efficiently.
Implementing heteroskedastic modeling in hetGP The typical GP formulation, utilizing a covariance structure based on Euclidean distance between inputs, emits a stationary process where features in input-output relationships are identical throughout the input space. Many data-generating mechanisms are at odds with the assumption of stationarity. A process can diverge from stationarity (i. e. , be nonstationary) in various ways, but few are well accommodated by computationally viable modeling methodology. Even fewer come with public software, such as tgp  (Gramacy 2007; Gramacy and Taddy 2010)  and laGP  (Gramacy 2016). Both of those packages offer a divide-and-conquer approach to accommodate disparate covariance structures throughout the input space. Such a mechanism is computationally thrifty, but it is not without drawbacks such as lack of continuity of the resulting predictive surfaces.
R> library("MASS") R> hom <-mleHomGP(mcycle$times, mcycle$accel, covtype = "Matern5_2") R> het <-mleHetGP(mcycle$times, mcycle$accel, covtype = "Matern5_2") The code below comprises of our first demonstration of the generic predict method provided by the hetGP package, utilizing a dense grid from the smallest to the largest values in the input space (times). R>  60, length = 301) , ncol = 1) R> p <-predict(x = Xgrid, object = hom) R> p2 <-predict(x = Xgrid, object = het) Figure  2  shows the resulting predictive surface overlaid on a scatter plot of the data. The thick lines are the predictive mean, and the thin lines trace out a 90% predictive interval. The output from predict separates variance in terms of mean (p$sd2) and residual (nugget-based p$nugs) estimates, which we combine in the figure to show the full predictive uncertainty of the response Y (x) | D N.
The overall method, which is effectively on the order of O(T N 4 ) for T MCMC samples, is totally impractical but works well on small problems. Several authors have economized on aspects of this framework  (Kersting, Plagemann, Pfaff, and Burgard 2007; Lazaro-Gredilla and Titsias 2011)  with approximations and simplifications of various sorts, but none (to our knowledge) have resulted in public R software. 1  The key ingredient in these works, of latent variance quantities smoothed by a GP, has merits and can be effective when handled gingerly. The methods implemented in hetGP are built around the following methodology. Let δ 1 , δ 2 ,. , δ n be latent variance variables (or equivalently latent nuggets), each corresponding to one of the n unique design sites xi under study. It is important to introduce latents only for the unique-n locations.
(2018)  describe how to stay within a (Woodbury) MLE framework, by defining a joint log-likelihood over both mean and variance GPs: l = c - N 2 log ν2 N - 1 2 n i=1 [(a i -1) log λ i + log a i ] - 1 2 log |K n | - n 2 log ν(g) - 1 2 log |K (g) |. (9) Maximization may be facilitated via closed-form derivatives with respect to all unknown quantities, all in O(n 3 ) time. For example, the derivative with respect to the latent ∆ n 1 A partial implementation is available for Python via GPy: https://sheffieldml. github. io/GPy/. values follows ∂ l ∂∆ n = ∂Λ n ∂∆ n ∂ l ∂Λ n = C (g) K -1 (g) ∂ℓ ∂Λ n - K -1 (g) ∆ n ν(g) where ∂ℓ ∂λ i = N 2 × (a i -1)σ 2 i λ 2 i + (K -1 n ȳn) 2 i a i νN - a i -1 2λ i - 1 2a i (K n ) -1 i,i Just like in Section 2, implementation is fairly straightforward once this likelihood and derivative are coded.