Erasure Coding in Distributed Storage: Mathematics, Implementation, and Trade-offs

# Reed-Solomon Encoding in Python using pyfinite

import numpy as np
from pyfinite import ffield

class ReedSolomonEncoder:
    def __init__(self, k: int, m: int):
        """
        k: number of data blocks
        m: number of parity blocks
        """
        self.k = k
        self.m = m
        self.n = k + m  # total blocks
        
        # Galois Field GF(2^8) for byte operations
        self.gf = ffield.FField(8)  # GF(256)
        
        # Generate encoding matrix
        self.encoding_matrix = self._generate_encoding_matrix()
    
    def _generate_encoding_matrix(self):
        """Generate Vandermonde matrix for encoding"""
        matrix = np.zeros((self.n, self.k), dtype=np.int32)
        
        # Identity matrix for data blocks (no transformation)
        for i in range(self.k):
            matrix[i, i] = 1
        
        # Vandermonde matrix for parity blocks
        for i in range(self.m):
            for j in range(self.k):
                # matrix[k+i, j] = (i+1)^j in GF(2^8)
                matrix[self.k + i, j] = self.gf.Multiply(i + 1, j)
        
        return matrix
    
    def encode(self, data_blocks: list[bytes]) -> list[bytes]:
        """
        Encode k data blocks into n = k+m total blocks
        
        Args:
            data_blocks: k blocks of equal size
            
        Returns:
            n blocks (k data + m parity)
        """
        assert len(data_blocks) == self.k, f"Expected {self.k} blocks"
        
        block_size = len(data_blocks[0])
        assert all(len(b) == block_size for b in data_blocks)
        
        # Convert blocks to numpy arrays
        data = np.array([list(block) for block in data_blocks], dtype=np.uint8)
        
        # Allocate output for all blocks
        encoded = np.zeros((self.n, block_size), dtype=np.uint8)
        
        # Copy data blocks (identity portion)
        encoded[:self.k] = data
        
        # Compute parity blocks
        for i in range(self.m):
            parity = np.zeros(block_size, dtype=np.uint8)
            for j in range(self.k):
                # Multiply data block by encoding matrix coefficient
                coeff = self.encoding_matrix[self.k + i, j]
                for byte_idx in range(block_size):
                    parity[byte_idx] ^= self.gf.Multiply(
                        coeff,
                        data[j, byte_idx]
                    )
            encoded[self.k + i] = parity
        
        # Convert back to bytes
        return [bytes(encoded[i]) for i in range(self.n)]
    
    def can_decode(self, available_blocks: list[int]) -> bool:
        """Check if we can decode from available blocks"""
        return len(available_blocks) >= self.k

# Example usage:
encoder = ReedSolomonEncoder(k=6, m=3)

# Original data: 6 blocks of 4KB each
data_blocks = [os.urandom(4096) for _ in range(6)]

# Encode: produces 9 blocks total (6 data + 3 parity)
all_blocks = encoder.encode(data_blocks)

# Simulate 3 block failures
surviving_blocks = [0, 1, 2, 4, 5, 8]  # Blocks 3, 6, 7 lost
print(f"Can recover: {encoder.can_decode(surviving_blocks)}")  # True!

# Reed-Solomon Decoding

class ReedSolomonDecoder:
    def __init__(self, encoder: ReedSolomonEncoder):
        self.encoder = encoder
        self.k = encoder.k
        self.m = encoder.m
        self.n = encoder.n
        self.gf = encoder.gf
    
    def decode(
        self,
        blocks: dict[int, bytes],  # block_index -> block_data
    ) -> list[bytes]:
        """
        Recover original k data blocks from any k surviving blocks
        
        Args:
            blocks: dict mapping block index to block data
                   (must have at least k entries)
        
        Returns:
            k recovered data blocks
        """
        assert len(blocks) >= self.k, "Not enough blocks to recover"
        
        available_indices = sorted(blocks.keys())[:self.k]
        block_size = len(next(iter(blocks.values())))
        
        # Extract submatrix for available blocks
        encoding_submatrix = self.encoder.encoding_matrix[available_indices]
        
        # Invert submatrix in GF(2^8)
        decoding_matrix = self._invert_matrix(encoding_submatrix)
        
        # Convert available blocks to array
        available_data = np.array(
            [list(blocks[i]) for i in available_indices],
            dtype=np.uint8
        )
        
        # Matrix multiplication to recover data blocks
        recovered = np.zeros((self.k, block_size), dtype=np.uint8)
        for i in range(self.k):
            for j in range(self.k):
                coeff = decoding_matrix[i, j]
                for byte_idx in range(block_size):
                    recovered[i, byte_idx] ^= self.gf.Multiply(
                        coeff,
                        available_data[j, byte_idx]
                    )
        
        return [bytes(recovered[i]) for i in range(self.k)]
    
    def _invert_matrix(self, matrix):
        """Invert matrix in Galois Field GF(2^8)"""
        n = len(matrix)
        # Create augmented matrix [A | I]
        augmented = np.hstack([matrix, np.eye(n, dtype=np.int32)])
        
        # Gaussian elimination in GF(2^8)
        for i in range(n):
            # Find pivot
            pivot = augmented[i, i]
            if pivot == 0:
                # Swap rows
                for j in range(i + 1, n):
                    if augmented[j, i] != 0:
                        augmented[[i, j]] = augmented[[j, i]]
                        pivot = augmented[i, i]
                        break
            
            # Scale row
            pivot_inv = self.gf.Inverse(pivot)
            for j in range(2 * n):
                augmented[i, j] = self.gf.Multiply(augmented[i, j], pivot_inv)
            
            # Eliminate column
            for j in range(n):
                if i != j and augmented[j, i] != 0:
                    factor = augmented[j, i]
                    for k in range(2 * n):
                        augmented[j, k] ^= self.gf.Multiply(
                            factor,
                            augmented[i, k]
                        )
        
        # Extract inverse from augmented matrix
        return augmented[:, n:]

