i/PXdZddlmZddlmZddlmZddZeZ dZ da dZ dddZ dZd Zd Zdad Zej& d dd Z d dd ZdZdadZej& d ddZ d ddZdZdadZddZy)a5 Custom fused Metal kernels for Demucs-MLX inference. Each kernel includes: - Metal source code - Python wrapper with shape handling - float32 accumulation for numerical stability When Metal is not available (e.g., MLX on Linux/CPU), all functions fall back to equivalent pure-MLX operations automatically. ) annotationsNc ttdd}|yt|dd}t|rt|Sy#t$rYywxYw)z-Check whether Metal GPU backend is available.metalNF is_availableT)getattrmxcallablebool Exception)ris_avails m/Volumes/fast/ai/experiments/voice-extract-mac/.venv/lib/python3.12/site-packages/demucs_mlx/metal_kernels.py _has_metalrsW GT* =5.$7 H  # # s?'? A  A a uint gid = thread_position_in_grid.x; uint total = params[0]; // N * C (total output elements) uint half_dim = params[1]; // C (half of last dimension) if (gid >= total) return; uint row = gid / half_dim; uint col = gid % half_dim; uint full_dim = half_dim * 2; float a = (float)x[row * full_dim + col]; float b = (float)x[row * full_dim + half_dim + col]; float sig_b = 1.0f / (1.0f + metal::exp(-b)); out[gid] = (T)(a * sig_b); cnt*tjjdddgdgtatS)N fused_gluxparamsoutname input_names output_namessource) _glu_kernelrfast metal_kernel _GLU_SOURCEr _get_glu_kernelrBs;gg**h + rc ts3tj|d|\}}|tj|zS|j}||z}t |j }||dzdk7rtd|d||||dz k7r4t t|}|d||c||<|d<|j|}|j d}|dz}tj|jd|} | j d} | |z} tj| |gtj} t| jd| gd |jfg| ddft!d | ddf| fg|jg d} | jg|j d d|}||dz k7r4t t|}|d||c||<|d<|j|}|S) zFused GLU: split x in half along axis, compute a * sigmoid(b). Equivalent to: a, b = mx.split(x, 2, axis=axis) return a * mx.sigmoid(b) axisrzAxis z size must be even, got dtypeTinputstemplategrid threadgroup output_shapes output_dtypesN) HAS_METALrsplitsigmoidndimlistshape ValueErrorrange transpose contiguousreshapearrayint32rr'min)rr#abr4in_shapepermlast_dimhalfx_2dNtotalr result_flatresults r rrNs xx14(12::a=  66D $;DAGG}H~Q5&>x~>NOPP taxE$K #BxdT DH AKK wwr{H q=D ==2x0 1D 1 A HE XXudm288 4F#/# R &).!Q]e_a+xjwwi  K![ 5!''#2, 5 5F taxE$K #BxdT DH!!!4( Mriai // Abramowitz & Stegun approximation of erf, max error ~1.5e-7 inline float erf_approx(float x) { // erf(-x) = -erf(x) float sign = (x >= 0.0f) ? 1.0f : -1.0f; float a = metal::abs(x); // A&S formula 7.1.26 float t = 1.0f / (1.0f + 0.3275911f * a); float t2 = t * t; float t3 = t2 * t; float t4 = t3 * t; float t5 = t4 * t; float poly = 0.254829592f * t - 0.284496736f * t2 + 1.421413741f * t3 - 1.453152027f * t4 + 1.061405429f * t5; float result = 1.0f - poly * metal::exp(-a * a); return sign * result; } a uint bg = threadgroup_position_in_grid.x; uint tid = thread_index_in_threadgroup; uint tg_size = threads_per_threadgroup.x; uint sid = thread_index_in_simdgroup; uint wid = simdgroup_index_in_threadgroup; uint num_simdgroups = tg_size / 32; uint num_groups = params[0]; uint channels_per_group = params[1]; uint spatial_size = params[2]; uint total_channels = params[3]; uint batch_idx = bg / num_groups; uint group_idx = bg % num_groups; uint elems_per_group = channels_per_group * spatial_size; uint base = batch_idx * total_channels * spatial_size + group_idx * channels_per_group * spatial_size; // Pass 1: Compute mean float local_sum = 0.0f; for (uint i = tid; i < elems_per_group; i += tg_size) { local_sum += (float)x[base + i]; } local_sum = simd_sum(local_sum); threadgroup float shared_sums[32]; if (sid == 0) shared_sums[wid] = local_sum; threadgroup_barrier(mem_flags::mem_threadgroup); if (wid == 0) { float val = (sid < num_simdgroups) ? shared_sums[sid] : 0.0f; val = simd_sum(val); if (sid == 0) shared_sums[0] = val; } threadgroup_barrier(mem_flags::mem_threadgroup); float mean = shared_sums[0] / (float)elems_per_group; // Pass 2: Compute variance float local_var = 0.0f; for (uint i = tid; i < elems_per_group; i += tg_size) { float diff = (float)x[base + i] - mean; local_var += diff * diff; } local_var = simd_sum(local_var); if (sid == 0) shared_sums[wid] = local_var; threadgroup_barrier(mem_flags::mem_threadgroup); if (wid == 0) { float val = (sid < num_simdgroups) ? shared_sums[sid] : 0.0f; val = simd_sum(val); if (sid == 0) shared_sums[0] = val; } threadgroup_barrier(mem_flags::mem_threadgroup); float var = shared_sums[0] / (float)elems_per_group; float inv_std = metal::rsqrt(var + eps[0]); // Pass 3: Normalize, apply affine, apply erf-based GELU float rsqrt2 = 0.7071067811865475f; // 1/sqrt(2) for (uint i = tid; i < elems_per_group; i += tg_size) { uint c_local = i / spatial_size; uint c_global = group_idx * channels_per_group + c_local; float val = ((float)x[base + i] - mean) * inv_std; val = val * (float)weight[c_global] + (float)bias[c_global]; // Exact GELU: 0.5 * x * (1 + erf(x / sqrt(2))) val = 0.5f * val * (1.0f + erf_approx(val * rsqrt2)); out[base + i] = (T)val; } cxt/tjjdgddgtt atS)Nfused_groupnorm_gelurweightbiasepsrr)rrrheaderr)_groupnorm_gelu_kernelrrr_GROUPNORM_GELU_HEADER_GROUPNORM_GELU_SOURCErrr _get_groupnorm_gelu_kernelrTs;%!#!5!5'@)) "6"  "!rc8|jd|jd}}||z}|j|||g|jdd}ttd|j} |j | d} t j|| d} || z t j| |zz} | j|j} d|gdg|jdz zz}| |j|z|j|z} tj| S)z.Pure-MLX GroupNorm + GELU (no Metal required).rr$r!NTr#keepdims) r6r;tupler8r4meanrvarrsqrtnngelu)rrMrN num_groupsrOBCcpgx_raxesrYrZx_normx_outw_shapes r _groupnorm_gelu_fallbackrgs 771:qwwqzqA z/C !))Az3 5 5C q#((# $D 888 -D &&4$ /CDjBHHS3Y/ /F NN177 #E!fsaffqj))G FNN7+ +dll7.C CE 775>rc tst|||||S|j}|jd}|jd}d}|jddD]} || z} ||zdk7rtd|d|||z} | |z} | tkDrt|||||St j |j|||} t j |jt j}t j |jt j}t j|gt j} t j|| ||gt j}||z}tdtd | d zd zd z}t| ||| |gd |jfg||zddf|ddf|||fg|jg d}|j|S) zFused GroupNorm + GELU for NCL or NCHW layout. Equivalent to: x = groupnorm(x, weight, bias, num_groups, eps) x = gelu(x) rr$r!N channels  not divisible by num_groups r& r(r*)r1rgr6r7_HYBRID_THRESHOLDrr:r;astypefloat32r<r=r>maxrTr')rrMrNr^rO orig_shaper_r` spatial_sizedchannels_per_groupelems_per_groupx_contigeps_arrr total_groupstgrIs r rKrK s '64SIIJ  A  AL WWQR[  :~9QC'DZLQRRj(<7O**'64SII}}QYYq!\:;H ]]6==4 5F ==RZZ0 1DhhuBJJ/G XXz#5|QGrxx XFz>L T3r_r1b8B>? @B ) ' )&$8.!RA &AJ1l+,wwi  F >>* %%ra uint bg = threadgroup_position_in_grid.x; uint tid = thread_index_in_threadgroup; uint tg_size = threads_per_threadgroup.x; uint sid = thread_index_in_simdgroup; uint wid = simdgroup_index_in_threadgroup; uint num_simdgroups = tg_size / 32; uint num_groups = params[0]; uint channels_per_group = params[1]; uint spatial_size = params[2]; uint total_channels = params[3]; uint half_channels = params[4]; uint batch_idx = bg / num_groups; uint group_idx = bg % num_groups; uint elems_per_group = channels_per_group * spatial_size; uint base = batch_idx * total_channels * spatial_size + group_idx * channels_per_group * spatial_size; // Pass 1: Compute mean float local_sum = 0.0f; for (uint i = tid; i < elems_per_group; i += tg_size) { local_sum += (float)x[base + i]; } local_sum = simd_sum(local_sum); threadgroup float shared_sums[32]; if (sid == 0) shared_sums[wid] = local_sum; threadgroup_barrier(mem_flags::mem_threadgroup); if (wid == 0) { float val = (sid < num_simdgroups) ? shared_sums[sid] : 0.0f; val = simd_sum(val); if (sid == 0) shared_sums[0] = val; } threadgroup_barrier(mem_flags::mem_threadgroup); float mean = shared_sums[0] / (float)elems_per_group; // Pass 2: Compute variance float local_var = 0.0f; for (uint i = tid; i < elems_per_group; i += tg_size) { float diff = (float)x[base + i] - mean; local_var += diff * diff; } local_var = simd_sum(local_var); if (sid == 0) shared_sums[wid] = local_var; threadgroup_barrier(mem_flags::mem_threadgroup); if (wid == 0) { float val = (sid < num_simdgroups) ? shared_sums[sid] : 0.0f; val = simd_sum(val); if (sid == 0) shared_sums[0] = val; } threadgroup_barrier(mem_flags::mem_threadgroup); float var = shared_sums[0] / (float)elems_per_group; float inv_std = metal::rsqrt(var + eps[0]); // Intermediate buffer for normalized values (stored in shared memory would be // too large for typical channel counts, so we write to a temp output and read back). // Instead, we normalize on-the-fly and store in the output buffer. // Pass 3: Normalize, apply affine, then apply GLU // GLU: output[c, s] = a[c, s] * sigmoid(b[c, s]) // where a = first half of channels, b = second half // After GroupNorm+affine, the normalized values are in (B, 2C, S) layout. // We need to pair channel c with channel c + C. // For GLU, we iterate over the OUTPUT elements (half the channels). // Each output element at (c_out, s) reads: // a = norm_affine(x[base_a + c_out * S + s]) // b = norm_affine(x[base_b + c_out * S + s]) // where base_a/base_b account for the group offset. // But groups span the full 2C channels. We need to handle the case // where the GLU split crosses group boundaries. // For standard usage (num_groups divides 2C and C), the first half of groups // correspond to 'a' channels and second half to 'b' channels. // Actually for simplicity and correctness, this kernel should handle the common // case where num_groups divides both 2C and C (e.g., groups=1 or groups=4 with // 2C divisible by 4). The GLU split is on the channel dimension, independent of groups. // We compute all normalized values first, then apply GLU. // Strategy: each thread processes output elements. For each output (c_out, s), // compute both a and b normalized values on the fly. // Half channels per group on the output side: uint half_cpg = channels_per_group / 2; uint out_epg = half_cpg * spatial_size; uint out_base = batch_idx * half_channels * spatial_size + group_idx * half_cpg * spatial_size; for (uint i = tid; i < out_epg; i += tg_size) { uint c_local = i / spatial_size; uint s = i % spatial_size; // 'a' channel: first half of the group uint c_a = c_local; uint c_a_global = group_idx * channels_per_group + c_a; float val_a = ((float)x[base + c_a * spatial_size + s] - mean) * inv_std; val_a = val_a * (float)weight[c_a_global] + (float)bias[c_a_global]; // 'b' channel: second half of the group uint c_b = c_local + half_cpg; uint c_b_global = group_idx * channels_per_group + c_b; float val_b = ((float)x[base + c_b * spatial_size + s] - mean) * inv_std; val_b = val_b * (float)weight[c_b_global] + (float)bias[c_b_global]; // GLU: a * sigmoid(b) float sig_b = 1.0f / (1.0f + metal::exp(-val_b)); out[out_base + i] = (T)(val_a * sig_b); } cnt*tjjdgddgtatS)Nfused_groupnorm_glurLrr)_groupnorm_glu_kernelrrr_GROUPNORM_GLU_SOURCErrr _get_groupnorm_glu_kernelrs8$ " 4 4&@( !5!  ! rct|jd|jd}}||z}|j|||g|jdd}ttd|j} |j | d} t j|| d} || z t j| |zz} | j|j} d|gdg|jdz zz}| |j|z|j|z}t j|dd\}}|t j|zS)z-Pure-MLX GroupNorm + GLU (no Metal required).rr$r!NTrVr") r6r;rXr8r4rYrrZr[r2r3)rrMrNr^rOr_C_fullrarbrcrYrZrdrerfnormedr?r@s r _groupnorm_glu_fallbackrs  AGGAJvA J C !))Az3 5 5C q#((# $D 888 -D &&4$ /CDjBHHS3Y/ /F NN177 #E&kQC166A:..G V^^G, ,t||G/D DF 88FAA &DAq rzz!} rc |j}|jd}|jd}|dz}d} |jddD]} | | z} ||zdk7rtd|d||dzdk7rtd|d||z} | | z} tr|dkDs | tkDrt |||||St j |j||| } t j |jt j}t j |jt j}t j|gt j}t j|| | ||gt j}||z}td td | d zd zd z}t|}||d<t| ||||gd |j fg||zddf|ddf||| fg|j g d}|j|S)a-Fused GroupNorm + GLU for NCL or NCHW layout. GroupNorm normalizes over 2C channels, then GLU splits in half on axis=1. Equivalent to: x = groupnorm(x, weight, bias, num_groups, eps) a, b = split(x, 2, axis=1) return a * sigmoid(b) Input shape: (B, 2C, ...) -> Output shape: (B, C, ...) Note: Only supported when num_groups=1 or when channels_per_group is even (so the GLU split aligns with group boundaries). When num_groups > 1 and channels_per_group is odd, falls back to separate GroupNorm + GLU. rr$r!Nrirjz must be even for GLUr&rkrlrmr(r*)r6r7r1rnrrr:r;rorpr<r=r>rqr5rr')rrMrNr^rOrrr_rC_halfrsrtrurvrwrxrryrz out_shaperIs r r|r|s,J  A WWQZF q[FL WWQR[  a9VH,I*VWW zQ9VH,ABCC:-(<7O  Q/L T3r_r1b8B>? @BZ IIaL ( & (&$8.!RA &AJ6<01wwi  F >>) $$ra uint i = thread_position_in_grid.x; uint total = params[3]; if (i >= total) return; uint T_dim = params[0]; uint Fr_dim = params[1]; uint C_dim = params[2]; uint FrT = Fr_dim * T_dim; uint CFrT = C_dim * FrT; uint b = i / (2 * CFrT); uint rem = i % (2 * CFrT); uint c2 = rem / FrT; uint c = c2 / 2; uint is_imag = c2 % 2; uint ft = rem % FrT; uint in_idx = (b * CFrT + c * FrT + ft) * 2 + is_imag; out[i] = x[in_idx]; cnt*tjjdddgdgtatS)Ncomplex_to_interleavedrrrr)_complex_to_interleaved_kernelrrr_COMPLEX_TO_INTERLEAVED_SOURCErrr "_get_complex_to_interleaved_kernelrEs<%-)+)=)=)h1 *>* & *)rc |j\}}}}tsXtj|}tj|}tj ||gdj ||dz||Stj|tj}tj|j d}|dz|z|z|z} tj|||| gtj} t|| gg| ddftd| ddf| fgtjgd} | j |d|z||S) zConvert complex (B, C, Fr, T) to interleaved real (B, 2*C, Fr, T). Equivalent to: real = mx.real(z) imag = mx.imag(z) m = mx.stack([real, imag], axis=2).reshape(B, C*2, Fr, T) r!r"r%r&r$r)r*r)r6r1rrealimagstackr;viewrpr:r<r=rr>) zr_r`Frr(rrz_realz_flatrGrrIs r fused_complex_to_interleavedrQs%''KAq"a wwqzwwqzxxt 1-55aQAFF WWQ #F ]]6 " * *2 .F EAINQ E XXq"a'rxx 8F 1 / 1Q]e_a+xjzzl  F >>!QUB **r)returnr )r$)rmx.arrayr#intrr)gh㈵>) rrrMrrNrr^rrOfloatrr)rrrr)__doc__ __future__rmlx.corecorermlx.nnr\rr1rrrrrnrRrSrQrTcompilergrKr~r}rrr|rrrrrrr rs #" L  "  0j,DL "#'!)1 . 3&3& 3& 3& 3&  3&  3&zqf !#'!)1 0 E%E% E% E% E%  E%  E%X","& *+r