Skip to content

Bases

The Reduced-Basis Method seeks to reproduce a high-fidelity calculation to a desired accuracy with a set of equations in a lower dimensional space. The main step of the reduction is approximating the high-fidelity wave function by

$$ \phi_{\rm HF} \approx \hat{\phi} = \phi_0 + \sum_i c_i \tilde{\phi}_i~. $$

These classes calculate and store the basis state $\phi_0$ and $\tilde{\phi}_i$.

Basis(solver, theta_train, rho_mesh, n_basis)

Base class / template

Builds a reduced basis.

Parameters:

Name Type Description Default
solver SchroedingerEquation

high-fidelity solver

 required
theta_train ndarray

training space

 required
rho_mesh ndarray

discrete $s=kr$ mesh points

 required
n_basis int

number of states in the expansion

 required
l int

orbital angular momentum

 required

Attributes:

Name Type Description
solver SchroedingerEquation

high-fidelity solver
theta_train ndarray

training space
rho_mesh ndarray

discrete $s=kr$ mesh points
n_basis int

number of states in the expansion
l int

orbital angular momentum
Source code in src/rose/basis.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def __init__(
    self,
    solver: SchroedingerEquation,
    theta_train: np.array,
    rho_mesh: np.array,
    n_basis: int,
):
    r"""Builds a reduced basis.

    Parameters:
        solver (SchroedingerEquation): high-fidelity solver
        theta_train (ndarray): training space
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): number of states in the expansion
        l (int): orbital angular momentum

    Attributes:
        solver (SchroedingerEquation): high-fidelity solver
        theta_train (ndarray): training space
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): number of states in the expansion
        l (int): orbital angular momentum

    """
    self.solver = solver
    self.l = solver.interaction.ell
    self.theta_train = theta_train
    self.rho_mesh = rho_mesh
    self.n_basis = n_basis

load(filename) classmethod

Loads a previously saved Basis.

Source code in src/rose/basis.py 
17
18
19
20
21
22
@classmethod
def load(cls, filename):
    r"""Loads a previously saved Basis."""
    with open(filename, "rb") as f:
        basis = pickle.load(f)
    return basis

phi_exact(theta)

Exact wave function.

Parameters:

Name Type Description Default
theta ndarray

parameters

 required
l int)

partial wave

 required

Returns:

Name Type Description
phi ndarray

wave function
Source code in src/rose/basis.py 
70
71
72
73
74
75
76
77
78
79
80
81
def phi_exact(self, theta: np.array):
    r"""Exact wave function.

    Parameters:
        theta (ndarray): parameters
        l (int) : partial wave

    Returns:
        phi (ndarray): wave function

    """
    return self.solver.phi(theta, self.rho_mesh)

phi_hat(coefficients)

Emulated wave function.

Every basis should know how to reconstruct hat{phi} from a set of coefficients. However, this is going to be different for each basis, so we will leave it up to the subclasses to implement this.

Parameters:

Name Type Description Default
coefficients ndarray

expansion coefficients

 required

Returns:

Name Type Description
phi_hat ndarray

approximate wave function
Source code in src/rose/basis.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def phi_hat(self, coefficients):
    r"""Emulated wave function.

    Every basis should know how to reconstruct hat{phi} from a set of
    coefficients. However, this is going to be different for each basis, so
    we will leave it up to the subclasses to implement this.

    Parameters:
        coefficients (ndarray): expansion coefficients

    Returns:
        phi_hat (ndarray): approximate wave function

    """
    raise NotImplementedError

save(filename)

Saves a basis to file.

Parameters:

Name Type Description Default
filename string

name of file

 required
Source code in src/rose/basis.py 
83
84
85
86
87
88
89
90
91
def save(self, filename):
    """Saves a basis to file.

    Parameters:
        filename (string): name of file

    """
    with open(filename, "wb") as f:
        pickle.dump(self, f)

CustomBasis(solutions, phi_0, rho_mesh, n_basis, expl_var_ratio_cutoff=None, solver=None, subtract_phi0=True, use_svd=None, center=None, scale=None)

Bases: Basis

Builds a custom basis. Allows the user to supply their own.

$$ \phi_{\rm HF} \approx \hat{\phi} = \phi_0 + \sum_i c_i \tilde{\phi}_i~. $$

Parameters:

Name Type Description Default
solutions ndarray

HF solutions

 required
phi_0 ndarray

free solution (no interaction)

 required
rho_mesh ndarray

discrete $s=kr$ mesh points

 required
n_basis int

min number of states in the expansion

 required
expl_var_ratio_cutoff float)

the cutoff in sv2/sum(sv2), sv being the singular values, at which the number of kept bases is chosen

 None
use_svd bool

Use principal components for $\tilde{\phi}$?

 None

Attributes:

Name Type Description
solver SchroedingerEquation

not specified or assumed at construction
theta_train ndarray

not specified or assumed at construction
rho_mesh ndarray

discrete $s=kr$ mesh points
n_basis int

number of states in the expansion
phi_0 ndarray

free solution (no interaction)
solutions ndarray