# Example: Recovery from failures
encoder = ReedSolomonEncoder(k=6, m=3)
decoder = ReedSolomonDecoder(encoder)

# Encode data
data = [f"Block {i}".encode().ljust(4096, b'\0') for i in range(6)]
all_blocks = encoder.encode(data)

# Simulate failures: lose blocks 2, 5, 7
surviving = {
    i: all_blocks[i] 
    for i in [0, 1, 3, 4, 6, 8]  # 6 of 9 blocks survived
}

# Recover original data
recovered = decoder.decode(surviving)

# Verify
assert recovered == data[:6], "Recovery failed!"
print("Successfully recovered from 3 block failures!")

# HDFS Erasure Coding Configuration

# hdfs-site.xml
<configuration>
  <!-- Enable erasure coding -->
  <property>
    <name>dfs.namenode.ec.policies.enabled</name>
    <value>true</value>
  </property>
  
  <!-- Default policy: RS-6-3-1024k -->
  <property>
    <name>dfs.namenode.ec.policies.default</name>
    <value>RS-6-3-1024k</value>
  </property>
  
  <!-- Stripe size: 1MB per block -->
  <property>
    <name>dfs.namenode.ec.stripe.size</name>
    <value>1048576</value>
  </property>
  
  <!-- Reconstruction threads -->
  <property>
    <name>dfs.datanode.ec.reconstruction.threads</name>
    <value>8</value>
  </property>
</configuration>

# Set EC policy on directory
hdfs ec -setPolicy -path /data/warehouse -policy RS-6-3-1024k

# View EC statistics
hdfs ec -listPolicies

# Output:
# Policy: RS-6-3-1024k
#   Data units: 6
#   Parity units: 3
#   Cell size: 1048576 bytes
#   Codec: rs
#   Status: ENABLED

# HDFS EC Architecture:

┌───────────────────────────────────────────────────────┐
│                     NameNode                          │
│  - Tracks block -> DataNode mapping                  │
│  - Triggers reconstruction on failures                │
└───────────────────────────────────────────────────────┘
                          │
          ┌───────────────┼───────────────┐
          │               │               │
┌─────────▼────┐ ┌────────▼────┐ ┌───────▼─────┐
│  DataNode 1  │ │ DataNode 2  │ │ DataNode 3  │
│  [D₀][D₃][P₀]│ │ [D₁][D₄][P₁]│ │ [D₂][D₅][P₂]│
└──────────────┘ └─────────────┘ └─────────────┘

# Client reads:
# - Reads data blocks D₀-D₅ directly (fast path)
# - On failure, reads any 6 of 9 blocks + reconstructs

# Storage savings:
# - 10PB raw data
# - Replication (3x): 30PB storage
# - EC (6+3): 15PB storage
# - Savings: 15PB (50% reduction!)

# Production-Ready Erasure Coding Library

from typing import List, Dict, Optional
import numpy as np
from concurrent.futures import ThreadPoolExecutor
import hashlib

