text:
Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Algorithmic decomposition for efficient multiple nuclear spin detection in diamond

## Abstract

Efficiently detecting and characterizing individual spins in solid-state hosts is an essential step to expand the fields of quantum sensing and quantum information processing. While selective detection and control of a few 13C nuclear spins in diamond have been demonstrated using the electron spin of nitrogen-vacancy (NV) centers, a reliable, efficient, and automatic characterization method is desired. Here, we develop an automated algorithmic method for decomposing spectral data to identify and characterize multiple nuclear spins in diamond. We demonstrate efficient nuclear spin identification and accurate reproduction of hyperfine interaction components for both virtual and experimental nuclear spectroscopy data. We conduct a systematic analysis of this methodology and discuss the range of hyperfine interaction components of each nuclear spin that the method can efficiently detect. The result demonstrates a systematic approach that automatically detects nuclear spins with the aid of computational methods, facilitating the future scalability of devices.

## Introduction

Diamond nitrogen-vacancy (NV) centers have recently attracted much attention in magnetic resonance and related fields of research such as quantum metrology1,2,3,4,5 and quantum information processing6,7,8,9 because of their good quantum coherence. The spin states of NV centers show long coherence times10,11 even at room temperature12,13, can be manipulated by microwaves14,15, can be connected optically over long distances16, and can be operated in a wide range of temperatures through the measurement of optically detected magnetic resonance (ODMR) using spin-dependent fluorescence17. The NV center electron spin is also used as a sensitive probe of surrounding spins (most commonly used are 13C nuclear spins interacting with NV center electron spin by hyperfine coupling) as schematically shown in Fig. 1a. Owing to the long coherence times of these nuclear spins18, detecting and manipulating naturally trapped nuclear spins in diamond through NV centers is emerging as a promising pathway to enable controlled quantum systems such as spin qubit registers19,20,21,22,23,24, quantum networks25,26, quantum simulation platforms18,27, and quantum memories28,29,30. Because naturally occurring 13C nuclear spins are randomly located around an NV center electron spin31,32, the efficient and automatic characterization of hyperfine interactions of individual nuclear spins is a prerequisite step towards building scalable quantum systems based on this platform.

Automatic methods that efficiently and objectively detect and characterize the underlying nuclear spins are required for further scalability of this platform. The detection of nuclear spins in diamond relies on the hyperfine coupling to NV center electron spin, and combined effects from each spin are reflected in the spectrum signal31,32. To identify each spin individually, the signal should be decomposed into multiple spectra, each resulting from a single nuclear spin. Therefore, as the spin system size increases, detecting individual nuclear spin signals from the observed magnetic resonance spectrum requires a more intricate analysis process. Applying for larger-scale spin systems, systematic and reliable approaches, currently lacking in the field of NV center-based nuclear spectrum analysis, are required.

In this study, we propose and demonstrate a computer-aided algorithm for analyzing nuclear spectra in diamond, resolving periodic dips automatically and determining hyperfine interaction tensor components of the nuclear spins weakly coupled to an NV center. Central to the algorithm, we adapted the Gaussian mixture model based on signal partitioning algorithms33 to identify each dip from the signal, enabling analysis of the hyperfine interaction tensor components of each nuclear spin. We demonstrate systematic and automatic decomposition of both virtual and experimental dynamical decoupling-based nuclear spectroscopy data. We also show that the algorithm is capable of detecting the magnitude of hyperfine interaction tensor within the uncertainty of 20%, by comparing the results to those from a more complex and demanding multidimensional spectroscopy method applied to the same nuclear spins23. Our result provides a useful, fast and general tool for analyzing a wide variety of spectroscopy data; this is a key tool for expanding the number of available coherent resources in a system using hyperfine coupling to NV centers and nuclear spins.

## Materials and methods

### The characteristic signal

In NV center-based nuclear magnetic resonance spectroscopy, dynamical decoupling pulse sequences are applied to detect a weak signature of individual nuclear spin signals. Periodic pulses are used to decouple the interaction of the NV center with a spin bath while amplifying specific nuclear Larmor precession resonant with the inter $$\pi$$-pulse period (Fig. 1b)31. In this study, we analyze Carr-Purcell-Meiboom-Gilbert (CPMG)34,35 type dynamical decoupling sequences. We choose this sequence because it forms the basic work horse of spin detection20,21,31,36, is of low experimental complexity and forms the starting point for more extensive characterization methods23,37,38.

The interaction of the NV center with an individual nuclear spin is reflected in sudden dips in the measured electron coherence occurring periodically as a function of total phase evolution time. The analytic expression for the NV center coherence in the presence of nuclear spin is given by31,

$$P_{x} = \frac{M + 1}{2}$$
(1)
$$M = 1 - (1 - \hat{n}_{0} \cdot \hat{n}_{1} )\sin^{2} \frac{N\phi }{2}$$
(2)

where $$P_{x}$$ is the probability that the initial state is preserved, N is the number of the repeated unit decoupling sequence, $$\hat{n}_{0}$$ and $$\hat{n}_{1}$$ are the rotation axes of the nuclear spin depending on the initial state of the electron spin with magnetic quantum number $$m_{s}$$ = 0, − 1, respectively, $$\phi$$ is the net rotation angle around the axis $$\hat{n}_{i}$$ (i = 0,1) following the relation $$cos\phi = \cos \alpha \cos \beta - m_{z} \sin \alpha \sin \beta$$, and $$1 - \hat{n}_{0} \cdot \hat{n}_{1} = m_{x}^{2} \frac{(1 - \cos \alpha )(1 - \cos \beta )}{{1 + \cos \alpha \cos \beta - m_{z} \sin \alpha \sin \beta }}$$, where $$\alpha = \tilde{\omega }\tau$$, $$\beta = \omega_{L} \tau$$, $$m_{z} = (A + \omega_{L} )/\tilde{\omega }$$, $$m_{x} = B/\tilde{\omega }$$, and $$\tilde{\omega } = \sqrt {(A + \omega_{L} )^{2} + B^{2} }$$. For example, Fig. 1c shows numerically simulated spectroscopy data with a single nuclear spin near an NV center. The spectra show sharp dips in the coherence of the electron spin, with a position and amplitude determined by the parallel (perpendicular) hyperfine interaction tensor components A (B).

In a typical dynamical decoupling experiment using a diamond chip with naturally abundant 13C nuclear spins, the spectroscopy signal shows a complicated spectrum stemming from each nuclear spin surrounding the electron spin superposed with a broad background signal. The random location of nuclear spins and the non-linearity of the signal makes stereotypical signal library construction impractical, and the peak detection and analysis should be performed for each given NV center and the uniquely associated nuclear spin bath31,32. These factors forbid direct application of conventional Nuclear Magnetic Resonance (NMR) spectrum analysis methods to NV center-based NMR signal processing.

The purpose of our computer-aided spectrum analysis is to accurately determine individual nuclear spins and obtain both longitudinal and transverse hyperfine parameters automatically. The signal analysis process can be divided into three steps: decomposition of a given signal into Gaussians using signal partitioning algorithm33; extraction of a single spin signal; and calculation of the hyperfine interaction of each spin. Each stage is fully automated and executed in sequential order, using user-defined parameters which are summarized in Table 1.

### Signal decomposition into Gaussians

The first step of this method is to apply the signal partitioning algorithm to decompose the given CPMG-based dynamical decoupling signal into a combination of Gaussian functions. The algorithm first automatically divides the spectroscopy data into fragments, then each fragment is determined by a distinct splitting dip in the data. The data is divided so that nominally only one coherence dip exists in one fragment. All fragments are then decomposed into Gaussians based on Expectation–Maximization (EM) iterations. All Gaussian outputs with an  amplitude smaller than the threshold are eliminated during the post-processing. MATLAB packages: Signal Processing Toolbox, Bioinformatics Toolbox, Curve Fitting Toolbox, and Parallel Computing Toolbox are used for the implementation of the algorithm33.

### Detection of single nuclear spins

Resonance dips stemming from nuclear spins appear in the CPMG signal when the condition $$\hat{n}_{0} \cdot \hat{n}_{1} = - 1$$ is met. In terms of phase evolution time in the kth fragment $$\tau_{k}$$, this condition can be approximated as

$$\tau_{k} = \frac{\pi }{{\tilde{\omega } + \omega_{L} }}\left( {\left( {2k - 1} \right) + \left( { - 1} \right)^{k - 1} \left( {\frac{1}{2\sqrt 2 }\left( {\frac{B}{{\omega_{L} + A}}} \right)^{2} - \frac{1}{4\sqrt 2 }\left( {\frac{B}{{\omega_{L} + A}}} \right)^{4} } \right)} \right) + O\left( {\left( {\frac{B}{{\omega_{L} + A}}} \right)^{6} } \right)$$
(3)

where $$\omega_{L}$$ is the Larmor frequency of the 13C nuclear spins, and $$\tilde{\omega } = \sqrt {(\omega_{L} + A)^{2} + B^{2} }$$. This approximation can be derived from finding the intersection between the curves $$\tan (\alpha /2)\tan (\beta /2) = 1/m_{z}$$ and $$\beta = \frac{{\omega_{L} }}{{\tilde{\omega }}}\alpha$$ in the βα plot, using Taylor expansion with $$\frac{B}{{\omega_{L} + A}}$$, where A(B) is the parallel(perpendicular) hyperfine interaction tensor component, $$\alpha = \tilde{\omega }\tau$$, $$\beta = \omega_{L} \tau$$, and $$m_{z} = (A + \omega_{L} )/\tilde{\omega }$$. When $$B < < \omega_{L}$$, further approximation yields a linear $$\tau_{k}$$ vs. k as follows:

$$\tau_{k} = \frac{(2k - 1)\pi }{{\tilde{\omega } + \omega_{L} }}$$
(4)

Therefore, the nuclear spin is identified by collecting the position of the kth coherence dip as a function of k, and fitting the result with a straight line. In this study, the dip positions from the signal decomposition process are partitioned into time windows, following the period $$T = \frac{\pi }{{\omega_{L} }}$$, and each dip position is marked with the relative time $$\Delta \tau_{k}$$ with respect to the center of the time window $$\Delta \tau_{k} = \tau_{k} - \left( {k - \frac{1}{2}} \right)T$$.

When collected, dip positions $$\Delta \tau_{k}$$ vs. k form a fan diagram, and the CPMG line fit method composed of the dip grouping and a line fitting process is applied to the plot. With the constraint that all lines start from $$(k,\Delta \tau_{k} ) = (0.5,0)$$, satisfying Eq. (4), the objective of this process is to find a set of line configurations that explains the positions of all dips in the $$\Delta \tau - k$$ plot. We use a sequential process to reduce the computational load during the process. First, a set of line candidate lists is built based on the individual dip positions. For each dip located at $$(k_{0} ,\Delta \tau_{0} )$$ in the $$\Delta \tau - k$$ plot, a list of line candidates starting from $$(k,\Delta \tau ) = (0.5,0)$$ and having slope $$\frac{d\tau }{{dk}} = \frac{{\Delta \tau_{0} }}{{\left( {k_{0} - 0.5} \right)}}$$ is generated.

