improved lowerpass

lowpass_oracle.cpp

@@ -0,0 +1,164 @@
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+using namespace std;
+
+typedef vector<vector<double>> Arr;
+typedef pair<Arr, double> Cut;
+
+Arr ndarray(int rows, int cols) {
+    return vector<vector<double>>(rows, vector<double>(cols));
+}
+
+class LowpassOracle {
+private:
+    bool more_alt;
+    Arr A;
+    int nwpass;
+    int nwstop;
+    double Lpsq;
+    double Upsq;
+
+public:
+    LowpassOracle(Arr A, int nwpass, int nwstop, double Lpsq, double Upsq) {
+        this->A = A;
+        this->nwpass = nwpass;
+        this->nwstop = nwstop;
+        this->Lpsq = Lpsq;
+        this->Upsq = Upsq;
+        this->more_alt = true;
+        this->i_A = 0;
+        this->count = 0;
+    }
+
+    pair<pair<Arr, double>, double> assess_optim(Arr x, double Spsq) {
+        int n = x.size();
+        this->more_alt = true;
+
+        // case 2, passband constraints
+        int N = A.size();
+        for (int k = 0; k < nwpass; k++) {
+            double v = 0;
+            for (int j = 0; j < n; j++) {
+                v += A[k][j] * x[j];
+            }
+            if (v > Upsq) {
+                Arr g = A[k];
+                pair<double, double> f = make_pair(v - Upsq, v - Lpsq);
+                return make_pair(make_pair(g, f), 0.0);
+            }
+            if (v < Lpsq) {
+                Arr g = A[k];
+                pair<double, double> f = make_pair(-v + Lpsq, -v + Upsq);
+                return make_pair(make_pair(g, f), 0.0);
+            }
+        }
+
+        // case 3, stopband constraint
+        double fmax = -INFINITY;
+        int imax = 0;
+        for (int k = nwstop; k < N; k++) {
+            double v = 0;
+            for (int j = 0; j < n; j++) {
+                v += A[k][j] * x[j];
+            }
+            if (v > Spsq) {
+                Arr g = A[k];
+                pair<double, double> f = make_pair(v - Spsq, v);
+                return make_pair(make_pair(g, f), 0.0);
+            }
+            if (v < 0) {
+                Arr g = A[k];
+                pair<double, double> f = make_pair(-v, -v + Spsq);
+                return make_pair(make_pair(g, f), 0.0);
+            }
+            if (v > fmax) {
+                fmax = v;
+                imax = k;
+            }
+        }
+
+        // case 4, nonnegative-real constraint on other frequencies
+        for (int k = nwpass; k < nwstop; k++) {
+            double v = 0;
+            for (int j = 0; j < n; j++) {
+                v += A[k][j] * x[j];
+            }
+            if (v < 0) {
+                double f = -v;
+                Arr g = A[k];
+                return make_pair(make_pair(g, f), 0.0); // single cut
+            }
+        }
+
+        this->more_alt = false;
+
+        // case 1 (unlikely)
+        if (x[0] < 0) {
+            Arr g(n, 0.0);
+            g[0] = -1.0;
+            double f = -x[0];
+            return make_pair(make_pair(g, f), 0.0);
+        }
+
+        // Begin objective function
+        Spsq = fmax;
+        pair<double, double> f = make_pair(0.0, fmax);
+        Arr g = A[imax];
+        return make_pair(make_pair(g, f), Spsq);
+    }
+};
+
+pair<LowpassOracle, double> create_lowpass_case(int N = 48) {
+    double delta0_wpass = 0.025;
+    double delta0_wstop = 0.125;
+    double delta1 = 20 * log10(1 + delta0_wpass);
+    double delta2 = 20 * log10(delta0_wstop);
+
+    int m = 15 * N;
+    vector<double> w(m);
+    for (int i = 0; i < m; i++) {
+        w[i] = i * M_PI / (m - 1);
+    }
+
+    Arr An(m, vector<double>(N));
+    for (int i = 0; i < m; i++) {
+        for (int j = 0; j < N; j++) {
+            An[i][j] = 2 * cos(w[i] * (j + 1));
+        }
+    }
+
+    Arr A(m, vector<double>(N + 1));
+    for (int i = 0; i < m; i++) {
+        A[i][0] = 1.0;
+        for (int j = 0; j < N; j++) {
+            A[i][j + 1] = An[i][j];
+        }
+    }
+
+    int nwpass = floor(0.12 * (m - 1)) + 1;
+    int nwstop = floor(0.20 * (m - 1)) + 1;
+
+    double Lp = pow(10, -delta1 / 20);
+    double Up = pow(10, delta1 / 20);
+    double Sp = pow(10, delta2 / 20);
+
+    double Lpsq = Lp * Lp;
+    double Upsq = Up * Up;
+    double Spsq = Sp * Sp;
+
+    LowpassOracle omega(A, nwpass, nwstop, Lpsq, Upsq);
+    return make_pair(omega, Spsq);
+}
+
+int main() {
+    pair<LowpassOracle, double> result = create_lowpass_case();
+    LowpassOracle omega = result.first;
+    double Spsq = result.second;
+
+    // Use the omega object and Spsq value in further calculations
+
+    return 0;
+}
+

