Author: jcouteau Date: 2009-04-07 15:09:00 +0000 (Tue, 07 Apr 2009) New Revision: 97 Added: trunk/sensitivity/SensitivityCalculatorRFast.java trunk/sensitivity/SensitivityCalculatorRFrF2.java trunk/sensitivity/SensitivityCalculatorRSobol.java Modified: trunk/sensitivity/SensitivityCalculatorRMorris.java trunk/sensitivity/SensitivityCalculatorROptimumLHS.java trunk/sensitivity/SensitivityCalculatorRRandomLHS.java Log: All the Sensitivity analysis functional, with documentation and default values. Added: trunk/sensitivity/SensitivityCalculatorRFast.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRFast.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorRFast.java 2009-04-07 15:09:00 UTC (rev 97) @@ -0,0 +1,341 @@ +/* *##% + * Copyright (C) 2006 - 2009 Code Lutin + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *##%*/ + +package sensitivity; + +import java.io.File; +import java.io.Serializable; +import java.util.List; +import java.util.Vector; + +import org.apache.commons.logging.Log; +import org.apache.commons.logging.LogFactory; +import org.codelutin.j2r.REngine; +import org.codelutin.j2r.RException; +import org.codelutin.j2r.RProxy; +import org.codelutin.math.matrix.MatrixFactory; +import org.codelutin.math.matrix.MatrixND; +import org.codelutin.util.FileUtil; +import org.rosuda.JRI.REXP; + +import fr.ifremer.isisfish.datastore.SimulationStorage; +import fr.ifremer.isisfish.simulator.sensitivity.DesignPlan; +import fr.ifremer.isisfish.simulator.sensitivity.Factor; +import fr.ifremer.isisfish.simulator.sensitivity.Scenario; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityCalculator; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityScenarios; +import fr.ifremer.isisfish.simulator.sensitivity.domain.ContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.DiscreteDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; +import fr.ifremer.isisfish.util.Doc; + +/** + * Implementation of Fast method using R. + * + * @author jcouteau + * @version $Revision: 89 $ + * + * Last update : $Date: 2009-03-25 13:45:16 +0100 (mer., 25 mars 2009) $ By : + * $Author: jcouteau $ + */ + +public class SensitivityCalculatorRFast implements SensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + private static Log log = LogFactory + .getLog(SensitivityCalculatorRFast.class); + + @Doc("an integer giving the sample size, i.e. the length of the discretization of the s-space (see Cukier et al.). (default=17)") + public int param_n; + + @Doc("an integer specifying the interference parameter, i.e. the number of harmonics to sum in the Fourier series decomposition (see Cukier et al.). (default=1)") + public int param_M; + + //public int[] param_omega; + + /** + * Retourne vrai si le calculateur sait gerer la cardinalité des facteurs + * continue. + * + * @return <tt>true</tt> s'il sait la gerer + */ + public boolean canManageCardinality() { + return true; + } + + public SensitivityScenarios compute(DesignPlan plan, File outputdirectory) { + + if(param_n==0){ + param_n=17; + } + + if(param_M==0){ + param_M=1; + } + + double[] dataframe = new double[0]; + int nbExperiments = 0; + int factorNumber = plan.getFactors().size(); + List<Factor<? extends Serializable>> factors = plan.getFactors(); + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + //Test all factors, if one is discrete, return null + for (int i = 0; i < factorNumber; i++) { + if (factors.get(i).getDomain() instanceof DiscreteDomain) { + return null; + } + } + + String rInstruction = "a<-fast99(model=NULL,factors=c(" + + factors.size(); + + /* // Creating the factors vector. + for (int i = 0; i < factorNumber; i++) { + if (i != (factorNumber - 1)) { + rInstruction = rInstruction + "\"" + factors.get(i).getName() + + "\","; + } else { + rInstruction = rInstruction + "\"" + factors.get(i).getName() + + "\""; + } + }*/ + + rInstruction = rInstruction + "), n=" + param_n + ",M=" + param_M; + + rInstruction = rInstruction + + ", q = \"qunif\", q.arg=list(min=0,max=1))"; + + REngine engine = new RProxy(); + try { + // Load sensitivity package into R (if package already loaded, + // nothing happens. + engine.voidEval("library(sensitivity)"); + log.info("Message sent to R" + "library(sensitivity)"); + + // Run sensitivity analysis + engine.voidEval(rInstruction); + log.info("Message sent to R" + rInstruction); + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R" + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + // Export the fast99 object for the second run in a .fast99 file + engine.voidEval("dput(a,file=\".fast99\")"); + log.info("Message sent to R" + "dput(a,file=\".fast99\")"); + + // Get back experiment plan + + Vector<REXP> dataframeVector = (Vector<REXP>) engine.eval("a$X"); + dataframe = new double[dataframeVector.size() + * dataframeVector.get(0).asDoubleArray().length]; + for (int i = 0; i < dataframeVector.size(); i++) { + for (int j = 0; j < dataframeVector.get(0).asDoubleArray().length; j++) { + dataframe[i * dataframeVector.get(0).asDoubleArray().length + + j] = dataframeVector.get(i).asDoubleArray()[j]; + } + } + log.info("Message sent to R" + "a$X"); + + if (log.isDebugEnabled()) { + log.debug("rInstruction = " + rInstruction); + } + + nbExperiments = dataframe.length / factorNumber; + + } catch (RException e) { + e.printStackTrace(); + // Error while retrieving scenario + } + + // Transform the result from R in a matrix + MatrixND fast = MatrixFactory.getInstance().create(dataframe, + new int[] { factorNumber, nbExperiments }); + + // Setting up the scenarios. + for (int j = 0; j < nbExperiments; j++) { + Scenario experimentScenario = new Scenario(); + for (int i = 0; i < factorNumber; i++) { + Factor<? extends Serializable> factor = plan.getFactors() + .get(i); + if ((factor.getDomain() instanceof MatrixContinuousDomain) + || (factor.getDomain() instanceof EquationContinuousDomain)) { + factor.setValueForIdentifier(Double.valueOf( + fast.getValue(new int[] { i, j })).toString()); + } else { + Double value = (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound() + + ((Double) ((ContinuousDomain<?>) factor + .getDomain()).getMaxBound() - (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound()) + * fast.getValue(new int[] { i, j }); + factor.setValueForIdentifier(value); + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + return thisExperiment; + + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) { + + if(param_n==0){ + param_n=17; + } + + if(param_M==0){ + param_M=1; + } + + REngine engine = new RProxy(); + try { + + // Call R + // Load sensitivity package into R (if package already loaded, + // nothing happens. + engine.voidEval("library(sensitivity)"); + log.info("Message sent to R : " + "library(sensitivity)"); + + //Set the working directory (for import and exports) + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : " + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + // Export the morris object for the second run in a .morris file + engine.voidEval("a<-dget(\".fast99\")"); + log.info("Message sent to R : " + "a<-dget(\".fast99\")"); + + /*int scenariosNumber = sensitivityScenarios.getScenarios().size();*/ + int scenariosNumber = (Integer)engine + .eval("length(a$X[,1])"); + log.info("Message sent to R : " + "length(a$X[,1])"); + + int sensitivityNumber = simulationStorages.get(0).getParameter() + .getSensitivityExport().size(); + + for (int k = 0; k < sensitivityNumber; k++) { + + // Creates the R expression to import results in R + String rInstruction = "results<-c("; + for (int l = 0; l < scenariosNumber; l++) { + File importFile = new File(simulationStorages.get(l) + .getDirectory().toString() + + File.separator + + SimulationStorage.RESULT_EXPORT_DIRECTORY, + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExtensionFilename()); + String simulResult = FileUtil.readAsString(importFile); + double simulationResult = Double.valueOf(simulResult); + if (l < simulationStorages.size() - 1) { + rInstruction = rInstruction + simulationResult + ","; + } else { + rInstruction = rInstruction + simulationResult; + } + } + rInstruction = rInstruction + ")"; + + log.info("Message sent to R : " + rInstruction); + + // Send the simulation results + engine.voidEval(rInstruction); + + //Compute results + engine.voidEval("tell(a,y=results)"); + log.info("Message sent to R : " + "tell(a,y=results)"); + + //Create the data.frame of scenarios and results for export purpose + engine.voidEval("dfresults=data.frame(a$X,results)"); + log.info("Message sent to R : " + + "dfresults=data.frame(a$X,results)"); + + //Set working directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Export V + engine.voidEval("write.csv(a$V,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + "_V.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$V,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_V.csv\")"); + + //Export D1 + engine.voidEval("write.csv(a$D1,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_D1.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$D1,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_D1.csv\")"); + + //Export Dt + engine.voidEval("write.csv(a$Dt,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Dt.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$Dt,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Dt.csv\")"); + + //Export results + engine.voidEval("write.csv(dfresults,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + log.info("Message sent to R : " + + "write.csv(results,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + //FIXME export through java to enable export when using Rserve + + } + } catch (Exception e) { + e.printStackTrace(); + } + } + + public String getDescription() { + return "Implementation of FAST method using R"; + } + +} Added: trunk/sensitivity/SensitivityCalculatorRFrF2.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRFrF2.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorRFrF2.java 2009-04-07 15:09:00 UTC (rev 97) @@ -0,0 +1,447 @@ +/* *##% + * Copyright (C) 2006 - 2009 Code Lutin + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *##%*/ + +package sensitivity; + +import java.io.File; +import java.io.Serializable; +import java.util.List; + +import org.apache.commons.logging.Log; +import org.apache.commons.logging.LogFactory; +import org.codelutin.j2r.REngine; +import org.codelutin.j2r.RException; +import org.codelutin.j2r.RProxy; +import org.codelutin.math.matrix.MatrixFactory; +import org.codelutin.math.matrix.MatrixND; +import org.codelutin.util.FileUtil; + +import fr.ifremer.isisfish.datastore.SimulationStorage; +import fr.ifremer.isisfish.simulator.sensitivity.DesignPlan; +import fr.ifremer.isisfish.simulator.sensitivity.Factor; +import fr.ifremer.isisfish.simulator.sensitivity.Scenario; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityCalculator; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityScenarios; +import fr.ifremer.isisfish.simulator.sensitivity.domain.ContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.DiscreteDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; +import fr.ifremer.isisfish.util.Doc; + +/** + * Implementation of FrF2 method using R. + * + * @author jcouteau + * @version $Revision: 94 $ + * + * Last update : $Date: 2009-04-03 13:13:35 +0200 (ven., 03 avr. 2009) $ By : + * $Author: chatellier $ + */ +@Doc("FrF2 method") +public class SensitivityCalculatorRFrF2 implements SensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + private static Log log = LogFactory + .getLog(SensitivityCalculatorRFrF2.class); + + @Doc("is the arabic numeral for the requested resolution of the design( 3 <= resolution <= 5). (if resolution=3, model = sum(Xi), if resolution = 4 or 5, model = sum(Xi)+sum(XiXj)") + public int param_resolution; + + /** + * Retourne vrai si le calculateur sait gerer la cardinalité des facteurs + * continue. + * + * @return <tt>true</tt> s'il sait la gerer + */ + public boolean canManageCardinality() { + return true; + } + + public SensitivityScenarios compute(DesignPlan plan, File outputdirectory) { + + double[] dataframe = new double[0]; + int nbExperiments = 0; + int factorNumber = plan.getFactors().size(); + List<Factor<? extends Serializable>> factors = plan.getFactors(); + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + //Test all factors, if one is discrete, return null + for (int i = 0; i < factorNumber; i++) { + if (factors.get(i).getDomain() instanceof DiscreteDomain) { + return null; + } + } + + REngine engine = new RProxy(); + + try { + engine.voidEval("library(FrF2)"); + log.info("Message sent to R : " + "library(FrF2)"); + + //Set the working directory (for import and exports) + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : " + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + String rInstruction = "a<-FrF2(nfactors=" + factorNumber + + ",resolution=" + param_resolution + ")"; + + engine.voidEval(rInstruction); + log.info("Message sent to R : " + rInstruction); + + // Export the FrF2 object for the second run in a .FrF2 file + engine.voidEval("dput(a,file=\".FrF2\")"); + log.info("Message sent to R" + "dput(a,file=\".FrF2\")"); + + // Get back experiment plan + dataframe = (double[]) engine.eval("a$desnum"); + log.info("Message sent to R : " + "a$desnum"); + + } catch (RException e) { + e.printStackTrace(); + } + + nbExperiments = dataframe.length / factorNumber; + + // Transform the result from R in a matrix + MatrixND frf2 = MatrixFactory.getInstance().create(dataframe, + new int[] { factorNumber, nbExperiments }); + + // Setting up the scenarios. + for (int j = 0; j < nbExperiments; j++) { + Scenario experimentScenario = new Scenario(); + for (int i = 0; i < factorNumber; i++) { + Factor<? extends Serializable> factor = plan.getFactors() + .get(i); + if ((factor.getDomain() instanceof MatrixContinuousDomain) + || (factor.getDomain() instanceof EquationContinuousDomain)) { + if (frf2.getValue(new int[] { i, j }) == -1) { + factor.setValueForIdentifier(Double.valueOf(0) + .toString()); + } + + if (frf2.getValue(new int[] { i, j }) == 1) { + factor.setValueForIdentifier(Double.valueOf(1) + .toString()); + } + } else { + if (frf2.getValue(new int[] { i, j }) == -1) { + factor + .setValueForIdentifier(((ContinuousDomain<?>) factor + .getDomain()).getMinBound()); + } + + if (frf2.getValue(new int[] { i, j }) == 1) { + factor + .setValueForIdentifier(((ContinuousDomain<?>) factor + .getDomain()).getMaxBound()); + } + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + try { + //Create the factors vectors + for (int j = 0; j < factorNumber; j++) { + Factor<? extends Serializable> factor = thisExperimentScenarios + .get(0).getFactors().get(j); + if (factor.getDomain() instanceof EquationContinuousDomain) { + + String vector = "factor" + j + "<-c("; + for (int i = 0; i < nbExperiments; i++) { + if (i < (nbExperiments - 1)) { + vector = vector + + ((EquationContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue() + ","; + } else { + vector = vector + + ((EquationContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue(); + } + + } + vector = vector + ")"; + engine.voidEval(vector); + log.info("Message sent to R : " + vector); + } else if (factor.getDomain() instanceof MatrixContinuousDomain) { + String vector = "factor" + j + "<-c("; + for (int i = 0; i < nbExperiments; i++) { + if (i < (nbExperiments - 1)) { + vector = vector + + ((MatrixContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue() + ","; + } else { + vector = vector + + ((MatrixContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue(); + } + + } + vector = vector + ")"; + engine.voidEval(vector); + log.info("Message sent to R : " + vector); + } else { + String vector = "factor" + j + "<-c("; + for (int i = 0; i < nbExperiments; i++) { + if (i < (nbExperiments - 1)) { + vector = vector + + thisExperimentScenarios.get(i) + .getFactors().get(j).getValue() + + ","; + } else { + vector = vector + + thisExperimentScenarios.get(i) + .getFactors().get(j).getValue(); + } + + } + vector = vector + ")"; + engine.voidEval(vector); + log.info("Message sent to R : " + vector); + } + } + + //Create the data data.frame from the factors + String data = "data<-data.frame("; + for (int j = 0; j < factorNumber; j++) { + if (j < factorNumber - 1) { + data = data + "factor" + j + "=factor(factor" + j + "),"; + } else { + data = data + "factor" + j + "=factor(factor" + j + "))"; + } + + } + engine.voidEval(data); + log.info("Message sent to R : " + data); + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : " + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + // Export the scenario matrix for the second run in a .FrF2.csv file + engine.voidEval("write.csv(data,file=\".FrF2.csv\")"); + log.info("Message sent to R : " + + "write.csv(data,file=\".FrF2.csv\")"); + + engine.terminate(); + + } catch (RException e) { + e.printStackTrace(); + } + return thisExperiment; + + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) { + + REngine engine = new RProxy(); + try { + + int sensitivityNumber = simulationStorages.get(0).getParameter() + .getSensitivityExport().size(); + + for (int k = 0; k < sensitivityNumber; k++) { + + engine.voidEval("library(FrF2)"); + log.info("Message sent to R : " + "library(FrF2)"); + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Get back the FrF2 object + engine.voidEval("a<-dget(\".FrF2\")"); + log.info("Message sent to R : a<-dget(\".FrF2\")"); + + //Get back the scenarios + engine.voidEval("factors<-read.csv(\".FrF2.csv\")"); + log + .info("Message sent to R : factors<-read.csv(\".FrF2.csv\")"); + + //Get back the factors number + int factorNumber = ((Double) engine.eval("length(factors)-1")) + .intValue(); + + //Create the results vectors + String result = "result<-c("; + for (int l = 0; l < simulationStorages.size(); l++) { + File importFile = new File(simulationStorages.get(l) + .getDirectory().toString() + + File.separator + + SimulationStorage.RESULT_EXPORT_DIRECTORY, + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExtensionFilename()); + String simulResult = FileUtil.readAsString(importFile); + double simulationResult = Double.valueOf(simulResult); + if (l < simulationStorages.size() - 1) { + result = result + simulationResult + ","; + } else { + result = result + simulationResult; + } + } + result = result + ")"; + engine.voidEval(result); + log.info("Message sent to R : " + result); + + //Create the dataforaov data.frame + String dataframe = "dataforaov<-data.frame(factors,result=result)"; + engine.voidEval(dataframe); + log.info("Message sent to R : " + dataframe); + + //get back the resolution + int resolution = (Integer) engine.eval("res.catlg(a$catentry)"); + log.info("Message sent to R : " + "res.catlg(a$catentry)"); + log.info(resolution); + + String aovCall; + + switch (resolution) { + case 3: + //Call aov() + aovCall = "aovresult<-aov(result~"; + for (int j = 0; j < factorNumber; j++) { + if (j < (factorNumber - 1)) { + aovCall = aovCall + "factor" + j + "+"; + } else { + aovCall = aovCall + "factor" + j + + ",data=dataforaov)"; + } + } + engine.voidEval(aovCall); + log.info("Message sent to R : " + aovCall); + break; + case 4: + aovCall = "aovresult<-aov(result~"; + for (int j = 0; j < factorNumber; j++) { + aovCall = aovCall + "factor" + j + "+"; + } + for (int i = 0; i < factorNumber; i++) { + for (int j = 0; j < factorNumber; j++) { + if (i < j) { + aovCall = aovCall + "factor" + i + "*factor" + + j + "+"; + } + } + } + aovCall = aovCall.substring(0, aovCall.length()-1); + aovCall = aovCall + ",data=dataforaov)"; + engine.voidEval(aovCall); + log.info("Message sent to R : " + aovCall); + break; + default: + log.info("pas bien"); + break; + } + + /*Export the results + *Export format is csv, data separated by ',' + *Results Export name is sensitivityExportName_Results.csv + *Sensitivity Indices export name is sensitivityExportName_SensitivityIndices.csv + */ + + //Compute Sum of Squares and Sensitivity indices + engine.voidEval("SoS<-summary(aovresult)[[1]][1:" + + factorNumber + ",2]"); + log.info("Message sent to R : SoS<-summary(aovresult)[[1]][1:" + + factorNumber + ",2]"); + engine + .voidEval("names(SoS)<-dimnames(summary(aovresult)[[1]])[[1]][1:" + + factorNumber + "]"); + log + .info("Message sent to R : names(SoS)<-dimnames(summary(aovresult)[[1]])[[1]][1:" + + factorNumber + "]"); + engine.voidEval("IndSensibilite<-SoS/sum(SoS)"); + log.info("Message sent to R : IndSensibilite<-SoS/sum(SoS)"); + + //Create a data.frame to export sensitivity important results in one file. + engine.voidEval("exportsensitivity=data.frame(SoS[1:" + + factorNumber + "],IndSensibilite[1:" + factorNumber + + "])"); + log + .info("Message sent to R : exportsensitivity=data.frame(SoS[1:" + + factorNumber + + "],IndSensibilite[1:" + + factorNumber + "])"); + engine + .voidEval("names(exportsensitivity)<-c(\"Sum Of Squares\",\"Sensitivity indices\")"); + log + .info("Message sent to R : names(exportsensitivity)<-c(\"Sum Of Squares\",\"Sensitivity indices\")"); + + + + /*Set the export directory + *Export directory is the first simulation export directory. + */ + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Save the results with the scenarios. + engine.voidEval("write.csv(dataforaov,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + log.info("Message sent to R : write.csv(dataforaov,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + + //Save the sensitivity indices + engine.voidEval("write.csv(exportsensitivity,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + "_SensitivityIndices.csv\")"); + log.info("Message sent to R : write.csv(exportsensitivity,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + "_SensitivityIndices.csv\")"); + //FIXME export through java to enable export when using Rserve (when distant Rserve). + engine.terminate(); + } + + } catch (Exception e) { + e.printStackTrace(); + // Error while processing + } + + } + + public String getDescription() { + return "Implementation of FrF2 method method using R"; + } +} Modified: trunk/sensitivity/SensitivityCalculatorRMorris.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRMorris.java 2009-04-07 12:44:34 UTC (rev 96) +++ trunk/sensitivity/SensitivityCalculatorRMorris.java 2009-04-07 15:09:00 UTC (rev 97) @@ -56,6 +56,9 @@ @Doc("Morris method") public class SensitivityCalculatorRMorris implements SensitivityCalculator { + @Doc("Integer giving the number of repetitions of the design, i.e. the number of elementary effect computed per factor. (Default value : 4)") + public int param_r; + /** to use log facility, just put in your code: log.info("..."); */ private static Log log = LogFactory .getLog(SensitivityCalculatorRMorris.class); @@ -77,6 +80,10 @@ int factorNumber = plan.getFactors().size(); List<Factor<? extends Serializable>> factors = plan.getFactors(); + if (((Integer) param_r == null) || ((Integer) param_r == 0)) { + param_r = 4; + } + //Test all factors, if one is discrete, return null for (int i = 0; i < factorNumber; i++) { if (factors.get(i).getDomain() instanceof DiscreteDomain) { @@ -99,8 +106,8 @@ // Adding the number of repetition parameter (r), the morris type // (type="oat") and the level vector - rInstruction = rInstruction - + "),r=4,design=list(type=\"oat\",levels=c("; + rInstruction = rInstruction + "),r=" + param_r + + ",design=list(type=\"oat\",levels=c("; // Creating the levels vector. for (int i = 0; i < factorNumber; i++) { @@ -129,8 +136,36 @@ } // Adding the grid.jump parameter - rInstruction = rInstruction + "),grid.jump=1),binf=c("; + rInstruction = rInstruction + "),grid.jump=c("; + // Creating the grid.jump vector. + for (int i = 0; i < factorNumber; i++) { + Domain<? extends Serializable> domain = factors.get(i).getDomain(); + if (i != (factorNumber - 1)) { + if (domain instanceof DiscreteDomain) { + rInstruction = rInstruction + + (Integer) (((DiscreteDomain<? extends Serializable>) domain) + .getValues().size() / 2) + ","; + } else if (domain instanceof ContinuousDomain) { + rInstruction = rInstruction + + (Integer) (((ContinuousDomain<? extends Serializable>) domain) + .getCardinality() / 2) + ","; + } + } else { + if (domain instanceof DiscreteDomain) { + rInstruction = rInstruction + + (Integer) (((DiscreteDomain<? extends Serializable>) domain) + .getValues().size() / 2); + } else if (domain instanceof ContinuousDomain) { + rInstruction = rInstruction + + (Integer) (((ContinuousDomain<? extends Serializable>) domain) + .getCardinality() / 2); + } + } + } + + rInstruction += ")),binf=c("; + // Adding the binf parameter for (int i = 0; i < factorNumber; i++) { Domain<? extends Serializable> domain = factors.get(i).getDomain(); Modified: trunk/sensitivity/SensitivityCalculatorROptimumLHS.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorROptimumLHS.java 2009-04-07 12:44:34 UTC (rev 96) +++ trunk/sensitivity/SensitivityCalculatorROptimumLHS.java 2009-04-07 15:09:00 UTC (rev 97) @@ -21,7 +21,6 @@ import java.io.File; import java.io.Serializable; import java.util.List; -import java.util.Vector; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; @@ -30,7 +29,6 @@ import org.codelutin.math.matrix.MatrixFactory; import org.codelutin.math.matrix.MatrixND; import org.codelutin.util.FileUtil; -import org.rosuda.JRI.REXP; import fr.ifremer.isisfish.datastore.SimulationStorage; import fr.ifremer.isisfish.simulator.sensitivity.DesignPlan; @@ -262,10 +260,10 @@ REngine engine = new RProxy(); try { - param_simulationNumber = simulationStorages.get(0).getParameter() + int sensitivityNumber = simulationStorages.get(0).getParameter() .getSensitivityExport().size(); - for (int k = 0; k < param_simulationNumber; k++) { + for (int k = 0; k < sensitivityNumber; k++) { // Set output directory engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() Modified: trunk/sensitivity/SensitivityCalculatorRRandomLHS.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRRandomLHS.java 2009-04-07 12:44:34 UTC (rev 96) +++ trunk/sensitivity/SensitivityCalculatorRRandomLHS.java 2009-04-07 15:09:00 UTC (rev 97) @@ -76,7 +76,7 @@ double[] dataframe = new double[0]; SensitivityScenarios thisExperiment = new SensitivityScenarios(); List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); - + if ((Integer) param_simulationNumber == 0) { param_simulationNumber = 10; } @@ -102,7 +102,6 @@ MatrixND morris = MatrixFactory.getInstance().create(dataframe, new int[] { factornumber, param_simulationNumber }); - // Setting up the scenarios. for (int j = 0; j < param_simulationNumber; j++) { Scenario experimentScenario = new Scenario(); @@ -129,25 +128,25 @@ thisExperimentScenarios.add(experimentScenario); thisExperiment.setScenarios(thisExperimentScenarios); } - + //Create the factors vectors for (int j = 0; j < factornumber; j++) { Factor<? extends Serializable> factor = thisExperimentScenarios - .get(0).getFactors().get(j); - if (factor - .getDomain() instanceof EquationContinuousDomain) { + .get(0).getFactors().get(j); + if (factor.getDomain() instanceof EquationContinuousDomain) { String vector = "factor" + j + "<-c("; for (int i = 0; i < param_simulationNumber; i++) { if (i < (param_simulationNumber - 1)) { vector = vector - + ((EquationContinuousDomain<?>)thisExperimentScenarios - .get(i).getFactors().get(j).getDomain()).getValue() - + ","; + + ((EquationContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue() + ","; } else { vector = vector - + ((EquationContinuousDomain<?>)thisExperimentScenarios - .get(i).getFactors().get(j).getDomain()).getValue(); + + ((EquationContinuousDomain<?>) thisExperimentScenarios + .get(i).getFactors().get(j) + .getDomain()).getValue(); } } @@ -198,18 +197,14 @@ String data = "data<-data.frame("; for (int j = 0; j < factornumber; j++) { if (j < factornumber - 1) { - data = data + "factor" + j + "=factor(factor" + j - + "),"; + data = data + "factor" + j + "=factor(factor" + j + "),"; } else { - data = data + "factor" + j + "=factor(factor" + j - + "))"; + data = data + "factor" + j + "=factor(factor" + j + "))"; } } engine.voidEval(data); log.info("Message sent to R : " + data); - - // Set output directory engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() @@ -234,16 +229,16 @@ public void analyzeResult(List<SimulationStorage> simulationStorages, File outputdirectory) throws SensitivityException { - + if ((Integer) param_simulationNumber == 0) { param_simulationNumber = 10; } - + REngine engine = new RProxy(); try { - + param_simulationNumber = simulationStorages.get(0).getParameter() - .getSensitivityExport().size(); + .getSensitivityExport().size(); for (int k = 0; k < param_simulationNumber; k++) { Added: trunk/sensitivity/SensitivityCalculatorRSobol.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRSobol.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorRSobol.java 2009-04-07 15:09:00 UTC (rev 97) @@ -0,0 +1,328 @@ +/* *##% + * Copyright (C) 2006 - 2009 Code Lutin + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + *##%*/ + +package sensitivity; + +import java.io.File; +import java.io.Serializable; +import java.util.List; +import java.util.Vector; + +import org.apache.commons.logging.Log; +import org.apache.commons.logging.LogFactory; +import org.codelutin.j2r.REngine; +import org.codelutin.j2r.RProxy; +import org.codelutin.math.matrix.MatrixFactory; +import org.codelutin.math.matrix.MatrixND; +import org.codelutin.util.FileUtil; +import org.rosuda.JRI.REXP; + +import fr.ifremer.isisfish.datastore.SimulationStorage; +import fr.ifremer.isisfish.simulator.sensitivity.DesignPlan; +import fr.ifremer.isisfish.simulator.sensitivity.Factor; +import fr.ifremer.isisfish.simulator.sensitivity.Scenario; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityCalculator; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityScenarios; +import fr.ifremer.isisfish.simulator.sensitivity.domain.ContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.DiscreteDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; +import fr.ifremer.isisfish.util.Doc; + +/** + * Implementation of Sobol method using R. + * + * @author jcouteau + * @version $Revision: 89 $ + * + * Last update : $Date: 2009-03-25 13:45:16 +0100 (mer., 25 mars 2009) $ By : + * $Author: jcouteau $ + */ + +public class SensitivityCalculatorRSobol implements SensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + private static Log log = LogFactory + .getLog(SensitivityCalculatorRFast.class); + + @Doc("doc") + public int param_n; + + @Doc("an integer, the maximum order in the ANOVA decomposition (all indices up to this order will be computed)") + public int param_order; + + @Doc("doc") + public int param_nboot; + + /** + * Retourne vrai si le calculateur sait gerer la cardinalité des facteurs + * continue. + * + * @return <tt>true</tt> s'il sait la gerer + */ + public boolean canManageCardinality() { + return true; + } + + public SensitivityScenarios compute(DesignPlan plan, File outputdirectory) { + + double[] dataframe = new double[0]; + int nbExperiments = 0; + int factorNumber = plan.getFactors().size(); + List<Factor<? extends Serializable>> factors = plan.getFactors(); + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + //Test all factors, if one is discrete, return null + for (int i = 0; i < factorNumber; i++) { + if (factors.get(i).getDomain() instanceof DiscreteDomain) { + return null; + } + } + + REngine engine = new RProxy(); + + try { + + engine.voidEval("library(sensitivity)"); + log.info("Message sent to R : " + "library(sensitivity)"); + + //Set the working directory (for import and exports) + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : " + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + engine.voidEval("X1<-data.frame(matrix(runif(" + factorNumber + "*" + + param_n + "),nrow=" + param_n + "))"); + log + .info("Message sent to R : " + + "X1<-data.frame(matrix(runif(" + factorNumber + + "*" + param_n + "),nrow=" + param_n + "))"); + engine.voidEval("X2<-data.frame(matrix(runif(" + factorNumber + "*" + + param_n + "),nrow=" + param_n + "))"); + log + .info("Message sent to R : " + + "X2<-data.frame(matrix(runif(" + factorNumber + + "*" + param_n + "),nrow=" + param_n + "))"); + + engine.voidEval("a<-sobol(model=NULL,X1=X1,X2=X2,order=" + + param_order + ",nboot=" + param_nboot + ")"); + log.info("Message sent to R : " + + "a<-sobol(model=NULL,X1=X1,X2=X2,order=" + param_order + + ",nboot=" + param_nboot + ")"); + + // Export the sobol object for the second run in a .sobol file + engine.voidEval("dput(a,file=\".sobol\")"); + log.info("Message sent to R" + "dput(a,file=\".sobol\")"); + // Export the X1 object for the second run in a .X1 file + engine.voidEval("dput(X1,file=\".X1\")"); + log.info("Message sent to R" + "dput(a,file=\".sobol\")"); + // Export the X2 object for the second run in a .X2 file + engine.voidEval("dput(X2,file=\".X2\")"); + log.info("Message sent to R" + "dput(a,file=\".sobol\")"); + + // Get back experiment plan + /*dataframe = (double[]) engine.eval("a$X");*/ + + Vector<REXP> dataframeVector = (Vector<REXP>) engine.eval("a$X"); + int dataframevectorsize = dataframeVector.size(); + int dataframevectorlength = dataframeVector.get(0).asDoubleArray().length; + dataframe = new double[dataframevectorsize * dataframevectorlength]; + for (int i = 0; i < dataframeVector.size(); i++) { + for (int j = 0; j < dataframeVector.get(0).asDoubleArray().length; j++) { + dataframe[i * dataframeVector.get(0).asDoubleArray().length + + j] = dataframeVector.get(i).asDoubleArray()[j]; + } + } + log.info("Message sent to R : " + "a$X"); + + } catch (Exception e) { + e.printStackTrace(); + } + + nbExperiments = dataframe.length / factorNumber; + + // Transform the result from R in a matrix + MatrixND fast = MatrixFactory.getInstance().create(dataframe, + new int[] { factorNumber, nbExperiments }); + + // Setting up the scenarios. + for (int j = 0; j < nbExperiments; j++) { + Scenario experimentScenario = new Scenario(); + for (int i = 0; i < factorNumber; i++) { + Factor<? extends Serializable> factor = plan.getFactors() + .get(i); + if ((factor.getDomain() instanceof MatrixContinuousDomain) + || (factor.getDomain() instanceof EquationContinuousDomain)) { + factor.setValueForIdentifier(Double.valueOf( + fast.getValue(new int[] { i, j })).toString()); + } else { + Double value = (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound() + + ((Double) ((ContinuousDomain<?>) factor + .getDomain()).getMaxBound() - (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound()) + * fast.getValue(new int[] { i, j }); + factor.setValueForIdentifier(value); + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + return thisExperiment; + + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) { + + REngine engine = new RProxy(); + try { + + // Call R + // Load sensitivity package into R (if package already loaded, + // nothing happens. + engine.voidEval("library(sensitivity)"); + log.info("Message sent to R : " + "library(sensitivity)"); + + //Set the working directory (for import and exports) + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : " + "setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + // Get the .X1 file + engine.voidEval("X1<-dget(\".X1\")"); + log.info("Message sent to R : " + "X1<-dget(\".X1\")"); + // Get the .X2 file + engine.voidEval("X2<-dget(\".X2\")"); + log.info("Message sent to R : " + "X2<-dget(\".X2\")"); + // Get the .sobol file + engine.voidEval("a<-dget(\".sobol\")"); + log.info("Message sent to R : " + "a<-dget(\".sobol\")"); + + int scenariosNumber = (Integer) engine.eval("length(a$X[,1])"); + log.info("Message sent to R : " + "length(a$X[,1])"); + + int sensitivityNumber = simulationStorages.get(0).getParameter() + .getSensitivityExport().size(); + + for (int k = 0; k < sensitivityNumber; k++) { + + // Creates the R expression to import results in R + String rInstruction = "results<-c("; + for (int l = 0; l < scenariosNumber; l++) { + File importFile = new File(simulationStorages.get(l) + .getDirectory().toString() + + File.separator + + SimulationStorage.RESULT_EXPORT_DIRECTORY, + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + simulationStorages.get(l).getParameter() + .getSensitivityExport().get(k) + .getExtensionFilename()); + String simulResult = FileUtil.readAsString(importFile); + double simulationResult = Double.valueOf(simulResult); + if (l < simulationStorages.size() - 1) { + rInstruction = rInstruction + simulationResult + ","; + } else { + rInstruction = rInstruction + simulationResult; + } + } + rInstruction = rInstruction + ")"; + + log.info("Message sent to R : " + rInstruction); + + // Send the simulation results + engine.voidEval(rInstruction); + + //Compute results + engine.voidEval("tell(a,y=results)"); + log.info("Message sent to R : " + "tell(a,y=results)"); + + //Create the data.frame of scenarios and results for export purpose + engine.voidEval("dfresults=data.frame(a$X,results)"); + log.info("Message sent to R : " + + "dfresults=data.frame(a$X,results)"); + + //Set working directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Export V + engine.voidEval("write.csv(a$V,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + + "_SensitivityIndices.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$V,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_V.csv\")"); + + //Export DD + engine.voidEval("write.csv(a$D,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_D.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$D,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_D.csv\")"); + + //Export S + engine.voidEval("write.csv(a$S,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_S.csv\")"); + log.info("Message sent to R : " + + "write.csv(a$S,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_S.csv\")"); + + //Export results + engine.voidEval("write.csv(dfresults,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + log.info("Message sent to R : " + + "write.csv(results,\"" + + simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k) + .getExportFilename() + "_Results.csv\")"); + //FIXME export through java to enable export when using Rserve + + } + } catch (Exception e) { + e.printStackTrace(); + } + } + + public String getDescription() { + return "Implementation of FAST method using R"; + } + +}