Wie wird der kubische Spline-Interpolationsalgorithmus beschleunigt?

Question

Ich merke, dass wenn ich "la_solve" verwende, bekomme ich den Malloc-Fehler aufgrund einer Begrenzung in der Anzahl der Proben, die "la_solve" verarbeiten kann.

Ich bemerkte, dass das Programm gut läuft, wenn ich etwa 900 Proben verwende.

Daher verändert ich meinen Code auf gesamten Daten in aufeinanderfolgenden Blöcken von 901 Proben zu suchen, wenn Proben> 901.

Allerdings habe ich ziemlich naiv Techniken und bin sicher, dass dies effizienter und robuster gemacht werden können: (1) Ich brauche eine Bedingung, dass, wenn die Proben < 900 sind, nicht klumpen; (2) passe auf Wiederholungen in x- und y-Werten auf, wenn ich Chunked-Daten angehängt habe; und (3) den Algorithmus effizienter, schneller und robuster machen.

Ich brauche vor allem Hilfe bei der „// Aufruf der Re-Sampler-Funktion“ am Ende des Codes

Dank !!!

import UIKit 
import Accelerate.vecLib.LinearAlgebra 

public class Resample { 

public class func resample(xi: [Double], xii: [Double], a: [Double]) ->[Double] 
{ 
    // xi - original time stamps (original x) 
    // xii - desired time stamps (desired x) 
    // a - orignal y values 
    // ->[Double] - interpolated y values 

    let ni = xii.count // Desired x count 

    let n = xi.count // Actual x count 

    var h: [Double] = Array(repeating:0, count:n-1) 

    for j in 0..<n-1 { 
     h[j] = (xi[j+1] - xi[j]) 
    } 

    var A: [Double] = Array(repeating:0, count:(n)*(n)) 

    A[0] = 1.0 

    A[(n*n)-1] = 1.0 

    for i in 1..<(n-1) { 
     A[(n+1)*i-1] = h[i-1] 
     A[(n+1)*i] = 2*(h[i-1] + h[i]) 
     A[(n+1)*i+1] = h[i] 
    } 

    var b: [Double] = Array(repeating:0, count:n) 

    for i in 1..<n-1 { 
     b[i] = (3/h[i])*(a[i+1]-a[i]) - (3/h[i-1])*(a[i]-a[i-1]) 
    } 

    let matA = la_matrix_from_double_buffer(A, la_count_t(n), la_count_t(n), la_count_t(n), la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES)) 

    let vecB = la_matrix_from_double_buffer(b, la_count_t(n), 1, 1, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES)) 

    let vecCj = la_solve(matA, vecB) 

    var cj: [Double] = Array(repeating: 0.0, count: n) 

    let status = la_matrix_to_double_buffer(&cj, 1, vecCj) 

    if status == la_status_t(LA_SUCCESS) { 
     //print(cj.count) 
    } 
    else { 
     print("Failure: \(status)") 
    } 

    var bj: [Double] = Array(repeating:0, count:n-1) 

    for i in 0..<n-1 { 
     bj[i] = (1/h[i])*(a[i+1]-a[i]) - (1/3*h[i])*(2*cj[i]+cj[i+1]) 
    } 

    var dj: [Double] = Array(repeating:0, count:n-1) 

    for i in 0..<n-1 { 
     dj[i] = (1/(3*h[i])) * (cj[i+1]-cj[i]) 
    } 

    var P: [Double] = Array(repeating: 0.0, count: (n-1)*4) 

    for i in 0..<n-1 { 
     P[(i*4)] = dj[i] 
     P[(i*4)+1] = cj[i] 
     P[(i*4)+2] = bj[i] 
     P[(i*4)+3] = a[i] 
    } 

    var ai: [Double] = Array(repeating: 0.0, count: ni) 


    var jl = Double(1) 

    for i in 0..<ni { 

     let inter = Double(xii[i]) 

     while (inter > xi[Int(jl)] && (Int(jl) < n)) { 
      jl += 1 

     } 

     let ind = 4*Int(jl)-4 

     var pp = Array(P[ind...ind+3]) 

     let xx = Double(xii[i]-xi[Int(jl)-1]) 

     ai[i] = pp[0]*pow(xx, 3) + pp[1]*pow(xx, 2) + pp[2]*xx + pp[3] 
    } 
    return(ai) 
} 
} 



// Calling the Re-sampler FUNCTION 

let x = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig X (Given, say) 

let y = Array(stride(from: 0, through: 16649, by: 1.0)) // Orig Y (Given, say) 

let f = 2.0 // Hz (Desired Sampling Freq., say) 

let g = Array(stride(from: 0, through: x.count-1, by: 900)) 
// Helper array for chunking (901 samples each time, say) 

var xxx = [Double]() // Results for appended chunks X 

var yyy = [Double]() // Results for appended chunks Y 

for i in 0..<g.count-1 { // loop through 

let xc = Array(x[Int(g[i])...Int(g[i+1])]) // x chunk (count = 901) 

let yc = Array(y[Int(g[i])...Int(g[i+1])]) // y chunk (count = 901) 

let xci = Array(stride(from: xc[0], through: xc[xc.count-1], by: 1.0/f)) 
// Desired time stamps for chunk of 901 samples 

let yci = Resample.resample(xi: xc, xii: xci, a: yc) // Interpolate chunk 

xxx += xci // Append interpolation X 

yyy += yci // Append interpolation X 

} 

if(Int(g[g.count-1])<x.count){ 
// If helper function last index < original no. of samples 

let glb = (Int(g[g.count-1])) // last chunk begin index 

let gle = (x.count) // last chunk end index 

let xce = Array(x[glb...gle-1]) // last x chunk 

let yce = Array(y[glb...gle-1]) // last y chunk 

let xcei = Array(stride(from: xce[0], through: xce[xce.count-1], by: 1.0/f)) 
// Desired time stamps for last chunk of samples 

let ycei = Resample.resample(xi: xce, xii: xcei, a: yce) // Interpolate last chunk 

xxx += xcei // Append interpolation X for last chunk 

yyy += ycei // Append interpolation X for last chunk 
} 

print(xxx) // would like to identify repeated x values and remove them 
print(yyy) // remove y values corresponsding to repated x values (as found above) 

// Calling the Re-sampler FUNCTION 

Answer 1

Beachten Sie, dass für binäre Suche alleind mit der Suche nach, während es Wert ist erhöht monoton. Sie können also einfach prüfen, ob der nächste x-Wert im selben Intervall bleibt oder zum nächsten springt. Pseudo-Code:

jl = 1 //before cycle 

//inside cycle 
    while (current_x > xi[jl]) and (jl < n) 
    jl++ 
    ind = 4 * jl - 4 

Ist la_solve allgemeine Methode zur linearen eq.system Lösung? Wenn ja, sollten Sie den dedizierten Ansatz für tridiagonale Systeme (kubische Komplexität vs. lineare Komplexität) verwenden. Classic text. Python implementation. Wiki. Es ist wichtig für lange Eingabe-Array.