From: miettinen Date: Thu, 27 Jun 2013 10:27:35 +0000 (+0000) Subject: Multiple varied parameters to sensitivity analysis. (refs #4238) X-Git-Tag: 1.8.1~283 X-Git-Url: https://gerrit.simantics.org/r/gitweb?a=commitdiff_plain;h=133bc5f7bd8afee0a73cb63e60519e66abdcb0b2;p=simantics%2Fsysdyn.git Multiple varied parameters to sensitivity analysis. (refs #4238) Now the parameters are defined outside of the distributions. This allows firstly to randomize the parameter values for each run so that no parameter values repeat between individual runs. Secondly this enables using more sophisticated techniques of selecting the sequence of parameter sets (not yet implemented). git-svn-id: https://www.simantics.org/svn/simantics/sysdyn/trunk@27664 ac1ea38d-2e2b-0410-8846-a27921b304fc --- diff --git a/org.simantics.modelica/src/org/simantics/modelica/data/SimulationResult.java b/org.simantics.modelica/src/org/simantics/modelica/data/SimulationResult.java index 369a2814..d4091285 100644 --- a/org.simantics.modelica/src/org/simantics/modelica/data/SimulationResult.java +++ b/org.simantics.modelica/src/org/simantics/modelica/data/SimulationResult.java @@ -21,6 +21,7 @@ import java.io.IOException; import java.io.InputStream; import java.io.LineNumberReader; import java.util.ArrayList; +import java.util.ConcurrentModificationException; import java.util.HashMap; import java.util.List; import java.util.regex.Matcher; @@ -458,12 +459,17 @@ public class SimulationResult { if(resArray.length == 2 && timesArray.length > 2) timesArray = new double[] {timesArray[0], timesArray[timesArray.length - 1]}; ds = new DataSet(variable, timesArray, resArray); - for (DataSet tempds : variables){ - if (tempds.name.equals(variable)) - // We should never need to go there unless some change in - // the logic of reading variables are made. - // Also helps in seeking memory leaks. - return null; + try { + for (DataSet tempds : variables){ + if (tempds.name.equals(variable)) + // We should never need to go there unless some change in + // the logic of reading variables are made. + // Also helps in seeking memory leaks. + System.err.println("Find me and add a comment that I was printed, please."); + // If it seems that in no case we need to go here, comment this try catch. + return null; + } + } catch (ConcurrentModificationException e) { } variables.add(ds); return ds; diff --git a/org.simantics.sysdyn.ontology/graph.tg b/org.simantics.sysdyn.ontology/graph.tg index d11b9bca..a01167be 100644 Binary files a/org.simantics.sysdyn.ontology/graph.tg and b/org.simantics.sysdyn.ontology/graph.tg differ diff --git a/org.simantics.sysdyn.ontology/graph/Sysdyn.pgraph b/org.simantics.sysdyn.ontology/graph/Sysdyn.pgraph index cfce20a5..9b887312 100644 --- a/org.simantics.sysdyn.ontology/graph/Sysdyn.pgraph +++ b/org.simantics.sysdyn.ontology/graph/Sysdyn.pgraph @@ -23,8 +23,8 @@ SYSDYN.ImportedOntologies : PROJ.NamespaceRequirement "http://www.simantics.org/Sysdyn-1.1" : L0.URI "http://www.simantics.org/Layer0-1.1" : L0.URI "http://www.simantics.org/SelectionView-1.2" : L0.URI - "http://www.simantics.org/Documentation-1.1" : L0.URI // Experimental documentation tool - "http://www.simantics.org/DocumentWorkbench-1.0" : L0.URI // Experimental documentation tool +// "http://www.simantics.org/Documentation-1.1" : L0.URI // Experimental documentation tool +// "http://www.simantics.org/DocumentWorkbench-1.0" : L0.URI // Experimental documentation tool SYSDYN.SharedFunctionOntology -- SYSDYN.SensitivityAnalysisExperiment.parameterList --> L0.List -- SYSDYN.SensitivityAnalysisExperiment.randomSeed --> L0.Integer -- SYSDYN.SensitivityAnalysisExperiment.resultRefreshRate --> L0.Integer -- SYSDYN.SensitivityAnalysisExperiment.numberOfValues --> L0.Integer -- SYSDYN.SensitivityAnalysisExperiment.Parameter.propabilityDistribution --> SYSDYN.ProbabilityDistribution -- SYSDYN.SensitivityAnalysisExperiment.Parameter.variable --> L0.String -- SYSDYN.SensitivityAnalysisExperiment.Parameter.indexes --> L0.StringArray -- SYSDYN.SensitivityAnalysisExperiment.Parameter.numberOfValues --> L0.Integer -- SYSDYN.UniformDistribution.minValue --> L0.Double ", - sr.SensitivityAnalysisExperiment_Parameter_numberOfValues, 10, L0.PartOf, experiment); ArrayList parameterList = new ArrayList(); diff --git a/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/SensitivityAnalysisExperimentTab.java b/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/SensitivityAnalysisExperimentTab.java index a88f633f..ae039710 100644 --- a/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/SensitivityAnalysisExperimentTab.java +++ b/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/SensitivityAnalysisExperimentTab.java @@ -1,3 +1,14 @@ +/******************************************************************************* + * Copyright (c) 2007, 2011 Association for Decentralized Information Management in + * Industry THTH ry. + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: + * VTT Technical Research Centre of Finland - initial API and implementation + *******************************************************************************/ package org.simantics.sysdyn.ui.properties; import java.util.Collection; @@ -117,7 +128,7 @@ public class SensitivityAnalysisExperimentTab extends LabelPropertyTabContributo // Label Composite labelComposite = new Composite(content, SWT.NONE); GridDataFactory.fillDefaults().grab(true, false).span(2, 1).applyTo(labelComposite); - GridLayoutFactory.fillDefaults().numColumns(4).applyTo(labelComposite); + GridLayoutFactory.fillDefaults().numColumns(6).applyTo(labelComposite); Label label = new Label(labelComposite, SWT.NONE); label.setText("Name"); @@ -128,6 +139,16 @@ public class SensitivityAnalysisExperimentTab extends LabelPropertyTabContributo name.setColorProvider(new SysdynBasicColorProvider(new LocalResourceManager(JFaceResources.getResources(), name.getWidget()))); GridDataFactory.fillDefaults().grab(true, false).applyTo(name.getWidget()); + label = new Label(labelComposite, SWT.NONE); + label.setText("Number of runs"); + + TrackedText n = new TrackedText(labelComposite, support, SWT.BORDER); + n.setTextFactory(new IntegerPropertyFactory(SysdynResource.URIs.SensitivityAnalysisExperiment_numberOfValues)); + n.addModifyListener(new IntegerPropertyModifier(context, SysdynResource.URIs.SensitivityAnalysisExperiment_numberOfValues)); + n.setInputValidator(new IntegerValidator()); + n.setColorProvider(new SysdynBasicColorProvider(new LocalResourceManager(JFaceResources.getResources(), n.getWidget()))); + GridDataFactory.fillDefaults().hint(80, SWT.DEFAULT).applyTo(n.getWidget()); + label = new Label(labelComposite, SWT.NONE); label.setText("Seed"); @@ -138,7 +159,6 @@ public class SensitivityAnalysisExperimentTab extends LabelPropertyTabContributo seed.setColorProvider(new SysdynBasicColorProvider(new LocalResourceManager(JFaceResources.getResources(), seed.getWidget()))); GridDataFactory.fillDefaults().hint(80, SWT.DEFAULT).applyTo(seed.getWidget()); - // (Ontology-based) GraphExplorer displaying range axis and variables mapped to those axis explorer = new AxisAndVariablesExplorerComposite(ArrayMap.keys( "displaySelectors", "displayFilter", "treeView").values(false, false, true), site, content, support, SWT.FULL_SELECTION | SWT.BORDER | SWT.SINGLE); @@ -242,18 +262,6 @@ public class SensitivityAnalysisExperimentTab extends LabelPropertyTabContributo // variable.setColorProvider(new SysdynBasicColorProvider(new LocalResourceManager(JFaceResources.getResources(), variable.getWidget()))); // GridDataFactory.fillDefaults().grab(true, false).applyTo(variable.getWidget()); - label = new Label(parameterProperties, SWT.NONE); - label.setText("Number of values:"); - GridDataFactory.fillDefaults().align(SWT.END, SWT.CENTER).applyTo(label); - - - TrackedText numValues = new TrackedText(parameterProperties, parameterSupport, SWT.BORDER); - numValues.setTextFactory(new IntegerPropertyFactory(SysdynResource.URIs.SensitivityAnalysisExperiment_Parameter_numberOfValues)); - numValues.addModifyListener(new IntegerPropertyModifier(context, SysdynResource.URIs.SensitivityAnalysisExperiment_Parameter_numberOfValues)); - numValues.setInputValidator(new IntegerValidator()); - numValues.setColorProvider(new SysdynBasicColorProvider(new LocalResourceManager(JFaceResources.getResources(), numValues.getWidget()))); - GridDataFactory.fillDefaults().align(SWT.BEGINNING, SWT.CENTER).hint(80, SWT.DEFAULT).applyTo(numValues.getWidget()); - label = new Label(parameterProperties, SWT.NONE); label.setText("Distribution:"); GridDataFactory.fillDefaults().align(SWT.END, SWT.CENTER).applyTo(label); @@ -363,7 +371,6 @@ public class SensitivityAnalysisExperimentTab extends LabelPropertyTabContributo Resource parameter = GraphUtils.create2(graph, sr.SensitivityAnalysisExperiment_Parameter, sr.SensitivityAnalysisExperiment_Parameter_propabilityDistribution, distribution, sr.SensitivityAnalysisExperiment_Parameter_variable, "", - sr.SensitivityAnalysisExperiment_Parameter_numberOfValues, 10, L0.PartOf, input); Resource parameterList = graph.getPossibleObject(input, sr.SensitivityAnalysisExperiment_parameterList); diff --git a/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/widgets/sensitivity/ParameterLabelRule.java b/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/widgets/sensitivity/ParameterLabelRule.java index 87984f7e..fd67dc8b 100644 --- a/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/widgets/sensitivity/ParameterLabelRule.java +++ b/org.simantics.sysdyn.ui/src/org/simantics/sysdyn/ui/properties/widgets/sensitivity/ParameterLabelRule.java @@ -8,6 +8,7 @@ * * Contributors: * Semantum Oy - initial API and implementation + * VTT Technical Research Centre of Finland *******************************************************************************/ package org.simantics.sysdyn.ui.properties.widgets.sensitivity; @@ -33,7 +34,6 @@ public class ParameterLabelRule implements LabelRule { HashMap result = new HashMap(); String variable = graph.getPossibleRelatedValue((Resource)content, SysdynResource.getInstance(graph).SensitivityAnalysisExperiment_Parameter_variable); - Integer n = graph.getPossibleRelatedValue((Resource)content, SysdynResource.getInstance(graph).SensitivityAnalysisExperiment_Parameter_numberOfValues); StringBuilder sb = new StringBuilder(); @@ -42,15 +42,6 @@ public class ParameterLabelRule implements LabelRule { else sb.append("No variable"); - sb.append(" (n="); - - if(n != null) - sb.append(n); - else - sb.append("undefined"); - - sb.append(")"); - result.put(ColumnKeys.SINGLE, sb.toString()) ; return result; diff --git a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/IDistribution.java b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/IDistribution.java index 8a6b8c31..760ed446 100644 --- a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/IDistribution.java +++ b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/IDistribution.java @@ -8,26 +8,20 @@ * * Contributors: * Semantum Oy - initial API and implementation + * VTT Technical Research Centre of Finland *******************************************************************************/ package org.simantics.sysdyn.adapter.distribution; public interface IDistribution { - - /** - * Return next value of the distribution - * @return - */ - public double next(); - + /** + * Map a probability to the random variable. + * When random numbers are used to generate the Monte Carlo simulation parameters, the probability + * is itself selected randomly (between (0,1); at the domain endpoints the value is undefined). * - * @return does the distribution have next value; + * @param randomVariable + * @return Inverse cumulative distribution function at probability, i.e. the random variable at the + * point where the cumulative distribution function yields probability. */ - public boolean hasNext(); - - /** - * Initialize distribution - */ - public void initialize(); - + public double inverseCDF(double probability); } diff --git a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/Interval.java b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/Interval.java index 3def0d33..39fd8760 100644 --- a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/Interval.java +++ b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/Interval.java @@ -8,11 +8,10 @@ * * Contributors: * Semantum Oy - initial API and implementation + * VTT Technical Research Centre of Finland *******************************************************************************/ package org.simantics.sysdyn.adapter.distribution; -import java.util.Random; - import org.simantics.databoard.Bindings; import org.simantics.db.ReadGraph; import org.simantics.db.Resource; @@ -22,13 +21,10 @@ import org.simantics.sysdyn.SysdynResource; public class Interval implements IDistribution { - private Random random; - private double min = 0; private double max = 10; private double intervalLength = 1; private int numberOfValues = 10; - private int seed = 123; private int index = 0; public Interval(ReadGraph graph, Resource distribution) { @@ -43,38 +39,20 @@ public class Interval implements IDistribution { Double maxValue = graph.getPossibleRelatedValue(distribution, SR.Interval_maxValue, Bindings.DOUBLE); if(maxValue != null) this.max = maxValue; - - + Resource parameter = graph.getPossibleObject(distribution, SR.SensitivityAnalysisExperiment_Parameter_propabilityDistribution_Inverse); + + Resource experiment = graph.getPossibleObject(parameter, Layer0.getInstance(graph).PartOf); - Integer numberOfValues = graph.getPossibleRelatedValue(parameter, SR.SensitivityAnalysisExperiment_Parameter_numberOfValues, Bindings.INTEGER); + Integer numberOfValues = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_numberOfValues, Bindings.INTEGER); if(numberOfValues != null) this.numberOfValues = numberOfValues; intervalLength = (maxValue - minValue) / (this.numberOfValues - 1); - Resource experiment = graph.getPossibleObject(parameter, Layer0.getInstance(graph).PartOf); - - Integer seed = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_randomSeed, Bindings.INTEGER); - if(seed != null) - this.seed = seed; - } catch (DatabaseException e) { } - - initialize(); - } - - @Override - public double next() { - double value = min + (intervalLength * index); - index++; - return value; - } - - public Random getRandom() { - return random; } public double getMin() { @@ -89,23 +67,18 @@ public class Interval implements IDistribution { return numberOfValues; } - public int getSeed() { - return seed; - } - public double getIntervalLength() { return intervalLength; } - @Override - public void initialize() { - random = new Random(this.seed); - index = 0; - } - - @Override - public boolean hasNext() { - return index < numberOfValues; - } + @Override + public double inverseCDF(double probability) { + // This is a bit awkward... + double value = min + (intervalLength * index); + index++; + if (index >= numberOfValues) + index = 0; + return value; + } } diff --git a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/NormalDistribution.java b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/NormalDistribution.java index efa2ad1b..7887a62a 100644 --- a/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/NormalDistribution.java +++ b/org.simantics.sysdyn/src/org/simantics/sysdyn/adapter/distribution/NormalDistribution.java @@ -8,37 +8,30 @@ * * Contributors: * Semantum Oy - initial API and implementation + * VTT Technical Research Centre of Finland *******************************************************************************/ package org.simantics.sysdyn.adapter.distribution; -import java.util.Random; - import org.simantics.databoard.Bindings; import org.simantics.db.ReadGraph; import org.simantics.db.Resource; import org.simantics.db.exception.DatabaseException; -import org.simantics.layer0.Layer0; import org.simantics.sysdyn.SysdynResource; public class NormalDistribution implements IDistribution { - - private Random random; - + private Double min = null; private Double max = null; + private double minProbability = 0; + private double maxProbability = 1; private double mean; private double stdDeviation; - private int numberOfValues = 10; - private int seed = 123; - private int index = 0; - - private static int MAX_TRIES = 100; - public NormalDistribution(ReadGraph graph, Resource distribution) { try { + SysdynResource SR = SysdynResource.getInstance(graph); min = graph.getPossibleRelatedValue(distribution, SR.NormalDistribution_minValue, Bindings.DOUBLE); @@ -46,66 +39,17 @@ public class NormalDistribution implements IDistribution { stdDeviation = graph.getPossibleRelatedValue(distribution, SR.NormalDistribution_stdDeviation, Bindings.DOUBLE); mean = graph.getPossibleRelatedValue(distribution, SR.NormalDistribution_mean, Bindings.DOUBLE); - - Resource parameter = graph.getPossibleObject(distribution, SR.SensitivityAnalysisExperiment_Parameter_propabilityDistribution_Inverse); + // Determine the max and min probabilities. + if (min != null) + minProbability = cdf(min, mean, stdDeviation, false); + if (max != null) + maxProbability = cdf(max, mean, stdDeviation, false); - Integer numberOfValues = graph.getPossibleRelatedValue(parameter, SR.SensitivityAnalysisExperiment_Parameter_numberOfValues, Bindings.INTEGER); - if(numberOfValues != null) - this.numberOfValues = numberOfValues; - - - - Resource experiment = graph.getPossibleObject(parameter, Layer0.getInstance(graph).PartOf); - - Integer seed = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_randomSeed, Bindings.INTEGER); - if(seed != null) - this.seed = seed; - - - random = new Random(this.seed); - } catch (DatabaseException e) { } } - @Override - public double next() { - - /* Try to find value between min and max, if they are defined. - * If value cannot be found after MAX_TRIES number of tries, - * return the last try. - */ - double value = 0; - for(int i = 0; i < MAX_TRIES; i++) { - value = random.nextGaussian() * stdDeviation + mean; - if(accept(value)) - break; - } - - index++; - return value; - } - - - private boolean accept(double value) { - if(min != null) { - if(value < min) - return false; - } - - if(max != null) { - if(value > max) - return false; - } - - return true; - } - - public Random getRandom() { - return random; - } - public double getMin() { return min; } @@ -113,25 +57,348 @@ public class NormalDistribution implements IDistribution { public double getMax() { return max; } + + @Override + public double inverseCDF(double probability) { + // Map probability to [min,max] + double mappedProbability = minProbability + (maxProbability - minProbability) * probability; + // Calculate the inverse CDF + return quantile(mappedProbability, mean, stdDeviation); + } - public int getNumberOfValues() { - return numberOfValues; - } - - public int getSeed() { - return seed; + //The following is from BEAST-MCMC (https://code.google.com/p/beast-mcmc/) + + /* + * NormalDistribution.java + * + * Copyright (c) 2002-2011 Alexei Drummond, Andrew Rambaut and Marc Suchard + * + * This file is part of BEAST. + * See the NOTICE file distributed with this work for additional + * information regarding copyright ownership and licensing. + * + * BEAST is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * BEAST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with BEAST; if not, write to the + * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, + * Boston, MA 02110-1301 USA + */ + + /** + * quantiles (=inverse cumulative density function) + * + * @param z argument + * @param m mean + * @param sd standard deviation + * @return icdf at z + */ + public static double quantile(double z, double m, double sd) { + return m + Math.sqrt(2.0) * sd * inverseErf(2.0 * z - 1.0); } + /** A more accurate and faster implementation of the cdf (taken from function pnorm in the R statistical language) + * This implementation has discrepancies depending on the programming language and system architecture + * In Java, returned values become zero once z reaches -37.5193 exactly on the machine tested + * In the other implementation, the returned value 0 at about z = -8 + * In C, this 0 value is reached approximately z = -37.51938 + * + * Will later need to be optimised for BEAST + * + * @param x argument + * @param mu mean + * @param sigma standard deviation + * @param log_p is p logged + * @return cdf at x + */ + public static double cdf(double x, double mu, double sigma, boolean log_p) { + boolean i_tail=false; + double p, cp = Double.NaN; - @Override - public void initialize() { - random = new Random(this.seed); - index = 0; - } + if(Double.isNaN(x) || Double.isNaN(mu) || Double.isNaN(sigma)) { + return Double.NaN; + } + if(Double.isInfinite(x) && mu == x) { /* x-mu is NaN */ + return Double.NaN; + } + if (sigma <= 0) { + if(sigma < 0) { + return Double.NaN; + } + return (x < mu) ? 0.0 : 1.0; + } + p = (x - mu) / sigma; + if(Double.isInfinite(p)) { + return (x < mu) ? 0.0 : 1.0; + } + x = p; + if(Double.isNaN(x)) { + return Double.NaN; + } + + double xden, xnum, temp, del, eps, xsq, y; + int i; + boolean lower, upper; + eps = DBL_EPSILON * 0.5; + lower = !i_tail; + upper = i_tail; + + y = Math.abs(x); + if (y <= 0.67448975) { /* Normal.quantile(3/4, 1, 0) = 0.67448975 */ + if (y > eps) { + xsq = x * x; + xnum = a[4] * xsq; + xden = xsq; + for (i = 0; i < 3; i++) { + xnum = (xnum + a[i]) * xsq; + xden = (xden + b[i]) * xsq; + } + } + else { + xnum = xden = 0.0; + } + temp = x * (xnum + a[3]) / (xden + b[3]); + if(lower) { + p = 0.5 + temp; + } + if(upper) { + cp = 0.5 - temp; + } + if(log_p) { + if(lower) { + p = Math.log(p); + } + if(upper) { + cp = Math.log(cp); + } + } + } + + + else if (y <= M_SQRT_32) { + /* Evaluate pnorm for 0.67448975 = Normal.quantile(3/4, 1, 0) < |x| <= sqrt(32) ~= 5.657 */ + + xnum = c[8] * y; + xden = y; + for (i = 0; i < 7; i++) { + xnum = (xnum + c[i]) * y; + xden = (xden + d[i]) * y; + } + temp = (xnum + c[7]) / (xden + d[7]); + + //do_del(y); + //swap_tail; + //#define do_del(X) \ + xsq = ((int) (y * CUTOFF)) * 1.0 / CUTOFF; + del = (y - xsq) * (y + xsq); + if(log_p) { + p = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp); + if((lower && x > 0.0) || (upper && x <= 0.0)) { + cp = Math.log(1.0-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp); + } + } + else { + p = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp; + cp = 1.0 - p; + } + //#define swap_tail \ + if (x > 0.0) { + temp = p; + if(lower) { + p = cp; + } + cp = temp; + } + } + /* else |x| > sqrt(32) = 5.657 : + * the next two case differentiations were really for lower=T, log=F + * Particularly *not* for log_p ! + * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50 + * Note that we do want symmetry(0), lower/upper -> hence use y + */ + else if(log_p || (lower && -37.5193 < x && x < 8.2924) + || (upper && -8.2924 < x && x < 37.5193)) { + + /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */ + xsq = 1.0 / (x * x); + xnum = p_[5] * xsq; + xden = xsq; + for (i = 0; i < 4; i++) { + xnum = (xnum + p_[i]) * xsq; + xden = (xden + q[i]) * xsq; + } + temp = xsq * (xnum + p_[4]) / (xden + q[4]); + temp = (M_1_SQRT_2PI - temp) / y; + + //do_del(x); + xsq = ((int) (x * CUTOFF)) * 1.0 / CUTOFF; + del = (x - xsq) * (x + xsq); + if(log_p) { + p = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp); + if((lower && x > 0.0) || (upper && x <= 0.0)) { + cp = Math.log(1.0-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp); + } + } + else { + p = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp; + cp = 1.0 - p; + } + //swap_tail; + if (x > 0.0) { + temp = p; + if(lower) { + p = cp; + } + cp = temp; + } + } + else { /* no log_p , large x such that probs are 0 or 1 */ + if(x > 0) { + p = 1.0; + cp = 0.0; + } + else { + p = 0.0; + cp = 1.0; + } + } + return p; - @Override - public boolean hasNext() { - return index < numberOfValues; } + // Private + + protected double m, sd; + + private static final double[] a = { + 2.2352520354606839287, + 161.02823106855587881, + 1067.6894854603709582, + 18154.981253343561249, + 0.065682337918207449113 + }; + private static final double[] b = { + 47.20258190468824187, + 976.09855173777669322, + 10260.932208618978205, + 45507.789335026729956 + }; + private static final double[] c = { + 0.39894151208813466764, + 8.8831497943883759412, + 93.506656132177855979, + 597.27027639480026226, + 2494.5375852903726711, + 6848.1904505362823326, + 11602.651437647350124, + 9842.7148383839780218, + 1.0765576773720192317e-8 + }; + private static final double[] d = { + 22.266688044328115691, + 235.38790178262499861, + 1519.377599407554805, + 6485.558298266760755, + 18615.571640885098091, + 34900.952721145977266, + 38912.003286093271411, + 19685.429676859990727 + }; + private static final double[] p_ = { + 0.21589853405795699, + 0.1274011611602473639, + 0.022235277870649807, + 0.001421619193227893466, + 2.9112874951168792e-5, + 0.02307344176494017303 + }; + private static final double[] q = { + 1.28426009614491121, + 0.468238212480865118, + 0.0659881378689285515, + 0.00378239633202758244, + 7.29751555083966205e-5 + }; + + private static final int CUTOFF = 16; /* Cutoff allowing exact "*" and "/" */ + + private static final double M_SQRT_32 = 5.656854249492380195206754896838; /* The square root of 32 */ + private static final double M_1_SQRT_2PI = 0.398942280401432677939946059934; + private static final double DBL_EPSILON = 2.2204460492503131e-016; + + /* + * ErrorFunction.java + * + * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut + * + * This file is part of BEAST. + * See the NOTICE file distributed with this work for additional + * information regarding copyright ownership and licensing. + * + * BEAST is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * BEAST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with BEAST; if not, write to the + * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, + * Boston, MA 02110-1301 USA + */ + + /** + * inverse error function + * + * @param z argument + * + * @return function value + */ + public static double inverseErf(double z) + { + return pointNormal(0.5*z+0.5)/Math.sqrt(2.0); + } + + // Private + + // Returns z so that Prob{x next() { + ArrayList randoms = new ArrayList(); + for (int i = 0; i < dimension; ++i) + randoms.add(random.nextDouble()); + return randoms; + } + + public void initialize() { + random = new Random(this.seed); + } + +} diff --git a/org.simantics.sysdyn/src/org/simantics/sysdyn/manager/SysdynSensitivityAnalysisExperiment.java b/org.simantics.sysdyn/src/org/simantics/sysdyn/manager/SysdynSensitivityAnalysisExperiment.java index 432ff044..0273ec9f 100644 --- a/org.simantics.sysdyn/src/org/simantics/sysdyn/manager/SysdynSensitivityAnalysisExperiment.java +++ b/org.simantics.sysdyn/src/org/simantics/sysdyn/manager/SysdynSensitivityAnalysisExperiment.java @@ -44,32 +44,57 @@ import org.simantics.sysdyn.adapter.SensitivityExperimentParameter; */ public class SysdynSensitivityAnalysisExperiment extends SysdynExperiment { - - private ArrayList results = null; - + private int seed = 124; + private int numberOfRuns = 0; + public SysdynSensitivityAnalysisExperiment(Resource experiment, Resource model) { super(experiment, model); } - private void findValuesAndRun(List parameters, HashMap values, HashMap experimentParameters) { - findValuesAndRun(parameters, 0, values, experimentParameters); - } - - private void findValuesAndRun(List parameters, int index, HashMap values, HashMap experimentParameters) { - SensitivityExperimentParameter p = parameters.get(index); + ArrayList> parameterMaps = new ArrayList>(); + int parametersSize = parameters.size(); + + // Get seed + try { + session.syncRequest(new ReadRequest() { + @Override + public void run(ReadGraph graph) throws DatabaseException { + SysdynResource SR = SysdynResource.getInstance(graph); + seed = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_randomSeed, Bindings.INTEGER); + numberOfRuns = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_numberOfValues, Bindings.INTEGER); + //refreshRate = numberOfRuns / 20 + 1; + refreshRate = 1; + } + }); + } catch (DatabaseException e) { + e.printStackTrace(); + } - while(p.getDistribution().hasNext()) { - values.put(p.getFullName(), Double.toString(p.getDistribution().next())); - if((index + 1) < parameters.size()) { - findValuesAndRun(parameters, index + 1, values, experimentParameters); - } else { - // Run with values - runSensitivityRun(values, experimentParameters); - } - } - p.getDistribution().initialize(); + // Generate values for parameters (randomly) + RandomGenerator randomGenerator = new RandomGenerator(seed, parametersSize); + + // Determine the changed paramters and their values for each individual run + for (int i = 0; i < numberOfRuns; ++i) { + ArrayList randoms = randomGenerator.next();// multi-dimensional random + HashMap parameterMap = new HashMap(parametersSize); + for (int j = 0; j < parametersSize; ++j) { + SensitivityExperimentParameter p = parameters.get(j); + // Get the random value + double value = p.getDistribution().inverseCDF(randoms.get(j)); + // Add the new parameter-value-pair + parameterMap.put(p.getFullName(), Double.toString(value)); + } + // Add the complete list of one individual run parameter-value-pairs. + parameterMaps.add(parameterMap); + } + + // Set the parameters for each individual experiment and run them. + for (HashMap parameterSet : parameterMaps) { + values.putAll(parameterSet); + runSensitivityRun(values, experimentParameters); + } } private void runSensitivityRun(HashMap values, HashMap experimentParameters) { @@ -97,7 +122,8 @@ public class SysdynSensitivityAnalysisExperiment extends SysdynExperiment { ModelicaManager.printProcessOutput(process, monitor); File resFile = new File(experimentParameters.get(ModelicaManager.RESULT_FILE_NAME)); - Thread resultThread = getResultThread(resFile, experimentParameters, monitor, progressMonitor, currentRun % refreshRate == 0); + boolean updateMonitors = currentRun % refreshRate == 0 || currentRun == numberOfRuns - 1; + Thread resultThread = getResultThread(resFile, experimentParameters, monitor, progressMonitor, updateMonitors); resultThread.run(); process = null; @@ -111,7 +137,6 @@ public class SysdynSensitivityAnalysisExperiment extends SysdynExperiment { private IModelicaMonitor monitor; private IProgressMonitor progressMonitor; // private HashMap experimentParameters; - private int numberOfRuns = 0; private int currentRun = 0; private int refreshRate = 1; private List parameters = new ArrayList(); @@ -293,21 +318,7 @@ public class SysdynSensitivityAnalysisExperiment extends SysdynExperiment { @Override public Integer perform(ReadGraph graph) throws DatabaseException { SysdynResource SR = SysdynResource.getInstance(graph); - - Resource parameterListResource = graph.getPossibleObject(experiment, SR.SensitivityAnalysisExperiment_parameterList); - List parameterResources = ListUtils.toList(graph, parameterListResource); - - Integer numberOfIterations = null; - - for(Resource parameter : parameterResources) { - Integer n = graph.getPossibleRelatedValue(parameter, SR.SensitivityAnalysisExperiment_Parameter_numberOfValues, Bindings.INTEGER); - if(n != null) { - if(numberOfIterations == null) - numberOfIterations = n; - else - numberOfIterations = numberOfIterations * n; - } - } + Integer numberOfIterations = graph.getPossibleRelatedValue(experiment, SR.SensitivityAnalysisExperiment_numberOfValues, Bindings.INTEGER); return numberOfIterations; }