Second, for all line candidates in the list, we find the nearest dips to the line for each k, and calculate the squared distance of τ between the dips and the line. The mean squared distance of each line was calculated using $$\overline{{d^{2} }} = \frac{{\sum\limits_{k = 1}^{{k_{\max } }} {\min \left( {\Delta \tau_{line} (k) - \Delta \tau_{dip} (k)_{\min } ,d_{\max } } \right)}^{2} }}{N}$$, where $$\Delta \tau_{line} (k) = \left( {\frac{d\tau }{{dk}}\left( {k - 0.5} \right) + 0.5} \right)$$ is the position of the line at period number k, and Δτdip(k)min is the position of the dip with the closest distance to the line in consideration. To prevent picking dips too far from the line, a maximum value dmax of the distance between the line and a dip is set. The line candidate having the least mean distance value among the k values is considered as the best fit, and the dip positions used in the calculation are grouped as signals of the same nuclear spin. Third, the examined line candidate as well as a determined group of dip positions are removed from the $$\Delta \tau - k$$ plot. This process is repeated until all the points in the $$\Delta \tau - k$$ plot are exhausted.

During the line selection process where kmax is the maximum time window index within the data aimed to be analyzed, lines with slopes larger than $$T/2k_{\max }$$ (smaller than $$- T/2k_{\max }$$) have k values with 0(2) dips within them. Furthermore, all of the following dips within time windows with k values larger than the window with 0(2) dips show the discrepancy between its time window index k and dip numbering index. To solve this discrepancy, additional dips with coordinates $$(k \pm m,\Delta \tau_{k} \mp mT), \, m < M{, }m = {\text{integer}}$$, where M is the number of additional layers of dips, are added for each dip in position $$(k,\Delta \tau_{k} )$$ in the $$\Delta \tau - k$$ plot. This procedure expands the period window to (2 M − 1)T, and the range of slope that can be fit into a straight line is expanded to $$( - (2M - 1)T/2k_{\max } ,(2M - 1)T/2k_{\max } )$$.

### Fitting and post selection

The full-width-half-maximum (FWHM) of each dip grouped in the CPMG line fit and the slope obtained in the CPMG line fit are used as the two constraints determining (A, B) of each nuclear spin. For the first constraint, the envelope of the CPMG signal represented by a Lorentzian is used,

$$1 - \hat{n}_{0} \cdot \hat{n}_{1} = \frac{2}{{1 + \left( {\frac{{(\tilde{\omega } + \omega_{L} )\tilde{\omega }(\tau - \tau_{k} )}}{B}} \right)^{2} }}$$
(5)

which is valid for $$\tau - \tau_{k} < < \frac{B}{{(\tilde{\omega } + \omega_{L} )\tilde{\omega }}}$$. From the Gaussian peak decomposition, we gain the standard deviation σ of the best Gaussian fit to the dip, giving the FWHM of the fit as $$2\sigma \sqrt {2\ln 2}$$. Using this value as the FWHM of the Lorentzian envelope, the first relation between A, B, and σ can be derived as follows:

$$\sigma = \frac{1}{{\sqrt {2\ln 2} }}\frac{B}{{(\tilde{\omega } + \omega_{L} )\tilde{\omega }}}$$
(6)

The second constraint comes from the slope obtained during the CPMG line fit process from the relation,

$$\frac{d\tau }{{dk}} = \left( {\frac{2\pi }{{\tilde{\omega } + \omega_{L} }} - T} \right)$$
(7)

Using the calculated A and B as the initial guess, the fit function implemented in MATLAB is used to obtain the final estimation of A and B. While using the fit function, we apply a filter on the root-mean-square error (RMSE) calculation to prevent unphysical fitting results due to interception of nearby dips originating from other nuclear spins. The filter is built based on the initial (A, B) values as follows:

