add backprop in elecsolver

README.md

@@ -10,29 +10,40 @@ This repository is **not** a general-purpose electrical system solver. Instead,
 - The graph-based description of an electric network
 - The corresponding sparse linear system to solve
 
+In a very simple way, ElecSolver takes as an input the Resistances $R$, Capacitances $C$, Inductances $L$, Mutuals $M$, Current sources $I$  and Voltage sources $V$ along with the connectivity graph $G$ of the system and outputs the linear system: matrix $S$ and vector $b$ to solve in order to get the solution of the electric problem.
+
+$$f_{\text{ElecSolver}}(R,C,L,M,I,V,G) = (S,b)$$
+
+It is then up to the user to choose the solver they want for solving the system $Sx=b$ and to manage the time iterations if needed for temporal problems.
+
+The main goal of ElecSolver is to provide a friendly Python interface for simulating and optimizing analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
+
 
-Its main goal is to provide a friendly Python interface for simulating analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
 
 
 > [!IMPORTANT]
 > ElecSolver has been designed with the following specifications in mind:
 > - The time needed for building the linear system must be negligible compared to the time needed for solving it.
-> - Handle natively inductive mutuals and resistive mutuals
-> - Handle as many coupled electric systems that one wants.
-> - Deal with lonely nodes and lonely edges in the electric graph: the problem can be well posed and thus solved.
+> - Handle natively inductive mutuals and resistive mutuals (skin effects).
+> - Handle as many coupled electric systems as one wants.
+> - Deal with lonely nodes and lonely edges in the electric graph: the problem can be well-posed and thus solved.
+> - Allow backpropagation of gradients through the system for optimization purposes.
 
 
 
 > [!NOTE]
 > Non-linear components are not supported. You must manage event detection and system updates yourself.
 
 