HF solutions provided by the user
pillars ndarray

$\tilde{\phi}_i$
singular_values ndarray

singular values from SVD
vectors ndarray

copy of pillars
phi_0_interp interp1d

interpolating function for $\phi_0$
vectors_interp interp1d

interpolating functions for vectors (basis states)
Source code in src/rose/basis.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def __init__(
    self,
    solutions: np.array,  # HF solutions, columns
    phi_0: np.array,  # "offset", generates inhomogeneous term
    rho_mesh: np.array,  # rho mesh; MUST BE EQUALLY SPACED POINTS!!!
    n_basis: int,
    expl_var_ratio_cutoff: float = None,
    solver: SchroedingerEquation = None,
    subtract_phi0=True,
    use_svd: bool = None,
    center: bool = None,
    scale: bool = None,
):
    r"""Builds a custom basis. Allows the user to supply their own.

    $$
    \phi_{\rm HF} \approx \hat{\phi} = \phi_0 + \sum_i c_i \tilde{\phi}_i~.
    $$

    Parameters:
        solutions (ndarray): HF solutions
        phi_0 (ndarray): free solution (no interaction)
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): min number of states in the expansion
        expl_var_ratio_cutoff (float) : the cutoff in sv**2/sum(sv**2), sv
            being the singular values, at which the number of kept bases is chosen
        use_svd (bool): Use principal components for $\tilde{\phi}$?

    Attributes:
        solver (SchroedingerEquation): not specified or assumed at construction
        theta_train (ndarray): not specified or assumed at construction
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): number of states in the expansion
        phi_0 (ndarray): free solution (no interaction)
        solutions (ndarray): HF solutions provided by the user
        pillars (ndarray): $\tilde{\phi}_i$
        singular_values (ndarray): singular values from SVD
        vectors (ndarray): copy of `pillars`
        phi_0_interp (interp1d): interpolating function for $\phi_0$
        vectors_interp (interp1d): interpolating functions for vectors (basis states)

    """

    super().__init__(solver, None, rho_mesh, n_basis)

    self.rho_mesh = rho_mesh
    self.n_basis = n_basis
    self.phi_0 = phi_0

    self.pillars, self.singular_values, self.phi_0 = pre_process_solutions(
        solutions, self.phi_0, self.rho_mesh, center, scale, use_svd, subtract_phi0
    )

    # keeping at min n_basis PC's, find cutoff
    if expl_var_ratio_cutoff is not None:
        expl_var = self.singular_values**2 / np.sum(self.singular_values**2)
        n_basis_svs = np.sum(expl_var > expl_var_ratio_cutoff)
        self.n_basis = max(n_basis_svs, self.n_basis)
    else:
        self.n_basis = n_basis

    self.vectors = self.pillars[:, : self.n_basis]

percent_explained_variance()

Returns:

Type Description

(float) : percent of variance explained in the training set by the first n_basis principal

components
Source code in src/rose/basis.py 
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
def percent_explained_variance(self):
    r"""
    Returns:
        (float) : percent of variance explained in the training set by the first n_basis principal
        components
    """
    if self.singular_values is None:
        return 100
    else:
        return np.array(
            [
                100
                * np.sum(self.singular_values[:i] ** 2)
                / np.sum(self.singular_values**2)
            ]
            for i in range(self.nbasis)
        )

phi_hat(coefficients)

Emulated wave function.

Parameters:

Name Type Description Default
coefficients ndarray

expansion coefficients

 required

Returns:

Name Type Description
phi_hat ndarray

approximate wave function
Source code in src/rose/basis.py 
282
283
284
285
286
287
288
289
290
291
292
def phi_hat(self, coefficients):
    r"""Emulated wave function.

    Parameters:
        coefficients (ndarray): expansion coefficients

    Returns:
        phi_hat (ndarray): approximate wave function

    """
    return self.phi_0 + np.sum(coefficients * self.vectors, axis=1)

project(x)

Return projection of x onto vectors

Source code in src/rose/basis.py 
294
295
296
297
298
299
300
301
302
303
def project(self, x):
    r"""
    Return projection of x onto vectors
    """
    x -= self.phi_0
    x /= np.trapz(np.absolute(x), self.rho_mesh)
    return [
        np.trapz(self.vectors[:, i].conj() * x, self.rho_mesh)
        for i in range(self.n_basis)
    ]

RelativeBasis(solver, theta_train, rho_mesh, n_basis, expl_var_ratio_cutoff=None, phi_0_energy=None, use_svd=True, center=None, scale=None)

Bases: Basis

Builds a "relative" reduced basis. This is the default choice.

$$ \phi_{\rm HF} \approx \hat{\phi} = \phi_0 + \sum_i c_i \tilde{\phi}_i~. $$

Parameters:

Name Type Description Default
solver SchroedingerEquation

high-fidelity solver

 required
theta_train ndarray

training space

 required
rho_mesh ndarray

discrete $s=kr$ mesh points

 required
n_basis int

number of states in the expansion

 required
l int

orbital angular momentum

 required
use_svd bool

Use principal components for $\tilde{\phi}$?

 True