$$filter(\tau ) = \left\{ {\begin{array}{*{20}c} {1,} & {\tau \in [\tau_{k} - \sigma ,\tau_{k} + \sigma ]} \\ {0,} & {otherwise} \\ \end{array} } \right.$$
(8)

The fitting process is applied to all lines that were selected in the CPMG line fit method. However, some of the chosen lines can be based on mixed dips that come from several different nuclear spins with differing interaction strength. To remove the incorrectly selected lines, a post-selection process is applied to the fitted (A, B) pairs. The (A, B) pairs that come out from the fitting process are regarded as candidates for a good nuclear spin interaction representation. The Beam Search39, a heuristic search strategy for memory saving based on the breadth-first search, is applied to find the optimal (A, B) pair configuration, which reconstructs the most similar signal compared to the given CPMG signal. The algorithm calculates the RMSE value with the filter for each (A, B) pair configuration similar to the fitting process.

## Results and discussion

We first build up our systematic approach with an example of virtually generated CPMG signals. We simulate the CPMG signal with randomly selected hyperfine interaction tensor components of ten individual nuclear spins within 2 nm radius from the NV center for a given $$\pi$$ pulse repetition number N = 32. Approximately 50,000 random coordinates, which represent nuclear spins more than 2 nm away from the electron spin, in the diamond lattice having interaction strength of A < 8 kHz, B < 0.25 kHz were chosen to reproduce the spin bath signal (Fig. 2a). Figure 2b shows the Gaussian decomposition33 of the input signal where the position, amplitude, and standard deviation value of each Gaussian are obtained using the algorithm described in the method section. Figure 2c shows the plot of the positions of the dips in the ∆τ- k fan diagram. As explained above, the automatic CPMG line fit method is applied to obtain the hyperfine interaction, which finds the best configuration of straight lines starting near origin covering most of the data points as shown in Fig. 2c, using solid lines. For each line found from the fit, A and B values are calculated by solving Eqs. (4), and (6). These hyperfine parameters are used as initial guesses to fit each dip iteratively based on Eq. (2) to reach the final hyperfine interaction tensor components. The procedure is applied to each detected nuclear spin, and the accuracy of the fit is determined by comparing the reconstructed and input CPMG signal. As shown in Fig. 2d, the reconstructed CPMG signal, using the converged hyperfine interaction tensor components, shows excellent agreement with the input configuration of the spins. This shows less than 5% error compared to the values used in the signal generation. It also demonstrates the capability of the algorithm to automatically search the given nuclear spectra. The error of hyperfine interaction was calculated as,

$$Error = \frac{{\sqrt {(A_{orig} - A_{obt} )^{2} + (B_{orig} - B_{obt} )^{2} } }}{{\sqrt {A_{orig}^{2} + B_{orig}^{2} } }}$$
(9)

where $$A_{orig}$$,$$B_{orig}$$ are the (A, B) value of the original hyperfine interaction, and $$A_{obt}$$,$$B_{obt}$$ are the (A, B) value of the obtained hyperfine interaction, the result of our program.

We turn to discuss more systematic analyses of the performance of the developed algorithm. We first test the performance with numerically generated CPMG spectrums of a single nuclear spin with varying A and B. Figure 3a shows the error of the hyperfine interaction strength from a single nuclear spin signal. The small A domain ($$\left| A \right| < {\text{5 kHz}}$$) was not detectable because of the overlap of the coherence dip signal with virtually generated spin bath and other nuclear spin signals with similar dip periods. As A becomes larger in the large A domain ($$\left| A \right| > 70{\text{ kHz}}$$), the slope of the line in ∆τ − k plot becomes steeper. In the spin detection stage, this program considers only two additional clone layers (M = 3) of the point set while applying the CPMG line fit method to the ∆τ − k plot to ensure a manageable computation time. As a result, the range of detectable slopes is limited, and a step jump of the error in the large A domain is observed. Additionally, the small B domain (B < 5 kHz) was not detectable because of the small size of the dips (smaller than the size of the threshold 0.05 given during the Signal decomposition process), and the large B domain (B > 120 kHz) was hardly detectable because the shape of each of the dips deviates from a simple Lorentzian form to a more complex shape with several ripples around the dip typical for a CPMG dynamical decoupling sequence. We expect that more advanced techniques such as pattern analysis based on machine learning combined with numerical fits can enhance the detectability in this strongly coupled regime.

Since the algorithm is based on a nuclear spectrum, stemming from all nuclear spins interacting with the NV center together, it is important to determine the minimum resolution of the hyperfine coupling parameter that the algorithm can reliably distinguish. We analyze simulated CPMG signal, originating from two different spins to estimate the resolution of the dip frequency $$f_{p} = \frac{{\tilde{\omega } + \omega_{L} }}{4\pi }$$. One nuclear spin with fixed hyperfine interaction tensor components was used as the reference signal while the hyperfine coupling of the second nuclear spin was varied. Figure 3b,c show the error of the obtained hyperfine parameters of the second nuclear spin as a function of parameter variation from the reference spin. The results show that the main factor significantly affecting the resolution of the program is A since it is directly related to the slope of each line in the spin detecting stage as well as the dip frequency of the nuclear spin. We empirically find that the obtained dip frequency resolution of the current algorithm is about 2 kHz, which showed an error less than 20%.

We further analyze the accuracy and applicability of this program by repeatedly testing the algorithm with virtual CPMG data sets containing multiple random nuclear spins, covering a wide range of A and B. The values are based on possible lattice points in diamond and calculating the hyperfine interaction from a purely dipole–dipole coupling—neglecting the contact hyperfine and extended wavefunctions for simplicity—between the 13C nuclear spin and the electron spin of the NV center using,

$$H = - \frac{{\mu_{0} \gamma_{e} \gamma_{n} \hbar^{2} }}{{4\pi \left| {\mathbf{r}} \right|^{3} }}(3({\mathbf{S}}_{e} \cdot {\hat{\mathbf{r}}})({\mathbf{I}}_{n} \cdot {\hat{\mathbf{r}}}) - {\mathbf{S}}_{e} \cdot {\mathbf{I}}_{n} )$$
(10)

where $$\mu_{0} = 4\pi \times 10^{ - 7}$$ H/m is the magnetic permeability of free space, $$\gamma_{e}$$ = 20,824 MHz/T ($$\gamma_{n}$$ = 10.708 MHz/T ) is the gyromagnetic ratio of the NV center electron (13C nuclear ) spin, $${\mathbf{S}}_{e}$$ is the vector of electron spin-1 operators, and $${\mathbf{I}}_{n}$$ is the nuclear spin-1/2 operators vector. We confine the lattice points to be within 2 nm from the NV center, and consider spins having hyperfine interaction values within the range − 100 kHz < A < 100 kHz (5 kHz < B < 100 kHz). Twenty lattice points in this range were randomly chosen to simulate the signal, and by repeating our algorithm, we record the error of the analysis result for each nuclear spin used to generate the signal.

Figure 3d shows the average error between the input hyperfine coefficients and the output of this method by accumulating 1,000 configurations of 20 random spins. From this result, we find that the region of 80% confidence is 5 kHz <|A|< 70 kHz and 15 kHz <|B|< 80 kHz. The region of both |A| and |B| is constrained from the limitation of detecting a single nuclear spin signal as shown in Fig. 3a. In addition, for |B|, the overlap of the dips coming from different nuclear spins appears in the many spin signals, narrowing the allowed range of this parameter. The loss in confidence can be ascribed to two main factors; (1) the accuracy of the period of the dips, and (2) the accuracy of estimating the amplitude of each dip. The period of each nuclear spin signal is estimated using a large collection of dip points. This is generally robust to the presence of other nuclear spin signals, leading to frequency resolution 2 kHz (Fig. 3e). However, the amplitude of each dip can be disturbed when overlapped with other spin signals. As a result, the error of the hyperfine parameters mainly comes from the error in fitting the amplitude of the dip, affecting the accuracy of the perpendicular hyperfine interaction element B. As shown in Fig. 3b,c, when the dip frequency coming from two different nuclear spins is similar within 2 kHz, the dips coming from different spins interfere with each other. This brings in a non-linear distortion to the dip amplitude. It is found that among spins in the 80% confidence region, an average of 20% of the spins has another nuclear spin whose dip frequency difference is within 2 kHz, making it difficult to fit the signal from each nuclear spin independently (Fig. 3f). The effect of the interfering nuclear spin signal can be approximated from Eq. (2). When some error in B comes in, represented as $$\delta B$$, the size of the dip is $$P = \sin^{2} \left( {\frac{N}{2}\frac{B}{{\sqrt {(\omega_{L} + A)^{2} + B^{2} } }}} \right)$$, and the error of the dip becomes $$\frac{\delta P}{P}\sim 2\frac{\delta B}{B}$$ up to the order of $$O\left( {\left( {\frac{B}{{\omega_{L} + A}}} \right)^{4} } \right)$$, for A, BωL. This result shows that the error in the amplitude of the dip P directly affects the error of B. From this relation, it can be interpreted that if an additional spin signal with its dip size of 0.2 is overlapped with the original signal, the output B value is expected to show 20% error.

We now show the application of the method to the experimental NV center-based nuclear spectroscopy data. The data was collected from a single NV center system at a low temperature with varying τ from 6 to longer than 50 μs11. The experimental data in ref. 11 were re-analyzed, and our algorithm reports 14 distinct nuclear spins (Fig. 4a). As an output of the algorithm, the hyperfine interaction tensor components of each spin were obtained. The simulated signal from these hyperfine parameters matched well with the experimental data, as shown in Fig. 4b. The obtained hyperfine interaction tensor components are given in Table 2. Based on the accuracy analysis performed above, spins 1, 2, 4, 10, 11, 12, and 14, marked with asterisk in the table, are in the reliable (A, B) region, and as expected, they show values less than 20% difference compared to the reported values from the experimental multidimensional spectroscopy on the same NV center and nuclear spin environment23 (Table 2).

This algorithm reported additional spins 3, 5, 6, 7, 8, 9 and 13, whose (A, B) values are out of the 80% confidence range. Spin 3 and 13 have their consistent match among the values of ref. 23, even though their B values were below the confidence threshold. Spin 5–9 show dip periods similar to the spin bath and overlap with each other as shown in Fig. 4c. For longer τ, the coherence dips for spin 5–9 become distinguishable from one another as shown in Fig. 4d. The reported hyperfine values of spin 5, 7, 8 and 9 show large deviations from all values in ref. 23, and spin 6 has a consistent match in ref. 23. While the algorithm treats each line/dip position group as coming from a single nuclear spin, due to limited spectral resolution of the CPMG data, it is highly likely that each spin label contains more than one nuclear spin. In the small A (< 5 kHz) regime, only the coherence dip period information, but not the individual A and B values, is reliable. Nevertheless, to explain the experimental coherence values significantly below 1 at the location of each of the dips coming from spin 5–9, at least one nuclear spin is necessary per coherence dip period. Overall, for the reported spins with (A, B) values out of the 80% confidence range, the dip period information of the spins can aid to set a spectral region to focus on when adapting more advanced spin spectroscopy methods with higher spectral resolution.

In this algorithm, the CPMG signal was simulated and analyzed based on Eqs. (1) and (2). For the purpose of extracting the hyperfine interaction coefficients A and B, these equations can be applied to spin-1/2 systems within dipolar hyperfine interaction approximation. Therefore, the application of this algorithm is not limited to the NV center-13C coupling in the natural diamond environment which we demonstrated. We expect that applying this algorithm to different spin environments, such as in silicon carbide, would be possible.

In these equations, the interaction between the 13C nuclear spins was neglected because of the small interaction strength compared to the one between the NV center electron spin and the 13C nuclear spin. However, the nuclear-nuclear couplings can be significant, especially for larger τ and N, and can cause broadening or splitting of the dips11,23. This could impair the analytical performance of this program. Additional consideration of the interaction between 13C nuclear spins in the fitting stage will increase the capability of this program even for results obtained from the large time region.

Additionally, under the same nuclear spin environment, the CPMG signal of larger N values shows more disordered results because the number of nuclear spins with detectable dips within the CPMG signal increase, elevating the probability of dip overlap between different nuclear spins in the resulting signal. Similar problems can arise in analyzing the CPMG signal obtained under an environment with13C concentration higher than natural abundance. In this algorithm, it is assumed that a single straight line in the spin detection stage represents a single nuclear spin. However, it is possible that multiple nuclear spins with different hyperfine interaction tensor components share their dip frequency, and form a single line as implied in Fig. 3b,c. In this case, because the multiple spin signal is considered as a signal originating from a single nuclear spin in this algorithm, the fitting accuracy can be reduced. Since treating signals from a single straight line as a multiple spin signal takes more complex calculations and a longer time, it was not considered in this algorithm. Considering multiple nuclear spins in the fitting stage will increase the accuracy of the analysis.

These limitations bring in difficulty in analyzing CPMG signals coming from sequences with a large pulse repetition number. Nevertheless, the CPMG signals with large N values are useful for detecting weakly coupled nuclear spins because the size of dips in the CPMG signal coming from a nuclear spin with small perpendicular hyperfine interaction strength becomes magnified, and the coherence time of the NV center increases. Further improvements in solving the non-linearity (for example, a combination of our numerical fit procedure with a more advanced machine learning approach, and including the nuclear spin-nuclear spin interactions to the model) should improve the general performance of our algorithm in the future.

## Conclusion

In this study, we have proposed and demonstrated an automated algorithmic method to analyze the CPMG dynamical decoupling nuclear spectrum to detect 13C nuclear spins near the NV center in diamond. Our algorithm decomposes spectroscopy signals into Gaussians using the signal partitioning algorithm33, detects nuclear spins automatically, and returns hyperfine interaction tensor components of each nuclear spin. We confirmed that our program successfully works for both simulated signals and an extended experimental data set. Moreover, we analyzed the range of applicability and limitations of this methodology, which can be improved with additional consideration of the interaction between nuclear spins, and signal processing regarding non-linearity due to multiple nuclear spins. With the possibility to be generalized to other spin-1/2 systems, our results show the first step of a systematic algorithmic approach, providing a useful and general tool for sensing complex spin structures, and furthermore, expanding quantum systems based on spin qubits.

## Data availability

All data generated within this study are available from the corresponding author on request.

## References

1. 1.

Unden, T. et al. Quantum metrology enhanced by repetitive quantum error correction. Phys. Rev. Lett. 116, 230502 (2016).

2. 2.

Zaiser, S. et al. Enhancing quantum sensing sensitivity by a quantum memory. Nat. Commun. 7, 12279 (2016).

3. 3.

Degen, C. L., Reinhard, F. & Cappellaro, P. Quantum sensing. Rev. Mod. Phys. 89, 035002 (2017).

4. 4.

Rosskopf, T., Zopes, J., Boss, J. M. & Degen, C. L. A quantum spectrum analyzer enhanced by a nuclear spin memory. NPJ Quantum Inf. 3, 33 (2017).

5. 5.

Kim, D. et al. A CMOS-integrated quantum sensor based on nitrogen–vacancy centres. Nat. Electron. 2, 284–289 (2019).

6. 6.

Wrachtrup, J. & Jelezko, F. Processing quantum information in diamond. J. Phys.: Condens. Matter 18, S807 (2006).

7. 7.

Doherty, M. W. et al. The nitrogen-vacancy colour centre in diamond. Phys. Rep. 528, 1–45 (2013).

8. 8.

Waldherr, G. et al. Quantum error correction in a solid-state hybrid spin register. Nature 506, 204 (2014).

9. 9.

Cramer, J. et al. Repeated quantum error correction on a continuously encoded qubit by real-time feedback. Nat. Commun. 7, 1–7 (2016).

10. 10.

Bar-Gill, N., Pham, L. M., Jarmola, A., Budker, D. & Walsworth, R. L. Solid-state electronic spin coherence time approaching one second. Nat. Commun. 4, 1–6 (2013).

11. 11.

Abobeih, M. H. et al. One-second coherence for a single electron spin coupled to a multi-qubit nuclear-spin environment. Nat. Commun. 9, 2552 (2018).

12. 12.

Balasubramanian, G. et al. Ultralong spin coherence time in isotopically engineered diamond. Nat. Mater. 8, 383 (2009).

13. 13.

Herbschleb, E. et al. Ultra-long coherence times amongst room-temperature solid-state spins. Nat. Commun. 10, 1–6 (2019).

14. 14.

Hanson, R., Dobrovitski, V., Feiguin, A., Gywat, O. & Awschalom, D. Coherent dynamics of a single spin interacting with an adjustable spin bath. Science 320, 352–355 (2008).

15. 15.

Tamarat, P. et al. Spin-flip and spin-conserving optical transitions of the nitrogen-vacancy centre in diamond. New J. Phys. 10, 045004 (2008).

16. 16.

Hensen, B. et al. Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres. Nature 526, 682–686 (2015).

17. 17.

Taylor, J. et al. High-sensitivity diamond magnetometer with nanoscale resolution. Nat. Phys. 4, 810 (2008).

18. 18.

Unden, T. et al. Coherent control of solid state nuclear spin nano-ensembles. NPJ Quantum Inf. 4, 39 (2018).

19. 19.

Dutt, M. G. et al. Quantum register based on individual electronic and nuclear spin qubits in diamond. Science 316, 1312–1316 (2007).

20. 20.

Zhao, N. et al. Sensing single remote nuclear spins. Nat. Nanotechnol. 7, 657 (2012).

21. 21.

Shi, F. et al. Sensing and atomic-scale structure analysis of single nuclear-spin clusters in diamond. Nat. Phys. 10, 21 (2014).

22. 22.

Liu, G.-Q. et al. Single-shot readout of a nuclear spin weakly coupled to a nitrogen-vacancy center at room temperature. Phys. Rev. Lett. 118, 150504 (2017).

23. 23.

Abobeih, M. et al. Atomic-scale imaging of a 27-nuclear-spin cluster using a quantum sensor. Nature 576, 411–415 (2019).

24. 24.

Bradley, C. et al. A ten-qubit solid-state spin register with quantum memory up to one minute. Phys. Rev. X 9, 031045 (2019).

25. 25.

Nemoto, K. et al. Photonic quantum networks formed from NV—centers. Sci. Rep. 6, 26284 (2016).

26. 26.

Humphreys, P. C. et al. Deterministic delivery of remote entanglement on a quantum network. Nature 558, 268 (2018).

27. 27.

Ju, C. et al. NV-center-based digital quantum simulation of a quantum phase transition in topological insulators. Phys. Rev. B 89, 045432 (2014).

28. 28.

Fuchs, G., Burkard, G., Klimov, P. & Awschalom, D. A quantum memory intrinsic to single nitrogen–vacancy centres in diamond. Nat. Phys. 7, 789 (2011).

29. 29.

Poem, E. et al. Broadband noise-free optical quantum memory with neutral nitrogen-vacancy centers in diamond. Phys. Rev. B 91, 205108 (2015).

30. 30.

Lai, Y.-Y., Lin, G.-D., Twamley, J. & Goan, H.-S. Single-nitrogen-vacancy-center quantum memory for a superconducting flux qubit mediated by a ferromagnet. Phys. Rev. A 97, 052303 (2018).

31. 31.

Taminiau, T. H. et al. Detection and control of individual nuclear spins using a weakly coupled electron spin. Phys. Rev. Lett. 109, 5 (2012).

32. 32.

Nizovtsev, A. P. et al. Non-flipping 13C spins near an NV center in diamond: hyperfine and spatial characteristics by density functional theory simulation of the C510 [NV] H252 cluster. New J. Phys. 20, 023022 (2018).

33. 33.

Polanski, A., Marczyk, M., Pietrowska, M., Widlak, P. & Polanska, J. Signal partitioning algorithm for highly efficient gaussian mixture modeling in mass spectrometry. PLoS ONE 10, e0134256 (2015).

34. 34.

Carr, H. Y. & Purcell, E. M. Effects of diffusion on free precession in nuclear magnetic resonance experiments. Phys. Rev. 94, 630 (1954).

35. 35.

Meiboom, S. & Gill, D. Modified spin-echo method for measuring nuclear relaxation times. Rev. Sci. Instrum. 29, 688–691 (1958).

36. 36.

Kolkowitz, S., Unterreithmeier, Q. P., Bennett, S. D. & Lukin, M. D. Sensing distant nuclear spins with a single electron spin. Phys. Rev. Lett. 109, 137601 (2012).

37. 37.

Zopes, J. et al. Three-dimensional localization spectroscopy of individual nuclear spins with sub-Angstrom resolution. Nat. Commun. 9, 1–8 (2018).

38. 38.

Zopes, J., Herb, K., Cujia, K. & Degen, C. L. Three-dimensional nuclear spin positioning using coherent radio-frequency control. Phys. Rev. Lett. 121, 170801 (2018).

39. 39.

Lowerre, B. T. The Harpy Speech Recognition System. Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA, USA (1976).

## Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIT) (No.2018R1A2A3075438, No.2019M3E4A1080144, No.2019M3E4A1080145) and the Creative-Pioneering Researchers Program through Seoul National University (SNU). This work was supported by the Netherlands Organisation for Scientific Research (NWO/OCW) through a Vidi grant, and as part of the Quantum Software Consortium programme (Project No. 024.003.037/3368). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 852410). This project (QIA) has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 820445.

## Author information

Authors

### Contributions

H.O. and J.Y. designed and implemented the program. K.J. and K.K. provided technical help. M.H.A. and T.H.T. provided theoretical background and experimental data. D.K. conceived and supervised the project. H.O, J.Y., and D.K. wrote the manuscript.

### Corresponding author

Correspondence to Dohun Kim.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

### Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Oh, H., Yun, J., Abobeih, M.H. et al. Algorithmic decomposition for efficient multiple nuclear spin detection in diamond. Sci Rep 10, 14884 (2020). https://doi.org/10.1038/s41598-020-71339-6

• Accepted:

• Published:

• ### Deep learning enhanced individual nuclear-spin detection

• Kyunghoon Jung
• , M. H. Abobeih
• , Jiwon Yun
• , Gyeonghun Kim
• , Hyunseok Oh
• , Ang Henry
• , T. H. Taminiau
•  & Dohun Kim

npj Quantum Information (2021)


text:
# Article

Full entry | Fulltext not available (moving wall 24 months)
Keywords:
multiplier ideal; monomial ideal; convex set
Summary:
Let $\mathfrak {a}\subseteq {\mathbb C}[x_1,\ldots ,x_n]$ be a monomial ideal and ${\mathcal J}(\mathfrak {a}^c)$ the multiplier ideal of $\mathfrak {a}$ with coefficient $c$. Then ${\mathcal J}(\mathfrak {a}^c)$ is also a monomial ideal of ${\mathbb C}[x_1,\ldots ,x_n]$, and the equality ${\mathcal J}(\mathfrak {a}^c)=\mathfrak {a}$ implies that $0<c<n+1$. We mainly discuss the problem when ${\mathcal J}(\mathfrak {a})=\mathfrak {a}$ or ${\mathcal J}(\mathfrak {a}^{n+1-\varepsilon })=\mathfrak {a}$ for all $0<\varepsilon <1$. It is proved that if ${\mathcal J}(\mathfrak {a})=\mathfrak {a}$ then $\mathfrak {a}$ is principal, and if ${\mathcal J}(\mathfrak {a}^{n+1-\varepsilon })=\mathfrak {a}$ holds for all $0<\varepsilon <1$ then $\mathfrak {a}=(x_1,\ldots ,x_n)$. One global result is also obtained. Let $\tilde {\frak {a}}$ be the ideal sheaf on ${\mathbb P}^{n-1}$ associated with $\frak {a}$. Then it is proved that the equality ${\mathcal J}(\tilde {\mathfrak {a}})=\tilde {\mathfrak {a}}$ implies that $\tilde {\mathfrak {a}}$ is principal.
References:
[1] Blickle, M.: Multiplier ideals and modules on toric varieties. Math. Z. 248 (2004), 113-121. DOI 10.1007/s00209-004-0655-y | MR 2092724 | Zbl 1061.14055
[2] Blickle, M., Lazarsfeld, R.: An informal introduction to multiplier ideals. Trends in Commutative Algebra. Mathematical Sciences Research Institute Publications 51 Cambridge University Press, Cambridge (2004), 87-114 L. L. Avramov et al. MR 2132649 | Zbl 1084.14015
[3] Demailly, J.-P., Ein, L., Lazarsfeld, R.: A subadditivity property of multiplier ideals. Mich. Math. J. 48 (2000), 137-156. DOI 10.1307/mmj/1030132712 | MR 1786484 | Zbl 1077.14516
[4] Eisenbud, D.: Commutative Algebra with a View Toward Algebraic Geometry. Graduate Texts in Mathematics 150 Springer, Berlin (1995). MR 1322960 | Zbl 0819.13001
[5] Esnault, H., Viehweg, E.: Lectures on Vanishing Theorems. DMV Seminar 20 Birkhäuser, Basel (1992). MR 1193913 | Zbl 0779.14003
[6] Fulton, W.: Introduction to Toric Varieties. Annals of Mathematics Studies 131 Princeton University Press, Princeton (1993). MR 1234037 | Zbl 0813.14039
[7] Hara, N., Yoshida, K.-I.: A generalization of tight closure and multiplier ideals. Trans. Am. Math. Soc. 355 (2003), 3143-3174. DOI 10.1090/S0002-9947-03-03285-9 | MR 1974679 | Zbl 1028.13003
[8] Hartshorne, R.: Algebraic Geometry. Graduate Texts in Mathematics 52 Springer, New York (1977). MR 0463157 | Zbl 0367.14001
[9] Howald, J. A.: Multiplier ideals of monomial ideals. Trans. Am. Math. Soc. 353 (2001), 2665-2671. DOI 10.1090/S0002-9947-01-02720-9 | MR 1828466 | Zbl 0979.13026
[10] Hübl, R., Swanson, I.: Adjoints of ideals. Mich. Math. J. 57 (2008), 447-462. DOI 10.1307/mmj/1220879418 | MR 2492462 | Zbl 1180.13005
[11] Lazarsfeld, R.: Positivity in Algebraic Geometry I\kern-1ptI. Positivity for Vector Bundles, and Multiplier Ideals. Ergebnisse der Mathematik und ihrer Grenzgebiete 3. Folge {\it 49} Springer, Berlin (2004). MR 2095471
[12] Lipman, J.: Adjoints and polars of simple complete ideals in two-dimensional regular local rings. Bull. Soc. Math. Belg., Sér. A 45 (1993), 223-244. MR 1316244 | Zbl 0796.13020
[13] Nadel, A. M.: Multiplier ideal sheaves and Kähler-Einstein metrics of positive scalar curvature. Ann. Math. (2) 132 (1990), 549-596. MR 1078269 | Zbl 0731.53063
[14] Siu, Y.-T.: Multiplier ideal sheaves in complex and algebraic geometry. Sci. China, Ser. A 48 (2005), 1-31. DOI 10.1007/BF02884693 | MR 2156488 | Zbl 1131.32010
[15] Swanson, I., Huneke, C.: Integral Closure of Ideals, Rings, and Modules. London Mathematical Society Lecture Note Series 336 Cambridge University Press, Cambridge (2006). MR 2266432 | Zbl 1117.13001

Partner of

format: plain


text:
Last time, we discussed the Fundamental Theorem of Dynamic Programming, which then lead to the efficient “value iteration” algorithm for finding the optimal value function. And then we could find the optimal policy by greedifying w.r.t. the optimal value function. In this lecture we will do two things:

1. Elaborate more on the the properties of value iteration as a way of obtaining near-optimal policies;
2. Discuss the computational complexity of planning in finite MDPs.

## Finding a Near-Optimal Policy using Value Iteration

In the previous lecture we found that the iterative computation that starts with some $v_0\in \mathbb{R}^{\mathrm{S}}$ and then obtains $v_{k+1}$ using the “Bellman update”

\begin{align} v_{k+1} = T v_k \label{eq:vi} \end{align}

leads to a sequence $$\{v_k\}_{k\ge 0}$$ whose $$k$$th term approaches $$v^*$$, the optimal value function, at a geometric rate:

\begin{align} \| v_k - v^* \|_\infty \le \gamma^k \| v_0 - v^* \|_\infty\,. \label{eq:vierror} \end{align}

While this is reassuring, our primary goal is to obtain an optimal, or at least a near-optimal policy. Since any policy that is greedy with respect to (w.r.t) $$v^*$$ is optimal, a natural idea is to stop the value iteration after some finite number of iteration steps and return a policy that is greedy w.r.t. the approximation of $$v^*$$ that was just obtained. If we stop the process after the $$k$$th step, this defines a policy $$\pi_k$$ such that $$\pi_k$$ is greedy w.r.t. $$v_k$$: $$T_{\pi_k} v_k = T v_k$$. The hope is that as $$v_k$$ approaches $$v^*$$, the policies $$\{\pi_k\}$$ will also get better in the sense that $$\|v^*-v^{\pi_k}\|_\infty$$ decreases.

The next theorem guarantees that this will indeed be the case.

Theorem (Policy Error Bound): Let $$v: \mathcal{S} \to \mathbb{R}$$ be arbitrary and $$\pi$$ be the greedy policy w.r.t. $$v$$: $$T_\pi v = T v$$. Then,

$v^\pi \geq v^* - \frac{2 \gamma \|v^*-v\|_\infty}{1 - \gamma} \boldsymbol{1}.$

In words, the theorem states that the policy error ($$\|v^*-v^{\pi}\|_\infty$$) of a policy that is greedy with respect to a function $$v$$ is controlled by the distance of $v$ to $$v^*$$. Let $\Gamma$ be an operator that maps functions over $\mathcal{S}$ to policies by assigning to a function $v$ the policy that is greedy w.r.t. $v$. With this, a shorter form of the theorem statement is that for any $v\in \mathbb{R}^{\mathcal{S}}$,

$\| v^*-v^{\Gamma(v)}\|_\infty \le \frac{2\gamma \|v^*-v\|_\infty}{1-\gamma}\,.$

Alternatively, the theorem can be seen as stating that the map $$v \mapsto \|v^*-v^{\Gamma(v)}\|_\infty$$ is “$$2\gamma/(1-\gamma)$$-smooth” with respect to the maximum norm distance at $$v=v^*$$.

The proof is an archetypical example of proofs of using contraction and monotonicity arguments to prove error bounds. We will see variations of this proof many times. Before the proof, let us introduce the notation $|x|$ for a vector $\mathbb{R}^d$ to mean the componentwise absolute value of the vector: $$|x|_i = |x_i|$$, $$i\in [d]$$.

As a way of using this notation, note that for any memoryless policy $\pi$, \begin{align} |P_\pi x |\le P_\pi |x| \le \|x\|_\infty P_\pi \boldsymbol{1} = \|x\|_\infty \boldsymbol{1}\,, \label{eq:ppineb} \end{align} and hence \begin{align} \|P_\pi x \|_\infty \le \|x\|_\infty\,. \label{eq:stochmxne} \end{align} In Eq. $$\eqref{eq:ppineb}$$ the first inequality follows because $P_\pi$ is monotone and $$x\le |x| \le \|x\|_\infty \boldsymbol{1}$$. For the proof it will also be useful to recall that we also have

\begin{align} T_\pi (v+c \boldsymbol{1}) &= T_\pi v \,\, + c \gamma \boldsymbol{1}\,, \label{eq:tpiadd1} \\ T (v+c \boldsymbol{1}) &= T v \,\, + c \gamma \boldsymbol{1}\,, \label{eq:tadd1} \end{align}

for any $$v\in \mathbb{R}^{\mathrm{S}}$$, $$c\in \mathbb{R}$$ and memoryless policy $$\pi$$. These two identities follow just by the definitions of $T$ and $T_\pi$, as the reader can easily verify them.

Proof: Let $$v,v^*,\pi$$ be as in the theorem statement and let $$\varepsilon = \|v^*-v\|_\infty$$. Let $$\delta = v^*-v^\pi$$. The result follows by algebra once we prove that $$\|\delta\|_\infty \le \gamma \|\delta\|_\infty + 2\gamma \varepsilon$$. Hence, we only need to prove this inequality.

By our assumptions on $$v$$ and $$v^*$$, $$-\varepsilon\boldsymbol{1}\le v^*-v \le \varepsilon\boldsymbol{1}$$. Now,

\begin{align*} \delta & = v^*-v^\pi \\ & = \textcolor{red}{T} v^* - \textcolor{red}{T_\pi} v^\pi & \text{(Fundamental Theorem, T_\pi v^\pi = v^\pi)}\\ & \le T(v+\textcolor{red}{\varepsilon\boldsymbol{1}})-T_\pi v^\pi & \text{(T monotone)}\\ & = Tv-T_\pi v^\pi +\textcolor{red}{\gamma\varepsilon\boldsymbol{1}} & \text{(Eq. \eqref{eq:tadd1})}\\ & = \textcolor{red}{T_\pi} v-T_\pi v^\pi +\gamma\varepsilon\boldsymbol{1} & \text{(\pi def.)}\\ & \le T_\pi(v^*+\textcolor{red}{\varepsilon\boldsymbol{1}})-T_\pi v^\pi + \gamma \varepsilon \boldsymbol{1} & \text{(T_\pi monotone)}\\ & = T_\pi v^* - T_\pi v^\pi + \textcolor{red}{2}\gamma \varepsilon\boldsymbol{1} & \text{(Eq. \eqref{eq:tpiadd1})}\\ & = \textcolor{red}{\gamma P_\pi}(v^*-v^\pi)+2\gamma \varepsilon\boldsymbol{1} & \text{(T_\pi def.)}\\ & = \gamma P_\pi \textcolor{red}{\delta}+2\gamma \varepsilon\boldsymbol{1}\,. & \text{(\delta def.)} \end{align*}

Taking the (pointwise) absolute value of both sides and using the triangle inequality, taking and Eq. $$\eqref{eq:stochmxne}$$ we find that \begin{align*} |\delta| \le \gamma \|\delta\|_\infty \boldsymbol{1} + 2\gamma \varepsilon\boldsymbol{1}\,. \end{align*} The proof is finished by taking the maximum over the components, noting that $$\max_s |\delta|_s = \|\delta\|_\infty$$. $$\qquad \blacksquare$$

An alternative way of finishing the proof is to note that from $\delta = \gamma P_\pi \delta + 2\gamma \varepsilon \boldsymbol{1}$, by reordering and using that $(I-\gamma P_\pi)^{-1} = \sum_{i\ge 0} \gamma^i P_\pi^i$ is a monotone operator, $\delta \le 2\gamma \varepsilon \sum_{i\ge 0} \gamma^i P_\pi \boldsymbol{1} = 2\gamma \varepsilon/(1-\gamma) \boldsymbol{1}$. Taking the max-norm of both sides, we get $$\|\delta\|_\infty \le 2\gamma \varepsilon/(1-\gamma)$$.

## Value Iteration as an Approximate Planning Algorithm

From Eq. $$\eqref{eq:vierror}$$ we see for $$k \geq H_{\gamma, \varepsilon} = \frac{\ln(1 / (\varepsilon (1 - \gamma)))}{1 - \gamma}$$, started with $$v_0 =0$$, value iteration yields $$v_k$$ such that $$\|v_k - v^*\|_\infty \leq \varepsilon$$ and consequently, for a policy $$\pi_k$$ that is greedy w.r.t. $$v_k$$, $$v^{\pi_k} \geq v^* - \frac{2 \gamma \varepsilon}{1 - \gamma} \boldsymbol{1}$$. Now, for a fixed $\delta>0$ setting $\varepsilon$ so that $$\delta = \frac{2 \gamma \varepsilon}{1 - \gamma}$$ holds, we see that after $$k \geq H_{\gamma, \frac{\delta(1 - \gamma)}{2\gamma}}$$ iterations, we get a $$\delta$$-optimal policy $$\pi_k$$: $$v^{\pi_k} \geq v^* - \delta \boldsymbol{1}$$. Computing $$v_{k+1}$$ using $$\eqref{eq:vi}$$ takes $$O(\mathrm{S}^2 \mathrm{A})$$ elementary arithmetic (and logic) operations. Putting things together we get the following result:

Theorem (Runtime of Approximate Planning with Value Iteration): Fix a finite discounted MDP and a target accuracy $\delta>0$. Then, after

$O \left(\mathrm{S}^2 \mathrm{A} H_{\gamma, \frac{\delta(1 - \gamma)}{2\gamma}} \right) = \tilde O\left( \frac{\mathrm{S}^2 \mathrm{A} }{1 - \gamma}\, \ln\left(\frac{1}{\delta}\right)\right)$

elementary arithmetic operations, value iteration produces a policy $\pi$ that is $\delta$-optimal: $v^\pi \ge v^* - \delta \boldsymbol{1}$, where the $\tilde{O}(\cdot)$ result holds when $\delta \le 1/e$ is fixed and $\tilde{O}(\cdot)$ hides $\log(2/(1-\gamma))$.

Note that the number of operations needed depends very mildly on the target accuracy. However, accuracy here means an additive error. While the optimal value could be as high as $$1/(1-\gamma)$$, it can easily happen that the best value that can be achieved, $$\|v^*\|_\infty$$, is significantly smaller than $$1/(1-\gamma)$$. It may be for example that $$\|v^*\|_\infty = 0.01$$, in which case a guarantee with $$\delta = 0.5$$ is vacuous.

By a careful inspection of $$\eqref{eq:vierror}$$ we can improve the previous result so that this problem is avoided:

Theorem (Runtime when Controlling for the Relative Error): Fix a finite discounted MDP and a target accuracy $$\delta_{\text{rel}}>0$$. Then, stopping value iteration after $k \ge H_{\gamma,\frac{\delta_{\text{rel}}}{2\gamma}}$ iterations, the policy $\pi$ produced satisfies the relative error bound

$v^\pi \ge v^* - \delta_{\text{rel}} \|v^*\|_\infty \boldsymbol{1}\,,$

while the total number of elementary arithmetic operations is

$O \left(\mathrm{S}^2 \mathrm{A} H_{\gamma, \frac{\delta_{\text{rel}}}{2\gamma}} \right) = \tilde O\left( \frac{\mathrm{S}^2 \mathrm{A} }{1 - \gamma}\, \ln\left(\frac{1}{\delta_{\text{rel}}}\right)\right)$

where $\tilde{O}(\cdot)$ hides $\log(1/(1-\gamma))$.

Notice that the runtime required to achieve a fixed relative accuracy appears to be the same as the runtime required to achieve the same level of absolute accuracy. In fact, the runtime slightly decreases. This should make sense: The worst-case for the fixed absolute accuracy is when $$\|v^*\|_\infty=1/(1-\gamma)$$, and in this case the relative accuracy is significantly less demanding: With $$\delta_{\text{rel}}=0.5$$, value iteration can stop after guaranteeing values of $$0.5/(1-\gamma)$$, which, as a value, is much smaller than $$1/(1-\gamma)-0.5$$, the target with the absolute accuracy level of $$\delta = 0.5$$.

Note that the relative error bound is not without problems either: It is possible that for some states $s$, $$v^*(s)-\delta_{\text{rel}} \|v^*\|_\infty$$ is negative, a vacuous guarantee. A reasonable stopping criteria would be top stop when the policy that we read out satisfies

$v^{\pi_k} \ge (1-\delta_{\text{rel}}) v^*\,.$

Since $$v^*$$ is not available, to arrive at a stopping condition that can be verified and which implies the above inequality, one can replace $v^*$ above with an upper bound on it, such as $$v_k +\gamma^k \|v_k\|_\infty/(1-\gamma^k) \boldsymbol{1}$$. In this imagined procedure, in each iteration, one also needs to compute the value function of policy $\pi_k$ to verify whether the stopping condition is met. If we do this much computation, we may as well replace $v_k$ with $v^{\pi_k}$ in the update equation $$\eqref{eq:vi}$$ hoping that this will further speed up convergence. This results in what is known as policy iteration, which is the subject of the next lecture.

## The Computational Complexity of Planning in MDPs

Now that we have our first results for the computation of approximately optimal policies, it is time to ask whether the algorithm we discovered is doing unnecessary work. That is, what is the minimax computational cost of calculating an optimal, or approximately optimal policy?

To precisely formulate this problem, we need to specify the inputs and the outputs of the algorithms considered. The simplest setting is when the inputs to the algorithms are arrays, describing the transition probabilities and the rewards for each state action pair with some ordering of state-action pairs (and next states in the case of transition probabilities). The output, by the Fundamental Theorem, can be a memoryless policy, either deterministic or stochastic. To describe such a policy, the algorithm could write a table. Clearly, the runtime of the algorithm will be at least the size of the table that needs to be written, so the shorter the output, the better the runtime can be. To be nice with the algorithms, we should allow them to output deterministic policies. After all, the Fundamental Theorem also guarantees that we can always find a deterministic memoryless policy which is optimal. Further, greedy policies can also be chosen to be deterministic, so the value-iteration algorithm would also satisfy this requirement. The shortest specification for a deterministic policy is an array of the size of the state space that has $$\mathrm{S}$$ entries.

Thus, the runtime of any algorithm that needs to “produce” a fully specified policy is at least $$\Omega(\mathrm{S})$$.

This is quite bad! As was noted before, $$\mathrm{S}$$, the number of states, in typical problems is expected to be gigantic. But by this easy argument we see that if we demand algorithms to produce fully specified policies then without any further help, they have to do as much work as the number of states. However, things are a bit even worse.

In Homework 0, we have seen that no algorithm can find a given value in an array without looking at all entries of the array (curiously, we saw that if we allow randomized computation, that on expectation it is enough to check half of the entries).

Based on this, it is not hard to show the following result:

Theorem (Computation Complexity of Planning in MDPs):

Let $0\le \delta < \gamma/(1-\gamma)$. Any algorithm that is guaranteed to produce $\delta$-optimal policies in any finite MDP described with tables, with a fixed discount factor $0\le \gamma <1$ and rewards in the $[0,1]$ interval needs at least $$\Omega(\mathrm{S}^2\mathrm{A})$$ elementary arithmetic operations on some MDP with the above properties and whose state space is of size $\mathrm{S}$ and action space is of size $\mathrm{A}$.

Proof sketch: We construct a family of MDPs such that no matter the algorithm, the algorithm will need to perform the said number of operations in at least one of the MDPs.

One-third of the states is reserved for “heaven”, one One-third is reserved for “hell” states. The remaining one-third set of states, call them $R$, is where the algorithms will need to make some nontrivial amount of work. The MDPs are going to be deterministic. In the tables given to the algorithms as input, we (conveniently for the algorithms) order the states so that the “hell” states come first, followed by the “heaven” states, followed by the states in $R$.

In the “heaven” class, all states self-loop under all actions and give a reward of one. The optimal value of any of these states is $1/(1-\gamma)$. In the “hell” class, states also self-loops under all actions but give a reward of zero. The optimal value of these states is $0$. For the remaining states, all actions except one lead to some hell state, while the chosen special action leads to some state in the heaven class.

The optimal value of all states in set $R$ have a value of $\gamma/(1-\gamma)$ and the value of a policy that in a state in $R$ does not choose the special optimal action gets the value of $0$ in that state. It follows that any algorithm that is guaranteed to be $\delta$ optimal needs to identify the unique optimal action at every state in $R$.

In particular, for every state $s\in R$ and action $a\in \mathcal{A}$, the algorithm needs to read $\Omega(\mathrm{S})$ entries of the transition probability vector $P_a(s)$ or it can’t find out whether $a$ leads to a state in the heaven class or the hell class: The probability vector $P_a(s)$ will have a single one at such an entry, either among the $\mathrm{S}/3$ entries representing the hell, or the $\mathrm{S}/3$ entries representing the heaven states. By the aforementioned homework problem, any algorithm that needs to find this “needle” requires to check $\Omega(\mathrm{S})$ entries. Since the number of states in $R$ is also $\Omega(\mathrm{S})$, we get that the algorithm needs to do $\Omega( \mathrm{S}\times \mathrm{A}) \mathrm{S}) = \Omega( \mathrm{S}^2 \mathrm{A})$ work. $$\qquad \blacksquare$$