-## Table of content
+## Table of contents
 
 - [ElecSolver](#elecsolver)
   - [Overview](#overview)
-  - [Table of content](#table-of-content)
+  - [Table of contents](#table-of-contents)
   - [How to install](#how-to-install)
+    - [Using pip](#using-pip)
+    - [Using conda/mamba](#using-condamamba)
+    - [From source](#from-source)
   - [Components](#components)
     - [FrequencySystemBuilder](#frequencysystembuilder)
       - [Features](#features)
@@ -44,20 +55,33 @@ Its main goal is to provide a friendly Python interface for simulating analog el
   - [Solver suggestions](#solver-suggestions)
   - [Extra uses: Hydraulic or Thermal system modeling](#extra-uses-hydraulic-or-thermal-system-modeling)
   - [Netlist import feature](#netlist-import-feature)
+  - [Gradient backpropagation](#gradient-backpropagation)
+    - [Example of backpropagation of TemporalSystemBuilder](#example-of-backpropagation-of-temporalsystembuilder)
+    - [Example of backpropagation of FrequencySystemBuilder](#example-of-backpropagation-of-frequencysystembuilder)
+    - [Conclusion on gradient backpropagation](#conclusion-on-gradient-backpropagation)
+
 
 ## How to install
-For now this package is distributed on pypi and can be installed using pip and conda/mamba
+### Using pip
+ElecSolver is distributed on PyPI and can be installed using pip or conda/mamba.
 ```
 pip install ElecSolver
 ```
-or
+### Using conda/mamba
+ElecSolver is also available on conda-forge and can be installed using conda or mamba:
 ```
 conda install elecsolver
 ```
-For solving the linear systems we advise using MUMPS through `python-mumps` available for linux, macOS and Windows. It can be installed via conda:
+For solving linear systems, we advise using MUMPS through `python-mumps`, available for Linux, macOS, and Windows. It can be installed via conda:
 ```
 conda install python-mumps
 ```
+### From source
+You can also install ElecSolver from source by cloning the repository and running:
+```
+pip install .
+```
+you will need to have `numpy`, `scipy` and `networkx` installed in your environment
 
 
 
@@ -74,6 +98,7 @@ This class handles **frequency-domain** analysis of linear electric systems.
 - Detects and couples multiple subsystems
 - Accepts arbitrary complex impedances and mutuals
 - Constructs sparse linear systems (COO format)
+- Allows backpropagation of gradients from the constructed linear system (`S` and `rhs`) to the parameters of the system for optimization purposes
 
 
 > [!TIP]
@@ -84,7 +109,7 @@ This class handles **frequency-domain** analysis of linear electric systems.
 We would like to study the following system:
 ![Multiple system](docs/img/schema.png)
 
-this can simply be defined in the following manner (We took R=1, L=1 and M=2):
+This can simply be defined in the following manner (we take $R=1$, $L=1$, and $M=2$):
 ```python
 import numpy as np
 from scipy.sparse.linalg import spsolve
@@ -112,7 +137,7 @@ electric_sys = FrequencySystemBuilder(
 # Add source (Current source here)
 electric_sys.add_current_source(intensity=10, input_node=2, output_node=0)
 # Set ground
-# 2 values because one for each subsystem
+# Two values because one is needed for each subsystem
 electric_sys.set_ground(0, 3)
 # Building system
 electric_sys.build_system()
@@ -127,18 +152,18 @@ frequencial_response = electric_sys.build_intensity_and_voltage_from_vector(sol)
 print(frequencial_response.potentials[3]-frequencial_response.potentials[4])
 ```
 #### Adding a Parallel Resistance
-We want to add components in parallel with existing components for instance inserting a resistor in parallel with the first inductance (between nodes 0 and 2)
+We want to add components in parallel with existing components, for instance by inserting a resistor in parallel with the first inductance (between nodes 0 and 2).
 ![Parallel system](docs/img/schema3.png)
 
-In python, simply add the resistance to the list of impedence in the very first lines of the script:
+In Python, simply add the resistance to the list of impedence in the very first lines of the script:
 
 ```python
 import numpy as np
 from scipy.sparse.linalg import spsolve
 from ElecSolver import FrequencySystemBuilder
 
 
-# We add an additionnal resistance between 0 and 2
+# We add an additional resistance between 0 and 2
 impedence_coords = np.array([[0,0,1,3,0],[1,2,2,4,2]], dtype=int)
 impedence_data = np.array([1, 1j,1, 1j,1], dtype=complex)
 
@@ -160,14 +185,15 @@ This class models **time-dependent** systems using resistors, capacitors, coils,
 - Detects and couples multiple subsystems
 - Accepts 3 dipole types: resistances, capacities and coils
 - Constructs sparse linear systems (COO format)
+- Allows backpropagation of gradients from the constructed linear systems (plural! gradients from `S_init`, `S1`, `S2` and `rhs`) to the parameters of the system for optimization purposes
 
 #### Example
 
 
 We would like to study the following system:
 ![Temporal system](docs/img/schema2.png)
 
-with R=1, L=0.1, C=2 this gives:
+With $R=1$, $L=0.1$, and $C=2$, this gives:
 
 ```python
 import numpy as np
@@ -243,13 +269,13 @@ This outputs the following graph that displays the intensity passing through the
 
 ## Extra uses: Hydraulic or Thermal system modeling
 
-This repository can be used as is in order to model the mass flow or thermal flux in respectively Hydraulic networks or Thermal networks where a difference of pressure or a difference of temperature can be assimilated to a tension source. Since electric potentials are always computed relatively to the ground node you might need to rescale the resulting potentials:
+This repository can be used as is to model mass flow or thermal flux in hydraulic or thermal networks, respectively, where a pressure difference or a temperature difference can be assimilated to a tension source. Since electric potentials are always computed relative to the ground node, you may need to rescale the resulting potentials:
 
 We are considering the following hydraulic problem:
 
 ![Hydraulic system](docs/img/hydraulic.png)
 
-Taking R=1 this gives
+Taking $R=1$ this gives
 
 ```python
 import numpy as np
@@ -304,7 +330,7 @@ print("Debit through the system",solution.intensities_sources[0])
 ```
 ## Netlist import feature
 
-A new class, named NetlistParser allows importing passive netlist and building a TemporalSystem instance.
+A new class named NetlistParser allows importing a passive netlist and building a TemporalSystemBuilder instance.
 Solving the system can then be performed like any other example above.
 
 ```python
@@ -376,4 +402,130 @@ plt.plot(arange(steps, dtype=float)*dt, voltage_src, label="V(1-0)") # equal to
 plt.legend()
 plt.savefig("intensities_res.png")
 ```
+## Gradient backpropagation
+
+As we said in the introduction, ElecSolver is roughly a function that takes as an input the graph-based description of an electric system and outputs the linear system to solve:
+
+$$f_{\text{ElecSolver}}(R,C,L,M,I,V,G) = (S,b)$$
+
+
+In order to allow optimization of the system, ElecSolver now gives the possibility to backpropagate gradients from the linear system to the parameters of the system (i.e., by giving $\frac{\partial R}{\partial S}$ , $\frac{\partial R}{\partial b}$ , $\frac{\partial C}{\partial S}$ , $\frac{\partial C}{\partial b}$ , etc.). This allows one to use ElecSolver in an optimization loop where one wants to optimize the parameters of the system with gradient descent, for instance.
+
+### Example of backpropagation of TemporalSystemBuilder
+In this example, we want to optimize the capacity values of the system in order to reach a measured target current for `t=0.8`.
+We can do this by backpropagating the gradients from the solution of the linear system to the capa_data array of the TemporalSystemBuilder instance and then performing gradient descent on capa_data.
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import TemporalSystemBuilder
+## Simple tetrahedron
+res_coords  = np.array([[0,2],[1,3]],dtype=int)
+res_data = np.array([1,1],dtype=float)
+
+coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+coil_data = np.array([1,1],dtype=float)
+
+capa_coords = np.array([[1,2],[3,0]],dtype=int)
+capa_data = np.array([1,1],dtype=float)
+## The target solution we want to reach is the solution of the system when capa_data = np.array([0.1,1],dtype=float)
+
+## mutuals
+mutuals_coords=np.array([[],[]],dtype=int)
+mutuals_data = np.array([],dtype=float)
+
+
+res_mutuals_coords=np.array([[],[]],dtype=int)
+res_mutuals_data = np.array([],dtype=float)
+
+## initializing system
+elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+elec_sys.add_current_source(10,1,0)
+elec_sys.set_ground(0)
+elec_sys.build_system()
+
+# Getting initial condition system
+S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+# Getting temporal systems
+S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+# Solving initial condition system
+sol_init =spsolve(S_i.tocsr(),rhs.todense())
+## Time iteration with euler implicit scheme for 1 timestep
+dt=0.8
+B = rhs*dt+S2@sol_init
+A = S2+dt*S1
+sol = spsolve(A,B)
+## Target solution (artificially made by setting capa_data = np.array([0.1,1],dtype=float))
+sol_target = [ 3.24786325, -1.1965812,  -6.16809117,  0.61253561,  0.58404558, -2.63532764, 0., 6.16809117,  2.10826211,  1.4957265 ]
+for i in range(1000):
+
+    ## computing gradients
+    dB = 2*spsolve(A.T, sol - sol_target)
+    ## chain rule for gradients of capa_data (S2 appears twice in the computation graph)
+    dS2 = -( dB[S2.row]*sol[S2.col]) ## Gradient from the solver (A\B) to S2
+    dS2 += (dB[S2.row]*sol_init[S2.col]) ## Gradient from the solver (A\B) to S2 through B
+
+    ## Backpropagate gradients from dS2 to capa_data
+    gradients = elec_sys.backpropagate_gradients(dS2=dS2)
+
+    ## change the values of capa_data using gradient descent
+    elec_sys.capa_data = elec_sys.capa_data - 0.01*gradients.capa_data
+    ## After the update of capa_data we need to rebuild the system to update S1, S2 and rhs
+    elec_sys.build_system()
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    ## recomputing solution with euler implicit scheme
+    B = rhs*dt+S2@sol_init
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(elec_sys.capa_data, np.array([0.1,1],dtype=float))
+```
+
+The function `backpropagate_gradients` of TemporalSystemBuilder allows to backpropagate gradients from `S_init`, `S1`, `S2` and `rhs`. See tests.test_gradients for more examples of backpropagation of gradients from different systems.
+
+### Example of backpropagation of FrequencySystemBuilder
+We want to optimize the value of a tension source in order to reach a target solution for the frequency response of the system.
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import FrequencySystemBuilder
+
+## sparse python res matrix
+impedence_coords = np.array([[0,0,1],[1,2,2]],dtype=int)
+impedence_data = np.array([1,1,1],dtype=complex)
+
+## mutuals
+mutuals_coords=np.array([[0],[1]],dtype=int)
+mutuals_data = np.array([2.j],dtype=complex)
+
+
+electric_sys = FrequencySystemBuilder(impedence_coords,impedence_data,mutuals_coords,mutuals_data)
+## target solution is the solution of the system when voltage=5
+electric_sys.add_voltage_source(voltage=10,input_node=1,output_node=0)
+# setting the ground
+electric_sys.set_ground(0)
+electric_sys.build_system()
+
+## getting system
+sys,b = electric_sys.get_system(sparse_rhs=True)
+sol = spsolve(sys.tocsr(),b.todense())
+## Target solution (Artificially made by setting voltage_source_data = np.array([5],dtype=complex))
+sol_target = np.array([-1.66666667+1.66666667j, -0.83333333+1.66666667j, 0.83333333-1.66666667j, 2.5-3.33333333j, 0.+0.j, 5+0.j, 4.16666667+1.66666667j])
+for i in range(3000):
+    ## computing gradients of the squared error between sol and sol_target with respect to the second member b of the system (the one that contains the voltage source data)
+    db = 2*spsolve(sys.tocsr().conj().T, sol - sol_target)
+    drhs = db[b.row]
+    ## Backpropagate gradients from drhs to voltage_source_data
+    gradients = electric_sys.backpropagate_gradients(drhs=drhs)
+    ## Performing gradient descent on voltage_source_data
+    electric_sys.voltage_source_data = electric_sys.voltage_source_data - 0.01*gradients.voltage_source_data
+    ## After the update of voltage_source_data we need to rebuild the system to update sys and b
+    electric_sys.build_system()
+    sys,b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(),b.todense())
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(electric_sys.voltage_source_data, np.array([5],dtype=complex))
+```
+The function `backpropagate_gradients` of FrequencySystemBuilder allows to backpropagate gradients from `S`, and `rhs`. See tests.test_gradients for more examples of backpropagation of gradients from different systems.
 
+### Conclusion on gradient backpropagation
+Although it provides the backpropagation feature, ElecSolver does not provide an automatic differentiation mechanism. Many libraries can be used to this extent, such as `autograd`, `jax`, `PyTorch` and many more, to avoid the hassle of computing gradients manually. ElecSolver aims to limit its dependencies to keep its versatility and remain a simple, user-friendly tool.

docs/components/frequency-system-builder.md

@@ -87,3 +87,76 @@ impedence_data = np.array([1, 1j, 1, 1j, 1], dtype=complex)
 mutuals_coords = np.array([[1], [3]], dtype=int)
 mutuals_data = np.array([2.0j], dtype=complex)
 ```
+
+## Gradient Backpropagation
+
+`FrequencySystemBuilder` can backpropagate gradients from `S` and `rhs` to model parameters.
+
+This enables gradient-based optimization loops directly on source values or component data.
+
+### Example: Optimize a Voltage Source
+
+In this example, we optimize the voltage source value so the system response matches a target solution.
+
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import FrequencySystemBuilder
+
+## sparse python res matrix
+impedence_coords = np.array([[0, 0, 1], [1, 2, 2]], dtype=int)
+impedence_data = np.array([1, 1, 1], dtype=complex)
+
+## mutuals
+mutuals_coords = np.array([[0], [1]], dtype=int)
+mutuals_data = np.array([2.0j], dtype=complex)
+
+electric_sys = FrequencySystemBuilder(
+    impedence_coords,
+    impedence_data,
+    mutuals_coords,
+    mutuals_data,
+)
+
+# Target solution is the solution of the system when voltage=5
+electric_sys.add_voltage_source(voltage=10, input_node=1, output_node=0)
+electric_sys.set_ground(0)
+electric_sys.build_system()
+
+## Getting system
+sys, b = electric_sys.get_system(sparse_rhs=True)
+sol = spsolve(sys.tocsr(), b.todense())
+
+# Target solution (artificially made by setting voltage_source_data = np.array([5], dtype=complex))
+sol_target = np.array([
+    -1.66666667 + 1.66666667j,
+    -0.83333333 + 1.66666667j,
+    0.83333333 - 1.66666667j,
+    2.5 - 3.33333333j,
+    0.0 + 0.0j,
+    5.0 + 0.0j,
+    4.16666667 + 1.66666667j,
+])
+
+for _ in range(3000):
+    ## Computing gradients of squared error with respect to b
+    db = 2 * spsolve(sys.tocsr().conj().T, sol - sol_target)
+    drhs = db[b.row]
+
+    ## Backpropagate gradients from drhs to voltage_source_data
+    gradients = electric_sys.backpropagate_gradients(drhs=drhs)
+    ## Performing gradient descent on voltage_source_data
+    electric_sys.voltage_source_data = (
+        electric_sys.voltage_source_data - 0.01 * gradients.voltage_source_data
+    )
+
+    ## After updating voltage_source_data, rebuild the system to update sys and b
+    electric_sys.build_system()
+    sys, b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(), b.todense())
+
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(electric_sys.voltage_source_data, np.array([5], dtype=complex))
+```
+
+For additional backpropagation examples, see `tests/test_gradients.py`.