gsl-shell.git - gsl-shell

index : gsl-shell.git
gsl-shell
summary refs log tree commit diff
diff options
context:
space:
mode:
authorfrancesco <francesco.bbt@gmail.com>2011年02月28日 23:14:27 +0100
committerfrancesco <francesco.bbt@gmail.com>2011年02月28日 23:14:27 +0100
commit9af3e54b51a4095dd1202ae213a35d4b385f2d69 (patch)
treef0b569c8448c43d803424f5dabde9abf3094a020
parent8ff6e8c84d8ef10561ad6a889a6333a1a473c21f (diff)
downloadgsl-shell-9af3e54b51a4095dd1202ae213a35d4b385f2d69.tar.gz
Implemented interface to QR decomposition
Only the most important QR-related functions are implemented. Other functions should be added to provide a more complete implementation.
Diffstat
-rw-r--r--Makefile 2
-rw-r--r--examples/qr-lssolve.lua 28
-rw-r--r--lua-gsl.c 2
-rw-r--r--qr_decomp.c 171
-rw-r--r--qr_decomp.h 8
5 files changed, 210 insertions, 1 deletions
diff --git a/Makefile b/Makefile
index 2eb3e33c..788b7e8b 100644
--- a/Makefile
+++ b/Makefile
@@ -52,7 +52,7 @@ C_SRC_FILES = gs-types.c lcomplex.c matrix.c matrix_arith.c nlinfit_helper.c \
pdf.c cdf.c sf.c fmultimin.c gradcheck.c fdfmultimin.c \
multimin.c eigen-systems.c mlinear.c bspline.c interp.c \
cmatrix.c cnlinfit.c code.c fft.c lua-graph.c lu_decomp_real.c \
- lu_decomp_complex.c lu_decomp.c lua-gsl.c
+ lu_decomp_complex.c lu_decomp.c qr_decomp.c lua-gsl.c
LUA_BASE_DIRS =
LUA_BASE_FILES = igsl.lua base.lua integ.lua csv.lua
diff --git a/examples/qr-lssolve.lua b/examples/qr-lssolve.lua
new file mode 100644
index 00000000..e707f627
--- /dev/null
+++ b/examples/qr-lssolve.lua
@@ -0,0 +1,28 @@
+
+function demo1()
+ local x0, x1, n = 0, 12.5, 32
+ local a, b = 0.55, -2.4
+ local xsmp = |i| (i-1)/(n-1) * x1
+
+ local r = rng()
+ local x = new(n, 1, xsmp)
+ local y = new(n, 1, |i| a*xsmp(i) + b + rnd.gaussian(r, 0.4))
+
+ X = new(n,2, |i,j| j == 1 and a or b * xsmp(i))
+
+ qr = QR(X)
+
+ xls, res = qr:lssolve(y)
+
+ print('Linear fit coefficients: ', tr(xls))
+
+ p = plot()
+ p:addline(xyline(x, X * xls))
+ p:add(xyline(x, y), 'blue', {{'stroke'}, {'marker', size=5}})
+ p.title = 'Linear Fit'
+ p:show()
+
+ return p
+end
+
+echo 'demo1() - examples of linear regression using QR decomposition'
diff --git a/lua-gsl.c b/lua-gsl.c
index 236feef8..694d9bf3 100644
--- a/lua-gsl.c
+++ b/lua-gsl.c
@@ -47,6 +47,7 @@
#include "bspline.h"
#include "interp.h"
#include "lu_decomp.h"
+#include "qr_decomp.h"
#ifdef AGG_PLOT_ENABLED
#include "lua-graph.h"
@@ -104,6 +105,7 @@ luaopen_gsl (lua_State *L)
bspline_register (L);
interp_register (L);
lu_decomp_register (L);
+ qr_decomp_register (L);
fft_register (L);
matrix_complex_register (L);
diff --git a/qr_decomp.c b/qr_decomp.c
new file mode 100644
index 00000000..c13035f3
--- /dev/null
+++ b/qr_decomp.c
@@ -0,0 +1,171 @@
+
+#include <lua.h>
+#include <lauxlib.h>
+
+#include <gsl/gsl_linalg.h>
+
+#include "gs-types.h"
+#include "matrix.h"
+#include "lua-utils.h"
+
+typedef struct {
+ gsl_matrix *m;
+ gsl_vector *tau;
+} qr_decomp;
+
+static int qr_decomp_new (lua_State *L);
+static int qr_decomp_free (lua_State *L);
+static int qr_decomp_index (lua_State *L);
+static int qr_decomp_solve (lua_State *L);
+static int qr_decomp_lssolve (lua_State *L);
+static int qr_decomp_unpack (lua_State *L);
+
+const struct luaL_Reg qr_decomp_metatable[] = {
+ {"__gc", qr_decomp_free},
+ {"__index", qr_decomp_index},
+ {NULL, NULL}
+};
+
+static const struct luaL_Reg qr_decomp_methods[] = {
+ {"solve", qr_decomp_solve},
+ {"lssolve", qr_decomp_lssolve},
+ {"unpack", qr_decomp_unpack},
+ {NULL, NULL}
+};
+
+
+static const struct luaL_Reg qr_decomp_functions[] = {
+ {"QR", qr_decomp_new},
+ {NULL, NULL}
+};
+
+int
+qr_decomp_new (lua_State *L)
+{
+ int status;
+ gsl_matrix *a = matrix_check (L, 1);
+ size_t m = a->size1, n = a->size2;
+ qr_decomp *qr = gs_new_object (sizeof(qr_decomp), L, GS_QR_DECOMP);
+ size_t tn = (m <= n ? m : n);
+
+ qr->m = gsl_matrix_alloc (m, n);
+ gsl_matrix_memcpy (qr->m, a);
+
+ qr->tau = gsl_vector_alloc (tn);
+ if (qr->tau == NULL)
+ return luaL_error (L, "out of memory");
+
+ status = gsl_linalg_QR_decomp (qr->m, qr->tau);
+
+ if (status != GSL_SUCCESS)
+ {
+ return luaL_error (L, "error during QR decomposition: %s",
+ gsl_strerror (status));
+ }
+
+ lua_newtable (L);
+ lua_setfenv (L, -2);
+
+ return 1;
+}
+
+int
+qr_decomp_free (lua_State *L)
+{
+ qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP);
+ gsl_vector_free (qr->tau);
+ gsl_matrix_free (qr->m);
+ return 0;
+}
+
+int
+qr_decomp_unpack (lua_State *L)
+{
+ qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP);
+ size_t m = qr->m->size1, n = qr->m->size2;
+ gsl_matrix *q = matrix_push_raw (L, m, m);
+ gsl_matrix *r = matrix_push_raw (L, m, n);
+
+ gsl_linalg_QR_unpack (qr->m, qr->tau, q, r);
+
+ return 2;
+}
+
+int
+qr_decomp_solve (lua_State *L)
+{
+ qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP);
+ size_t m = qr->m->size1, n = qr->m->size2;
+ gsl_matrix *b = matrix_check (L, 2);
+ gsl_vector_view b1 = gsl_matrix_column (b, 0);
+
+ if (m != n)
+ return luaL_error (L, "\"solve\" requires a square matrix");
+
+ if (b->size1 != n || b->size2 != 1)
+ return luaL_error (L, "matrix dimensions mismatch");
+
+ {
+ int status;
+ gsl_matrix *x = matrix_push_raw (L, n, 1);
+ gsl_vector_view x1 = gsl_matrix_column (x, 0);
+ status = gsl_linalg_QR_solve (qr->m, qr->tau, &b1.vector, &x1.vector);
+ if (status != GSL_SUCCESS)
+ {
+ return luaL_error (L, "cannot solve the system: %s",
+ gsl_strerror (status));
+ }
+ }
+
+ return 1;
+}
+
+int
+qr_decomp_lssolve (lua_State *L)
+{
+ qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP);
+ size_t m = qr->m->size1, n = qr->m->size2;
+ gsl_matrix *b = matrix_check (L, 2);
+ gsl_vector_view b1 = gsl_matrix_column (b, 0);
+
+ if (m <= n)
+ return luaL_error (L, "too few rows");
+
+ if (b->size1 != m || b->size2 != 1)
+ return luaL_error (L, "matrix dimensions mismatch");
+
+ {
+ int status;
+ gsl_matrix *x = matrix_push_raw (L, n, 1);
+ gsl_matrix *r = matrix_push_raw (L, m, 1);
+ gsl_vector_view x1 = gsl_matrix_column (x, 0);
+ gsl_vector_view r1 = gsl_matrix_column (r, 0);
+
+ status = gsl_linalg_QR_lssolve (qr->m, qr->tau, &b1.vector, &x1.vector,
+ &r1.vector);
+
+ if (status != GSL_SUCCESS)
+ {
+ return luaL_error (L, "cannot solve the system: %s",
+ gsl_strerror (status));
+ }
+ }
+
+ return 2;
+}
+
+int
+qr_decomp_index (lua_State *L)
+{
+ return mlua_index_methods (L, qr_decomp_methods);
+}
+
+void
+qr_decomp_register (lua_State *L)
+{
+ luaL_newmetatable (L, GS_METATABLE(GS_QR_DECOMP));
+ luaL_register (L, NULL, qr_decomp_metatable);
+ lua_pop (L, 1);
+
+ luaL_register (L, NULL, qr_decomp_functions);
+}
diff --git a/qr_decomp.h b/qr_decomp.h
new file mode 100644
index 00000000..79ea8a84
--- /dev/null
+++ b/qr_decomp.h
@@ -0,0 +1,8 @@
+#ifndef QR_DECOMP_H
+#define QR_DECOMP_H
+
+#include <lua.h>
+
+extern void qr_decomp_register (lua_State *L);
+
+#endif
generated by cgit v1.2.3 (git 2.25.1) at 2025年10月03日 23:13:25 +0000

AltStyle によって変換されたページ (->オリジナル) /