lowpass_oracle_ai.rs

@@ -0,0 +1,129 @@
+use ndarray::Array2;
+use ndarray::Array1;
+use ndarray::Axis;
+use ndarray::s;
+use ndarray::stack;
+use ndarray::Array;
+use ndarray::ArrayView;
+use ndarray::ArrayViewMut;
+use ndarray::Array2;
+use ndarray::Array1;
+use ndarray::Axis;
+use ndarray::s;
+use ndarray::stack;
+use ndarray::Array;
+use ndarray::ArrayView;
+use ndarray::ArrayViewMut;
+
+type Arr = Array2<f64>;
+type Cut = (Arr, f64);
+
+struct LowpassOracle {
+    A: Arr,
+    nwpass: usize,
+    nwstop: usize,
+    Lpsq: f64,
+    Upsq: f64,
+    more_alt: bool,
+}
+
+impl LowpassOracle {
+    fn new(A: Arr, nwpass: usize, nwstop: usize, Lpsq: f64, Upsq: f64) -> Self {
+        LowpassOracle {
+            A,
+            nwpass,
+            nwstop,
+            Lpsq,
+            Upsq,
+            more_alt: true,
+        }
+    }
+
+    fn assess_optim(&mut self, x: ArrayView<f64, ndarray::Ix1>, Spsq: f64) -> (Cut, Option<f64>) {
+        let n = x.len();
+        self.more_alt = true;
+
+        for k in 0..self.nwpass {
+            let v = self.A.slice(s![k, ..]).dot(&x);
+            if v > self.Upsq {
+                let g = self.A.slice(s![k, ..]);
+                let f = (v - self.Upsq, v - self.Lpsq);
+                return ((g, f), None);
+            }
+
+            if v < self.Lpsq {
+                let g = -self.A.slice(s![k, ..]);
+                let f = (-v + self.Lpsq, -v + self.Upsq);
+                return ((g, f), None);
+            }
+        }
+
+        let mut fmax = f64::NEG_INFINITY;
+        let mut imax = 0;
+        for k in self.nwstop..self.A.shape()[0] {
+            let v = self.A.slice(s![k, ..]).dot(&x);
+            if v > Spsq {
+                let g = self.A.slice(s![k, ..]);
+                let f = (v - Spsq, v);
+                return ((g, f), None);
+            }
+
+            if v < 0.0 {
+                let g = -self.A.slice(s![k, ..]);
+                let f = (-v, -v + Spsq);
+                return ((g, f), None);
+            }
+
+            if v > fmax {
+                fmax = v;
+                imax = k;
+            }
+        }
+
+        for k in self.nwpass..self.nwstop {
+            let v = self.A.slice(s![k, ..]).dot(&x);
+            if v < 0.0 {
+                let f = -v;
+                let g = -self.A.slice(s![k, ..]);
+                return ((g, f), None);
+            }
+        }
+
+        self.more_alt = false;
+
+        if x[0] < 0.0 {
+            let g = Array1::zeros(n);
+            g[[0]] = -1.0;
+            let f = -x[0];
+            return ((g, f), None);
+        }
+
+        let Spsq = fmax;
+        let f = (0.0, fmax);
+        let g = self.A.slice(s![imax, ..]);
+        ((g, f), Some(Spsq))
+    }
+}
+
+fn create_lowpass_case(N: usize) -> (LowpassOracle, f64) {
+    let delta0_wpass = 0.025;
+    let delta0_wstop = 0.125;
+    let delta1 = 20.0 * (delta0_wpass).log10();
+    let delta2 = 20.0 * (delta0_wstop).log10();
+
+    let m = 15 * N;
+    let w = Array::linspace(0.0, std::f64::consts::PI, m);
+    let An = Array::from_shape_fn((m, N - 1), |(i, j)| 2.0 * (w[i] * (j + 1) as f64).cos());
+    let A = stack![Axis(1), Array::ones(m).insert_axis(Axis(1)), An];
+    let nwpass = (0.12 * (m - 1) as f64).floor() as usize + 1;
+    let nwstop = (0.20 * (m - 1) as f64).floor() as usize + 1;
+    let Lp = (10.0f64).powf(-delta1 / 20.0);
+    let Up = (10.0f64).powf(delta1 / 20.0);
+    let Sp = (10.0f64).powf(delta2 / 20.0);
+    let Lpsq = Lp * Lp;
+    let Upsq = Up * Up;
+    let Spsq = Sp * Sp;
+
+    let omega = LowpassOracle::new(A, nwpass, nwstop, Lpsq, Upsq);
+    (omega, Spsq)
+}