Homework 2

Overview

Three numerical computing problems: spectral peak estimation with splines, Gibbs sampling of the 2-D Ising model, and a command-line secant root finder with input validation.

Getting Started

Download the starter code: hw2-starter.zip

unzip hw2-starter.zip
cd hw2-starter
python generate_datasets.py --seed 541

This creates data/ with a placeholder raman.rod in the same layout as the real file but with random peaks. Download the real file from the Data link into data/ to replace it.

Each problem directory contains the script named in the problem, with a docstring stating how it is run and what it prints or writes. q1/ and q3/ also contain a test_interfaces.py that runs the script and checks its output. Run it from inside the problem directory:

python -m pytest test_interfaces.py

Problem 1: Peak Finding

Raman spectroscopy is a technique that uses inelastic scattering of light to identify unknown chemical substances. Spectral “peaks” indicate vibrational and rotational modes and are of special importance because they act like a chemical fingerprint. Raman spectroscopy measures photon intensity vs. Raman shift. The Raman shift relates the frequencies of the exciting laser and the scattered photons and is often reported as a wavenumber – the frequency difference in wavelengths per cm (i.e., cm\(^{-1}\)).

Generate a molecular fingerprint using the spectroscopic data in raman.rod. The file contains intensity vs. wavenumber data for an unknown chemical sample. A Raman Open Database (ROD) file includes content in addition to the raw intensity data:

# content
more content
_raman_spectrum.intensity
wavenumber1 intensity1
wavenumber2 intensity2
...
wavenumbern intensityn

Use string matching to ignore all lines before _raman_spectrum.intensity. Load valid (wavenumber, intensity) pairs until the first invalid intensity line (or upon reaching the end of file).

Use the method below to estimate the wavenumbers of all spectral peaks.

First detect peaks in the raw spectral data. Use the peak locations to focus on regions of interest within the spectrum. For instance: if you detect peaks at \(x_1\) cm\(^{-1}\) and \(x_2\) cm\(^{-1}\) use regions of interest: \([x_1 - n_1, x_1 + n_1]\) and \([x_2 - n_2, x_2 + n_2]\). Experiment to find “good” widths \(n_1\), \(n_2\), etc. Then use a spline to interpolate intensity within each region of interest. Calculate zero-crossings of the derivative to estimate wavenumbers with maximum intensity.

a. Print the wavenumber estimates for the eight largest spectral peak to STDOUT sorted by magnitude (largest first).

b. Create a figure that shows the Raman data (intensity vs. wavenumber) and mark each of the maximum intensity values.

c. Produce a “zoomed-in” figure for the “regions of interest” corresponding to the four largest peaks. Plot the raw spectral data and overlay your interpolating function. Use a marker to show the wavenumber with maximal intensity.

Deliverables

See Submission. peaks.py is your implementation. q1.pdf contains the eight wavenumber estimates, the intensity vs. wavenumber figure, and the four zoomed-in figures.

Problem 2: Gibbs Sampler

A Gibbs sampler generates samples from a joint probability distribution by iteratively sampling each variable conditioned on the others. In this problem, you will implement Gibbs sampling to analyze the 2-dimensional Ising model, which represents a lattice of binary spins interacting with nearest neighbors. This model is used in statistical physics to study phase transitions and collective behavior arising from local spin interactions.

Part A: Sampler

Implement a Gibbs sampler for a \(N = 20\) grid (i.e. \(20 \times 20\)) where each site has state \(s_i \in \{-1, +1\}\). The system energy is:

where the sum is over nearest-neighbor pairs. The probability of a configuration is:

Initialize the grid randomly and store the initial state.

Part B: Energy

Run the sampler for \(\beta = \{0.2, 0.4, 0.6\}\). For each site update:

where \(\Delta E_i = \text{energy change if site } i \text{ flips}\). Record the energy after every 10 iterations for 1000 total iterations. Create a figure showing \(E(\mathbf{s})\) versus iteration number for each \(\beta\).

Part C: Magnetization

Calculate the average magnetization:

Plot \(M\) versus iteration number. Create visualizations of the grid state at iterations \(\{0, 100, 500, 1000\}\) using:

plt.imshow(grid, cmap='binary')
plt.colorbar()
plt.savefig(f'state_{beta}_{iter}.png')

Deliverables

See Submission. gibbs.py is your implementation. q2.pdf contains the energy and magnetization plots for each \(\beta\) and the grid state visualizations.

Problem 3: Root Finding

The secant method is an iterative root-finding algorithm. It uses a sequence of secant line roots to approximate \(c\) such that \(f(c) = 0\) for a continuous function \(f\). Unlike Newton’s method it does not require knowledge or evaluation of the derivative \(f'\). The secant method is defined by the recurrence:

Write a python script that uses the secant method to approximate the root of a continuous function \(f\) in the interval \([a, b]\). You may assume that \(f\) has at most one root in \([a,b]\). Use \(|x_{k+1} - x_{k}| < 10^{-10}\) as the convergence criterion. Let \(N\) be the number of iterations to reach convergence. Output \(N\) followed by the three root approximations \(x_{N-2}\), \(x_{N-1}\), \(x_N\). Output each number to its own line and use precision sufficient to show convergence.

Import the function \(f\) from a file named func.py in the same directory as your script – i.e. from func import f – example func.py (included in starter code). You may assume that \(f\) is continuous on \([a,b]\) and that func.f(x) returns a scalar float for all \(x \in [a,b]\).

Your script should accept a and b as two numeric command line arguments, i.e. python secant.py "1.1" "1.4". Your script must validate that a and b are numeric, verify that a < b, and check that f(a)f(b) < 0 – see Bolzano’s Theorem. Write “Range error” to STDERR (standard error) if any of these three conditions fail and immediately terminate.

Deliverables

See Submission. secant.py is your script. Include the func.py you tested against.

README.md
.gitignore
requirements.txt
generate_datasets.py
q1/
├── peaks.py
├── test_interfaces.py
└── q1.pdf
q2/
├── gibbs.py
└── q2.pdf
q3/
├── secant.py
├── func.py
└── test_interfaces.py