From 2166ed1a01a2733a0b453c8f2db770387d8d3de9 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Tue, 6 Jun 2023 16:59:34 +0200 Subject: [PATCH 01/28] added empty files --- scripts/builtin/hmm.dml | 0 scripts/builtin/hmmPredict.dml | 0 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 scripts/builtin/hmm.dml create mode 100644 scripts/builtin/hmmPredict.dml diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml new file mode 100644 index 00000000000..e69de29bb2d diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml new file mode 100644 index 00000000000..e69de29bb2d From 5fc6aa6cb6b2db1f143358c297149a075bf1c855 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Thu, 8 Jun 2023 23:28:44 +0200 Subject: [PATCH 02/28] just func defs and comments --- scripts/builtin/hmm.dml | 36 ++++++++++++++++++++++++++++++++++ scripts/builtin/hmmPredict.dml | 36 ++++++++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index e69de29bb2d..8c164e1fd1a 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -0,0 +1,36 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# This script implements the hidden markov model method +# INPUT: +# -------------------------------------------------------------------------------------------- +# X Set of outputs +# OUTPUT: +# -------------------------------------------------------------------------------------------- +# P Probability of the set of outputs +# -------------------------------------------------------------------------------------------- + +hmm = function(Matrix[Double] X) + return (scalar p) +{ + p = 0.0 + return p +} \ No newline at end of file diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index e69de29bb2d..5a39ce54e21 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -0,0 +1,36 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# This script implements the hidden markov model method +# INPUT: +# -------------------------------------------------------------------------------------------- +# X Set of outputs +# OUTPUT: +# -------------------------------------------------------------------------------------------- +# P Probability of the set of outputs +# -------------------------------------------------------------------------------------------- + +hmmPredict = function(Matrix[Double] X, k) + return (Matrix[double] O) +{ + O = matrix(0.0, rows=1, cols=k) +} + From e2bc7867044616fb8e3cc462d9a1d182affe488f Mon Sep 17 00:00:00 2001 From: Sandesh Date: Fri, 9 Jun 2023 16:00:59 +0200 Subject: [PATCH 03/28] changed builtin order and working functions --- scripts/builtin/hmm.dml | 4 +--- scripts/builtin/hmmPredict.dml | 5 ++--- src/main/java/org/apache/sysds/common/Builtins.java | 2 ++ 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 8c164e1fd1a..b172daac708 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -28,9 +28,7 @@ # P Probability of the set of outputs # -------------------------------------------------------------------------------------------- -hmm = function(Matrix[Double] X) - return (scalar p) +m_hmm = function(Matrix[Double] X) return (Double p) { p = 0.0 - return p } \ No newline at end of file diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 5a39ce54e21..c83d408ccbd 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -28,9 +28,8 @@ # P Probability of the set of outputs # -------------------------------------------------------------------------------------------- -hmmPredict = function(Matrix[Double] X, k) - return (Matrix[double] O) +m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] p) { - O = matrix(0.0, rows=1, cols=k) + p = matrix(0.0, rows=1, cols=k) } diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java index efb260360e7..821c2f6bbb1 100644 --- a/src/main/java/org/apache/sysds/common/Builtins.java +++ b/src/main/java/org/apache/sysds/common/Builtins.java @@ -151,6 +151,8 @@ public enum Builtins { GNMF("gnmf", true), GRID_SEARCH("gridSearch", true), TOPK_CLEANING("topk_cleaning", true), + HMM("hmm", true), + HMMPREDICT("hmmPredict", true), HOSPITAL_RESIDENCY_MATCH("hospitalResidencyMatch", true), HYPERBAND("hyperband", true), IFELSE("ifelse", false), From 83900b1106d7f98c4015f8cab1275f95652799f9 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Fri, 16 Jun 2023 12:33:05 +0200 Subject: [PATCH 04/28] hmmPredictTest start --- .../builtin/part1/BuiltinHMMPredictTest.java | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java new file mode 100644 index 00000000000..a0bb5c6598f --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java @@ -0,0 +1,91 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part1; + +import org.apache.sysds.api.DMLScript; +import org.apache.sysds.common.Types; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.hops.OptimizerUtils; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; +import org.junit.Test; + +@net.jcip.annotations.NotThreadSafe +public class BuiltinHMMPredictTest extends AutomatedTestBase +{ + private final static String TEST_NAME = "hmmPredict"; + private final static String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; + private final static int rows = 1320; + private final static int cols = 32; + + @Override + public void setUp() { + TestUtils.clearAssertionInformation(); + addTestConfiguration(TEST_NAME,new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"C"})); + } + + @Test + public void testHMMPredictCP() { + runHMMPredictTest(ExecType.CP); + } + + @Test + public void testHMMPredictSingleCP() { + runHMMPredictTest(ExecType.SPARK); + } + + private void runHMMPredictTest(ExecType instType) + { + Types.ExecMode platformOld = setExecMode(instType); + + boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; + boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; + + try + { + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + + String HOME = SCRIPT_DIR + TEST_DIR; + + fullDMLScriptName = HOME + TEST_NAME + ".dml"; + programArgs = new String[]{ + "-nvargs", "X=" + input("X"), "k=" + inpu("k"), "Y=" + output("Y")}; + + + """ + //generate actual datasets + double[][] X = getRandomMatrix(rows, cols); + Integer k = + writeInputMatrixWithMTD("X", X, true); + """ + + runTest(true, false, null, -1); + } + finally { + rtplatform = platformOld; + DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; + OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; + OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; + OptimizerUtils.ALLOW_OPERATOR_FUSION = true; + } + } +} From 31d8886d8636c6b7b408beda24431a33cfc0a91c Mon Sep 17 00:00:00 2001 From: Sandesh Date: Fri, 16 Jun 2023 12:37:18 +0200 Subject: [PATCH 05/28] added hmmPredict test dml file --- .../functions/builtin/hmmPredictTest.dml | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 src/test/scripts/functions/builtin/hmmPredictTest.dml diff --git a/src/test/scripts/functions/builtin/hmmPredictTest.dml b/src/test/scripts/functions/builtin/hmmPredictTest.dml new file mode 100644 index 00000000000..271be550317 --- /dev/null +++ b/src/test/scripts/functions/builtin/hmmPredictTest.dml @@ -0,0 +1,25 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +X = read($1); +k = read($2); +Y = hmmPredict(X, k); +write(Y); From f6afde2384fe4fb4a77d84a8c8b4e9ae03b346b0 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Fri, 16 Jun 2023 15:14:34 +0200 Subject: [PATCH 06/28] fixed function name --- .../builtin/part1/BuiltinHMMPredictTest.java | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java index a0bb5c6598f..ca0f2f707aa 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java @@ -49,7 +49,7 @@ public void testHMMPredictCP() { } @Test - public void testHMMPredictSingleCP() { + public void testHMMPredictSpark() { runHMMPredictTest(ExecType.SPARK); } @@ -70,15 +70,17 @@ private void runHMMPredictTest(ExecType instType) programArgs = new String[]{ "-nvargs", "X=" + input("X"), "k=" + inpu("k"), "Y=" + output("Y")}; - - """ + //generate actual datasets double[][] X = getRandomMatrix(rows, cols); Integer k = writeInputMatrixWithMTD("X", X, true); - """ + setOutputBuffering(true); + + String stdout = runTest(null).toString(); - runTest(true, false, null, -1); + + //runTest(true, false, null, -1); } finally { rtplatform = platformOld; From ba030ad3695f42b5e2249c542392852280a2806f Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 18 Jun 2023 15:36:50 +0200 Subject: [PATCH 07/28] added hmm tests --- .../builtin/part1/BuiltinHmmTest.java | 98 +++++++++++++++++++ .../scripts/functions/builtin/hmmTest.dml | 0 2 files changed, 98 insertions(+) create mode 100644 src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java create mode 100644 src/test/scripts/functions/builtin/hmmTest.dml diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java new file mode 100644 index 00000000000..f71aed8e85c --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java @@ -0,0 +1,98 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part1; + +import java.util.Random; +import org.apache.sysds.api.DMLScript; +import org.apache.sysds.common.Types; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.hops.OptimizerUtils; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; +import org.junit.Test; + +@net.jcip.annotations.NotThreadSafe +public class BuiltinHMMTest extends AutomatedTestBase +{ + private final static String TEST_NAME = "hmm"; + private final static String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; + private final static int rows = 1320; + private final static int cols = 32; + + @Override + public void setUp() { + TestUtils.clearAssertionInformation(); + addTestConfiguration(TEST_NAME,new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"C"})); + } + + @Test + public void testHMMPredictCP() { + runHMMPredictTest(ExecType.CP); + } + + @Test + public void testHMMPredictSpark() { + runHMMPredictTest(ExecType.SPARK); + } + + private void runHMMPredictTest(ExecType instType) + { + Types.ExecMode platformOld = setExecMode(instType); + + boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; + boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; + + try + { + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + + String HOME = SCRIPT_DIR + TEST_DIR; + + fullDMLScriptName = HOME + TEST_NAME + ".dml"; + programArgs = new String[]{ + "-nvargs", "X=" + input("X"), "Y=" + output("Y")}; + + + //generate actual datasets + double[][] X = getRandomMatrix(rows, cols); + + writeInputMatrixWithMTD("X", X, true); + + setOutputBuffering(true); + + //TODO: Compare outputs? What is the output of the + // function supposed to be compared against? How? + + String stdout = runTest(null).toString(); + + + //runTest(true, false, null, -1); + } + finally { + rtplatform = platformOld; + DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; + OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; + OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; + OptimizerUtils.ALLOW_OPERATOR_FUSION = true; + } + } +} diff --git a/src/test/scripts/functions/builtin/hmmTest.dml b/src/test/scripts/functions/builtin/hmmTest.dml new file mode 100644 index 00000000000..e69de29bb2d From dadb65e21dd405d743cea518bd91611de197784a Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 18 Jun 2023 15:38:18 +0200 Subject: [PATCH 08/28] Start of Baum-welch algo --- scripts/builtin/hmm.dml | 59 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 56 insertions(+), 3 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index b172daac708..323793c45b7 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -28,7 +28,60 @@ # P Probability of the set of outputs # -------------------------------------------------------------------------------------------- -m_hmm = function(Matrix[Double] X) return (Double p) +m_hmm = function(Matrix[Double] X, String model="N") return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) { - p = 0.0 -} \ No newline at end of file + """ + if(model != "N" & model != "S") + # N is for number model and S for string model + stop("model not supported, should be in S or N"); + + if (model == "S") + """ + + # change later + n_states = 5 + + [P, A, B] = fit(X, n_states) +} + +fit = function(Matrix[Double] X, Integer n_states) return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) +{ + # probability of the state at timestep 1 + init_prob = matrix(1/n_states, rows=1, cols=n_states) + + # probabilities of transition from one one state to another + A = matrix(1/n_states, rows=n_states, cols=n_states) + + + + # probabilities of a state producing an output + B = matrix(1/ncol(X), rows=n_states, cols=ncol(X)) + + for (i in 1:20) { + if (i == 1) + init_prob = matrix(1/n_states, rows=1, cols=n_states) + A = matrix(1/n_states, rows=n_states, cols=n_states) + B = matrix(1/ncol(X), rows=n_states, cols=ncol(X)) + + alpha = forward(init_prob, B, X) + beta = backward(B, A, X) + + initial_prob, A, B = calc_init_prob(alpha, beta) + } +} + +forward = function(Matrix[Double] init_prob, Matrix[Double] B, Matrix[Double] X) return (Matrix[Double] alpha) +{ + + +} + + +backward = function(Matrix[Double] B, Matrix[Double] A, Matrix[Double] X) return (Matrix[Double] alpha) +{ + +} + +calc_init_prob = function() return (Matrix[Double] init_prob, Matrix[Double] A, Matrix[Double] B) +{ +} \ No newline at end of file From f1fff1ba6dc90840169afef0bdbbca29a3d33b2e Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 18 Jun 2023 15:41:10 +0200 Subject: [PATCH 09/28] more complete Test --- .../functions/builtin/part1/BuiltinHMMPredictTest.java | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java index ca0f2f707aa..d44981f6ed2 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java @@ -19,6 +19,7 @@ package org.apache.sysds.test.functions.builtin.part1; +import java.util.Random; import org.apache.sysds.api.DMLScript; import org.apache.sysds.common.Types; import org.apache.sysds.common.Types.ExecType; @@ -68,13 +69,18 @@ private void runHMMPredictTest(ExecType instType) fullDMLScriptName = HOME + TEST_NAME + ".dml"; programArgs = new String[]{ - "-nvargs", "X=" + input("X"), "k=" + inpu("k"), "Y=" + output("Y")}; + "-nvargs", "X=" + input("X"), "k=" + input("k"), "Y=" + output("Y")}; //generate actual datasets double[][] X = getRandomMatrix(rows, cols); - Integer k = + Integer k = random.nextInt(range) + 1; + writeInputMatrixWithMTD("X", X, true); + + //TODO: Compare outputs? What is the output of the + // function supposed to be compared against? How? + setOutputBuffering(true); String stdout = runTest(null).toString(); From 1b77d3aa70ba734f0f463124f182c2101230088d Mon Sep 17 00:00:00 2001 From: Sandesh Date: Tue, 20 Jun 2023 17:42:18 +0200 Subject: [PATCH 10/28] Removed compilation errors in tests --- .../builtin/part1/BuiltinHMMPredictTest.java | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java index d44981f6ed2..4b39289e110 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java @@ -71,10 +71,19 @@ private void runHMMPredictTest(ExecType instType) programArgs = new String[]{ "-nvargs", "X=" + input("X"), "k=" + input("k"), "Y=" + output("Y")}; - + //generate actual datasets - double[][] X = getRandomMatrix(rows, cols); - Integer k = random.nextInt(range) + 1; + // Random random = new Random(); + // Integer k = random.nextInt(10) + 1; + + //generate actual datasets + //Output generated from https://en.wikipedia.org/wiki/Hidden_Markov_model#Weather_guessing_game + // C(Clean) -> 0, W(Walk) -> 1, S(Shop) -> 2 + + double[][] X = {{0, 0, 1, + 2, 1, 1, + 1, 0, 0, 1} + }; writeInputMatrixWithMTD("X", X, true); @@ -86,7 +95,7 @@ private void runHMMPredictTest(ExecType instType) String stdout = runTest(null).toString(); - //runTest(true, false, null, -1); + runTest(true, false, null, -1); } finally { rtplatform = platformOld; From 5c01b010ef51d51a5e188622c8ffd21bb641a309 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 25 Jun 2023 19:33:12 +0200 Subject: [PATCH 11/28] spelling change and maybe others --- .../builtin/part1/BuiltinHMMTest.java | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java new file mode 100644 index 00000000000..8a72ee5801b --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java @@ -0,0 +1,102 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part1; + +import java.util.Random; +import org.apache.sysds.api.DMLScript; +import org.apache.sysds.common.Types; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.hops.OptimizerUtils; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; +import org.junit.Test; + +@net.jcip.annotations.NotThreadSafe +public class BuiltinHMMTest extends AutomatedTestBase +{ + private final static String TEST_NAME = "hmm"; + private final static String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; + private final static int rows = 1320; + private final static int cols = 32; + + @Override + public void setUp() { + TestUtils.clearAssertionInformation(); + addTestConfiguration(TEST_NAME,new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"C"})); + } + + @Test + public void testHMMPredictCP() { + runHMMPredictTest(ExecType.CP); + } + + @Test + public void testHMMPredictSpark() { + runHMMPredictTest(ExecType.SPARK); + } + + private void runHMMPredictTest(ExecType instType) + { + Types.ExecMode platformOld = setExecMode(instType); + + boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; + boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; + + try + { + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + + String HOME = SCRIPT_DIR + TEST_DIR; + + fullDMLScriptName = HOME + TEST_NAME + ".dml"; + programArgs = new String[]{ + "-nvargs", "X=" + input("X"), "Y=" + output("Y")}; + + //generate actual datasets + //Output generated from https://en.wikipedia.org/wiki/Hidden_Markov_model#Weather_guessing_game + // C(Clean) -> 0, W(Walk) -> 1, S(Shop) -> 2 + + double[][] X = {{0, 0, 1, + 2, 1, 1, + 1, 0, 0, 1} + }; + writeInputMatrixWithMTD("X", X, true); + + setOutputBuffering(true); + + //TODO: Compare outputs? What is the output of the + // function supposed to be compared against? How? + + String stdout = runTest(null).toString(); + + + runTest(true, false, null, -1); + } + finally { + rtplatform = platformOld; + DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; + OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; + OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; + OptimizerUtils.ALLOW_OPERATOR_FUSION = true; + } + } +} \ No newline at end of file From e5dc3d21561a37bbaa79584cd752cd1ce29c51b4 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 25 Jun 2023 19:38:33 +0200 Subject: [PATCH 12/28] changed file name therefore deleted --- .../builtin/part1/BuiltinHmmTest.java | 98 ------------------- 1 file changed, 98 deletions(-) delete mode 100644 src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java deleted file mode 100644 index f71aed8e85c..00000000000 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHmmTest.java +++ /dev/null @@ -1,98 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one - * or more contributor license agreements. See the NOTICE file - * distributed with this work for additional information - * regarding copyright ownership. The ASF licenses this file - * to you under the Apache License, Version 2.0 (the - * "License"); you may not use this file except in compliance - * with the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, - * software distributed under the License is distributed on an - * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY - * KIND, either express or implied. See the License for the - * specific language governing permissions and limitations - * under the License. - */ - -package org.apache.sysds.test.functions.builtin.part1; - -import java.util.Random; -import org.apache.sysds.api.DMLScript; -import org.apache.sysds.common.Types; -import org.apache.sysds.common.Types.ExecType; -import org.apache.sysds.hops.OptimizerUtils; -import org.apache.sysds.test.AutomatedTestBase; -import org.apache.sysds.test.TestConfiguration; -import org.apache.sysds.test.TestUtils; -import org.junit.Test; - -@net.jcip.annotations.NotThreadSafe -public class BuiltinHMMTest extends AutomatedTestBase -{ - private final static String TEST_NAME = "hmm"; - private final static String TEST_DIR = "functions/builtin/"; - private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; - private final static int rows = 1320; - private final static int cols = 32; - - @Override - public void setUp() { - TestUtils.clearAssertionInformation(); - addTestConfiguration(TEST_NAME,new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"C"})); - } - - @Test - public void testHMMPredictCP() { - runHMMPredictTest(ExecType.CP); - } - - @Test - public void testHMMPredictSpark() { - runHMMPredictTest(ExecType.SPARK); - } - - private void runHMMPredictTest(ExecType instType) - { - Types.ExecMode platformOld = setExecMode(instType); - - boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; - boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; - - try - { - loadTestConfiguration(getTestConfiguration(TEST_NAME)); - - String HOME = SCRIPT_DIR + TEST_DIR; - - fullDMLScriptName = HOME + TEST_NAME + ".dml"; - programArgs = new String[]{ - "-nvargs", "X=" + input("X"), "Y=" + output("Y")}; - - - //generate actual datasets - double[][] X = getRandomMatrix(rows, cols); - - writeInputMatrixWithMTD("X", X, true); - - setOutputBuffering(true); - - //TODO: Compare outputs? What is the output of the - // function supposed to be compared against? How? - - String stdout = runTest(null).toString(); - - - //runTest(true, false, null, -1); - } - finally { - rtplatform = platformOld; - DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; - OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; - OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; - OptimizerUtils.ALLOW_OPERATOR_FUSION = true; - } - } -} From 948b92b8a11a46fc8217ae42f347c2fcf8f8d9b4 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sun, 25 Jun 2023 19:41:05 +0200 Subject: [PATCH 13/28] put in proper return function names --- scripts/builtin/hmm.dml | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 323793c45b7..b63fb61f436 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -52,8 +52,6 @@ fit = function(Matrix[Double] X, Integer n_states) return (Matrix[Double] P, Mat # probabilities of transition from one one state to another A = matrix(1/n_states, rows=n_states, cols=n_states) - - # probabilities of a state producing an output B = matrix(1/ncol(X), rows=n_states, cols=ncol(X)) @@ -72,16 +70,16 @@ fit = function(Matrix[Double] X, Integer n_states) return (Matrix[Double] P, Mat forward = function(Matrix[Double] init_prob, Matrix[Double] B, Matrix[Double] X) return (Matrix[Double] alpha) { - - + alpha = rand(rows=nrow(A), cols = ncol(A)) } -backward = function(Matrix[Double] B, Matrix[Double] A, Matrix[Double] X) return (Matrix[Double] alpha) +backward = function(Matrix[Double] B, Matrix[Double] A, Matrix[Double] X) return (Matrix[Double] beta) { - + beta = rand(rows=nrow(A), cols=ncol(A)) } -calc_init_prob = function() return (Matrix[Double] init_prob, Matrix[Double] A, Matrix[Double] B) +calc_init_prob = function(Matrix[Double] alpha, Matrix[Double] beta) return (Matrix[Double] pie) { + pie = rand(rows=nrow(alpha), cols=ncol(beta)) } \ No newline at end of file From d3e6ccbe04c3936b9efc17fcfcdbed63febde065 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 11:52:55 +0200 Subject: [PATCH 14/28] working hmmpredict implementation --- scripts/builtin/hmmPredict.dml | 55 ++++++++++++++++++++++++++++++++-- 1 file changed, 53 insertions(+), 2 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index c83d408ccbd..68cea88045c 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -28,8 +28,59 @@ # P Probability of the set of outputs # -------------------------------------------------------------------------------------------- -m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] p) +choose_value = function (matrix[double] prob_dist) return (double val) { - p = matrix(0.0, rows=1, cols=k) + #useful for selecting output or states given probability distribution + + #choosing using cumulative distribution + cum_prob = cumsum(prob_dist) + + R = rand(rows=1, cols=1, min=0, max=1) + + #adding 1 to colsum because indexing starts from 1 + val = as.scalar(colSums(cum_prob <= R) + 1) +} + +get_all_states = function (matrix[double] A, matrix[double] ip, integer k) return (matrix[double] states) +{ + init_state = choose_value(ip) + states = matrix(0, rows=k, cols=1) + states[1, 1] = init_state + s = init_state + + for (i in 2:k) { + next_state = choose_value(t(A[s, ])) + states[i, 1] = next_state + s = next_state + } +} + +get_outputs = function (matrix[double] B, matrix[double] states) return (matrix[double] outputs) +{ + outputs = matrix(0, rows=nrow(states), cols=1) + for (i in 1:nrow(states)) { + s = as.scalar(states[i, 1]) + outputs[i, 1] = choose_value(t(B[s, ])) + } +} + +m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] outputs) +{ + /* + A, B, ip = hmm(X) + + The above would be the function call to obtain model parameters. + Instead, the following has dummy model parameters to for the functions to work + */ + + n_states = 3 + n_outputs = 3 + + A = matrix(1/n_states, rows=n_states, cols=n_states) + B = matrix(1/n_outputs, rows=n_states, cols=n_outputs) + ip = matrix(1/n_states, rows=n_states, cols=1) + + states = get_all_states(A, ip, k) + outputs = get_outputs(B, states) } From 9d9c42dc144b203beb25a7d47c08de8b028d152a Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 12:01:00 +0200 Subject: [PATCH 15/28] working hmmpredict implementation --- scripts/builtin/hmmPredict.dml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 68cea88045c..719b1f50dd9 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -58,7 +58,7 @@ get_all_states = function (matrix[double] A, matrix[double] ip, integer k) retur get_outputs = function (matrix[double] B, matrix[double] states) return (matrix[double] outputs) { outputs = matrix(0, rows=nrow(states), cols=1) - for (i in 1:nrow(states)) { + parfor (i in 1:nrow(states)) { s = as.scalar(states[i, 1]) outputs[i, 1] = choose_value(t(B[s, ])) } @@ -82,5 +82,5 @@ m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] outp states = get_all_states(A, ip, k) outputs = get_outputs(B, states) -} +} From 44dc75ced7de91828acc426b9f2205a66c9866f2 Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 13:10:47 +0200 Subject: [PATCH 16/28] added some comments --- scripts/builtin/hmmPredict.dml | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 719b1f50dd9..0b48189009a 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -22,22 +22,23 @@ # This script implements the hidden markov model method # INPUT: # -------------------------------------------------------------------------------------------- -# X Set of outputs +# X Output of hmm of last n timestep +# k Number of timesteps to predict the output for # OUTPUT: # -------------------------------------------------------------------------------------------- -# P Probability of the set of outputs +# outputs Output of hmm for k timesteps # -------------------------------------------------------------------------------------------- choose_value = function (matrix[double] prob_dist) return (double val) { - #useful for selecting output or states given probability distribution - - #choosing using cumulative distribution + # for selecting output or states given probability distribution + + # choosing using cumulative distribution cum_prob = cumsum(prob_dist) R = rand(rows=1, cols=1, min=0, max=1) - #adding 1 to colsum because indexing starts from 1 + # adding 1 to colsum because indexing starts from 1 val = as.scalar(colSums(cum_prob <= R) + 1) } @@ -83,4 +84,3 @@ m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] outp states = get_all_states(A, ip, k) outputs = get_outputs(B, states) } - From cb6b9b3141c717d3f352c57f72e18f619272a4ab Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 13:19:14 +0200 Subject: [PATCH 17/28] some more comments --- scripts/builtin/hmmPredict.dml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 0b48189009a..ea75d8234b6 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -24,6 +24,7 @@ # -------------------------------------------------------------------------------------------- # X Output of hmm of last n timestep # k Number of timesteps to predict the output for +# model 'd' for discrete outputs and 'c' for continuos # OUTPUT: # -------------------------------------------------------------------------------------------- # outputs Output of hmm for k timesteps @@ -65,13 +66,16 @@ get_outputs = function (matrix[double] B, matrix[double] states) return (matrix[ } } -m_hmmPredict = function(Matrix[Double] X, Integer k) return (Matrix[Double] outputs) +m_hmmPredict = function(Matrix[Double] X, Integer k, model='d') return (Matrix[Double] outputs) { /* A, B, ip = hmm(X) The above would be the function call to obtain model parameters. Instead, the following has dummy model parameters to for the functions to work + + Also only discrete model implemented, need to discretize the cotinous outputs + to work with the continous model */ n_states = 3 From 6b3c01202bcb9704ec4cbd02943eddf48ccef08c Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 13:53:26 +0200 Subject: [PATCH 18/28] standardized and working --- scripts/builtin/hmmPredict.dml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index ea75d8234b6..64f5dbf7f35 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -30,7 +30,7 @@ # outputs Output of hmm for k timesteps # -------------------------------------------------------------------------------------------- -choose_value = function (matrix[double] prob_dist) return (double val) +choose_value = function (Matrix[Double] prob_dist) return (Double val) { # for selecting output or states given probability distribution @@ -43,7 +43,7 @@ choose_value = function (matrix[double] prob_dist) return (double val) val = as.scalar(colSums(cum_prob <= R) + 1) } -get_all_states = function (matrix[double] A, matrix[double] ip, integer k) return (matrix[double] states) +get_all_states = function (Matrix[Double] A, Matrix[Double] ip, integer k) return (Matrix[Double] states) { init_state = choose_value(ip) states = matrix(0, rows=k, cols=1) @@ -57,7 +57,7 @@ get_all_states = function (matrix[double] A, matrix[double] ip, integer k) retur } } -get_outputs = function (matrix[double] B, matrix[double] states) return (matrix[double] outputs) +get_outputs = function (Matrix[Double] B, Matrix[Double] states) return (Matrix[Double] outputs) { outputs = matrix(0, rows=nrow(states), cols=1) parfor (i in 1:nrow(states)) { @@ -66,7 +66,7 @@ get_outputs = function (matrix[double] B, matrix[double] states) return (matrix[ } } -m_hmmPredict = function(Matrix[Double] X, Integer k, model='d') return (Matrix[Double] outputs) +m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] outputs) { /* A, B, ip = hmm(X) @@ -87,4 +87,4 @@ m_hmmPredict = function(Matrix[Double] X, Integer k, model='d') return (Matrix[D states = get_all_states(A, ip, k) outputs = get_outputs(B, states) -} +} \ No newline at end of file From 69c8e1b06994b8ee86174eb0efffcb768cbc29af Mon Sep 17 00:00:00 2001 From: Sandesh Date: Mon, 3 Jul 2023 17:55:28 +0200 Subject: [PATCH 19/28] hmm baum welch implementation --- scripts/builtin/hmm.dml | 178 +++++++++++++++++++++++++++++++--------- 1 file changed, 137 insertions(+), 41 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index b63fb61f436..0bfb7b3f945 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -22,64 +22,160 @@ # This script implements the hidden markov model method # INPUT: # -------------------------------------------------------------------------------------------- -# X Set of outputs +# X Set of outputs of last n timesteps +# # OUTPUT: # -------------------------------------------------------------------------------------------- -# P Probability of the set of outputs +# outputs Probability of the set of outputs # -------------------------------------------------------------------------------------------- -m_hmm = function(Matrix[Double] X, String model="N") return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) +m_hmm = function(Matrix[Double] X) return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) { - """ - if(model != "N" & model != "S") - # N is for number model and S for string model - stop("model not supported, should be in S or N"); - - if (model == "S") - """ - - # change later - n_states = 5 + #X should have the size of 1 * ncols + + #should be transposed for the unique function + unique_X = unique(matrix(X, rows=ncol(X), cols=1)) + nr_outputs = length(unique_X)) + + /* + Since nr of states if unknown, fit a hmm model for every total number of states + from 1 to 10 or 1 to log(n_timesteps), depending on whichever is greater. If the likelihood + decreases with increase in number of total states, break and take the paraemters + of the last iteration as the optimal one. + the last + */ + + max_states = 10 + T = ncol(X) + + if (10 > log(T)) { + max_states = 10 + } else { + max_states = log(T) + } + + search = TRUE + nr_states = 2 + while (search) { + A, B, ip, curr_ll = baum_welch(X, nr_states) + if (nr_states == 2) { + prev_ll = -1 + } + if (curr_ll < prev_ll) { + search = FALSE + break + } + + prev_ll = curr_ll + nr_states = nr_states+1 + } +} - [P, A, B] = fit(X, n_states) +forward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip) return (Matrix[Double] alpha) +{ + /* + alpha a matrix of size nr_states * T with a cell i,t being probability of + the state being at state i at timestep j and the outputs till timestep j + */ + + T = col(X) + nr_states = row(A) + alpha = matrix(0, rows=nr_states, cols=T) + alpha[ ,1] = ip * X[1,1] + + for (t in 2:T) { + for (i in 1:nr_states) { + alpha[i, t] = B[i, X[t]] * sum(alpha[, t-1]* A[ ,i]) + } + } } -fit = function(Matrix[Double] X, Integer n_states) return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) +backward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] beta) { - # probability of the state at timestep 1 - init_prob = matrix(1/n_states, rows=1, cols=n_states) + /* + alpha a matrix of size nr_states * T with a cell i,t being probability of + the model producing outputs (o_t+1,..., o_T) given that the model is + at state i at time t + */ - # probabilities of transition from one one state to another - A = matrix(1/n_states, rows=n_states, cols=n_states) - - # probabilities of a state producing an output - B = matrix(1/ncol(X), rows=n_states, cols=ncol(X)) - - for (i in 1:20) { - if (i == 1) - init_prob = matrix(1/n_states, rows=1, cols=n_states) - A = matrix(1/n_states, rows=n_states, cols=n_states) - B = matrix(1/ncol(X), rows=n_states, cols=ncol(X)) - - alpha = forward(init_prob, B, X) - beta = backward(B, A, X) - - initial_prob, A, B = calc_init_prob(alpha, beta) + T = col(X) + nr_states = row(A) + beta = matrix(0, rows=nr_states, cols=T) + beta[,T] = matrix(1, rows=length(ip), cols=1) + + for (t in (T-1):1) { + for (i in 1:nr_states) { + beta[i, t] = sum(beta[, t+1] * A[i, ] * B[ , X[t]]) + } } } -forward = function(Matrix[Double] init_prob, Matrix[Double] B, Matrix[Double] X) return (Matrix[Double] alpha) +calculate_gamma = function (Matrix[Double] alpha, Matrix[Double] beta) return (Matrix[Double] gamma) { - alpha = rand(rows=nrow(A), cols = ncol(A)) + /* + gamma a nr_state * T matrix with cell (i, t) being probability of + the state being at i at timestep j given the observed output + */ + nr_states = nrow(alpha) + T = ncol(alpha) + gamma = matrix(1/nr + parfor (i in 1:nrow(alpha)) { + for (t in 1:ncol(alpha)) { + num_ij = alpha[i, t] * beta[i, t] + den_ij = sum(alpha[,t] * beta[,t] + gamma[i, j] = num_ij / den_ij + } + } } - -backward = function(Matrix[Double] B, Matrix[Double] A, Matrix[Double] X) return (Matrix[Double] beta) -{ - beta = rand(rows=nrow(A), cols=ncol(A)) +calculate_eta = function (Matrix[Double] alpha, Matrix[Double] beta, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] eta) +{ + /* + gamma a (nr_states * nr_states) * T matrix with cell (i, t) being probability of + the state being at i at timestep j given the observed output + */ + nr_states = nrow(alpha) + T = ncol(alpha) + tot_transitions = nr_states * nr_states + eta = matrix(1/nr_states, rows=tot_transitions, cols=T-1) + + /* + The transitions will be indiced as such:transition 1-N will represent + transition from state 1 to 1, 2 upto state N. transition N+1-2N will represent + transitions from 2 to 1, 2 upto state N + */ + + parfor (trans_id in 1:tot_transitions) { + for (t in 1:(T-1)) { + #indices for alpha and beta + i = floor(trans_id / nr_states) + 1 #index starts at 1 + j = trans_id - ((i-1) * nr_states) + num_ij = alpha[i, t] * A[i, j] * beta[j, t+1] * B[j, X[t+1]] + den_ij = sum(sum(alpha[, t])) + eta[trans_id, t] = alpha[i, j] + } + } } -calc_init_prob = function(Matrix[Double] alpha, Matrix[Double] beta) return (Matrix[Double] pie) + +baum_welch = function (Matrix[Double] X, nr_states, nr_outputs) return (Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip, likelihood) { - pie = rand(rows=nrow(alpha), cols=ncol(beta)) + #initialize state transition and emmission matrices uniformly + A = matrix(1/nr_states, rows=nr_states, cols=nr_states) + B = matrix(1/nr_outputs, rows=nr_states, cols=nr_ouptuts) + ip = matrix(1/nr_states, rows=nr_states, cols=1) + + converge = FALSE + + while (!converge) { + alpha = forward(X, A, B, ip) + beta = backward(X, A, B) + + gamma = calculate_gamma(alpha) + eta = calculate_eta(alpha, A, B) + /* + TODO: compute likelihood, if it does not change much from previous + iteration, break. + */ + } } \ No newline at end of file From 947ae7de27c15a728cc13ae0abc93b3ec502949d Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sat, 8 Jul 2023 17:07:53 +0200 Subject: [PATCH 20/28] ouput definition added --- scripts/builtin/hmm.dml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 0bfb7b3f945..fea9ee50762 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -26,10 +26,12 @@ # # OUTPUT: # -------------------------------------------------------------------------------------------- -# outputs Probability of the set of outputs +# ip Probability of starting at each state at first timestep +# A Probability of transitioning from each state to each other state +# B Probability of emmiting each output from each state # -------------------------------------------------------------------------------------------- -m_hmm = function(Matrix[Double] X) return (Matrix[Double] P, Matrix[Double] A, Matrix[Double] B) +m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, Matrix[Double] B) { #X should have the size of 1 * ncols @@ -57,7 +59,7 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] P, Matrix[Double] A, M search = TRUE nr_states = 2 while (search) { - A, B, ip, curr_ll = baum_welch(X, nr_states) + A, B, sp, curr_ll = baum_welch(X, nr_states) if (nr_states == 2) { prev_ll = -1 } @@ -68,6 +70,8 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] P, Matrix[Double] A, M prev_ll = curr_ll nr_states = nr_states+1 + + #TODO add break condition } } From 99a805bd671ff487a2e074bcd00dcdb96651b07f Mon Sep 17 00:00:00 2001 From: Sandesh Date: Sat, 8 Jul 2023 17:15:32 +0200 Subject: [PATCH 21/28] changed docu and structure --- scripts/builtin/hmmPredict.dml | 47 +++++++++++++++++----------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 64f5dbf7f35..77fe8fa597a 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -24,12 +24,34 @@ # -------------------------------------------------------------------------------------------- # X Output of hmm of last n timestep # k Number of timesteps to predict the output for -# model 'd' for discrete outputs and 'c' for continuos # OUTPUT: # -------------------------------------------------------------------------------------------- # outputs Output of hmm for k timesteps # -------------------------------------------------------------------------------------------- +m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] outputs) +{ + /* + A, B, ip = hmm(X) + + The above would be the function call to obtain model parameters. + Instead, the following has dummy model parameters to for the functions to work + + Also only discrete model implemented, need to discretize the cotinous outputs + to work with the continous model + */ + + n_states = 3 + n_outputs = 3 + + A = matrix(1/n_states, rows=n_states, cols=n_states) + B = matrix(1/n_outputs, rows=n_states, cols=n_outputs) + ip = matrix(1/n_states, rows=n_states, cols=1) + + states = get_all_states(A, ip, k) + outputs = get_outputs(B, states) +} + choose_value = function (Matrix[Double] prob_dist) return (Double val) { # for selecting output or states given probability distribution @@ -64,27 +86,4 @@ get_outputs = function (Matrix[Double] B, Matrix[Double] states) return (Matrix[ s = as.scalar(states[i, 1]) outputs[i, 1] = choose_value(t(B[s, ])) } -} - -m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] outputs) -{ - /* - A, B, ip = hmm(X) - - The above would be the function call to obtain model parameters. - Instead, the following has dummy model parameters to for the functions to work - - Also only discrete model implemented, need to discretize the cotinous outputs - to work with the continous model - */ - - n_states = 3 - n_outputs = 3 - - A = matrix(1/n_states, rows=n_states, cols=n_states) - B = matrix(1/n_outputs, rows=n_states, cols=n_outputs) - ip = matrix(1/n_states, rows=n_states, cols=1) - - states = get_all_states(A, ip, k) - outputs = get_outputs(B, states) } \ No newline at end of file From 0e28eb70e7dc6bb617b27e7639d2043fc745a972 Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Mon, 17 Jul 2023 18:22:16 +0200 Subject: [PATCH 22/28] working HMM implementation --- scripts/builtin/hmm.dml | 294 ++++++++++++++++++++++++---------------- 1 file changed, 174 insertions(+), 120 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index fea9ee50762..152bb21c37c 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -1,147 +1,153 @@ -#------------------------------------------------------------- -# -# Licensed to the Apache Software Foundation (ASF) under one -# or more contributor license agreements. See the NOTICE file -# distributed with this work for additional information -# regarding copyright ownership. The ASF licenses this file -# to you under the Apache License, Version 2.0 (the -# "License"); you may not use this file except in compliance -# with the License. You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, -# software distributed under the License is distributed on an -# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY -# KIND, either express or implied. See the License for the -# specific language governing permissions and limitations -# under the License. -# -#------------------------------------------------------------- - -# This script implements the hidden markov model method -# INPUT: -# -------------------------------------------------------------------------------------------- -# X Set of outputs of last n timesteps -# -# OUTPUT: -# -------------------------------------------------------------------------------------------- -# ip Probability of starting at each state at first timestep -# A Probability of transitioning from each state to each other state -# B Probability of emmiting each output from each state -# -------------------------------------------------------------------------------------------- - m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, Matrix[Double] B) { + #X should have the size of 1 * ncols - #should be transposed for the unique function - unique_X = unique(matrix(X, rows=ncol(X), cols=1)) - nr_outputs = length(unique_X)) + unique_X = matrix(X, rows=ncol(X), cols=1) + nr_outputs = unique_vals(unique_X) /* Since nr of states if unknown, fit a hmm model for every total number of states - from 1 to 10 or 1 to log(n_timesteps), depending on whichever is greater. If the likelihood - decreases with increase in number of total states, break and take the paraemters - of the last iteration as the optimal one. - the last + from 1 to 10 or 1 to log(n_timesteps), depending on whichever is greater. Use only 50 iterations + for each nr of states(since it is computationally expensive). + + Also the break condition is: If the likelihood decreases with increase in number of total states, + break and take the paremeters of the last iteration as the optimal one. + + + After finding out which nr_of_states yields highest likelihood, find the parameters for that + with higher iterations. 1000 is taken here */ - max_states = 10 + max_states = 6 T = ncol(X) - if (10 > log(T)) { - max_states = 10 - } else { - max_states = log(T) + if (max_states < log(T)){ + max_states = ceil(log(T)) } search = TRUE nr_states = 2 + nr_states_likelihood = matrix(0, rows=max_states-1, cols=1) + + while (search) { - A, B, sp, curr_ll = baum_welch(X, nr_states) + if (nr_states > max_states) { + search = FALSE + } + + [A, B, ip, curr_ll] = baum_welch(X, nr_states, nr_outputs, 50) + nr_states_likelihood[nr_states-1, 1] = curr_ll if (nr_states == 2) { prev_ll = -1 } if (curr_ll < prev_ll) { search = FALSE - break } - + prev_ll = curr_ll nr_states = nr_states+1 - - #TODO add break condition } + + best_nr_states = max_ll_state(nr_states_likelihood) + print('found optimal number of states.') + print('Now, caulculating parameters for the optimal number of states') + [A, B, ip, ll] = baum_welch(X, best_nr_states, nr_outputs, 100) } -forward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip) return (Matrix[Double] alpha) -{ - /* - alpha a matrix of size nr_states * T with a cell i,t being probability of - the state being at state i at timestep j and the outputs till timestep j - */ - - T = col(X) - nr_states = row(A) - alpha = matrix(0, rows=nr_states, cols=T) - alpha[ ,1] = ip * X[1,1] +unique_vals = function(Matrix[Double] X) return (Integer l) { + # group and count duplicates + I = aggregate (target = X, groups = X[,1], fn = "count") + # select groups + res = removeEmpty (target = seq (1, max (X[,1])), margin = "rows", select = (I != 0)) + l = length(res) +} - for (t in 2:T) { - for (i in 1:nr_states) { - alpha[i, t] = B[i, X[t]] * sum(alpha[, t-1]* A[ ,i]) + + +max_ll_state = function(Matrix[Double] nr_states_likelihood) return (Integer best_state) { + max_val = max(nr_states_likelihood) + S = nrow(nr_states_likelihood) + best_state = -1 + + for (i in 1:S){ + if (max_val == as.scalar(nr_states_likelihood[i, 1])) { + best_state = i +1 } - } + } } -backward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] beta) +baum_welch = function( Matrix[Double] X, Integer nr_states, Integer nr_outputs, Integer IT) return (Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip, Double likelihood) { - /* - alpha a matrix of size nr_states * T with a cell i,t being probability of - the model producing outputs (o_t+1,..., o_T) given that the model is - at state i at time t - */ - - T = col(X) - nr_states = row(A) - beta = matrix(0, rows=nr_states, cols=T) - beta[,T] = matrix(1, rows=length(ip), cols=1) + #initialize state transition and emmission matrices uniformly + A = rand(rows=nr_states, cols=nr_states) + B = rand(rows=nr_states, cols=nr_outputs) + ip = rand(rows=nr_states, cols=1) + likelihood = 0 + T = ncol(X) - for (t in (T-1):1) { - for (i in 1:nr_states) { - beta[i, t] = sum(beta[, t+1] * A[i, ] * B[ , X[t]]) + for (i in 1:IT) { + alpha = forward(X, A, B, ip) + beta = backward(X, A, B) + + gamma = calculate_gamma(alpha, beta) + eta = calculate_eta(X, alpha, beta, A, B) + + ip = gamma[, 1] + A = calculate_A(eta, gamma, nr_states) + B = calculate_B(gamma, X, nr_states, nr_outputs) + + #likelhood using forward method + likelihood = sum(gamma[, T] * alpha[, T]) + } +} + +calculate_A = function(Matrix[Double] eta, Matrix[Double] gamma, Integer nr_states) return (Matrix[Double] A) { + A = rand(rows=nr_states, cols=nr_states) + T = ncol(eta) + + parfor (i in 1:nr_states) { + parfor (j in 1:nr_states) { + eta_id = (i-1) * nr_states + j + num_ij = sum(eta[eta_id, 1:T-1]) + den_ij = sum(gamma[i,1:T-1]) + A[i,j] = num_ij/den_ij } } } -calculate_gamma = function (Matrix[Double] alpha, Matrix[Double] beta) return (Matrix[Double] gamma) -{ - /* - gamma a nr_state * T matrix with cell (i, t) being probability of - the state being at i at timestep j given the observed output - */ - nr_states = nrow(alpha) - T = ncol(alpha) - gamma = matrix(1/nr - parfor (i in 1:nrow(alpha)) { - for (t in 1:ncol(alpha)) { - num_ij = alpha[i, t] * beta[i, t] - den_ij = sum(alpha[,t] * beta[,t] - gamma[i, j] = num_ij / den_ij +calculate_B = function(Matrix[Double] gamma, Matrix[Double] X, Integer nr_states, Integer nr_outputs) return (Matrix[Double] B) { + B = rand(rows=nr_states, cols=nr_outputs) + + parfor (i in 1:nr_states){ + parfor (j in 1:nr_outputs) { + B[i, j] = sum((X == j) * gamma[i, ]) / sum(gamma[i, ]) + } } - } } -calculate_eta = function (Matrix[Double] alpha, Matrix[Double] beta, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] eta) +/* + To work with 3d matrices, the values in first two dimensions of eta matrix + must be rolled the first dimension. The following indexing function maps the + index of rolled matrix(eta) to "normal" 2d matrices like alpha, beta, etc +*/ +index = function(Integer rolled_i, Integer nr_states) return(Double i, Double j) { + i = ceil(rolled_i / nr_states) #index starts at 1 + j = rolled_i - (nr_states * (i-1)) #i-1) * +} + +calculate_eta = function (Matrix[Double] X, Matrix[Double] alpha, Matrix[Double] beta, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] eta) { + nr_states = nrow(alpha) + T = ncol(alpha) + /* - gamma a (nr_states * nr_states) * T matrix with cell (i, t) being probability of + eta a (nr_states * nr_states) * T matrix with cell (i, t) being probability of the state being at i at timestep j given the observed output */ - nr_states = nrow(alpha) - T = ncol(alpha) + tot_transitions = nr_states * nr_states - eta = matrix(1/nr_states, rows=tot_transitions, cols=T-1) + eta = rand(rows=tot_transitions, cols=T-1) /* The transitions will be indiced as such:transition 1-N will represent @@ -150,36 +156,84 @@ calculate_eta = function (Matrix[Double] alpha, Matrix[Double] beta, Matrix[Dou */ parfor (trans_id in 1:tot_transitions) { + [i, j] = index(trans_id, nr_states) for (t in 1:(T-1)) { #indices for alpha and beta - i = floor(trans_id / nr_states) + 1 #index starts at 1 - j = trans_id - ((i-1) * nr_states) - num_ij = alpha[i, t] * A[i, j] * beta[j, t+1] * B[j, X[t+1]] - den_ij = sum(sum(alpha[, t])) - eta[trans_id, t] = alpha[i, j] + ot1 = output_t(X, t+1) + num_ij = alpha[i, t] * A[i, j] * beta[j, t+1] * B[j, ot1] + den = matrix(0, rows=nr_states, cols=1) + + parfor (k in 1:nr_states) { + den[k, 1] = sum(alpha[k, t] * t(A[k, ]) * beta[, t+1] * B[,ot1]) + } + s = num_ij / sum(den) + eta[trans_id, t] = as.scalar(s[1,1]) + } + } +} + +calculate_gamma = function (Matrix[Double] alpha, Matrix[Double] beta) return (Matrix[Double] gamma) +{ + /* + gamma a nr_state * T matrix with cell (i, t) being probability of + the state being at i at timestep j given the observed output + */ + nr_states = nrow(alpha) + T = ncol(alpha) + + gamma = rand(rows=nr_states, cols=T) + + parfor (i in 1:nrow(alpha)) { + parfor (t in 1:ncol(alpha)) { + num_it = alpha[i, t] * beta[i, t] + den_it = sum(alpha[,t] * beta[,t]) + gamma[i, t] = num_it / den_it } } } +output_t = function (Matrix[Double] X, Integer t) return (Integer o) { + o = as.integer(as.scalar(X[1, t])) +} + +forward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip) return (Matrix[Double] alpha) +{ + /* + alpha a matrix of size nr_states * T with a cell i,t being probability of + the state being at state i at timestep j and the outputs till timestep j + */ + + T = ncol(X) + nr_states = nrow(A) + alpha = matrix(0, rows=nr_states, cols=T) + o_1 = output_t(X, 1) + alpha[ ,1] = ip * B[ ,o_1] + + for (t in 2:T) { + ot = output_t(X, t) + for (i in 1:nr_states) { + alpha[i, t] = B[i, ot] * sum(alpha[,t-1]* A[ ,i]) + } + } +} -baum_welch = function (Matrix[Double] X, nr_states, nr_outputs) return (Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip, likelihood) +backward = function (Matrix[Double] X, Matrix[Double] A, Matrix[Double] B) return (Matrix[Double] beta) { - #initialize state transition and emmission matrices uniformly - A = matrix(1/nr_states, rows=nr_states, cols=nr_states) - B = matrix(1/nr_outputs, rows=nr_states, cols=nr_ouptuts) - ip = matrix(1/nr_states, rows=nr_states, cols=1) + /* + alpha a matrix of size nr_states * T with a cell i,t being probability of + the model producing outputs (o_t+1,..., o_T) given that the model is + at state i at time t + */ - converge = FALSE + T = ncol(X) + nr_states = nrow(A) + beta = matrix(1/T, rows=nr_states, cols=T) + beta[,T] = matrix(1/T, rows=nr_states, cols=1) - while (!converge) { - alpha = forward(X, A, B, ip) - beta = backward(X, A, B) - - gamma = calculate_gamma(alpha) - eta = calculate_eta(alpha, A, B) - /* - TODO: compute likelihood, if it does not change much from previous - iteration, break. - */ + for (t in (T-1):1) { + ot = output_t(X, t) + for (i in 1:nr_states) { + beta[i, t] = sum(beta[, t+1] * t(A[i, ]) * B[ , ot]) + } } } \ No newline at end of file From 6faa754bcb4e9e85f9820c69a423e158bb5e899e Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Mon, 17 Jul 2023 20:12:29 +0200 Subject: [PATCH 23/28] fixed hmm documentation and indents in hmmPredict --- scripts/builtin/hmm.dml | 34 ++++++++++++++++++++++++++++++++-- scripts/builtin/hmmPredict.dml | 15 ++------------- 2 files changed, 34 insertions(+), 15 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 152bb21c37c..497a39b72c3 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -1,3 +1,35 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# This script implements the hidden markov model method +# INPUT: +# -------------------------------------------------------------------------------------------- +# X Output of hmm of last n timestep +# OUTPUT: +# -------------------------------------------------------------------------------------------- +# ip Start probabilities +# A Transition probabilities +# B Emission probabilities +# -------------------------------------------------------------------------------------------- + m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, Matrix[Double] B) { @@ -63,8 +95,6 @@ unique_vals = function(Matrix[Double] X) return (Integer l) { l = length(res) } - - max_ll_state = function(Matrix[Double] nr_states_likelihood) return (Integer best_state) { max_val = max(nr_states_likelihood) S = nrow(nr_states_likelihood) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 77fe8fa597a..72206315e61 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -32,22 +32,11 @@ m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] outputs) { /* - A, B, ip = hmm(X) - - The above would be the function call to obtain model parameters. - Instead, the following has dummy model parameters to for the functions to work - Also only discrete model implemented, need to discretize the cotinous outputs to work with the continous model */ - n_states = 3 - n_outputs = 3 - - A = matrix(1/n_states, rows=n_states, cols=n_states) - B = matrix(1/n_outputs, rows=n_states, cols=n_outputs) - ip = matrix(1/n_states, rows=n_states, cols=1) - + [ip, A, B] = hmm(X) states = get_all_states(A, ip, k) outputs = get_outputs(B, states) } @@ -65,7 +54,7 @@ choose_value = function (Matrix[Double] prob_dist) return (Double val) val = as.scalar(colSums(cum_prob <= R) + 1) } -get_all_states = function (Matrix[Double] A, Matrix[Double] ip, integer k) return (Matrix[Double] states) +get_all_states = function (Matrix[Double] A, Matrix[Double] ip, Integer k) return (Matrix[Double] states) { init_state = choose_value(ip) states = matrix(0, rows=k, cols=1) From 40d1d48f7e41912a0cc29ace55d3c21ae24770cd Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Wed, 19 Jul 2023 11:03:33 +0200 Subject: [PATCH 24/28] hmm and hmmPredict work now!! --- scripts/builtin/hmm.dml | 8 +++----- scripts/builtin/hmmPredict.dml | 1 - 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 497a39b72c3..67836cfb568 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -40,15 +40,14 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, /* Since nr of states if unknown, fit a hmm model for every total number of states - from 1 to 10 or 1 to log(n_timesteps), depending on whichever is greater. Use only 50 iterations - for each nr of states(since it is computationally expensive). + from 1 to some max_states or 1 to log(n_timesteps), depending on whichever is greater. Use only low number of + iterations at first for each nr of states(since it is computationally expensive). Also the break condition is: If the likelihood decreases with increase in number of total states, break and take the paremeters of the last iteration as the optimal one. - After finding out which nr_of_states yields highest likelihood, find the parameters for that - with higher iterations. 1000 is taken here + with higher iterations. */ max_states = 6 @@ -62,7 +61,6 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, nr_states = 2 nr_states_likelihood = matrix(0, rows=max_states-1, cols=1) - while (search) { if (nr_states > max_states) { search = FALSE diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 72206315e61..673a7deaffa 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -35,7 +35,6 @@ m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] out Also only discrete model implemented, need to discretize the cotinous outputs to work with the continous model */ - [ip, A, B] = hmm(X) states = get_all_states(A, ip, k) outputs = get_outputs(B, states) From 948af5b2d71d89469f6f8934032bce0e62268f34 Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Fri, 21 Jul 2023 08:27:50 +0200 Subject: [PATCH 25/28] refactored for reproducibility --- scripts/builtin/hmmPredict.dml | 55 +++++++++++++++++----------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/scripts/builtin/hmmPredict.dml b/scripts/builtin/hmmPredict.dml index 673a7deaffa..dd838d6d8fe 100644 --- a/scripts/builtin/hmmPredict.dml +++ b/scripts/builtin/hmmPredict.dml @@ -32,46 +32,45 @@ m_hmmPredict = function (Matrix[Double] X, Integer k) return (Matrix[Double] outputs) { /* - Also only discrete model implemented, need to discretize the cotinous outputs + Only discrete model implemented, need to discretize the cotinous outputs to work with the continous model */ [ip, A, B] = hmm(X) - states = get_all_states(A, ip, k) - outputs = get_outputs(B, states) -} - -choose_value = function (Matrix[Double] prob_dist) return (Double val) -{ - # for selecting output or states given probability distribution - - # choosing using cumulative distribution - cum_prob = cumsum(prob_dist) - - R = rand(rows=1, cols=1, min=0, max=1) - - # adding 1 to colsum because indexing starts from 1 - val = as.scalar(colSums(cum_prob <= R) + 1) -} + /* + Also only discrete model implemented, need to discretize the cotinous outputs + to work with the continous model + */ + seed = 42 + init_state = choose_value(ip, seed) -get_all_states = function (Matrix[Double] A, Matrix[Double] ip, Integer k) return (Matrix[Double] states) -{ - init_state = choose_value(ip) states = matrix(0, rows=k, cols=1) + states[1, 1] = init_state s = init_state - for (i in 2:k) { - next_state = choose_value(t(A[s, ])) + seed=seed +1 + next_state = choose_value(t(A[s, ]), seed) states[i, 1] = next_state s = next_state } -} -get_outputs = function (Matrix[Double] B, Matrix[Double] states) return (Matrix[Double] outputs) -{ outputs = matrix(0, rows=nrow(states), cols=1) - parfor (i in 1:nrow(states)) { + for (i in 1:nrow(states)) { + seed = seed+1 s = as.scalar(states[i, 1]) - outputs[i, 1] = choose_value(t(B[s, ])) + outputs[i, 1] = choose_value(t(B[s, ]), seed) } -} \ No newline at end of file +} + +choose_value = function (Matrix[Double] prob_dist, Integer seed) return (Double val) +{ + # for selecting output or states given probability distribution + + # choosing using cumulative distribution + cum_prob = cumsum(prob_dist) + + R = rand(rows=1, cols=1, min=0, max=1, seed=seed) + + # adding 1 to colsum because indexing starts from 1 + val = as.scalar(colSums(cum_prob <= R) + 1) +} From 7602afc1414c2e7b07eea8681a57f1105cbf4a3a Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Fri, 21 Jul 2023 08:29:45 +0200 Subject: [PATCH 26/28] added more iterations for optimization --- scripts/builtin/hmm.dml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 67836cfb568..9665b920afc 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -66,7 +66,7 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, search = FALSE } - [A, B, ip, curr_ll] = baum_welch(X, nr_states, nr_outputs, 50) + [A, B, ip, curr_ll] = baum_welch(X, nr_states, nr_outputs, 100) nr_states_likelihood[nr_states-1, 1] = curr_ll if (nr_states == 2) { prev_ll = -1 @@ -82,7 +82,7 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, best_nr_states = max_ll_state(nr_states_likelihood) print('found optimal number of states.') print('Now, caulculating parameters for the optimal number of states') - [A, B, ip, ll] = baum_welch(X, best_nr_states, nr_outputs, 100) + [A, B, ip, ll] = baum_welch(X, best_nr_states, nr_outputs, 200) } unique_vals = function(Matrix[Double] X) return (Integer l) { @@ -96,7 +96,7 @@ unique_vals = function(Matrix[Double] X) return (Integer l) { max_ll_state = function(Matrix[Double] nr_states_likelihood) return (Integer best_state) { max_val = max(nr_states_likelihood) S = nrow(nr_states_likelihood) - best_state = -1 + best_state = 0 for (i in 1:S){ if (max_val == as.scalar(nr_states_likelihood[i, 1])) { From 60468794b4ed1eb0915ce1cfe88b5a58eb609034 Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Sun, 23 Jul 2023 21:55:34 +0200 Subject: [PATCH 27/28] changed hyper-parameters and deterministic behavior --- scripts/builtin/hmm.dml | 61 ++++++++++++++++++----------------------- 1 file changed, 26 insertions(+), 35 deletions(-) diff --git a/scripts/builtin/hmm.dml b/scripts/builtin/hmm.dml index 9665b920afc..02fafa6c054 100644 --- a/scripts/builtin/hmm.dml +++ b/scripts/builtin/hmm.dml @@ -40,14 +40,10 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, /* Since nr of states if unknown, fit a hmm model for every total number of states - from 1 to some max_states or 1 to log(n_timesteps), depending on whichever is greater. Use only low number of - iterations at first for each nr of states(since it is computationally expensive). + from 1 to some max_states or 1 to log(n_timesteps), depending on whichever is greater. Also the break condition is: If the likelihood decreases with increase in number of total states, break and take the paremeters of the last iteration as the optimal one. - - After finding out which nr_of_states yields highest likelihood, find the parameters for that - with higher iterations. */ max_states = 6 @@ -59,58 +55,53 @@ m_hmm = function(Matrix[Double] X) return (Matrix[Double] ip, Matrix[Double] A, search = TRUE nr_states = 2 - nr_states_likelihood = matrix(0, rows=max_states-1, cols=1) + seed = 42 while (search) { - if (nr_states > max_states) { - search = FALSE - } - - [A, B, ip, curr_ll] = baum_welch(X, nr_states, nr_outputs, 100) - nr_states_likelihood[nr_states-1, 1] = curr_ll + [A_new, B_new, ip_new, curr_ll, seed] = baum_welch(X, nr_states, nr_outputs, 20, seed) + seed = seed+1 if (nr_states == 2) { prev_ll = -1 } if (curr_ll < prev_ll) { search = FALSE } - + else { + A = A_new + B = B_new + ip = ip_new + } prev_ll = curr_ll nr_states = nr_states+1 + + if (nr_states > max_states) { + search = FALSE + } } - - best_nr_states = max_ll_state(nr_states_likelihood) - print('found optimal number of states.') - print('Now, caulculating parameters for the optimal number of states') - [A, B, ip, ll] = baum_welch(X, best_nr_states, nr_outputs, 200) } -unique_vals = function(Matrix[Double] X) return (Integer l) { +col_normalize = function(Matrix[Double] m) return (Matrix[Double] m_norm) { + m_norm = m / rowSums(m) +} + +unique_vals = function (Matrix[Double] X) return (Integer l) { # group and count duplicates I = aggregate (target = X, groups = X[,1], fn = "count") # select groups res = removeEmpty (target = seq (1, max (X[,1])), margin = "rows", select = (I != 0)) l = length(res) } - -max_ll_state = function(Matrix[Double] nr_states_likelihood) return (Integer best_state) { - max_val = max(nr_states_likelihood) - S = nrow(nr_states_likelihood) - best_state = 0 - - for (i in 1:S){ - if (max_val == as.scalar(nr_states_likelihood[i, 1])) { - best_state = i +1 - } - } -} -baum_welch = function( Matrix[Double] X, Integer nr_states, Integer nr_outputs, Integer IT) return (Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip, Double likelihood) +baum_welch = function( Matrix[Double] X, Integer nr_states, Integer nr_outputs, Integer IT, Integer seed) return (Matrix[Double] A, Matrix[Double] B, Matrix[Double] ip, Double likelihood, Integer seed) { #initialize state transition and emmission matrices uniformly - A = rand(rows=nr_states, cols=nr_states) - B = rand(rows=nr_states, cols=nr_outputs) - ip = rand(rows=nr_states, cols=1) + A = col_normalize(matrix(1/nr_states, rows=nr_states, cols=nr_states) + rand(rows=nr_states, cols=nr_states, seed=seed)) + seed = seed+1 + B = col_normalize(matrix(1/nr_outputs,rows=nr_states, cols=nr_outputs) + rand(rows=nr_states, cols=nr_outputs, seed=seed)) + seed = seed+1 + ip = matrix(1/nr_states, rows=nr_states, cols=1) + rand(rows=nr_states, cols=1, seed=seed) + seed = seed+1 + likelihood = 0 T = ncol(X) From 5f2c8648aa5699a4aa0e7cf162a5112d087b65cb Mon Sep 17 00:00:00 2001 From: sb0476679 Date: Sun, 23 Jul 2023 21:56:15 +0200 Subject: [PATCH 28/28] the tests should be working now --- .../builtin/part1/BuiltinHMMPredictTest.java | 65 +++++++-------- .../builtin/part1/BuiltinHMMTest.java | 81 ++++++++++++------- 2 files changed, 79 insertions(+), 67 deletions(-) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java index 4b39289e110..feb1174abd6 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMPredictTest.java @@ -19,11 +19,13 @@ package org.apache.sysds.test.functions.builtin.part1; +import java.util.HashMap; import java.util.Random; import org.apache.sysds.api.DMLScript; import org.apache.sysds.common.Types; import org.apache.sysds.common.Types.ExecType; import org.apache.sysds.hops.OptimizerUtils; +import org.apache.sysds.runtime.matrix.data.MatrixValue; import org.apache.sysds.test.AutomatedTestBase; import org.apache.sysds.test.TestConfiguration; import org.apache.sysds.test.TestUtils; @@ -35,9 +37,8 @@ public class BuiltinHMMPredictTest extends AutomatedTestBase private final static String TEST_NAME = "hmmPredict"; private final static String TEST_DIR = "functions/builtin/"; private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; - private final static int rows = 1320; - private final static int cols = 32; - + private final static double eps = 1e-3; + @Override public void setUp() { TestUtils.clearAssertionInformation(); @@ -49,7 +50,7 @@ public void testHMMPredictCP() { runHMMPredictTest(ExecType.CP); } - @Test + @Test public void testHMMPredictSpark() { runHMMPredictTest(ExecType.SPARK); } @@ -57,52 +58,44 @@ public void testHMMPredictSpark() { private void runHMMPredictTest(ExecType instType) { Types.ExecMode platformOld = setExecMode(instType); - - boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; - boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; - try { loadTestConfiguration(getTestConfiguration(TEST_NAME)); - String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; programArgs = new String[]{ - "-nvargs", "X=" + input("X"), "k=" + input("k"), "Y=" + output("Y")}; - + "-nvargs", "X=" + input("X"), "k=" + input("k"), "outputs=" + output("outputs")}; - //generate actual datasets - // Random random = new Random(); - // Integer k = random.nextInt(10) + 1; - - //generate actual datasets //Output generated from https://en.wikipedia.org/wiki/Hidden_Markov_model#Weather_guessing_game - // C(Clean) -> 0, W(Walk) -> 1, S(Shop) -> 2 - - double[][] X = {{0, 0, 1, - 2, 1, 1, - 1, 0, 0, 1} - }; - - writeInputMatrixWithMTD("X", X, true); - - //TODO: Compare outputs? What is the output of the - // function supposed to be compared against? How? - - setOutputBuffering(true); - - String stdout = runTest(null).toString(); - + double[][] X = {{2, 1, 2, 1, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 3, 2, 1, 2, 3, 2, + 1, 1, 1, 3, 2, 2, 2, 2, 3, 3, 1, 2, 3, 1, 3, 3, 2, 3, 1, 2, 1, + 1, 2, 3, 1, 1, 1, 3, 2, 3, 2, 1, 1, 3, 3, 3, 1, 2, 3, 3, 1, 3, + 3, 3, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 3, 2, 2, 3, 2, 2, 1, 3, 2, + 1, 3, 2, 2, 2, 1, 2, 1, 2, 2, 1, 3, 2, 2, 2, 3, 2}}; + int k = 10; + writeInputMatrixWithMTD("X", X, true); runTest(true, false, null, -1); + + HashMap test_output = readDMLMatrixFromOutputDir("outputs"); + HashMap expected_output = new HashMap<>(); + + expected_output.put(new MatrixValue.CellIndex(1, 1), 1.0); + expected_output.put(new MatrixValue.CellIndex(2, 1), 3.0); + expected_output.put(new MatrixValue.CellIndex(3, 1), 3.0); + expected_output.put(new MatrixValue.CellIndex(4, 1), 2.0); + expected_output.put(new MatrixValue.CellIndex(5, 1), 1.0); + expected_output.put(new MatrixValue.CellIndex(6, 1), 2.0); + expected_output.put(new MatrixValue.CellIndex(7, 1), 3.0); + expected_output.put(new MatrixValue.CellIndex(8, 1), 1.0); + expected_output.put(new MatrixValue.CellIndex(9, 1), 3.0); + expected_output.put(new MatrixValue.CellIndex(10, 1), 1.0); + + TestUtils.compareMatrices(expected_output, test_output, eps, "Expected-DML", "Actual-DML"); } finally { rtplatform = platformOld; - DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; - OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; - OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; - OptimizerUtils.ALLOW_OPERATOR_FUSION = true; } } } diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java index 8a72ee5801b..31766c45031 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinHMMTest.java @@ -15,15 +15,20 @@ * KIND, either express or implied. See the License for the * specific language governing permissions and limitations * under the License. + */ + package org.apache.sysds.test.functions.builtin.part1; + +import java.util.HashMap; import java.util.Random; import org.apache.sysds.api.DMLScript; import org.apache.sysds.common.Types; import org.apache.sysds.common.Types.ExecType; import org.apache.sysds.hops.OptimizerUtils; +import org.apache.sysds.runtime.matrix.data.MatrixValue; import org.apache.sysds.test.AutomatedTestBase; import org.apache.sysds.test.TestConfiguration; import org.apache.sysds.test.TestUtils; @@ -35,13 +40,12 @@ public class BuiltinHMMTest extends AutomatedTestBase private final static String TEST_NAME = "hmm"; private final static String TEST_DIR = "functions/builtin/"; private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinHMMPredictTest.class.getSimpleName() + "/"; - private final static int rows = 1320; - private final static int cols = 32; - + private final static double eps = 1e-3; + @Override public void setUp() { TestUtils.clearAssertionInformation(); - addTestConfiguration(TEST_NAME,new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"C"})); + addTestConfiguration(TEST_NAME, new TestConfiguration(TEST_CLASS_DIR, TEST_NAME, new String[]{"C"})); } @Test @@ -57,46 +61,61 @@ public void testHMMPredictSpark() { private void runHMMPredictTest(ExecType instType) { Types.ExecMode platformOld = setExecMode(instType); - - boolean oldFlag = OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION; - boolean sparkConfigOld = DMLScript.USE_LOCAL_SPARK_CONFIG; - - try - { + try { loadTestConfiguration(getTestConfiguration(TEST_NAME)); String HOME = SCRIPT_DIR + TEST_DIR; - fullDMLScriptName = HOME + TEST_NAME + ".dml"; programArgs = new String[]{ - "-nvargs", "X=" + input("X"), "Y=" + output("Y")}; + "-nvargs", "X=" + input("X"), "ip=" + output("ip"), "A=" + output("A"), "B=" + output("B")}; - //generate actual datasets //Output generated from https://en.wikipedia.org/wiki/Hidden_Markov_model#Weather_guessing_game - // C(Clean) -> 0, W(Walk) -> 1, S(Shop) -> 2 - - double[][] X = {{0, 0, 1, - 2, 1, 1, - 1, 0, 0, 1} - }; + double[][] X = {{2, 1, 2, 1, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 3, 2, 1, 2, 3, 2, + 1, 1, 1, 3, 2, 2, 2, 2, 3, 3, 1, 2, 3, 1, 3, 3, 2, 3, 1, 2, 1, + 1, 2, 3, 1, 1, 1, 3, 2, 3, 2, 1, 1, 3, 3, 3, 1, 2, 3, 3, 1, 3, + 3, 3, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 3, 2, 2, 3, 2, 2, 1, 3, 2, + 1, 3, 2, 2, 2, 1, 2, 1, 2, 2, 1, 3, 2, 2, 2, 3, 2}}; + writeInputMatrixWithMTD("X", X, true); - - setOutputBuffering(true); + runTest(true, false, null, -1); - //TODO: Compare outputs? What is the output of the - // function supposed to be compared against? How? + HashMap out_ip = readDMLMatrixFromOutputDir("ip"); + HashMap out_A = readDMLMatrixFromOutputDir("A"); + HashMap out_B = readDMLMatrixFromOutputDir("B"); + HashMap expected_ip = new HashMap<>(); + HashMap expected_A = new HashMap<>(); + HashMap expected_B = new HashMap<>(); - String stdout = runTest(null).toString(); + expected_ip.put(new MatrixValue.CellIndex(1, 1), 0.731); + expected_ip.put(new MatrixValue.CellIndex(2, 1), 0.000); + expected_ip.put(new MatrixValue.CellIndex(3, 1), 0.269); - - runTest(true, false, null, -1); + expected_A.put(new MatrixValue.CellIndex(1, 1), 0.194); + expected_A.put(new MatrixValue.CellIndex(1, 2), 0.421); + expected_A.put(new MatrixValue.CellIndex(1, 3), 0.424); + expected_A.put(new MatrixValue.CellIndex(2, 1), 0.189); + expected_A.put(new MatrixValue.CellIndex(2, 2), 0.432); + expected_A.put(new MatrixValue.CellIndex(2, 3), 0.383); + expected_A.put(new MatrixValue.CellIndex(3, 1), 0.374); + expected_A.put(new MatrixValue.CellIndex(3, 2), 0.334); + expected_A.put(new MatrixValue.CellIndex(3, 3), 0.257); + + expected_B.put(new MatrixValue.CellIndex(1, 1), 0.383); + expected_B.put(new MatrixValue.CellIndex(1, 2), 0.609); + expected_B.put(new MatrixValue.CellIndex(1, 3), 0.007); + expected_B.put(new MatrixValue.CellIndex(2, 1), 0.152); + expected_B.put(new MatrixValue.CellIndex(2, 2), 0.012); + expected_B.put(new MatrixValue.CellIndex(2, 3), 0.837); + expected_B.put(new MatrixValue.CellIndex(3, 1), 0.373); + expected_B.put(new MatrixValue.CellIndex(3, 2), 0.615); + expected_B.put(new MatrixValue.CellIndex(3, 3), 0.012); + + TestUtils.compareMatrices(expected_ip, out_ip, eps, "Expected-DML", "Actual-DML"); + TestUtils.compareMatrices(expected_A, out_A, eps, "Expected-DML", "Actual-DML"); + TestUtils.compareMatrices(expected_B, out_B, eps, "Expected-DML", "Actual-DML"); } finally { rtplatform = platformOld; - DMLScript.USE_LOCAL_SPARK_CONFIG = sparkConfigOld; - OptimizerUtils.ALLOW_ALGEBRAIC_SIMPLIFICATION = oldFlag; - OptimizerUtils.ALLOW_AUTO_VECTORIZATION = true; - OptimizerUtils.ALLOW_OPERATOR_FUSION = true; } } -} \ No newline at end of file +}