We immediately see two differences between the lower bound and our previous upper bound(s): In the lower bound there is no dependence on $1/(1-\gamma)$ (the effective horizon at a constant precision). Furthermore, there is no dependence on $1/\delta$, the inverse accuracy.

As it turns out, the dependence on $1/\delta$ of value-iteration is superfluous and can be removed. The algorithm that achieves this is policy iteration, which was mentioned earlier. However, this result is saved for the next lecture. After this, the only remaining gap will be the order of the polynomials and the dependence on $1/(1-\gamma)$, which is closely related to the said polynomial order.

And of course, we save for later the most pressing issue that we need to somehow be able to avoid the situation when the runtime depends on the size of the state space (forgetting about the action space for a moment). By the lower bound just presented we already know that this will require changing the problem setting. Just how to do this will be the core question that we will keep returning to in the class.

## Notes

### Value iteration

The idea of value iteration is probably due to Richard Bellman.

### Error bound for greedification

This theorem is due to Singh & Yee, 1994.

### Computational complexity lower bound

The last theorem is due to Chen and Wang (2017), but the construction is also (unsurprisingly) similar to one that appeared in an earlier paper that studied query complexity in the setting when the access to the MDP is provided by a simulation model. In fact, we will present this lower bound later in a lecture where we study batch RL. According to this result, the query-complexity (also known as sample-complexity) of finding a $\delta$-optimal policy with constant probability in discounted MDPs accessible through a random access simulator, apart from logarithmic factors, is $SA H^3/\delta^2$, where $H=1/(1-\gamma)$.

