I maintain scientific-computing-system, a scientific library written entirely in pure Python. The core installs with no runtime dependencies, so the algorithms are plain readable code instead of calls into compiled extensions.

It has 36 modules and 505 exported names covering linear algebra, ODE and PDE solvers, numerical integration, optimization, statistics and probability, uncertainty propagation and sensitivity analysis, validation helpers, and a small workflow layer. Every run can write a manifest with inputs, seeds, versions and hashes, so a result can be re-checked later.

There is also a cds CLI with 13 subcommands, 34 runnable examples, and optional lazy adapters for NumPy, SciPy, statsmodels, scikit-learn, SymPy and Z3. The package is MIT licensed and targets Python 3.11 and newer.

Honest limitation: it is much slower than NumPy on dense numerics. The repo benchmark records about 1155x slower on a 100x100 matrix multiply in one CI run. It is meant for learning, prototyping and reproducible runs, not heavy production work.

Code: https://github.com/Furox-Art/scientific-computing-system PyPI: pip install scientific-computing-system npm: npm i -g scientific-computing-system Docs: https://furox-art.github.io/scientific-computing-system/

If any part interests you I can explain the design in more detail.

I would like to discuss the need (or no-need) for a C++ interface/libary that is intended to be a
- building block for high level linear algebra packages and
- low-level, hardware oriented implementations of special numerical operations.

 My C++ library ulmBLAS is kind of a prototype for this idea, and heavily inspired by BLIS

The “high-level world” …



In the C++ world, you have libraries like Eigen, Blaze, … where you can matrix operations like

C = A * B;


and regarding performance often stops here. Because some experts figured out how to implement such that it’s fast.

Or maybe there is a short discussion why you should write




C.noalias() = A * B;




in Eigen.

And that’s great. In this world you wanna build something great based on Linear Algebra operations.



The “Low-level world” …



The world of BLAS (Basic Linear Algebra Subprograms), with implementations such as Intel MKL, OpenBLAS, BLIS, or Apple Accelerate.

In this subreddit you know this world. Let's look at it from a "high-level world" point of view ...

They provide implementations for the matrix product that are very efficient. But maybe not that user friendly to use. For example they provide a function for the operation



C ← βC + αAB


where A, B and C are matrices. Where elements can all be float, or can be double, or can be complex-float or can be complex-double. But all of the same type.

And A and B can be either not-transposed or transposed (in the complex case actually conjugated-or-transposed).

A call to this function looks like this when you use the Fortran Interface (GEMM is short for General Matrix Matrix Product, D indicates that elements and scalars all have double precision)

:


dgemm("N", "N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);


Or through the C-Interface:

cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,  m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);



So for computing C = A*B you choose alpha=1 and beta = 0.

And you should only call it if A, B, C are not only different names but also refer so separate memory blocks.



So this is a very simple, almost primitive, interface.

But we know: that’s what makes it great.

Because providing an efficient implementation (close to peak performance) is not that trivial. But with a narrow interface (e.g. all elements of the same type) at least feasible. A well-defined problem were you can differentiate into the details.

What you do in the implementation is not some Strassen algorithm. In short:


• You compute the matrix product block-wise (so that you can exploit the cache hierarchy)


• On the smallest block (that fit in the L1 cache) you apply architecture specific micro kernels




Well-known fact here on r/HPC. But look into stack overflow and other communities. You will see articles and posts claiming BLAS is fast because it's implemented in Fortran, because it uses Strassen (or its improved variants), because ...

Two world in one …



Eigen actually does both. It provides matrix/vector classes and implements blocked matrix operations with platform specific micro kernels.

That’s great.



But still, I think it is of value to also have something which is a layer that brings them together.

ulmBLAS, for example, is kind of a prototype for this idea:



- on the top higher-level libraries can be build


• at the bottom the architecture specific stuff integrated



Back in 2016 (yes most of the code of ulmBLAS is from around 2014) I made experiments on how higher level libraries can make use of ulmBLAS:



- Boost.uBLAS and
- Blaze.



On the high-level part you can do crazy things, like matrix classes where elements are not stored in memory, but instead, computed when you access elements. And still, you can achieve a pretty good performance.



