Scientific Programming
Best practices for research software and reproducible computation.
Project Structure
project/
├── README.md # Project overview, how to reproduce
├── LICENSE # MIT, Apache 2.0, or GPL
├── requirements.txt # or environment.yml (conda)
├── setup.py / pyproject.toml
├── data/
│ ├── raw/ # Never modify raw data
│ ├── processed/ # Cleaned/transformed data
│ └── external/ # Third-party data
├── src/ or scripts/
│ ├── data_processing.py
│ ├── analysis.py
│ ├── models.py
│ └── visualization.py
├── notebooks/ # Exploratory analysis
│ ├── 01_eda.ipynb
│ ├── 02_modeling.ipynb
│ └── 03_figures.ipynb
├── results/
│ ├── figures/
│ └── tables/
├── tests/
└── docs/
Reproducibility Checklist
Environment: Pin all dependencies with versions
pip freeze > requirements.txt
# or conda
conda env export > environment.yml
Random seeds: Set and document all random seeds
import numpy as np
import random
SEED = 42
np.random.seed(SEED)
random.seed(SEED)
# torch.manual_seed(SEED)
# tf.random.set_seed(SEED)
Data versioning: Use DVC or git-lfs for large data
dvc init
dvc add data/raw/dataset.csv
git add data/raw/dataset.csv.dvc
Configuration: Separate config from code
# config.yaml
# experiment:
# learning_rate: 0.001
# batch_size: 32
# epochs: 100
import yaml
with open('config.yaml') as f:
config = yaml.safe_load(f)
Logging: Record all experiments
import logging
logging.basicConfig(level=logging.INFO,
format='%(asctime)s %(levelname)s: %(message)s',
filename='experiment.log')
Parallel Computing
# Multiprocessing (CPU-bound)
from multiprocessing import Pool
import numpy as np
def process_chunk(data):
return heavy_computation(data)
with Pool(processes=8) as pool:
results = pool.map(process_chunk, data_chunks)
# Concurrent futures (simpler API)
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
with ProcessPoolExecutor(max_workers=8) as executor:
results = list(executor.map(process_func, items))
# For I/O-bound tasks (API calls, file reading)
with ThreadPoolExecutor(max_workers=20) as executor:
results = list(executor.map(fetch_data, urls))
Performance Optimization
# Profiling
import cProfile
cProfile.run('my_function()', sort='cumulative')
# Line profiling
# pip install line_profiler
# @profile decorator, then: kernprof -l -v script.py
# NumPy vectorization (avoid loops)
# Bad:
result = [x**2 + 2*x + 1 for x in data]
# Good:
result = data**2 + 2*data + 1
# Memory profiling
# pip install memory_profiler
# @profile decorator, then: python -m memory_profiler script.py
Data Management
FAIR Principles
- Findable: Persistent identifiers (DOI), rich metadata
- Accessible: Open protocols, authentication when needed
- Interoperable: Standard formats (CSV, JSON, HDF5, NetCDF)
- Reusable: Clear license, provenance, community standards
File Formats for Science
| Format |
Best For |
Size |
Speed |
| CSV |
Small tabular, universal |
Large |
Slow |
| Parquet |
Large tabular, columnar |
Small |
Fast |
| HDF5 |
Multidimensional arrays |
Small |
Fast |
| NetCDF |
Climate/geospatial |
Small |
Fast |
| FITS |
Astronomy |
Medium |
Fast |
| Feather |
DataFrame interchange |
Small |
Very fast |
# Parquet (recommended for large datasets)
df.to_parquet('data.parquet', compression='snappy')
df = pd.read_parquet('data.parquet')
# HDF5 (for arrays)
import h5py
with h5py.File('data.h5', 'w') as f:
f.create_dataset('experiment1', data=array)
Testing Scientific Code
import numpy as np
import pytest
def test_conservation_law():
"""Physical quantities should be conserved"""
initial_energy = compute_energy(initial_state)
final_energy = compute_energy(simulate(initial_state))
np.testing.assert_allclose(initial_energy, final_energy, rtol=1e-6)
def test_known_solution():
"""Compare against analytical solution"""
numerical = solve_numerically(params)
analytical = analytical_solution(params)
np.testing.assert_allclose(numerical, analytical, atol=1e-4)
def test_symmetry():
"""Result should be symmetric under transformation"""
result1 = compute(data)
result2 = compute(transform(data))
np.testing.assert_array_equal(result1, result2)
Tips
- Raw data is sacred — never modify it, only create processed copies
- Use version control (git) from day one
- Write README before writing code
- Automate the full pipeline (Makefile or Snakemake)
- Document assumptions and decisions in code comments
- Use type hints for clarity in scientific code
- Publish code alongside papers (GitHub + Zenodo for DOI)
1---2name: code-science3description: Scientific programming best practices including reproducible research, computational notebooks, version control for research, data management, HPC/parallel computing, and research software engineering. Use when user needs help with research code organization, reproducibility, scientific Python/R workflows, or computational infrastructure. Triggers on "reproducible research", "research code", "scientific computing", "HPC", "parallel computing", "Jupyter", "notebook", "data management plan", "research software", "code review for science".4---56# Scientific Programming78Best practices for research software and reproducible computation.910## Project Structure1112```13project/14├── README.md # Project overview, how to reproduce15├── LICENSE # MIT, Apache 2.0, or GPL16├── requirements.txt # or environment.yml (conda)17├── setup.py / pyproject.toml18├── data/19│ ├── raw/ # Never modify raw data20│ ├── processed/ # Cleaned/transformed data21│ └── external/ # Third-party data22├── src/ or scripts/23│ ├── data_processing.py24│ ├── analysis.py25│ ├── models.py26│ └── visualization.py27├── notebooks/ # Exploratory analysis28│ ├── 01_eda.ipynb29│ ├── 02_modeling.ipynb30│ └── 03_figures.ipynb31├── results/32│ ├── figures/33│ └── tables/34├── tests/35└── docs/36```3738## Reproducibility Checklist39401. **Environment**: Pin all dependencies with versions41 ```bash42 pip freeze > requirements.txt43 # or conda44 conda env export > environment.yml45 ```46472. **Random seeds**: Set and document all random seeds48 ```python49 import numpy as np50 import random51 SEED = 4252 np.random.seed(SEED)53 random.seed(SEED)54 # torch.manual_seed(SEED)55 # tf.random.set_seed(SEED)56 ```57583. **Data versioning**: Use DVC or git-lfs for large data59 ```bash60 dvc init61 dvc add data/raw/dataset.csv62 git add data/raw/dataset.csv.dvc63 ```64654. **Configuration**: Separate config from code66 ```python67 # config.yaml68 # experiment:69 # learning_rate: 0.00170 # batch_size: 3271 # epochs: 10072 import yaml73 with open('config.yaml') as f:74 config = yaml.safe_load(f)75 ```76775. **Logging**: Record all experiments78 ```python79 import logging80 logging.basicConfig(level=logging.INFO, 81 format='%(asctime)s %(levelname)s: %(message)s',82 filename='experiment.log')83 ```8485## Parallel Computing8687```python88# Multiprocessing (CPU-bound)89from multiprocessing import Pool90import numpy as np9192def process_chunk(data):93 return heavy_computation(data)9495with Pool(processes=8) as pool:96 results = pool.map(process_chunk, data_chunks)9798# Concurrent futures (simpler API)99from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor100101with ProcessPoolExecutor(max_workers=8) as executor:102 results = list(executor.map(process_func, items))103104# For I/O-bound tasks (API calls, file reading)105with ThreadPoolExecutor(max_workers=20) as executor:106 results = list(executor.map(fetch_data, urls))107```108109## Performance Optimization110111```python112# Profiling113import cProfile114cProfile.run('my_function()', sort='cumulative')115116# Line profiling117# pip install line_profiler118# @profile decorator, then: kernprof -l -v script.py119120# NumPy vectorization (avoid loops)121# Bad:122result = [x**2 + 2*x + 1 for x in data]123# Good:124result = data**2 + 2*data + 1125126# Memory profiling127# pip install memory_profiler128# @profile decorator, then: python -m memory_profiler script.py129```130131## Data Management132133### FAIR Principles134- **Findable**: Persistent identifiers (DOI), rich metadata135- **Accessible**: Open protocols, authentication when needed136- **Interoperable**: Standard formats (CSV, JSON, HDF5, NetCDF)137- **Reusable**: Clear license, provenance, community standards138139### File Formats for Science140| Format | Best For | Size | Speed |141|--------|----------|------|-------|142| CSV | Small tabular, universal | Large | Slow |143| Parquet | Large tabular, columnar | Small | Fast |144| HDF5 | Multidimensional arrays | Small | Fast |145| NetCDF | Climate/geospatial | Small | Fast |146| FITS | Astronomy | Medium | Fast |147| Feather | DataFrame interchange | Small | Very fast |148149```python150# Parquet (recommended for large datasets)151df.to_parquet('data.parquet', compression='snappy')152df = pd.read_parquet('data.parquet')153154# HDF5 (for arrays)155import h5py156with h5py.File('data.h5', 'w') as f:157 f.create_dataset('experiment1', data=array)158```159160## Testing Scientific Code161162```python163import numpy as np164import pytest165166def test_conservation_law():167 """Physical quantities should be conserved"""168 initial_energy = compute_energy(initial_state)169 final_energy = compute_energy(simulate(initial_state))170 np.testing.assert_allclose(initial_energy, final_energy, rtol=1e-6)171172def test_known_solution():173 """Compare against analytical solution"""174 numerical = solve_numerically(params)175 analytical = analytical_solution(params)176 np.testing.assert_allclose(numerical, analytical, atol=1e-4)177178def test_symmetry():179 """Result should be symmetric under transformation"""180 result1 = compute(data)181 result2 = compute(transform(data))182 np.testing.assert_array_equal(result1, result2)183```184185## Tips186- Raw data is sacred — never modify it, only create processed copies187- Use version control (git) from day one188- Write README before writing code189- Automate the full pipeline (Makefile or Snakemake)190- Document assumptions and decisions in code comments191- Use type hints for clarity in scientific code192- Publish code alongside papers (GitHub + Zenodo for DOI)