### Representations matter

We already saw that in order to just clearly define the computational problems (which is necessary for being able to talk about lower bounds), we need to be clear about the inputs (and the outputs). The table representation of MDPs is far from being the only possibility. We just mentioned the “simulation model”. Here the algorithm “learns” about the MDP by issuing next state and reward queries to the simulator at some state-action pair $(s,a)$ of its choice to which the simulator responds with a random next state (drawn fresh) and the $r_a(s)$. Interestingly, this can provably reduce the number of queries compared to the table representation.

Another alternative, which still keeps tables, is to give the algorithm a cumulative probability representation. In this representation, the states are identified with ${1,\dots,\mathrm{S}}$ as before but instead of giving the algorithm the tables $[P_a(s,1), \dots, P_a(s,\mathrm{S})]$ for fixed $(s,a)$, the algorithm is given

$[P_a(s,1), P_a(s,1)+P_a(s,2), \dots, 1]$

(the last entry could be saved, because it is always equal to one, but in the grand scheme of things, of course, this does not matter). Now, it is not hard to see that if the original probability vector had a single one and zeroes everywhere else, the “needle in the haystack problem” used in the lower bound, with the integral representation above, a clever algorithm can find the entry with the one with at most $O(\log( \mathrm{S}))$ queries. As it turns out, with this representation, the query complexity (number of queries required) of producing a good policy can indeed be reduced from the quadratic dependence on the size of the state-space to a log-linear dependence. Hence, we see that the input representation crucially matters. Chen and Wang (2017) also make this point and they discuss yet another, “tree” representation, which leads to a similar speedup.

