From e1e5f87fcfc453737760b25eceff051119df64ef Mon Sep 17 00:00:00 2001 From: scottsievert Date: Tue, 15 Jul 2014 17:18:16 -0500 Subject: [PATCH] svd, matrix+*double sped up, pointerToMatrix fixed --- swix/swix.xcodeproj/project.pbxproj | 10 +++- .../xcdebugger/Breakpoints_v2.xcbkptlist | 16 +++++ swix/swix/main.swift | 2 + swix/swix/swix/objc/conversion.swift | 2 +- swix/swix/swix/objc/math.m | 4 ++ swix/swix/swix/objc/svd.m | 39 ++++++++++++ swix/swix/swix/objc/swix-Bridging-Header.h | 4 +- swix/swix/swix/oneD/oneD-operators.swift | 59 +++++++++++-------- swix/swix/swix/twoD/twoD-complex-math.swift | 36 +++++++++++ swix/swix/swix/twoD/twoD-simple-math.swift | 25 +++----- 10 files changed, 153 insertions(+), 44 deletions(-) create mode 100644 swix/swix/swix/objc/svd.m create mode 100644 swix/swix/swix/twoD/twoD-complex-math.swift diff --git a/swix/swix.xcodeproj/project.pbxproj b/swix/swix.xcodeproj/project.pbxproj index dcb19a8..8149b17 100644 --- a/swix/swix.xcodeproj/project.pbxproj +++ b/swix/swix.xcodeproj/project.pbxproj @@ -12,6 +12,8 @@ D20087B4196E186600AB26AE /* twoD-operators.swift in Sources */ = {isa = PBXBuildFile; fileRef = D20087B3196E186600AB26AE /* twoD-operators.swift */; }; D206A2A119746886009232B7 /* opencv.mm in Sources */ = {isa = PBXBuildFile; fileRef = D206A2A019746886009232B7 /* opencv.mm */; }; D2333746197430B40027DFD5 /* ScalarArithmetic.swift in Sources */ = {isa = PBXBuildFile; fileRef = D2333745197430B40027DFD5 /* ScalarArithmetic.swift */; }; + D23648CB1975CE290020FB95 /* svd.m in Sources */ = {isa = PBXBuildFile; fileRef = D23648CA1975CE290020FB95 /* svd.m */; }; + D23648CD1975CED70020FB95 /* twoD-complex-math.swift in Sources */ = {isa = PBXBuildFile; fileRef = D23648CC1975CED70020FB95 /* twoD-complex-math.swift */; }; D24A8555196D7D73009C18AC /* main.swift in Sources */ = {isa = PBXBuildFile; fileRef = D24A8554196D7D73009C18AC /* main.swift */; }; D24A855C196D7E86009C18AC /* matrix.swift in Sources */ = {isa = PBXBuildFile; fileRef = D24A855B196D7E86009C18AC /* matrix.swift */; }; D271C9BC197031C1008E577A /* math.swift in Sources */ = {isa = PBXBuildFile; fileRef = D271C9BB197031C1008E577A /* math.swift */; }; @@ -49,6 +51,8 @@ D206A29F1974680A009232B7 /* OpenCV.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = OpenCV.h; sourceTree = ""; }; D206A2A019746886009232B7 /* opencv.mm */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.objcpp; path = opencv.mm; sourceTree = ""; }; D2333745197430B40027DFD5 /* ScalarArithmetic.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; name = ScalarArithmetic.swift; path = "swix/ScalarArithmetic-1.1.1/ScalarArithmetic/ScalarArithmetic.swift"; sourceTree = ""; }; + D23648CA1975CE290020FB95 /* svd.m */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.objc; path = svd.m; sourceTree = ""; }; + D23648CC1975CED70020FB95 /* twoD-complex-math.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; name = "twoD-complex-math.swift"; path = "swix/swix/twoD/twoD-complex-math.swift"; sourceTree = SOURCE_ROOT; }; D24A8551196D7D73009C18AC /* swix */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = swix; sourceTree = BUILT_PRODUCTS_DIR; }; D24A8554196D7D73009C18AC /* main.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = main.swift; sourceTree = ""; }; D24A855B196D7E86009C18AC /* matrix.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = matrix.swift; sourceTree = ""; }; @@ -87,10 +91,11 @@ D20087B0196E14FB00AB26AE /* twoD */ = { isa = PBXGroup; children = ( + D20087AE196E14EF00AB26AE /* matrix2d.swift */, D20087B1196E180C00AB26AE /* twoD-initing.swift */, D2BAB4C3196E656E005020F6 /* twoD-simple-math.swift */, - D20087AE196E14EF00AB26AE /* matrix2d.swift */, D20087B3196E186600AB26AE /* twoD-operators.swift */, + D23648CC1975CED70020FB95 /* twoD-complex-math.swift */, ); name = twoD; sourceTree = ""; @@ -142,6 +147,7 @@ D2BBBDAB196DD693004707F3 /* math.m */, D2BBBDAC196DD693004707F3 /* swix-Bridging-Header.h */, D206A29F1974680A009232B7 /* OpenCV.h */, + D23648CA1975CE290020FB95 /* svd.m */, D206A2A019746886009232B7 /* opencv.mm */, D2BD92AD19746E33001F9EDD /* opencv2.framework */, D2FF0F0A197308440036CD3B /* indexing.m */, @@ -231,11 +237,13 @@ D2A5CAD8196D89C500749F52 /* oneD-simple-math.swift in Sources */, D2FF0F0B197308440036CD3B /* indexing.m in Sources */, D24A855C196D7E86009C18AC /* matrix.swift in Sources */, + D23648CD1975CED70020FB95 /* twoD-complex-math.swift in Sources */, D2333746197430B40027DFD5 /* ScalarArithmetic.swift in Sources */, D24A8555196D7D73009C18AC /* main.swift in Sources */, D2BAB4C4196E656E005020F6 /* twoD-simple-math.swift in Sources */, D2A5CAD4196D805A00749F52 /* oneD-initing.swift in Sources */, D20087AF196E14EF00AB26AE /* matrix2d.swift in Sources */, + D23648CB1975CE290020FB95 /* svd.m in Sources */, D20087B4196E186600AB26AE /* twoD-operators.swift in Sources */, D271C9BE1970328C008E577A /* conversion.swift in Sources */, D2FF0F1319736AEC0036CD3B /* README.md in Sources */, diff --git a/swix/swix.xcodeproj/xcuserdata/scott.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist b/swix/swix.xcodeproj/xcuserdata/scott.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist index 841da13..df05a8d 100644 --- a/swix/swix.xcodeproj/xcuserdata/scott.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist +++ b/swix/swix.xcodeproj/xcuserdata/scott.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist @@ -51,5 +51,21 @@ landmarkType = "7"> + + + + diff --git a/swix/swix/main.swift b/swix/swix/main.swift index b5bc65b..713bec6 100644 --- a/swix/swix/main.swift +++ b/swix/swix/main.swift @@ -104,6 +104,8 @@ matrix2d_indexing_matrix_test() fft_test() dot_test() +var x = array("1 2; 4 8; 3 5") +var (u, s, v) = svd(x) diff --git a/swix/swix/swix/objc/conversion.swift b/swix/swix/swix/objc/conversion.swift index 6a2b7f3..e89d8d2 100644 --- a/swix/swix/swix/objc/conversion.swift +++ b/swix/swix/swix/objc/conversion.swift @@ -19,6 +19,6 @@ func matrixToPointer(x: [Int])->UnsafePointer{ func pointerTo2DMatrix(xPC: UnsafePointer, N: CInt, M:CInt) -> matrix2d{ var x = zeros((N.int, M.int)) var xP = matrixToPointer(x.flat) - xP = xPC; + copy_objc(xPC, xP, N*M); return x } \ No newline at end of file diff --git a/swix/swix/swix/objc/math.m b/swix/swix/swix/objc/math.m index 069f304..a396443 100644 --- a/swix/swix/swix/objc/math.m +++ b/swix/swix/swix/objc/math.m @@ -95,6 +95,10 @@ void index_xa_b_objc(double* x, double* a, double* b, int N){ void copy_objc(double*x, double*y, int N){ cblas_dcopy(N, x, 1, y, 1); } +void mul_scalar_objc(double* x, double A, double* y, int N){ + double C = 0; + vDSP_vsmsaD(x, 1, &A, &C, y, 1, N); +} diff --git a/swix/swix/swix/objc/svd.m b/swix/swix/swix/objc/svd.m new file mode 100644 index 0000000..e7e9afb --- /dev/null +++ b/swix/swix/swix/objc/svd.m @@ -0,0 +1,39 @@ +// +// svd.m +// swix +// +// Created by Scott Sievert on 7/15/14. +// Copyright (c) 2014 com.scott. All rights reserved. +// + +#import +#import +#import +#import "swix-Bridging-Header.h" // for copy_objc + +void svd_objc(double * xx, int m, int n, double* sigma, double* vt, double* u){ + // adapted from http://stackoverflow.com/questions/5047503/lapack-svd-singular-value-decomposition + // Setup a buffer to hold the singular values: + int lda = m; + int numberOfSingularValues = m < n ? m : n; + + // Workspace and status variables: + double workSize; + double *work = &workSize; + int lwork = -1; + int *iwork = malloc(8*numberOfSingularValues); + int info = 0; + + // Call dgesdd_ with lwork = -1 to query optimal workspace size: + dgesdd_("A", &m, &n, xx, &lda, sigma, u, &m, vt, &n, work, &lwork, iwork, &info); + + // Optimal workspace size is returned in work[0]. + lwork = workSize; + work = malloc(lwork * sizeof *work); + + // Call dgesdd_ to do the actual computation: + dgesdd_("A", &m, &n, xx, &lda, sigma, u, &m, vt, &n, work, &lwork, iwork, &info); + free(work); +} + + diff --git a/swix/swix/swix/objc/swix-Bridging-Header.h b/swix/swix/swix/objc/swix-Bridging-Header.h index cf9b1f3..395dd32 100644 --- a/swix/swix/swix/objc/swix-Bridging-Header.h +++ b/swix/swix/swix/objc/swix-Bridging-Header.h @@ -15,4 +15,6 @@ double min_objc(double* x, int N); double max_objc(double* x, int N); void mod_objc(double * x, double mod, double * y, int N); void index_xa_b_objc(double * x, double*a, double*b, int N); -void copy_objc(double*x, double*y, int N); \ No newline at end of file +void copy_objc(double*x, double*y, int N); +void mul_scalar_objc(double* x, double A, double* y, int N); +void svd_objc(double * xx, int m, int n, double* sigma, double* vt, double* u); diff --git a/swix/swix/swix/oneD/oneD-operators.swift b/swix/swix/swix/oneD/oneD-operators.swift index d018c01..980a3a0 100644 --- a/swix/swix/swix/oneD/oneD-operators.swift +++ b/swix/swix/swix/oneD/oneD-operators.swift @@ -64,11 +64,13 @@ func make_operator(lhs:matrix, operator:String, rhs:matrix) -> matrix{ } func make_operator(lhs:matrix, operator:String, rhs:Double) -> matrix{ var array = zeros(lhs.n) - if operator == "%"{ + if operator == "%" || operator=="*" { var xP = matrixToPointer(lhs) var arrayP = matrixToPointer(array) - // mod_objc is a for-loop in C - mod_objc(xP, rhs, arrayP, lhs.n.cint); + if operator == "%" + {mod_objc(xP, rhs, arrayP, lhs.n.cint); + } else if operator == "*" + {mul_scalar_objc(xP, rhs.cdouble, arrayP, lhs.n.cint)} } else{ for i in 0.. matrix{ } func make_operator(lhs:Double, operator:String, rhs:matrix) -> matrix{ var array = zeros(rhs.n) // lhs[i], rhs[i] - for i in 0.."{ - if lhs > rhs[i]{ array[i] = 1 } - } else if operator == "<"{ - if lhs < rhs[i]{ array[i] = 1 } - } else if operator == "~=="{ - if abs(lhs - rhs[i])<1e-9 { array[i] = 1 } - } else if operator == "<="{ - if lhs <= rhs[i]{ array[i] = 1 } - } else if operator == ">="{ - if lhs >= rhs[i]{ array[i] = 1 } - } else if operator == "+"{ - array[i] = lhs + rhs[i] - } else if operator == "-"{ - array[i] = lhs - rhs[i] - } else if operator == "*"{ - array[i] = lhs * rhs[i] - } else if operator == "/"{ - array[i] = lhs / rhs[i] - } else { assert(false, "Operator not reconginzed!") } + if operator=="*"{ + var xP = matrixToPointer(rhs) + var arrayP = matrixToPointer(array) + if operator == "*" + {mul_scalar_objc(xP, lhs.cdouble, arrayP, rhs.n.cint)} + } else{ + for i in 0.."{ + if lhs > rhs[i]{ array[i] = 1 } + } else if operator == "<"{ + if lhs < rhs[i]{ array[i] = 1 } + } else if operator == "~=="{ + if abs(lhs - rhs[i])<1e-9 { array[i] = 1 } + } else if operator == "<="{ + if lhs <= rhs[i]{ array[i] = 1 } + } else if operator == ">="{ + if lhs >= rhs[i]{ array[i] = 1 } + } else if operator == "+"{ + array[i] = lhs + rhs[i] + } else if operator == "-"{ + array[i] = lhs - rhs[i] + } else if operator == "*"{ + array[i] = lhs * rhs[i] + } else if operator == "/"{ + array[i] = lhs / rhs[i] + } else { assert(false, "Operator not reconginzed!") } + } } return array } diff --git a/swix/swix/swix/twoD/twoD-complex-math.swift b/swix/swix/swix/twoD/twoD-complex-math.swift new file mode 100644 index 0000000..f223174 --- /dev/null +++ b/swix/swix/swix/twoD/twoD-complex-math.swift @@ -0,0 +1,36 @@ +// +// twoD-complex-math.swift +// swix +// +// Created by Scott Sievert on 7/15/14. +// Copyright (c) 2014 com.scott. All rights reserved. +// + +import Foundation + +func svd(x: matrix2d) -> (matrix2d, matrix, matrix2d){ + // 2014-7-15: almost have this working. to change by tomorrow. + var (m, n) = x.shape + var nS = m < n ? m : n + var sigma = zeros(nS) + var vt = zeros((n,n)) + var u = zeros((m,m)) + + var xx = zeros_like(x) + xx.flat = x.flat + xx = transpose(xx) + var xP = matrixToPointer(xx.flat) + var sP = matrixToPointer(sigma) + var vP = matrixToPointer(vt.flat) + var uP = matrixToPointer(u.flat) + + println(x) + svd_objc(xP, m.cint, n.cint, sP, vP, uP); + if m > n {vt = transpose(vt)} + + println(sigma) + println(vt) + println(u) + + return (zeros((2,2)), zeros(2), zeros((2,2))) +} \ No newline at end of file diff --git a/swix/swix/swix/twoD/twoD-simple-math.swift b/swix/swix/swix/twoD/twoD-simple-math.swift index af83dd9..7987ee8 100644 --- a/swix/swix/swix/twoD/twoD-simple-math.swift +++ b/swix/swix/swix/twoD/twoD-simple-math.swift @@ -69,6 +69,14 @@ func min(x: matrix2d, absValue:Bool=false)-> Double{ func max(x: matrix2d, absValue:Bool=false)-> Double{ return max(x.flat, absValue:absValue) } +func norm(x: matrix2d, type:String="l2") -> Double{ + if type=="l0"{ return norm(x.flat, type:"l0")} + if type=="l1"{ return norm(x.flat, type:"l1")} + if type=="l2"{ return norm(x.flat, type:"l2")} + + assert(false, "type of norm unrecongnized") + return -1.0 +} //func pow(x: matrix, power: Double) -> matrix{ // var y = zeros(x.count) @@ -101,22 +109,7 @@ func max(x: matrix2d, absValue:Bool=false)-> Double{ // var z = x - y // return sum(pow(z, 2) / x.count.double) //} -//func norm(x: matrix, type:String="l2") -> Double{ -// if type=="l2"{ return sqrt(sum(pow(x, 2)))} -// if type=="l1"{ return sum(abs(x))} -// if type=="l0"{ -// var count = 0.0 -// for i in 0.. matrix{ // let N = x.count // var y = zeros(N)