Skip to content
This repository has been archived by the owner on Jul 26, 2020. It is now read-only.

Commit

Permalink
svd, matrix+*double sped up, pointerToMatrix fixed
Browse files Browse the repository at this point in the history
  • Loading branch information
scottsievert committed Jul 15, 2014
1 parent eb82f4f commit e1e5f87
Show file tree
Hide file tree
Showing 10 changed files with 153 additions and 44 deletions.
10 changes: 9 additions & 1 deletion swix/swix.xcodeproj/project.pbxproj
Original file line number Diff line number Diff line change
Expand Up @@ -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 */; };
Expand Down Expand Up @@ -49,6 +51,8 @@
D206A29F1974680A009232B7 /* OpenCV.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = OpenCV.h; sourceTree = "<group>"; };
D206A2A019746886009232B7 /* opencv.mm */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.objcpp; path = opencv.mm; sourceTree = "<group>"; };
D2333745197430B40027DFD5 /* ScalarArithmetic.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; name = ScalarArithmetic.swift; path = "swix/ScalarArithmetic-1.1.1/ScalarArithmetic/ScalarArithmetic.swift"; sourceTree = "<group>"; };
D23648CA1975CE290020FB95 /* svd.m */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.objc; path = svd.m; sourceTree = "<group>"; };
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 = "<group>"; };
D24A855B196D7E86009C18AC /* matrix.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = matrix.swift; sourceTree = "<group>"; };
Expand Down Expand Up @@ -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 = "<group>";
Expand Down Expand Up @@ -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 */,
Expand Down Expand Up @@ -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 */,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -51,5 +51,21 @@
landmarkType = "7">
</BreakpointContent>
</BreakpointProxy>
<BreakpointProxy
BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint">
<BreakpointContent
shouldBeEnabled = "No"
ignoreCount = "0"
continueAfterRunningActions = "No"
filePath = "swix/swix/oneD/oneD-operators.swift"
timestampString = "427148551.443173"
startingColumnNumber = "9223372036854775807"
endingColumnNumber = "9223372036854775807"
startingLineNumber = "66"
endingLineNumber = "66"
landmarkName = "make_operator(_:_:_:)"
landmarkType = "7">
</BreakpointContent>
</BreakpointProxy>
</Breakpoints>
</Bucket>
2 changes: 2 additions & 0 deletions swix/swix/main.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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)



Expand Down
2 changes: 1 addition & 1 deletion swix/swix/swix/objc/conversion.swift
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,6 @@ func matrixToPointer(x: [Int])->UnsafePointer<Int>{
func pointerTo2DMatrix(xPC: UnsafePointer<Double>, 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
}
4 changes: 4 additions & 0 deletions swix/swix/swix/objc/math.m
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}



Expand Down
39 changes: 39 additions & 0 deletions swix/swix/swix/objc/svd.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
//
// svd.m
// swix
//
// Created by Scott Sievert on 7/15/14.
// Copyright (c) 2014 com.scott. All rights reserved.
//

#import <Foundation/Foundation.h>
#import <Accelerate/Accelerate.h>
#import <stdint.h>
#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);
}


4 changes: 3 additions & 1 deletion swix/swix/swix/objc/swix-Bridging-Header.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
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);
59 changes: 34 additions & 25 deletions swix/swix/swix/oneD/oneD-operators.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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..<lhs.n{
if operator == "<"{
Expand Down Expand Up @@ -104,28 +106,35 @@ func make_operator(lhs:matrix, operator:String, rhs:Double) -> matrix{
}
func make_operator(lhs:Double, operator:String, rhs:matrix) -> matrix{
var array = zeros(rhs.n) // lhs[i], rhs[i]
for i in 0..<rhs.n{
if operator == "<"{
if lhs < rhs[i]{ 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 == "~=="{
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..<rhs.n{
if operator == "<"{
if lhs < rhs[i]{ 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 == "~=="{
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
}
Expand Down
36 changes: 36 additions & 0 deletions swix/swix/swix/twoD/twoD-complex-math.swift
Original file line number Diff line number Diff line change
@@ -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)))
}
25 changes: 9 additions & 16 deletions swix/swix/swix/twoD/twoD-simple-math.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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..<x.n{
// if x[i] != 0{
// count += 1
// }
// }
// return count
// }
//
// assert(false, "type of norm unrecongnized")
// return -1.0
//}

//func cumsum(x: matrix) -> matrix{
// let N = x.count
// var y = zeros(N)
Expand Down

0 comments on commit e1e5f87

Please sign in to comment.