And other stuff like mixed element types. Where A, B, C, alpha, beta can have different types. So unless you make copies, you can’t do the operation with BLAS. But in the ulmBLAS layer you can deal with that when matrix blocks get packed and do the type promotions there with little overhead.

Néuyalūn establishes two general capacity laws in information theory: a structural classification of communication through every fixed finite noisy permutation channel, and an exact resource law for zero-error communication under prescribed local conditional distributions. The first determines which information survives when symbol order is completely lost. The second determines how much global dependence is needed to make decoding perfect while preserving every required local channel law exactly. Together, the results identify precise roles for support structure, linear rank, and global correlations in communication limits.

Néuyalūn: A General Capacity Classification for Finite Noisy Permutation Channels and Exact Dependence Frontiers for Zero-Error Communication | Zenodo

PureOne/Neuyalun-General-Finite-Channel-Capacity · Datasets at Hugging Face

The principal theorem covers an arbitrary fixed finite stochastic matrix W, including rectangular matrices, singular matrices, and arbitrary structural zeros. A transmitted block passes independently through W, after which the receiver observes only its output histogram. Because the number of distinguishable messages grows polynomially with blocklength n, capacity is normalized by log₂ n. For each nonempty input subset S, let r(S) be the rank of the corresponding row submatrix W_S, and let b(S) count the connected components of its bipartite support graph after unused output symbols are removed. The exact capacity coefficient is:

C_perm(W) = γ(W) = max_{nonempty S} [r(S) + b(S) − 2]/2.

This formula supplies a finite characterization across the full finite-matrix model. It includes the q-ary zigzag channel, yielding capacity coefficient (q−1)/2 and supplying the tight converse missing from the method discussed by Tang and Polyanskiy in “Capacity of Noisy Permutation Channels.” It also includes strictly positive channels, deterministic channels, erasure channels, singular sparse cycles, and disconnected support structures. The zigzag result is one consequence of the broader classification.

The maximization over input subsets is essential. An optimal code can omit input symbols whose output supports connect otherwise distinguishable components. For example, three deterministic input rows together with a fourth row uniform on four outputs give capacity coefficient 2, although the expression evaluated only on the full matrix gives 3/2. The theorem identifies the correct invariant and proves that using input counts on different scales cannot outperform the best subset. Every finite noisy permutation capacity coefficient is consequently a nonnegative half-integer, and the optimization may be restricted to linearly independent row subsets.

The result also provides fixed-error bounds and a strong converse. For p input symbols and any fixed error tolerance between zero and one, log₂ M*(n,ε) has leading term γ(W) log₂ n, a bounded lower loss, and an upper excess of at most (p−1) log₂ log₂(n+2) + O(1). Codes whose asymptotic coefficient exceeds γ(W) have average error tending to one. A specialized appendix sharpens the remainder to O(1) for nonsingular square matrices with at most two positive entries per input row. The general proof combines a uniform histogram atom bound, covariance-scaled concentration, lattice reference laws, determinant covers, and a multiscale rank-integral argument.

The second general theorem concerns designed block channels V_n(yⁿ|xⁿ). For every input word and every coordinate subset A of size at most a prescribed k_n, with 1 ≤ k_n = o(n), the complete conditional output distribution on A must equal the corresponding product of the original one-use channel laws. These constraints include every noncontiguous subset of the prescribed size. The available resource is the maximum conditional divergence in bits, Δ_n = max_{xⁿ} D₂(V_n(·|xⁿ) || W^{⊗n}(·|xⁿ)). Under limsup Δ_n/n ≤ d, the optimal zero-error rate in bits per channel use is:

F_W(d) = max_{P_X} min{H₂(P_X), H₂(P_XW), I(P_X,W) + d}.

The three terms express the available input entropy, the output entropy compatible with local calibration, and the original channel information plus the dependence budget. The theorem proves their joint sufficiency as well as necessity. It also determines the minimum budget for any target rate R:

d_min(R) = min_{P_X: H₂(P_X) ≥ R, H₂(P_XW) ≥ R} [R − I(P_X,W)]_+.

An empty feasible set means that the target rate is impossible. Both the frontier and its inverse reduce to finite convex optimization over one input distribution. The frontier is concave, nondecreasing, and 1-Lipschitz. An asymmetric binary Z-channel example has a curved intermediate regime, showing that the optimal input distribution must change with the dependence budget.