class ProductionErasureCoder:
    def __init__(
        self,
        k: int,
        m: int,
        block_size: int = 1024 * 1024,  # 1MB blocks
        checksum: bool = True
    ):
        self.k = k
        self.m = m
        self.block_size = block_size
        self.checksum = checksum
        
        # Use Intel ISA-L for hardware acceleration
        try:
            import pyeclib
            self.backend = pyeclib.ECDriver(
                k=k,
                m=m,
                ec_type='liberasurecode_rs_vand',  # Vandermonde RS
                chksum_type='crc32' if checksum else None
            )
        except ImportError:
            # Fallback to pure Python
            self.backend = ReedSolomonEncoder(k, m)
    
    def encode_file(
        self,
        file_path: str,
        output_dir: str
    ) -> Dict[int, str]:
        """
        Encode file into EC blocks
        
        Returns:
            Dict mapping block index to output file path
        """
        # Read file in chunks
        with open(file_path, 'rb') as f:
            data = f.read()
        
        # Pad to multiple of block_size * k
        stripe_size = self.block_size * self.k
        padded_size = (len(data) + stripe_size - 1) // stripe_size * stripe_size
        data = data.ljust(padded_size, b'\0')
        
        # Split into stripes
        num_stripes = len(data) // stripe_size
        output_files = {}
        
        with ThreadPoolExecutor(max_workers=4) as executor:
            futures = []
            
            for stripe_idx in range(num_stripes):
                stripe_start = stripe_idx * stripe_size
                stripe_end = stripe_start + stripe_size
                stripe_data = data[stripe_start:stripe_end]
                
                # Split stripe into k data blocks
                data_blocks = [
                    stripe_data[i * self.block_size:(i + 1) * self.block_size]
                    for i in range(self.k)
                ]
                
                # Encode stripe
                future = executor.submit(self._encode_stripe, stripe_idx, data_blocks, output_dir)
                futures.append(future)
            
            # Wait for all stripes
            for future in futures:
                stripe_files = future.result()
                output_files.update(stripe_files)
        
        return output_files
    
    def _encode_stripe(
        self,
        stripe_idx: int,
        data_blocks: List[bytes],
        output_dir: str
    ) -> Dict[int, str]:
        """Encode single stripe"""
        # Encode with EC
        if hasattr(self.backend, 'encode'):
            # pyeclib backend
            all_blocks = self.backend.encode(b''.join(data_blocks))
        else:
            # Pure Python backend
            all_blocks = self.backend.encode(data_blocks)
        
        # Write blocks to disk
        output_files = {}
        for block_idx, block_data in enumerate(all_blocks):
            filename = f"stripe_{stripe_idx}_block_{block_idx}.ec"
            filepath = os.path.join(output_dir, filename)
            
            with open(filepath, 'wb') as f:
                if self.checksum:
                    # Add checksum
                    checksum = hashlib.sha256(block_data).digest()
                    f.write(checksum)
                f.write(block_data)
            
            output_files[block_idx] = filepath
        
        return output_files
    
    def reconstruct_file(
        self,
        block_files: Dict[int, str],
        output_path: str,
        original_size: int
    ) -> None:
        """Reconstruct original file from EC blocks"""
        # Group by stripes
        stripes = self._group_by_stripe(block_files)
        
        with open(output_path, 'wb') as out_f:
            for stripe_idx in sorted(stripes.keys()):
                stripe_blocks = stripes[stripe_idx]
                
                # Verify we have enough blocks
                if len(stripe_blocks) < self.k:
                    raise ValueError(f"Insufficient blocks for stripe {stripe_idx}")
                
                # Read block data
                blocks_data = {}
                for block_idx, filepath in stripe_blocks.items():
                    with open(filepath, 'rb') as f:
                        if self.checksum:
                            expected_checksum = f.read(32)
                            block_data = f.read()
                            actual_checksum = hashlib.sha256(block_data).digest()
                            if actual_checksum != expected_checksum:
                                raise ValueError(f"Checksum failed for block {block_idx}")
                        else:
                            block_data = f.read()
                        blocks_data[block_idx] = block_data
                
                # Decode stripe
                if hasattr(self.backend, 'decode'):
                    # pyeclib backend
                    recovered = self.backend.decode(
                        list(blocks_data.values()),
                        list(blocks_data.keys())
                    )
                else:
                    # Pure Python backend
                    recovered = self.backend.decode(blocks_data)
                
                # Write recovered data
                out_f.write(b''.join(recovered))
        
        # Truncate to original size
        with open(output_path, 'r+b') as f:
            f.truncate(original_size)

# Usage example:
ec = ProductionErasureCoder(k=6, m=3, block_size=1024*1024)

# Encode large file
output_files = ec.encode_file('/data/input.dat', '/data/ec_blocks/')

# Simulate failures (remove 3 blocks)
surviving_files = {k: v for k, v in list(output_files.items())[:-3]}

# Reconstruct from remaining blocks
ec.reconstruct_file(surviving_files, '/data/recovered.dat', original_size=...)