phi_0_energy float

energy at which $\phi_0$ is calculated

 None

Attributes:

Name Type Description
solver SchroedingerEquation

high-fidelity solver
theta_train ndarray

training space
rho_mesh ndarray

discrete $s=kr$ mesh points
n_basis int

number of states in the expansion
l int

orbital angular momentum
phi_0 ndarray

free solution (no interaction)
pillars ndarray

$\tilde{\phi}_i$
singular_values ndarray

singular values from SVD
vectors ndarray

copy of pillars
Source code in src/rose/basis.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def __init__(
    self,
    solver: SchroedingerEquation,
    theta_train: np.array,
    rho_mesh: np.array,
    n_basis: int,
    expl_var_ratio_cutoff: float = None,
    phi_0_energy: float = None,
    use_svd: bool = True,
    center: bool = None,
    scale: bool = None,
):
    r"""Builds a "relative" reduced basis. This is the default choice.

    $$
    \phi_{\rm HF} \approx \hat{\phi} = \phi_0 + \sum_i c_i \tilde{\phi}_i~.
    $$

    Parameters:
        solver (SchroedingerEquation): high-fidelity solver
        theta_train (ndarray): training space
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): number of states in the expansion
        l (int): orbital angular momentum
        use_svd (bool): Use principal components for $\tilde{\phi}$?
        phi_0_energy (float): energy at which $\phi_0$ is calculated

    Attributes:
        solver (SchroedingerEquation): high-fidelity solver
        theta_train (ndarray): training space
        rho_mesh (ndarray): discrete $s=kr$ mesh points
        n_basis (int): number of states in the expansion
        l (int): orbital angular momentum
        phi_0 (ndarray): free solution (no interaction)
        pillars (ndarray): $\tilde{\phi}_i$
        singular_values (ndarray): singular values from SVD
        vectors (ndarray): copy of `pillars`

    """

    super().__init__(
        solver,
        theta_train,
        rho_mesh,
        n_basis,
    )

    if phi_0_energy is not None:
        k = np.sqrt(2 * self.solver.interaction.mu * phi_0_energy / HBARC**2)
        eta = self.solver.interaction.k_c / k
    else:
        if isinstance(self.solver.interaction, EnergizedInteractionEIM):
            # k_mean = np.sqrt(2*self.solver.interaction.mu*np.mean(theta_train[:, 0])/HBARC**2)
            # eta = self.solver.interaction.k_c / k_mean
            # Does not support Coulomb (yet).
            eta = 0.0
        else:
            # If the phi_0 energy is not specified, we're only going to work
            # with non-Coulombic systems (for now).
            eta = 0.0

    # Returns Bessel functions when eta = 0.
    self.phi_0 = np.array(
        [coulombf(self.l, eta, rho) for rho in self.rho_mesh], dtype=np.complex128
    )
    self.solutions = np.array([self.phi_exact(theta) for theta in theta_train]).T

    self.pillars, self.singular_values, self.phi_0 = pre_process_solutions(
        self.solutions, self.phi_0, self.rho_mesh, center, scale, use_svd
    )

    # keeping at min n_basis PC's, find cutoff
    if expl_var_ratio_cutoff is not None:
        expl_var = self.singular_values**2 / np.sum(self.singular_values**2)
        n_basis_svs = np.sum(expl_var > expl_var_ratio_cutoff)
        self.n_basis = max(n_basis_svs, self.n_basis)
    else:
        self.n_basis = n_basis

    self.vectors = self.pillars[:, : self.n_basis].copy()

percent_explained_variance()

Returns:

Type Description

(float) : percent of variance explained in the training set by the first n_basis principal

components
Source code in src/rose/basis.py 
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def percent_explained_variance(self):
    r"""
    Returns:
        (float) : percent of variance explained in the training set by the first n_basis principal
        components
    """
    if self.singular_values is None:
        return 100
    else:
        return np.array(
            [
                100
                * np.sum(self.singular_values[:i] ** 2)
                / np.sum(self.singular_values**2)
            ]
            for i in range(self.nbasis)
        )

phi_hat(coefficients)

Emulated wave function.

Parameters:

Name Type Description Default
coefficients ndarray

expansion coefficients

 required

Returns:

Name Type Description
phi_hat ndarray

approximate wave function
Source code in src/rose/basis.py 
176
177
178
179
180
181
182
183
184
185
186
def phi_hat(self, coefficients):
    r"""Emulated wave function.

    Parameters:
        coefficients (ndarray): expansion coefficients

    Returns:
        phi_hat (ndarray): approximate wave function

    """
    return self.phi_0 + np.sum(coefficients * self.vectors, axis=1)

project(x)

Return projection of x onto vectors

Source code in src/rose/basis.py 
188
189
190
191
192
193
194
195
196
197
def project(self, x):
    r"""
    Return projection of x onto vectors
    """
    x -= self.phi_0
    x /= np.trapz(np.absolute(x), self.rho_mesh)
    return [
        np.trapz(self.vectors[:, i].conj() * x, self.rho_mesh)
        for i in range(self.n_basis)
    ]