### MDPs with short descriptions

The simulator model assumption addresses the problem that just reading the input may be the bottleneck. This is not the only possibility. One can imagine various classes of MDPs that have a short description, which may raise the hope that one can find out a good policy in them without touching each state-action pair. There are many examples of classes of MDPs that belong to this category. These include

• factored MDPs: The transition dynamics have a short, structured (factored) representation, and the same applies to the reward
• parametric MDPs: The transition dynamics and the rewards have a short, parametric representation. Examples include linear-quadratic regulation (linear dynamics, quadratic reward, Euclidean state and action spaces, Gaussian noise in the transition dynamics), robotic systems, various operations research problems.

For factored MDPs one is out of luck: In these, planning is provably “very hard” (computationally). For linear-quadratic regulation, on the other hand, planning is “easy”; once the data is read, all one has to do is to solve some algebraic equations, for which efficient solution methods have been worked out.

### Query vs. computational complexity

The key idea of the lower bound crucially hinges upon that good algorithms need to “learn” about their inputs: The number of arithmetic and logic operations of any algorithm is at least as large as the number of “read” operations it issues. The minimum number of required read operations to produce an input of some desired property is often called the problems query complexity and by the above reasoning we see that the computational complexity is lower bounded by the query complexity. As it happens, query complexity is much easier to bound than computational complexity in the sense that it is rare to see computational complexity lower bounds strictly larger than the query complexity (the exceptions to this come when a “compact” representation of the MDP is available, such as in the case of factored MDPs). At the heart of query complexity lower bounds is often the needle in the haystack problem. This seems to be generally true when the inputs are “deterministic”. When querying results in stochastic (random) outcomes, multiple queries may be necessary to “reject”, “reduce”, or “filter out” the noise and then new considerations appear.

