Author: echatellier Date: 2013-11-23 22:50:19 +0100 (Sat, 23 Nov 2013) New Revision: 3844 Url: http://forge.codelutin.com/projects/isis-fish/repository/revisions/3844 Log: Manage distribution for continuous domains Modified: branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/Distribution.java branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/domain/ContinuousDomain.java branches/4.0.1/src/test/java/fr/ifremer/isisfish/simulator/sensitivity/SensitivityAnalysisTest.java Modified: branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/Distribution.java =================================================================== --- branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/Distribution.java 2013-11-23 21:48:39 UTC (rev 3843) +++ branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/Distribution.java 2013-11-23 21:50:19 UTC (rev 3844) @@ -24,6 +24,16 @@ package fr.ifremer.isisfish.simulator.sensitivity; import static org.nuiton.i18n.I18n._; +import jdistlib.Binomial; +import jdistlib.Cauchy; +import jdistlib.ChiSquare; +import jdistlib.Exponential; +import jdistlib.F; +import jdistlib.Gamma; +import jdistlib.Geometric; +import jdistlib.HyperGeometric; +import jdistlib.Normal; +import jdistlib.Uniform; /** * Distribution used in sensitivity analysis (for continuous domain). @@ -39,59 +49,73 @@ QUNIFPC(_("QUnif %"), "runif", + Uniform.class, new DistributionParam("reference", _("Reference value"), true), - new DistributionParam("coefficient", _("Coefficient"))), + new DistributionParam("coefficient", _("Coefficient"), false)), QUNIFMM(_("QUnif Min/Max"), "runif", + Uniform.class, new DistributionParam("min", _("Minimum value"), true), new DistributionParam("max", _("Maximum value"), true)), DBINOM(_("Dbimom"), "rbinom", + Binomial.class, new DistributionParam("size", _("number of trials (zero or more)")), new DistributionParam("prob", _("probability of success on each trial"))), DCAUCHY(_("dcauchy"), "rcauchy", + Cauchy.class, new DistributionParam("location", _("location and scale parameters. (location = 0, scale = 1)")), new DistributionParam("scale", _("location and scale parameters. (location = 0, scale = 1)"))), DCHISQ(_("dchisq"), "rchisq", + ChiSquare.class, new DistributionParam("df", _("degrees of freedom (non-negative, but can be non-integer)")), new DistributionParam("ncp", _("non-centrality parameter (non-negative)."))), DEXP(_("dexp"), "rexp", + Exponential.class, new DistributionParam("rate", _("vector of rates"))), DF(_("df"), "rf", + F.class, new DistributionParam("df1", _("degrees of freedom. 'Inf' is allowed")), new DistributionParam("df2", _("degrees of freedom. 'Inf' is allowed")), new DistributionParam("ncp", _("non-centrality parameter. If omitted the central F is assumed"))), DGAMMA(_("dgamma"), "rgamma", + Gamma.class, new DistributionParam("shape", _("shape and scale parameters. Must be positive, 'scale' strictly")), new DistributionParam("rate", _("an alternative way to specify the scale.")), new DistributionParam("scale", _("shape and scale parameters. Must be positive, 'scale' strictly"))), DGEOM(_("dgeom"), "rgeom", + Geometric.class, new DistributionParam("prob", _("probability of success in each trial. '0 < prob <= 1'"))), DHYPER(_("dhyper"), "rhyper", + HyperGeometric.class, new DistributionParam("m", _("the number of white balls in the urn.")), new DistributionParam("n", _("the number of black balls in the urn.")), new DistributionParam("k", _("the number of balls drawn from the urn."))), DLNORM(_("dlnorm"), "rlnorm", - new DistributionParam("meanlog", _("mean and standard deviation of the distribution on the log scale with default values of '0' and '1' respectively.")), - new DistributionParam("sdlog", _("mean and standard deviation of the distribution on the log scale with default values of '0' and '1' respectively."))); + Normal.class, + new DistributionParam("meanlog", _("mean and standard deviation of the distribution on the log scale with default values of '0' and '1' respectively."), false), + new DistributionParam("sdlog", _("mean and standard deviation of the distribution on the log scale with default values of '0' and '1' respectively."), false)); protected String description; protected String instruction; + protected Class distClass; + protected DistributionParam[] params; - private Distribution(String description, String instruction, DistributionParam... params) { + private Distribution(String description, String instruction, Class distClass, DistributionParam... params) { this.description = description; this.instruction = instruction; + this.distClass = distClass; this.params = params; } @@ -103,6 +127,10 @@ return instruction; } + public Class getDistClass() { + return distClass; + } + public DistributionParam[] getDistibutionParams() { return params; } @@ -120,7 +148,7 @@ this.description = description; this.originalValue = originalValue; } - + public DistributionParam(String name, String description) { this(name, description, false); } Modified: branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/domain/ContinuousDomain.java =================================================================== --- branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/domain/ContinuousDomain.java 2013-11-23 21:48:39 UTC (rev 3843) +++ branches/4.0.1/src/main/java/fr/ifremer/isisfish/simulator/sensitivity/domain/ContinuousDomain.java 2013-11-23 21:50:19 UTC (rev 3844) @@ -25,19 +25,24 @@ package fr.ifremer.isisfish.simulator.sensitivity.domain; +import java.lang.reflect.InvocationTargetException; +import java.util.ArrayList; import java.util.HashMap; +import java.util.List; import java.util.Map; import java.util.Map.Entry; import java.util.SortedMap; +import org.apache.commons.beanutils.MethodUtils; +import org.apache.commons.logging.Log; +import org.apache.commons.logging.LogFactory; +import org.nuiton.math.matrix.MatrixND; + import fr.ifremer.isisfish.IsisFishRuntimeException; import fr.ifremer.isisfish.simulator.sensitivity.Distribution; import fr.ifremer.isisfish.simulator.sensitivity.Distribution.DistributionParam; import fr.ifremer.isisfish.simulator.sensitivity.Domain; import fr.ifremer.isisfish.simulator.sensitivity.visitor.DomainVisitor; -import jdistlib.Uniform; -import jdistlib.rng.QMersenneTwister; -import jdistlib.rng.QRandomEngine; /** * All the continuous domains are based on distribution definition. @@ -50,6 +55,8 @@ */ public class ContinuousDomain implements Domain { + private static Log log = LogFactory.getLog(ContinuousDomain.class); + /** serialVersionUID. */ private static final long serialVersionUID = -2037768174807839046L; @@ -168,8 +175,58 @@ @Override public Object getValueForIdentifier(Object identifier) { - Double result = Uniform.quantile((Double)identifier,0,1, true,false); + // first quantite args is always identifier value + List<Object> args = new ArrayList<Object>(); + args.add(identifier); + // build distribution param args list + Class distClass = distribution.getDistClass(); + if (distribution == Distribution.QUNIFPC) { + // special management : % to min/max + Object reference = distributionParameters.get(Distribution.QUNIFPC.getDistibutionParams()[0].getName()); + if (reference instanceof MatrixND) { + args.add(0); // min + args.add(1); // max + } else { + double ref = (Double)reference; + double coef = (Double)distributionParameters.get(Distribution.QUNIFPC.getDistibutionParams()[1].getName()); + args.add(ref * (1 - coef)); // min + args.add(ref * (1 + coef)); // max + } + } else if (distribution == Distribution.QUNIFMM) { + Object min = distributionParameters.get(Distribution.QUNIFMM.getDistibutionParams()[0].getName()); + if (min instanceof MatrixND) { + args.add(0); // min + } else { + args.add(min); + } + Object max = distributionParameters.get(Distribution.QUNIFMM.getDistibutionParams()[0].getName()); + if (max instanceof MatrixND) { + args.add(1); // max + } else { + args.add(max); + } + } else { + for (DistributionParam param : distribution.getDistibutionParams()) { + String name = param.getName(); + Object value = distributionParameters.get(name); + args.add(value); + } + } + + // add default value for lower.tail = TRUE, log.p = FALSE + args.add(true); + args.add(false); + + // invoke quantile method + Object result = null; + try { + result = MethodUtils.invokeStaticMethod(distClass, "quantile", args.toArray()); + } catch (NoSuchMethodException | IllegalAccessException | InvocationTargetException ex) { + if (log.isErrorEnabled()) { + log.error("Can't invoke quantile method", ex); + } + } return result; } Modified: branches/4.0.1/src/test/java/fr/ifremer/isisfish/simulator/sensitivity/SensitivityAnalysisTest.java =================================================================== --- branches/4.0.1/src/test/java/fr/ifremer/isisfish/simulator/sensitivity/SensitivityAnalysisTest.java 2013-11-23 21:48:39 UTC (rev 3843) +++ branches/4.0.1/src/test/java/fr/ifremer/isisfish/simulator/sensitivity/SensitivityAnalysisTest.java 2013-11-23 21:50:19 UTC (rev 3844) @@ -28,7 +28,6 @@ import java.io.IOException; import fr.ifremer.isisfish.IsisFish; -import jdistlib.Uniform; import org.apache.commons.io.FileUtils; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; @@ -82,14 +81,16 @@ protected DesignPlan getDesignPlan(Factor... factors) { DesignPlan result = new DesignPlan(); + int factorIndex = 1; for (Factor factor : factors) { + factor.setName("f" + factorIndex++); result.addFactor(factor); } return result; } - protected FactorGroup getFactorGroup(String name, boolean continuous, Factor... factors) { - FactorGroup result = new FactorGroup(name, continuous); + protected FactorGroup getFactorGroup(boolean continuous, Factor... factors) { + FactorGroup result = new FactorGroup("", continuous); for (Factor factor : factors) { result.addFactor(factor); } @@ -175,6 +176,18 @@ return factor1; } + protected void displayScenarios(SensitivityScenarios scenarii) { + if (log.isDebugEnabled()) { + int scIndex = 1; + for (Scenario scenario : scenarii.getScenarios()) { + log.debug("Scenario " + scIndex++); + for (Factor factor : scenario.getFactors()) { + log.debug("Factor " + factor.getName() + " = " + factor.getValue()); + } + } + } + } + /** * Test tmp files/dir. * @@ -199,13 +212,36 @@ // uniform min/max DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); - Assert.assertEquals(9, scenarii.getScenarios().size()); + Assert.assertTrue(scenarii.getScenarios().size() >= 3); // random, can be 9, 12 + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } // uniform pc designPlan = getDesignPlan(getContinuousDoubleUniformPcFactor(), getContinuousDoubleUniformPcFactor()); scenarii = script.compute(designPlan, simulationsDir); - Assert.assertEquals(9, scenarii.getScenarios().size()); + Assert.assertTrue(scenarii.getScenarios().size() >= 3); // random, can be 9, 12 + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 30.0 && value <= 50.0); + } + } + // matrix uniform + designPlan = getDesignPlan(getContinuousMatrixUniformMMFactor(), getContinuousMatrixUniformMMFactor()); + scenarii = script.compute(designPlan, simulationsDir); + Assert.assertTrue(scenarii.getScenarios().size() >= 3); // random, can be 9, 12 + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } + // check fail is discrete designPlan = getDesignPlan(getDiscreteIntFactor(),getContinuousDoubleUniformMMFactor()); try { @@ -230,7 +266,12 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); - + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } } /** @@ -248,11 +289,18 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(40, scenarii.getScenarios().size()); + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } // uniform pc designPlan = getDesignPlan(getContinuousDoubleUniformPcFactor(), getContinuousDoubleUniformPcFactor()); scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(40, scenarii.getScenarios().size()); + displayScenarios(scenarii); // check fail is discrete designPlan = getDesignPlan(getDiscreteIntFactor(),getContinuousDoubleUniformMMFactor()); @@ -280,11 +328,18 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(10, scenarii.getScenarios().size()); + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } // uniform pc designPlan = getDesignPlan(getContinuousDoubleUniformPcFactor(), getContinuousDoubleUniformPcFactor()); scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(10, scenarii.getScenarios().size()); + displayScenarios(scenarii); // check fail is discrete designPlan = getDesignPlan(getDiscreteIntFactor(),getContinuousDoubleUniformMMFactor()); @@ -312,11 +367,18 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(10, scenarii.getScenarios().size()); + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } // uniform pc designPlan = getDesignPlan(getContinuousDoubleUniformPcFactor(), getContinuousDoubleUniformPcFactor()); scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(10, scenarii.getScenarios().size()); + displayScenarios(scenarii); // check fail is discrete designPlan = getDesignPlan(getDiscreteIntFactor(),getContinuousDoubleUniformMMFactor()); @@ -342,7 +404,13 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); - + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } + } /** @@ -354,13 +422,15 @@ @Test public void testRegularFractions() throws IsisFishException, SensitivityException { SensitivityAnalysis script = SensitivityAnalysisStorage.getSensitivityAnalysis("RegularFractions").getNewSensitivityAnalysisInstance(); - SensitivityAnalysisStorage.setParameterValue(script,"pathToFunction", IsisFish.config.getDatabaseDirectory().getAbsolutePath() + "/sensitivityanalysis"); + // this script need pathToFunction parameter to find R script in current directory + SensitivityAnalysisStorage.setParameterValue(script,"pathToFunction", + IsisFish.config.getDatabaseDirectory().getAbsolutePath() + File.separator + "sensitivityanalysis"); + Assert.assertNotNull("RegularFractions script not found in test data", script); DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); - - System.out.println(); + displayScenarios(scenarii); } @@ -379,11 +449,23 @@ DesignPlan designPlan = getDesignPlan(getContinuousDoubleUniformMMFactor(), getContinuousDoubleUniformMMFactor()); SensitivityScenarios scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(80, scenarii.getScenarios().size()); + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 2.0); + } + } // sobol avec cauchy/chisq designPlan = getDesignPlan(getContinuousDoubleCauchyFactor(), getContinuousDoubleChisqFactor()); scenarii = script.compute(designPlan, simulationsDir); Assert.assertEquals(80, scenarii.getScenarios().size()); + for (Scenario scenario : scenarii.getScenarios()) { + for (Factor factor : scenario.getFactors()) { + double value = (Double)factor.getValue(); + Assert.assertTrue(value >= 0.0 && value <= 100.0); + } + } // check fail is discrete designPlan = getDesignPlan(getDiscreteIntFactor(),getContinuousDoubleUniformMMFactor());