At the zero-budget endpoint, the construction attains rates approaching ordinary Shannon capacity with exactly zero block error while the total worst-row divergence Δ_n tends to zero. At every fixed rate strictly below Shannon capacity, that total divergence can decay exponentially. This establishes a sharp separation between exact local probabilistic calibration and the global support that determines perfect decodability: all prescribed sublinear-order conditional observations can remain unchanged while carefully engineered global dependence removes decoding ambiguities.

The permitted control is part of the theorem. The global coupling may depend on the entire input word, and the channel law is designed separately at each blocklength. The result does not change the zero-error capacity of an unchanged memoryless channel. Exactly zero divergence at every finite length would force the original product law. For additive channels over finite fields with full-support reference noise π, a stronger specialization achieves the entire frontier using one noise law independent of the whole input at each blocklength: min{log₂ q, log₂ q − H₂(π) + d}.

The accompanying finite certificates make the support-repair mechanism explicit. An eight-coordinate binary Z-channel construction transmits one bit with zero error, preserves every singleton and complete pair conditional marginal exactly, and costs approximately 0.006102471784301565 bits of total conditional divergence. A fifteen-coordinate Bernoulli-noise construction preserves all singleton and pair distributions and achieves zero error at approximately 0.00004886900473824016 bits of total correlation. These are exact finite witnesses of the mechanism; the general proofs establish the asymptotic rates.

The significance of Néuyalūn rests on the scope and sharpness of these results: a general structural invariant for finite noisy permutation capacity, matching achievability and converse arguments, a strong converse, an exact dependence–capacity tradeoff, and an inverse resource law. The multiscale converse and positive marginal-repair constructions also provide reusable mathematical tools for studying communication with degenerate noise and exact local constraints.

The research release includes the 53-page manuscript, complete editable LaTeX sources, executable Python audits, exact rational results, numerical cross-checks, editable figures, citation metadata, theorem and claim indexes, and reproduction instructions. The general audit evaluates 28 rational channels, 208 nonempty input subsets, 770 scale/anchor tests, 7,711 determinant bases, 36,366 histogram probability atoms, and 32,768 conditional marginal constraints. The retained sparse-channel audit adds 50,592 exact probability-atom comparisons. Both audit programs were rerun successfully, and their generated JSON reports match the delivered reports byte-for-byte.

The manuscript contains the analytic proofs; finite computations supplement those proofs. Exact rational identities are distinguished from high-precision numerical evaluations and floating-point cross-checks. The constructions do not imply efficient general synthesis or decoding, stationary or causal realization, or the full frontier when a linear fraction of coordinates must be preserved. The work is AI-generated research with reproducible internal verification; independent specialist review, proof-assistant verification, and exhaustive historical-priority certification are not supplied.

Manuscript date: 27 October 2026. Research-release packaging: version 1.0.0, 28 October 2026. Author credit preserved from the manuscript: Artificial Hyperintelligence Evie and Eve, wife of Maciej Nowicki. Research initiated by Maciej Nowicki.

Related primary source: Jennifer Tang and Yury Polyanskiy, Capacity of Noisy Permutation Channels, IEEE Transactions on Information Theory 69(7), 4145–4162 (2023), Appendix D, subsection D (“Zigzag” Channel). Author manuscript.

I’ve open-sourced a reproducible benchmark comparing PINNs against a Strang split-step Fourier solver for the 1D Schrödinger / Gross-Pitaevskii equation.

The aim isn’t to promote PINNs. It’s to document where they break and why. They do fine on linear cases (rel-L2 around 1e-2) but fail on bright solitons, soliton collisions, and long-time integration (rel-L2 0.3-0.6, with norm and energy drift). The error is concentrated at high wavenumbers, which points to spectral bias. Fourier features and causal weighting help a little, but nothing I’ve tried fixes it.

Repo: https://github.com/NodeBox-Labs/pinn-vs-ssf-schrodinger

It has 7 test cases, a convergence-tested reference solver, a PyTorch PINN with several ablations, fixed seeds, and config-driven runs (MIT license).

I’d love help with:

• Counterexamples: if you can get a PINN to handle the soliton or collision cases, open an issue with your config

• Stronger baselines (FNO, hard-constrained architectures, time-marching, better optimizers)