In any case, query complexity is a question about quickly determining the information crucial to arrive at a good decision early and is in a way about “learning”: Before a table is read, the algorithm does not know which MDP it faces. Hence, query complexity is essentially an “information” question and is also sometimes called information complexity and we can think of query complexity as the most basic information theory question. This is a bit different though than mainstream information theory, which is somehow tied up in dealing with reducing the effect of random responses (random “corruptions” of the clean information).

### Query complexity everywhere

Query complexity is widely studied in a number of communities which, sadly, are almost entirely disjoint. Information-theory, mentioned above is one of them, though as was noted, here the problems are often tied to studying the speed of gaining information in the presence of noise. Besides information theory, there is the whole field of information-based complexity, which has its own journal, multiple books and more. Also notable is the theory community that studies the complexity of evolutionary algorithms. Besides these, of course, query complexity made appearances in the optimization literature (with or without noise), operations research, and of course in the machine learning and statistics community. In particular, in the machine learning and statistics community, when the algorithm is just handed over noisy data, “the sample”, one can ask how large this sample needs to be to achieve some good outcome (e.g., good predictions on unseen data). This leads to the notion of sample complexity, which is the same as our query complexity except that the queries are of the “dull”, “passive” nature of “give me the next datapoint”. As opposed to this, “active learning” refers to the case when the algorithms themselves control some aspects of how the data is collected.

### Free lunches, needles and a bit of philosophy

Everyone after going to a few machine learning conferences or reading their first book, or blog posts would have heard about David Wolpert’s “no-free lunch theorems”. Yet, I find that to most people the exact nature (or significance) of these theorems remain elusive. Everyone heard that these theorem essentially state that “in the lack of bias, all algorithms are equal” (and therefore there is no free lunch), from which we should conclude that the only way to choose between algorithms is by introducing bias.

But what does bias means? If one reads these results carefully (and the theory community of evolutionary computation made a good job of making them accessible) one finds that the results are nothing more that describing some corollaries that to find a needle in a haystack (the special entry in a long array), one needs to search the whole haystack (query almost all entries of the array).

Believers of the power of data like to dismiss the significance of the no-free lunch result by claiming that it is ridiculous in that it assumes no structure at all. I find these arguments weak. The main problem is that they are evasive. The evasiveness comes from the reluctance to be clear about what we expect the algorithms to achieve. The claim is that once we are clear about this, that is, clear about the goals, or just the problem specification, we can always hunt for the “needle in the haystack” subproblems within the problem class. This is about figuring out the symmetries (as symmetry equals no structure) that sneakily appear in pretty much any reasonable problem we think of worth studying. The only problems that do not have “needle in the haystack” situations embedded into them are the ones that are not specified at all.

What is the upshot of all this? In a way, the real problem is to be clear about what the problem we want to solve is. This is the problem that most theoreticians in my field struggle with every day. Just because this is hard, we cannot give up on this before even starting, or this will just lead to chaos.

As we shall see in this class, how to specify the problem is also at the very heart of reinforcement learning theory research. We constantly experiment with various problem definitions, tweaking them in various ways, trying to separate hopelessly hard problems from the easy, but reasonably general ones. Theoreticians like to build a library of various problem settings that they can classify in various ways, including relating the problem settings to each other. While algorithm design is the constructive side of RL (and computer science, more generally), understanding the relationship between the various problem settings is just as equally important.

## References

• Chen, Y., & Wang, M. (2017). Lower bound on the computational complexity of discounted markov decision problems. arXiv preprint arXiv:1705.07312. [link]
• Singh, S. P., & Yee, R. C. (1994). An upper bound on the loss from approximate optimal-value functions. Machine Learning, 16(3), 227-233. [link]


text:
# Math Help - I need help with definite integrals.

1. ## I need help with definite integrals.

In, class when we were attempting questions in which we used to find area bounded by curves through definite integrals, we used to look for the interval in which the function would be negative, then we'd Divide the whole GIVEN INTERVAL into Sub-intervals in which the function was either positive OR negative. We then used to calculate the areas in these subintervals and from here the confusing part starts...

We put a negative sign with the areas under the x-axis, and then added them to the rest of the areas...fair enough and we understood what it meant....

HOWEVER, when we were solving simple definite integrals (where we were not asked about area) we did not care abt the values for which the function was positive...Instead, we just solved it by putting in the limits of the integral and we got the answer......

Im now abit confused?...Why didnt we care for the "negative values of the function" in solving SIMPLE DEFINITE INTEGRALS, just like we did in solving INTEGRALS FOR AREAS?

By the way....Nevermind.......

2. A definite integral may or may not represent area. It depends on the context of the question.

If you're asked to evaluate the area, then you will need to add all the absolute values of the definite integrals (which is equivalent to subtracting the negative definite integrals) because there is no such thing as negative area.

However, if you are simply asked to find a definite integral, it may be for some other purpose, so you would just evaluate the definite integral as is.


text:
# Avogadro’s number from thermal conductivity in an ideal gas

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 1.65.

We can get a rough estimate of Avogadro’s number from the thermal conductivity and other macroscopic quantities measurable for an ideal gas. Schroeder’s derivation of the thermal conductivity formula gives the result

$\displaystyle k_{t}=\frac{C_{V}}{2V}\ell\bar{v} \ \ \ \ \ (1)$

Suppose we set up the system so that we have a box of volume ${V}$ and cross-sectional area ${A}$. We know that the average speed ${\bar{v}}$ is approximately the rms speed, which is

$\displaystyle \bar{v}\approx v_{rms}=\sqrt{\frac{3kT}{m}}=\sqrt{\frac{3NkT}{Nm}}=\sqrt{\frac{3PV}{M}} \ \ \ \ \ (2)$

where ${M}$ is the total mass of the gas. If we measure the thermal conductivity and heat capacity, we can then get the mean free path:

$\displaystyle \ell=\frac{2Vk_{t}}{C_{V}\bar{v}}=\frac{2k_{t}}{C_{V}}\sqrt{\frac{MV}{3P}} \ \ \ \ \ (3)$

Schroeder’s expression for ${\ell}$ is based on the idea that the mean path length is equal to the length of a cylinder of radius equal to the molecule’s diameter and volume equal to the average volume per molecule ${V/N}$, so that

$\displaystyle \ell=\frac{1}{4\pi r^{2}}\frac{V}{N} \ \ \ \ \ (4)$

where ${r}$ is the radius of the molecule. To get ${N}$ from this formula, we need to know ${r}$, but this is a microscopic quantity which we’re assuming we don’t know. I can’t see any way of progressing from here unless we take a different value for ${\ell}$. Since we’re after only a rough approximation of Avogadro’s number, we can take ${\ell}$ to be the average distance between molecules, rather than the average distance between collisions. That is

$\displaystyle \ell\approx\left(\frac{V}{N}\right)^{1/3} \ \ \ \ \ (5)$

We can now combine this with 3 to get an estimate for the number of molecules ${N}$:

$\displaystyle N\approx\left(\frac{C_{V}}{2k_{t}}\right)^{3}\left(\frac{3P}{MV}\right)^{3/2}V=\frac{1}{\sqrt{V}}\left(\frac{C_{V}}{2k_{t}}\right)^{3}\left(\frac{3P}{M}\right)^{3/2} \ \ \ \ \ (6)$

Assuming we know the gas constant ${R}$, we can get the number of moles from the ideal gas law

$\displaystyle n=\frac{PV}{RT} \ \ \ \ \ (7)$

$\displaystyle N_{A}=\frac{N}{n}\approx RT\sqrt{P}\left(\frac{C_{V}}{2k_{t}}\right)^{3}\left(\frac{3}{MV}\right)^{3/2} \ \ \ \ \ (8)$

As a check on this formula, we can verify that it’s dimensionless. ${RT}$ has the dimensions of ${PV}$, so dimensionally, the quantity is equivalent to

 $\displaystyle \left(\frac{C_{V}}{k_{t}}\right)^{3}\left(\frac{P}{M}\right)^{3/2}\frac{1}{\sqrt{V}}$ $\displaystyle \rightarrow$ $\displaystyle \left[\frac{\mbox{energy}/\mbox{K}}{\mbox{energy}/\mbox{s m K}}\right]^{3}\left[\frac{\mbox{force}/\mbox{m}^{2}}{\mbox{kg}}\right]^{3/2}\frac{1}{\mbox{m}^{3/2}}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left[\mbox{s}^{3}\mbox{m}^{3}\right]\left[\mbox{m}^{-1}\mbox{s}^{-2}\right]^{3/2}\mbox{m}^{-3/2}\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 1 \ \ \ \ \ (11)$

Thus the dimensions check out.

I don’t know if this is the solution Schroeder had in mind, or whether it’s possible to get ${N_{A}}$ using the formula 4 for the mean free path. If we’re not allowed to know a priori any microscopic quantities, that means we don’t know ${k}$ (Boltzmann’s constant), ${m}$ (the mass of a single molecule) or ${r}$ (the radius of a molecule). Together with ${N}$, that makes 4 unknowns, so we need 4 independent equations to find them. Comments welcome.

## 3 thoughts on “Avogadro’s number from thermal conductivity in an ideal gas”

1. snailmecream.ru

Today, Avogadro’s number is formally defined to be the number of carbon-12 atoms in 12 grams of unbound carbon-12 in its rest-energy electronic state. The current state of the art estimates the value of N

1. Jesse

When I worked through this problem, I had the same issue, not knowing of a way to eliminate the microscopic parameters. I’m not sure if this is what Schroeder had in mind, but I ended up deciding that the ‘b’ parameter in the Van der Waal equation of state was suitably macroscopic and measurable. This equation was discussed in problem 1.17 so I figured it’s fair game.

You can do a quick “proof” that b = V/A, where V is the Van der Waal volume associated with the Van der Was radius and A is Avogadro’s number. Just show that this is the limit of the volume per molecule at high pressure. This lets you eliminate r in the mean free path and then you can solve the ideal gas conductivity equation for Avogadro’s number.

Plugging in values that I looked up for helium at 20°C and atmospheric pressure, I calculated A=5.23×10^23 using only macroscopic quantities, which seems pretty good for all the approximations involved!


text:
Creating a hostile Mob

So I'm currently developing a 2D top down view java game. I have a mob called "Mutant" that currently does nothing. It generates a random number and depending on that number moves a few pixels in that direction. I"m trying to figure out how to make the mutant follow the player's x & y coordinates when the player comes near.

Here is what my mutant does currently:

    int xa = 0;
