From 6722fafe50c4548eb745c8c941da9b60720c97dc Mon Sep 17 00:00:00 2001 From: Guillaume <sagitta1618@gmail.com> Date: Thu, 4 Jul 2024 11:37:41 +0200 Subject: [PATCH] Move sequence generator to sequence.py and refactor the code --- ohmpi/ohmpi.py | 42 ++-- ohmpi/sequence.py | 541 ++++++++++++++-------------------------------- 2 files changed, 172 insertions(+), 411 deletions(-) diff --git a/ohmpi/ohmpi.py b/ohmpi/ohmpi.py index 8deca508..a7e510bb 100644 --- a/ohmpi/ohmpi.py +++ b/ohmpi/ohmpi.py @@ -251,47 +251,31 @@ class OhmPi(object): Each tuple is the form (<array_name>, param1, param2, ...) Dipole spacing is specified in terms of "number of electrode spacing". Dipole spacing is often referred to 'a'. Number of levels is a multiplier - of 'a', often refere to 'n'. For multigradient array, an additional parameter + of 'a', often referred to 'n'. For multigradient array, an additional parameter 's' is needed. Types of sequences available are : - ('wenner', a) - ('dpdp', a, n) - ('schlum', a, n) - ('multigrad', a, n, s) - ireciprocal : bool, optional + By default, if an integer is provided for a, n and s, the parameter + will be considered varying from 1 to this value. For instance, for + ('wenner', 3), the sequence will be generated for a = 1, a = 2 and a = 3. + If only some levels are desired, the user can use a list instead of an int. + For instance ('wenner', [3]) will only generate quadrupole for a = 3. + include_reciprocal : bool, optional If True, will add reciprocal quadrupoles (so MNAB) to the sequence. + ip_optimization: bool, optional + If True, will optimize for induced polarization measurement (i.e. will + try to put as much time possible between injection and measurement at + the same electrode). fpath : str, optional Path where to save the sequence (including filename and extension). By default, sequence is saved in ohmpi/sequences/sequence.txt. """ - # dictionnary of function to create sequence - fdico = { - 'dpdp': dpdp1, - 'wenner': wenner, - 'schlum': schlum1, - 'multigrad': multigrad, - } - # check arguments - if fpath is None: - fpath = os.path.join(os.path.dirname(__file__), '../sequences/sequence.txt') - qs = [] - - # create sequence - for p in params: - pok = [int(p[i]) for i in np.arange(1, len(p))] # make sure all are int - qs.append(fdico[p[0]](nelec, *pok).values.astype(int)) - quad = np.vstack(qs) - if len(quad.shape) == 1: # only one quadrupole - quad = quad[None, :] - - # add reciprocal - if ireciprocal: - quad = np.r_[quad, quad[:, [2, 3, 0, 1]]] - self.sequence = quad - - # save sequence + dfseq = create_sequence(nelec, params=params, include_reciprocal=include_reciprocal) + self.sequence = dfseq.astype(int).values np.savetxt(fpath, self.sequence, delimiter=' ', fmt='%d') - print('{:d} quadrupoles generated.'.format(self.sequence.shape[0])) @staticmethod def _find_identical_in_line(quads): diff --git a/ohmpi/sequence.py b/ohmpi/sequence.py index 50dd89d7..a5066cd4 100644 --- a/ohmpi/sequence.py +++ b/ohmpi/sequence.py @@ -1,440 +1,217 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Nov 14 16:05:59 2018 - -@author: pmclachlan -originally produced for the ResIPy project (https://gitlab.com/hkex/resipy) -""" import numpy as np import pandas as pd +import os + +def create_sequence(nelec, params=[('dpdp', 1, 8)], include_reciprocal=False, ip_optimization=False): + """Creates a sequence of quadrupole.Several type of + sequence or sequence with different parameters can be combined together. -def dpdp1(elec_num, a, n): - ''' Generates quadrupole matrix for dipole-dipole survey. - Parameters ---------- - elec_num : int - Number of electrodes - a : int or list of int - Spacing between C and V pairs in electrode spacing - (a = 1 is the same as skip 0). - n : int or list of int - Quadrupole separation in electrode spacing. - - Notes - ----- - Length of `a` should match `n`. - ''' - elec_id = np.arange(elec_num)+1 - abmn = [] - if isinstance(a, list) is False: - if a < 1: - raise ValueError('a must be >= 1, it\'s the electrode spacing between C and V pairs (a = 1 is the same as skip 0)') - n_ = np.array(range(0,n))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - B = A + a - M = B + n_[j] * a - N = M + a - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - if a[i] < 1: - raise ValueError('a must be >= 1, it\'s the electrode spacing between C and V pairs (a = 1 is the same as skip 0)') - n_ = np.array(range(0,n[i]))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - B = A + a[i] - M = B + n_[j] * a[i] - N = M + a[i] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,3] <= elec_num] - return proto_mtrx - - -def dpdp2(elec_num, a, n): + nelec : int + Number of electrodes. + params : list of tuple, optional + Each tuple is the form (<array_name>, param1, param2, ...) + Dipole spacing is specified in terms of "number of electrode spacing". + Dipole spacing is often referred to 'a'. Number of levels is a multiplier + of 'a', often referred to 'n'. For multigradient array, an additional parameter + 's' is needed. + Types of sequences available are : + - ('wenner', a) + - ('dpdp', a, n) + - ('schlum', a, n) + - ('multigrad', a, n, s) + By default, if an integer is provided for a, n and s, the parameter + will be considered varying from 1 to this value. For instance, for + ('wenner', 3), the sequence will be generated for a = 1, a = 2 and a = 3. + If only some levels are desired, the user can use a list instead of an int. + For instance ('wenner', [3]) will only generate quadrupole for a = 3. + include_reciprocal : bool, optional + If True, will add reciprocal quadrupoles (so MNAB) to the sequence. + ip_optimization: bool, optional + If True, will optimize for induced polarization measurement (i.e. will + try to put as much time possible between injection and measurement at + the same electrode). + """ + # dictionary of function to create sequence + fdico = { + 'dpdp': dpdp, + 'wenner': wenner, + 'schlum': schlum, + 'multigrad': multigrad, + } + + # check parameters (all int and > 0) + for param in params: + for i, p in enumerate(param[1:]): + if isinstance(p, float) or isinstance(p, int): + if p <= 0: + raise ValueError('parameters of sequence (a, n or s) needs to be > 0') + else: + for pp in p: + if pp <= 0: + raise ValueError('parameters of sequence (a, n or s) needs to be > 0') + + # create sequence + dfs = [] + for param in params: + dfs.append(fdico[param[0]](nelec, *param[1:])) + dfseq = pd.concat(dfs, axis=0) + + # add reciprocal + if include_reciprocal: + dfseqr = dfseq.copy().rename(columns={'a': 'm', 'b': 'n', 'm': 'a', 'n': 'b'}) + dfseqr = dfseqr.sort_values(['a', 'b', 'm', 'n']) + dfseq = pd.concat([dfseq, dfseqr]).reset_index(drop=True) + + # optimize for IP + if ip_optimization: + print('NOT IMPLEMENTED YET') + + # stats + print('{:d} quadrupoles generated.'.format(dfseq.shape[0])) + return dfseq + +def dpdp(nelec, a, n): ''' Generates quadrupole matrix for dipole-dipole survey. Parameters ---------- - elec_num : int + nelec : int Number of electrodes a : int or list of int - Spacing between C and V pairs in electrode spacing + Spacing between AB (current) and MN (voltage) pairs in electrode spacing (a = 1 is the same as skip 0). n : int or list of int - Quadrupole seperation in electrode spacing. - - Notes - ----- - Length of `a` should match `n`. - ''' - elec_id = np.arange(elec_num)+1 - abmn = [] - if isinstance(a, list) is False: - if a < 1: - raise ValueError('a must be >= 1, it\'s the electrode spacing between C and V pairs (a = 1 is the same as skip 0)') - n_ = np.array(range(0,n))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - B = A + a - M = B + n_[j] - N = M + a - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - if a[i] < 1: - raise ValueError('a must be >= 1, it\'s the electrode spacing between C and V pairs (a = 1 is the same as skip 0)') - n_ = np.array(range(0,n[i]))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - B = A + a[i] - M = B + n_[j] - N = M + a[i] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,3] <= elec_num] - return proto_mtrx - - -def wenner(elec_num, n): - ''' Generate quadrupoles matrix for Wenner alpha survey for a = 1 ... n. + Quadrupole separation in electrode spacing. ''' - seq = [] - for a in np.arange(n) + 1: - seq.append(wenner_alpha(elec_num, a)) - return pd.concat(seq).reset_index(drop=True) - - -def wenner_alpha(elec_num, a): - ''' Generates quadrupole matrix for Wenner alpha survey. - - Parameters - ---------- - elec_num : int - Number of electrodes - a : int or list of int - Spacing between electrodes (in electrode spacing). - ''' + # check parameters + if isinstance(a, int): + a = np.arange(a) + 1 + if isinstance(n, int): + n = np.arange(n) + 1 - elec_id = np.arange(elec_num)+1 + elec = np.arange(nelec) + 1 abmn = [] - if isinstance(a, list) is False: - A = elec_id - M = A + a - N = M + a - B = N + a - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - A = elec_id - M = A + a[i] - N = M + a[i] - B = N + a[i] + for aa in a: + if aa < 1: + raise ValueError('a must be >= 1, it is the electrode spacing between AB and MN pairs (a = 1 is the same as skip 0)') + for nn in n: + if nn < 1: + raise ValueError('n must be >= 1, it is the number of level between AB and MN') + A = elec + B = A + aa + M = B + nn * aa + N = M + aa abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return(proto_mtrx) - - -def wenner_beta(elec_num, a): - ''' Generates quadrupole matrix for Wenner beta survey. + abmn = np.vstack(abmn) + abmn = abmn[(abmn <= nelec).all(1), :] + df = pd.DataFrame(abmn, columns=['a', 'b', 'm', 'n']) + df = df.sort_values(['a', 'b', 'm', 'n']).reset_index(drop=True) + return df + +def wenner(nelec, a): + '''Generates quadrupole matrix for Wenner alpha survey. Parameters ---------- - elec_num : int + nelec : int Number of electrodes a : int or list of int Spacing between electrodes (in electrode spacing). ''' - - elec_id = np.arange(elec_num)+1 + if isinstance(a, int): + a = np.arange(a) + 1 + elec = np.arange(nelec) + 11 abmn = [] - - if type(a) == int: - B = elec_id - A = B + a - M = A + a - N = M + a + for aa in a: + A = elec + M = A + aa + N = M + aa + B = N + aa abmn.append(np.vstack([A, B, M, N]).T) - - if type(a) == list: - for i in np.array(range(0,len(a))): - B = elec_id - A = B + a[i] - M = A + a[i] - N = M + a[i] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return(proto_mtrx) + abmn = np.vstack(abmn) + abmn = abmn[(abmn <= nelec).all(1), :] + df = pd.DataFrame(abmn, columns=['a', 'b', 'm', 'n']) + df = df.sort_values(['a', 'b', 'm', 'n']).reset_index(drop=True) + return df -def wenner_gamma(elec_num, a): - ''' Generates quadrupole matrix for Wenner gamma survey. - - Parameters - ---------- - elec_num : int - Number of electrodes - a : int or list of int - Spacing between electrodes (in electrode spacing). - ''' - - elec_id = np.arange(elec_num)+1 - abmn = [] - if isinstance(a, list) is False: - A = elec_id - M = A + a - B = M + a - N = B + a - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - A = elec_id - M = A + a[i] - B = M + a[i] - N = B + a[i] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return(proto_mtrx) - - -def schlum1(elec_num, a, n): +def schlum(nelec, a, n): ''' Generates quadrupole matrix for Schlumberger survey. Parameters ---------- - elec_num : int + nelec : int Number of electrodes a : int or list of int Spacing between electrodes (in electrode spacing). n : int or list of int - Quadrupole seperation in electrode spacing. + Quadrupole separation in electrode spacing. ''' - elec_id = np.arange(elec_num)+1 - abmn = [] - if isinstance(a, list) is False: - n_ = np.array(range(0,n))+1 - for j in np.array(range(0,len(n_))): - if n_[j] >= a: - #print(n_[j],a) - A = elec_id - M = A + n_[j] - N = M + a - B = N + n_[j] - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - n_ = np.array(range(0,n[i]))+1 - for j in np.array(range(0,len(n_))): - if n_[j] >= a[i]-1: #minus 1 here to permit measurements on level 2 of pseudosection - A = elec_id - M = A + n_[j] - N = M + a[i] - B = N + n_[j] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return(proto_mtrx) - - -#def schlum1(elec_num, a, n): -# ''' Generates quadrupole matrix for Schlumberger survey. -# Parameters -# ---------- -# elec_num : int -# Number of electrodes -# a : int or list of int -# Spacing between electrodes (in electrode spacing). -# n : int or list of int -# Quadrupole seperation in electrode spacing. -# ''' -# elec_id = np.arange(elec_num)+1 -# proto_mtrx = pd.DataFrame([]) -# if isinstance(a, list) is False: -# n_ = np.array(range(0,n))+1 -# for j in np.array(range(0,len(n_))): -# if n_[j] >= a: -# #print(n_[j],a) -# A = elec_id -# M = A + n_[j] -# N = M + a -# B = N + n_[j] -# # if j == 0: -# # proto_mtrx = pd.DataFrame(np.column_stack((A, B, M, N))) -# # if j > 0: -# new_rows = pd.DataFrame(np.column_stack((A, B, M, N))) -# proto_mtrx = proto_mtrx.append(new_rows) - -# else: -# for i in np.array(range(0,len(a))): -# n_ = np.array(range(0,n[i]))+1 -# for j in np.array(range(0,len(n_))): -# if n_[j] >= a[i]: -# A = elec_id -# M = A + n_[j] -# N = M + a[i] -# B = N + n_[j] -# -# new_rows = pd.DataFrame(np.column_stack((A, B, M, N))) -# proto_mtrx = proto_mtrx.append(new_rows) - -# proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] -# return(proto_mtrx) + # check parameters + if isinstance(a, int): + a = np.arange(a) + 1 + if isinstance(n, int): + n = np.arange(n) + 1 - -#def schlum1(elec_num, a, n): -# ''' Generates quadrupole matrix for Schlumberger survey. - -# Parameters -# ---------- -# elec_num : int -# Number of electrodes -# a : int or list of int -# Spacing between electrodes (in electrode spacing). -# n : int or list of int -# Quadrupole seperation in electrode spacing. -# ''' -# elec_id = np.arange(elec_num)+1 - -# if isinstance(a, list) is False: -# n_ = np.array(range(0,n))+1 -# for j in np.array(range(0,len(n_))): -# A = elec_id -# M = A + n_[j] * a -# N = M + a -# B = N + n_[j] * a -# if j == 0: -# proto_mtrx = pd.DataFrame(np.column_stack((A, B, M, N))) -# if j > 0: -# new_rows = pd.DataFrame(np.column_stack((A, B, M, N))) -# proto_mtrx = proto_mtrx.append(new_rows) -# -# else: -# for i in np.array(range(0,len(a))): -# n_ = np.array(range(0,n[i]))+1 -# for j in np.array(range(0,len(n_))): -# A = elec_id -# M = A + n_[j] * a[i] -# N = M + a[i] -# B = N + n_[j] * a[i] -# if (i + j) == 0: -# proto_mtrx = pd.DataFrame(np.column_stack((A, B, M, N))) -# if (i + j) > 0: -# new_rows = pd.DataFrame(np.column_stack((A, B, M, N))) -# proto_mtrx = proto_mtrx.append(new_rows) -# -# proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] -# return(proto_mtrx) - - -def schlum2(elec_num, a, n): - ''' Generates quadrupole matrix for Schlumberger survey. - - Parameters - ---------- - elec_num : int - Number of electrodes - a : int or list of int - Spacing between electrodes (in electrode spacing). - n : int or list of int - Quadrupole seperation in electrode spacing. - ''' - - elec_id = np.arange(elec_num)+1 + elec = np.arange(nelec) + 1 abmn = [] - if isinstance(a, list) is False: - n_ = np.array(range(0,n))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - M = A + n_[j] * a - N = M + a - B = N + n_[j] * a - abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - n_ = np.array(range(0,n[i]))+1 - for j in np.array(range(0,len(n_))): - A = elec_id - M = A + n_[j] - N = M + a[i] - B = N + n_[j] + for aa in a: + for nn in n: + if nn >= aa-1: #minus 1 here to permit measurements on level 2 of pseudosection + A = elec + M = A + nn + N = M + aa + B = N + nn abmn.append(np.vstack([A, B, M, N]).T) + abmn = np.vstack(abmn) + abmn = abmn[(abmn <= nelec).all(1), :] + df = pd.DataFrame(abmn, columns=['a', 'b', 'm', 'n']) + df = df.sort_values(['a', 'b', 'm', 'n']).reset_index(drop=True) + return df - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return proto_mtrx - - -def multigrad(elec_num, a, n, s): - ''' Genetrate measurement matrix for multigradient array Torleif Dahlin. +def multigrad(nelec, a, n, s): + ''' Generate measurement matrix for multigradient array. Parameters ---------- - elec_num : int + nelec : int Number of electrodes - a : int + a : int or list of int Spacing between potential electrodes (in electrode spacing). - n : int + n : int or list of int Multiplier for `a` to determine spacing from A to M. - s : int - Seperation factor for current electrodes, should be the intermediate + s : int or list of int + Separation factor for current electrodes, should be the intermediate numbers. ''' - elec_id = np.arange(elec_num)+1 + # check parameters + if isinstance(a, int): + a = np.arange(a) + 1 + if isinstance(n, int): + n = np.arange(n) + 1 + if isinstance(s, int): + n = np.arange(s) + 1 + elec = np.arange(nelec) + 1 abmn = [] - - if isinstance(a, list) is False: - n_ = np.array(range(0, n))+1 - s_ = np.array(range(0, s))+1 - for j in np.array(range(0, len(n_))): - for k in np.array(range(0, len(s_))): - A = elec_id - B = A + s_[k] + 2 - M = A + n_[j] * a - N = M + a + for aa in a: + for nn in n: + for ss in s: + A = elec + B = A + ss + 2 + M = A + nn * aa + N = M + aa abmn.append(np.vstack([A, B, M, N]).T) - - else: - for i in np.array(range(0,len(a))): - n_ = np.array(range(0, n[i]))+1 - s_ = np.array(range(0, s[i]))+1 - for j in np.array(range(0, len(n_))): - for k in np.array(range(0, len(s_))): - A = elec_id - B = A + s_[k] + 2 - M = A + n_[j] * a[i] - N = M + a[i] - abmn.append(np.vstack([A, B, M, N]).T) - - proto_mtrx = pd.DataFrame(np.vstack(abmn)) - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] > proto_mtrx.iloc[:,3]] - proto_mtrx = proto_mtrx[proto_mtrx.iloc[:,1] <= elec_num] - return proto_mtrx + abmn = np.vstack(abmn) + abmn = abmn[(abmn <= nelec).all(1), :] + df = pd.DataFrame(abmn, columns=['a', 'b', 'm', 'n']) + df = df.sort_values(['a', 'b', 'm', 'n']).reset_index(drop=True) + return df # test code # x1 = dpdp1(24, 2, 8) -# x2 = dpdp2(24, 2, 8) # x3 = wenner_alpha(24, 1) -# x4 = wenner_beta(24, 1) -# x5 = wenner_gamma(24, 1) # x6 = schlum1(24, 1, 10) -# x7 = schlum2(24, 1, 10) # x8 = multigrad(24, 1, 10, 2) -- GitLab