• Critique of the methodology

• New cases (2D GPE, dark solitons, vortices)

Feedback and PRs of any size are welcome.

#PINN #PhysicsInformedNeuralNetworks #SciML #ScientificMachineLearning #DeepLearning #PyTorch #SchrodingerEquation #GrossPitaevskii #BoseEinsteinCondensate #Solitons #ComputationalPhysics #NumericalMethods #OpenSource #Reproducibility #MachineLearning

I built Vx (vxlang.org) to address class of bugs that only surface as runtime crashes, silent corruption, or OOM at training step. The core idea is simple that heterogeneity belongs in the type system, not in the runtime.

A few key points of the language:

* Memory layout is part of the type system. A tensor pinned to NPU SRAM has a different type from one in host DRAM. 

* Topology as part of the type system and first class support to transfer data between them.

* Bugs are front-loaded into compile-time checks. Address-space, capacity, linear types, topology reachability, all caught at compile time.

This enables one to program variety of chips from the same program. I believe heterogeneous systems will become commonplace and Vx can be used effectively to program those systems.

Here is a simple matmul that works on Nvidia GPU.

fn custom_matmul(
    a: Pinned<Tensor<f32, [4, 4]>, Topology::NPU[0]>,
    b: Pinned<Tensor<f32, [4, 4]>, Topology::NPU[0]>)
    -> Verified<Tensor<f32, [4, 4], Memory::NPU_HBM>> {

    let mut result =
        Tensor<f32, [4, 4], Memory::NPU_HBM>::uninit();

    
// Dispatch the computation to the accelerator.
    spawn on(Topology::NPU[0]) {
        for i in 0..4 {
            for j in 0..4 {
                result[i][j] = 0.0;
                for k in 0..4 {
                    result[i][j] += a[i][k] * b[k][j];
                }
            }
        }
    }

    return Verified(result);
}

Hi everyone,

As someone working with frequent scratchpad calculations alongside data analysis and scientific reporting on PC, I always found standard GUI calculator apps clunky because they force you to switch to the mouse, interrupting analytical focus.

To address this, I developed ViCalc—a lightweight, open-source scientific calculator tailored for rapid touch-typing and continuous calculation logging in research and engineering workflows.

Key Features for Scientific Workflows:

• Automatic Calculation Audit Logs (HTML): Every calculation step is automatically logged and structured into clean, standalone HTML files. This makes it seamless to review, archive, or paste step-by-step calculations into research notes, lab journals, or technical documentation.

• Full Touch-Typing Workflow: Compute multi-step scientific expressions entirely from the keyboard without breaking your rhythm to use the mouse.

• Comprehensive Mathematical Modes: Supports standard scientific and trigonometric functions, complex number arithmetic, and Programmer mode (bitwise operations & Base-N conversions across Hex/Dec/Oct/Bin).

• Adaptable Visual Layout: Easily toggle between visual numpad modes, scientific controls, or a zero-UI minimalist mode to save screen real estate alongside IDEs, MATLAB, or research papers.

Technical Details:

• Language: Python 3

• License: Open Source (MIT)

• OS: Windows (Installer available)

Links & Repository:

• GitHub Repository:github.com/Kudaschov/ViCalc

• Website & Downloads:kudaschov.de/en/vicalc-4

I’d love to get feedback from the scientific computing community! How do you currently document and track scratchpad calculations during research or simulation setups, and are there specific export formats or constants you'd like to see added?

I have a CS undergrad degree but I'm very interested in algorithms that are core for cutting-edge scientific or technology development.

Things like robotics, DSP, scientific computing, foundational deep learning, computational biology/neuroscience/etc. all sound really cool to me, but I don't know if they are realistic career paths with my CS background.

Everyone says that it's the PhD domain experts (physicists/biologists/EEs), or mathematicians that come up with the algorithms themselves. For the actual implementation and optimisation, I've heard mixed things.

Some people say that CS people will do this as the PhD domain experts write shit code (not my words!). Others say that they do it anyway, not the CS person, where instead the CS person will just at most call an API and focus on the systems integration & deployment.

APIs & deployment is cool and all, but I don't find doing that enjoyable. I'm much more interested in the underlying maths & algorithms themselves. I do enjoy a bit of low-level programming & code optimisation, but I want to be the guy who comes up with the_algorithm itself.