int ya = 0;
move = generator.nextInt(50) + 1;

if(move == (50))
ya = (random.nextInt(3) - 1) * random.nextInt(2);
if (!hasCollided(xa, ya))

if(move == 40)
ya = (random.nextInt(3) - 1) * random.nextInt(2);
if (!hasCollided(xa, ya))

if(move == 30)
xa = (random.nextInt(3) - 1) * random.nextInt(2);
if (!hasCollided(xa, ya))

if(move == 20)
xa = (random.nextInt(3) - 1) * random.nextInt(2);
if (!hasCollided(xa, ya))

if(xa != 0 || ya != 0) {
move(xa, ya);
isMoving = true;
} else {
isMoving = false;
}

Does anyone know how to make the mob simulate hostility towards the player? I haven't implemented a health system yet, I've only worked on the mob's A.I.

A very simple solution might be to compare the mutant's x and y position to the x and y position of the player.

Example: Player at {10, 10}, Mutant at {0, 0}

if (mutant.x < player.x)
mutant.x += amountToMove;
else if (mutant.x > player.x)
mutant.x -= amountToMove;

if (mutant.y < player.y)
mutant.y += amountToMove;
else if (mutant.y > player.y)
mutant.y -= amountToMove;

This would always move the monster towards the player, but you may need to split the amount to move between the x and y axis so it doesn't move more quickly in a diagonal direction than in a pure vertical or horizontal motion.

totalDistance = xDistance + yDistance;
xSpeed = amountToMove * (xDistance / totalDistance); //uses proportion of xDistance to set proportion of speed.
ySpeed = amountToMove * (yDistance / totalDistance); //same for yDistance

The mutant can acquire a reference to the player at a certain range. For example, every frame (or less frequently) the range between the mutant and the player can be checked. When the mutant is within a specific range, the mutant can acquire the player as target.

How to represent that is up to you, but a short version could be:

bool isHostile = true
Entity target = player

or maybe a temporary solution:

int targetX = player.x  // needs updated every frame, as player moves
int targetY = player.y

After the mutant has the player as a target, the mutant can movement can head towards the player. A piece of math called a "normalized vector" is useful for figuring out how much to move in a given direction. Without a normalized vector, you'll run into the classic problem of having a game where diagonal movement is significantly faster than left/right or up/down movement. Here is some pseudocode for determining the direction to move, and how much to move:

  // xDiff and yDiff would be the length of the line between mutant and player
var xDiff = target.x - npc.x
var yDiff = target.y - npc.y

// unitVector is a line which is 1.0 units long and points from mutant to player
var unitVector = math.normalizeVector(xDiff, yDiff)

// move along the line from mutant to player
npc.x += unitVector.x * npc.speed * delta
npc.y += unitVector.y * npc.speed * delta

In the above code, target would be the player, and npc would be your mutant. You can ignore/remove the delta for now, and also type in some test numbers for speed -- but it may be something for later, a topic called "frame independent movement."

The approximate math for normalizeVector is:

  length = Math.sqrt((x * x) + (y * y))
var xComponent = x / length
var yComponent = y / length

This will convert an x,y movement of [+1, +1] into something like [0.707106781186547, 0.707106781186547] which is an ugly set of numbers but carries the correct amount of movement. I'll note that this particular equation will hit an error if x and y are both 0 (it'll divide by zero). If you're using a math or game library it may already have a function like this.

As DoubleDouble said, then you can layer more complex behaviours on top, for example:

• Limited sense range, e.g. only attacks the player if they are nearby or ahead of the mob.
• Play a sound, animation and/or change the movement style when the player is seen to give feedback to the player.
• Predict the player's future position based on their current velocity.
• A simple state machine, e.g. monster patrols or acts randomly until the player is seen.


text:
Errors on a Minecraft server can be extremely scary and frustrating; are you about to lose 100+ hours of work, or is it just something you can fix quick and get right back into it? One of the scariest are chunk file location errors – what to do?!

Thankfully there’s an app for that.

So today was a routine day for Windows 7 – the good ol’ Blue Screen of Death first thing in the morning. A million Windows versions and updates and still the same old rigamarole… UGH.

Normally BSODs are generally horrible, but when they rip your computer into a powerless void sometimes you can lose data that your computer was in the middle of writing. Without fail my Tekkit Lite server seems to re-write chunk files into the wrong location.

So, how do we fix this?

## Diagnosing

The first thing you want to do is look in your server log – if you see a ton of spam that looks like this, then you have a chunk location error and you’re gonna have to get your hands dirty (a little – I’ll try and make it painless)

2013-02-13 09:40:09 [INFO] [STDOUT] Chunk file at -37,3 is in the wrong location; relocating. (Expected -37, 3, got -37, 9)
...
... (and then sometime later)
...
2013-02-13 10:07:20 [INFO] [STDOUT] Wrong location! qn['Skeleton'/23872, l='world', x=-582.50, y=4.00, z=56.50]
2013-02-13 10:07:20 [INFO] [STDERR] java.lang.Exception: Stack trace
2013-02-13 10:07:20 [INFO] [STDERR] at java.lang.Thread.dumpStack(Unknown Source)
...
... (repeat one or two BAJILLION times... like, seriously?)

So if you see something like “Chunk file at xx,xx is in the wrong location; relocating. (Expected xx, xx, got yy, yy)” you’re going to have really two options:

• Delete the Chunk
• Replace the Chunk from a server backup

Your best option is always going to be to replace the chunk with a backup; so hopefully you’ve been regularly backing up your server. I tend to back ours up daily because I know how reliable Windows is NOT.

## Triage (Assessing the Damage)

So what chunk is actually wrong? What is in it? Will deleting it ruin everything you worked on? To answer this question you get do some math. This step is optional and either confirms the doom and gloom, or gives sweet relief. For us it was relief

First, take a look at the Chunk location the error is mentioning (see above); for this example it was the Chunk located -37, 3. Now those numbers are in Chunk coordinates so to get to player coordinates we are going to have to multiply by 16 (use a calc if your math skills rival ChiefDiesel’s electrical prowess).

To help figure out which direction the chunk actually lies, and what it contains, you can use the “Wrong Location! … zomg mobs gone wild” coordinates. The server was complaining about a skeleton at -582, 56 so we need to verify that our chunk math contains the point -582, 56 since we know that one is bad. The image below shows sweet little png of the math I did for our chunk location error.

So, what’s inside the chunk?

• Meh, pretty much nothing ‘cept some sweet shrubs I planted
• Lucky for you. You can either delete the dumb thing, or try and replace it with an older version from the backup
• OH NOES, IT’s MAH BASE! THE WHOLE THING!
• Srsly? Who has a base that small? Well you better have been backing up your server/world… (e.g. copying the folder every once in a while)

## Fixing the Chunk

Ok so at this point you should have an idea of what’s in the bad chunk (in case you need to go back and add your machines in again). Even more hopefully you’ve BACKED UP the server prior to the crash.

For us there was nothing important in the chunk, but I’d still rather replace it with chunks from the backup just in case.

1. You probably already did this, but just make sure your server is NOT running…
2. BACK UP YOUR WHOLE SERVER FOLDER
1. Yeah I know it’s got issues just now, but at least the whole thing isn’t corrupt. What’s easier, rebuilding a 16×16 square of your base, or your entire world?
2. I typically ZIP the server folder and rename it with a timestamp (e.g. “Tekkit Lite Server –  v0.5.7 – 2013 02 13 at 1031am.ZIP” )
1. The link above is for a stable version compatible with MC 1.4.7. If you’re running a later version, or this post is way out of date, or you want to oogle the source code – click here and get a different version
4. Open the command prompt
1. Windows -> Start -> Run -> cmd.exe
2. (alternatively) <windows-key>+<r>, then type cmd and hit <enter>
5. In the command prompt (henceforth, cmd), navigate (change directory) to the Region Fixer folder (the one you just unzipped) with the cd command
1. If you’re lazy like me then you can do a sweet shortcut by first typing cd then dragging and dropping the Region Fixer folder into the cmd window to have Windows type the whole path for you – woot!
6. Next type region-fixer.exe –help to see the available options
7. Type region-fixer.exe “C:Pathtoyourserver folderworld” to have Region Fixer check for issues. The results should match with the errors you were seeing in your server log; i.e. 1 chunk in wrong location.
8. Now it’s go-time.
1. If you have a backup (e.g. a folder copy of your server or world from before the crash) you can call Region Fixer again, passing it the –rw parameter (to make it replace wrong chunk) and the  –backups=”C:PathtotheBack Upserverworld” parameter to make it use a previous version for the relocation.
1. CAUTION: This is the safer of the two options, but if your server backup is really old you might lose stuff you built inside the chunk.
2. If you don’t have a backup anywhere, you’re going to have to use the –dw parameter, ouch :/
1. CAUTION: This will DELETE wrong located chunks (they will regenerate as fresh chunks when you log in again), and you will LOSE ERRYYYTHING you had built inside the chunk.
9. DONE! At this point you should be all set. Region Fixer should have either deleted or replaced the chunk(s) and you can go ahead and re-launch your server. You should see a nice clean, error-free, server log (ignoring the normal mod loading spam) – all the chunks happily in the right places!

## Note on Region-Fixer parameters

All the parameters typed out in this tutorial use a double dash denotation (two dashes in a row) but for some reason the font, or WordPress itself, make it look like a single line. It’s not.

To dramatize this… the –rw parameter from before should contain two dashes, like this:

--rw

I would have used the preformatted font style everywhere but WordPress doesn’t play nice.

## Credit

The Region Fixer tool used in this tutorial is clearly awesome and thanks and donations should go to Fenixin

The tutorial here was put together entirely by us (The Mindcrafters) so if you found it helpful, or thought some points could use clarification please let us know!

Happy Chunking.

### 7 Responses to “Minecraft Chunk File in the Wrong Location”

I’m having trouble with step #5. When i type in “cd” nothing happens. I tried typing “change directory but it just says CHANGE PLEASE help me! I have a Windows 8 computer now but it’s not one of those new ones. I have a normal laptop with windows 8 installed. I really put a lot of work into this server! Help!

OK let’s get you back up and running!

Make sure to open the command prompt first (Start->Run-> “cmd”) then when the window comes up you need to to change the directory. To get it to change you have to issue the “cd” command like so
 cd C:\path\to\the\folder

If your username or any folders along the way have spaces in them you need to do the following
cd "C:\path\to\the\folder with\spaces\"

Give that a try and let me know how it goes!

So I’m confused on step #8. I don’t have a back up so I have to go the –dw route. What do I do exactly from that point on? What do I have to type in?

Step 8 without a back up:

(sorry my english, I hope you understand me)

the command is: region-fixer.exe –rw (and here the direcction of your server folder)