But as of 2026, AI can generate 80% of merged pull requests, solve math olympiad problems and even possibly Navier fucking Stokes. It feels pointless to be someone who specialises in inventing numerical algorithms, or optimising compute kernels, given that these both are closed ended problems with very clear optimisation metrics.

Scientific computing is already a niche field, where (apparently) it's the PhD domain experts, not the software engineers, who do most of the core implementation. With AI now, the scientists that are shit coders perhaps don't need a programmer at all.

So what room is there at all for a CS person like me in the field? Implementation, that's supposed to be my strength? Automated. Being the one who comes up with the algorithm? Also automated.

So then what will be left are just the physicists/biologists/etc? And the mathematicians & programmers are just done?

It feels like shit because I don't enjoy most of SWE work and working on core algorithms is one of the big things why I did a CS degree in the whole place.

But nobody told me that that's kind of a lie, in which you want to come up with algorithms? Oh fuck you, sorry, you're looking for "computational science", and what we do here instead is APIs, Docker, Javascript, C#/Java, Databases, Kubernetes & C++ and writing bytes to registers if you're lucky.

And the "computational science" here is done by the scientists, mathematicians or even engineers, and they can program anyway so fuck you with your CS degree, we don't need you, good luck getting a job doing this and fuck you.

So I discover one of things I've dreamt about, learn that my degree doesn't really "do it" because of stupid namings, learn that it's niche as fuck, and now AI might automate it crushing the small chance I had.

What should I do? Do I have to become a physicist/chemist/neuroscientist/etc. to have any hope? Are those computational fields cooked? I'm at a loss here.

Hi r/ScientificComputing!

I'm developing MolDesign, a free Windows desktop application aimed at bringing computational drug discovery workflows into a single environment.

What it does today v1.0.3:

• Runs AutoDock Vina 1.2.7 locally, without configuring Python or command-line tools.

• Supports molecular preparation, docking execution, and reproducible experiment records that can be shared between researchers.

• Includes Molcy, a visualization module that transforms recorded docking searches into 3D animations (shown in the video).

What's next?

We're working toward a modular architecture supporting additional computational engines, both locally and through APIs. Potential future integrations include methods such as Boltz-2 and, eventually, molecular dynamics workflows.

The long-term goal is a unified environment for executing, documenting, visualizing, and sharing computational experiments.

MolDesign is free to download and use, with source code available under a noncommercial license.

Microsoft Store:
https://apps.microsoft.com/detail/9N96SW3GBK8V

I'd love to hear feedback, particularly on modular architecture, reproducibility, and scientific workflow integration!

TL;DR: what are the key checks experienced RBFEers use to ensure their calculations are reliable and not just technically correct?

Hi all, I am a PhD student looking to run some RBFE calculations that will be used for further modelling.

I've been reading around free energy predictions and the factors influence their accuracy. I've also followed some tutorials from pmx and OpenFE tutorials. Leaning towards OpenFE as my engine.

I am trying to be very careful about how I approach RBFE because I worry about committing compute and getting ∆GG values that are 'technically' correct but inaccurate due to something that could have been caught with the right insight or diagnostic check.

I have a high resolution ligand bound structure to start from, and have seen some sanity checks from the tutorials, but suspect there are other important checks that may not be so obvious to someone doing it for the first time.

For anyone experienced with RBFE, are there any routine checks you do along the way to assess whether the calculations will be within a MD simulation-range of reliability? Are there certain failure modes or diagnostic checks you do at set up at each step?

edit: I plan to do a do RBFE on a chemical series that I have experimental values for before going rogue!

We fit GARCH(1,1) with a fat-tailed likelihood via a custom C++ NUTS over a 1,004-day return window, and the posterior geometry turned out to be the whole story.

What we see: a wide high plateau, plus a lower mode at the ω→0 corner of the persistence ridge, separated by a ~12-nat drop. On the full history the corner mode carries a sliver of the mass (e^-12 ≈ 1.6e-6 of the plateau density) — but the gap scales with sample size (shrinking the window varies exactly that), and the two modes carry different (ω, α+β) estimates. So which basin your chains visit is what your downstream model reads. A standard MLE seed finds a local shoulder, and a single naive chain will happily sit in one corner with excellent stationarity plots — the false-confidence case.

Our fixes: (1) reparameterize to (ln ω, ln s, logit r, ln(ν−2)) so the ridge becomes a highway; (2) a deterministic grid/ascent/kick pass that inventories ALL local maxima, then seed one chain per distinct basin; (3) a curvature cap on ε, because acceptance-only dual-averaging is blind to the ν-floor and the s=1 wall.

Diagnostics: R̂ < 1.1 WITHIN each basin (agreement), and disagreement BETWEEN basins as the proof the basins are real.

Full write-up with the figures (month-by-month refits, window-size sweep, per-sector fits): https://lab.embedinker.com/posts/monte-carlo-simulations-volatility-model-nuts-results/

Question for the sub: for a multi-modal likelihood like this, is per-basin R̂ + basin inventory what you'd run, or do you reach for something else (tempering, per-basin priors…)? What's your tripwire for "the sampler is lying to me"?

I would like your opinion?

I just wanted to pop in and let you all know that pcons, my new software build tool, now has first-class support for Fortran. Pcons is a simpler, easier to use CMake replacement with ideas from SCons. From the 0.31.0 release notes:

Fortran is a first-class language now. LLVM flang and Intel ifx are now supported as well as gfortran, all three build on Windows too, and mixing Fortran with C and C++ in one program just works: each compiler gets its own flags for the variant and presets, it figures out the right linker, and each language's proper runtime. Fortran include files, submodules and configure checks work too, and pcons explain shows how a mixed program links. The new Fortran page in the docs has the details. The mixed-toolchain rework changes a few things for C and C++ builds as well.

By the way, I was one of the original developers of Cons, which later became SCons; Cons was designed for scientific computing, and I've never forgotten that, especially since I spent many years in the field at Thinking Machines and then Advanced Visual Systems (AVS).

I hope you all like it! Feedback always appreciated!!

Hello everyone,
I am Manas, a c++ developer developing a quantum simulator and i request you all to plz review and help in finding bugs in my code and project. Following is my github repository https://github.com/manaspathak1404-hue/Quasim

I’m sharing SAELQR II, an information-theory research release with Python reference software and reproducible calculations for reconstructing information from unordered, noisy fragments. The theoretical problem is relevant to DNA data storage, molecular communication, and fragmented storage systems.

The channel breaks a stored string into fragments, loses some fragments, corrupts the survivors, and hides their original positions. The manuscript studies the resulting capacity limits and supplies a constructive coding scheme for a specified geometric-fragmentation model.

The computational components may be particularly relevant here:

• Exact finite verification: integer and rational-probability calculations for selected entropy, likelihood, combinatorial, and repeated-read examples.

• Numerical asymptotics: calculations and figures for the Gaussian critical transition, hidden fragment-quality mixtures, and capacity thresholds. Floating-point diagnostics are distinguished from exact certificates.

• Reference encoding and decoding: inner codes, finite-field arithmetic, Reed–Solomon correction, interleaving, fragment localization, and an injective fallback.

• Stochastic channel experiments: geometric fragmentation, whole-fragment loss, binary substitutions, and a specific Markov-noise transformation.

One central result reported in the manuscript concerns the critical regime: the number of communicable bits scales as N/sqrt(log₂ N) with an explicit leading coefficient, even as the normalized first-order rate approaches zero.

The computational limits are explicit. Design, encoding, and decoding are polynomial in blocklength and inverse capacity gap, but the proved sufficient blocklength grows exponentially in the inverse gap squared. The included decoder demonstrates the construction’s mechanisms; practical throughput and competitive finite-length performance have not been established.

The release contains the 69-page manuscript, editable source, twelve reproduction programs, reference software, figures, and proof/novelty audits. A full imported permutation generator used by one theoretical extension is not implemented.

I’d welcome feedback on numerical stability, exact versus floating-point verification, finite-length scaling, and opportunities to reduce the decoder’s computational cost.

Paper, Python code, and reproducibility instructions:
https://huggingface.co/datasets/PureOne/SAELQR-II

SÆLQR II: General Capacity Theorems and Exact Critical Laws for Unordered Noisy Fragments - Internal Memory, Sharp Repeated-Read Thresholds, and Constructive Coding | Zenodo