Author: jcouteau Date: 2009-03-25 13:53:22 +0000 (Wed, 25 Mar 2009) New Revision: 90 Added: trunk/sensitivity/SensitivityCalculatorROptimumLHS.java trunk/sensitivity/SensitivityCalculatorRRandomLHS.java Log: First implementation of Random and Optimum Latin Hypercube Sensitivity calculators Added: trunk/sensitivity/SensitivityCalculatorROptimumLHS.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorROptimumLHS.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorROptimumLHS.java 2009-03-25 13:53:22 UTC (rev 90) @@ -0,0 +1,333 @@ +/* *##% + * 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.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.SensitivityException; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityScenarios; +import fr.ifremer.isisfish.simulator.sensitivity.domain.ContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; + +/** + * Implementation of Random Latin Hypercube method using R. + * + * @author jcouteau + * @version $Revision: 1.0 $ + * + */ + +public class SensitivityCalculatorROptimumLHS implements SensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + static private Log log = LogFactory + .getLog(SensitivityCalculatorROptimumLHS.class); + + public int param_simulationNumber; + + public int param_MaxSweeps; + + public int param_eps; + + /** + * 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) + throws SensitivityException { + + int factornumber = plan.getFactors().size(); + double[] dataframe = new double[0]; + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + REngine engine = new RProxy(); + try { + + //Load the lhs library + engine.voidEval("library(lhs)"); + log.info("Message sent to R : " + "library(lhs)"); + + //Create the scenarios + engine.voidEval("x<-optimumLHS(" + param_simulationNumber + "," + + factornumber + "," + param_MaxSweeps + "," + param_eps + + ")"); + log.info("Message sent to R : " + "x<-optimumLHS(" + + param_simulationNumber + "," + factornumber + "," + + param_MaxSweeps + "," + param_eps + ")"); + + // Get back experiment plan + dataframe = (double[]) engine.eval("x"); + log.info("Message sent to R" + "x"); + + // Transform the result from R in a matrix + 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(); + 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( + morris.getValue(new int[] { i, j })) + .toString()); + } else { + Double value = (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound() + + ((Double) ((ContinuousDomain<?>) factor + .getDomain()).getMaxBound() - (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound()) + * morris.getValue(new int[] { i, j }); + factor.setValueForIdentifier(value); + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + + //Create the factors vectors + for (int j = 0; j < factornumber; j++) { + String vector = "factor" + j + "<-c("; + for (int i = 0; i < param_simulationNumber; i++) { + if (i < (param_simulationNumber - 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 .randomlhs.csv file + engine.voidEval("write.csv(data,file=\".optimumlhs.csv\")"); + log.info("Message sent to R : " + + "write.csv(data,file=\".optimumlhs.csv\")"); + + engine.terminate(); + + } catch (Exception e) { + e.printStackTrace(); + // Error while processing + } + + return thisExperiment; + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) throws SensitivityException { + + REngine engine = new RProxy(); + try { + + for (int k = 0; k < param_simulationNumber; k++) { + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Get back the scenarios + //engine.voidEval("factors<-read.csv(\".expandgrid.csv\",row.names=1,col.names=1)"); + engine.voidEval("factors<-read.csv(\".optimumlhs.csv\")"); + log + .info("Message sent to R : factors<-read.csv(\".optimumlhs.csv\")"); + + //Get back the factors number + int factorNumber = ((Double) engine.eval("length(factors)-1")) + .intValue(); + //factorNumber=factorNumber-1; + + //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); + + //Call aov() + String 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); + + /*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(); + } + log.info("end"); + + } catch (Exception e) { + e.printStackTrace(); + // Error while processing + } + + } + + public String getDescription() { + return "Implementation of Random Latin Hypercube method method using R"; + } + +} Added: trunk/sensitivity/SensitivityCalculatorRRandomLHS.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRRandomLHS.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorRRandomLHS.java 2009-03-25 13:53:22 UTC (rev 90) @@ -0,0 +1,332 @@ +/* *##% + * 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.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.SensitivityException; +import fr.ifremer.isisfish.simulator.sensitivity.SensitivityScenarios; +import fr.ifremer.isisfish.simulator.sensitivity.domain.ContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; + +/** + * Implementation of Random Latin Hypercube method using R. + * + * @author jcouteau + * @version $Revision: 1.0 $ + * + */ + +public class SensitivityCalculatorRRandomLHS implements SensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + static private Log log = LogFactory + .getLog(SensitivityCalculatorRRandomLHS.class); + + public int param_simulationNumber; + + /** + * 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) + throws SensitivityException { + + int factornumber = plan.getFactors().size(); + double[] dataframe = new double[0]; + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + REngine engine = new RProxy(); + try { + + //Load the lhs library + engine.voidEval("library(lhs)"); + log.info("Message sent to R : " + "library(lhs)"); + + //Create the scenarios + engine.voidEval("x<-randomLHS(" + param_simulationNumber + "," + + factornumber + ")"); + log.info("Message sent to R : " + "x<-randomLHS(" + + param_simulationNumber + "," + factornumber + ")"); + + // Get back experiment plan + dataframe = (double[]) engine.eval("x"); + log.info("Message sent to R" + "x"); + + // Transform the result from R in a matrix + 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(); + 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( + morris.getValue(new int[] { i, j })) + .toString()); + } else { + Double value = (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound() + + ((Double) ((ContinuousDomain<?>) factor + .getDomain()).getMaxBound() - (Double) ((ContinuousDomain<?>) factor + .getDomain()).getMinBound()) + * morris.getValue(new int[] { i, j }); + factor.setValueForIdentifier(value); + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + + //Create the factors vectors + for (int j = 0; j < factornumber; j++) { + String vector = "factor" + j + "<-c("; + for (int i = 0; i < param_simulationNumber; i++) { + if (i < (param_simulationNumber - 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 .randomlhs.csv file + engine.voidEval("write.csv(data,file=\".randomlhs.csv\")"); + log.info("Message sent to R : " + + "write.csv(data,file=\".randomlhs.csv\")"); + + engine.terminate(); + + } catch (Exception e) { + e.printStackTrace(); + // Error while processing + } + + return thisExperiment; + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) throws SensitivityException { + + REngine engine = new RProxy(); + try { + + for (int k = 0; k < param_simulationNumber; k++) { + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + //Get back the scenarios + //engine.voidEval("factors<-read.csv(\".expandgrid.csv\",row.names=1,col.names=1)"); + engine.voidEval("factors<-read.csv(\".randomlhs.csv\")"); + log + .info("Message sent to R : factors<-read.csv(\".randomlhs.csv\")"); + + //Get back the factors number + int factorNumber = ((Double) engine.eval("length(factors)-1")) + .intValue(); + //factorNumber=factorNumber-1; + + //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); + + //Call aov() + String 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); + + /*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(); + } + log.info("end"); + + } catch (Exception e) { + e.printStackTrace(); + // Error while processing + } + + } + + public String getDescription() { + return "Implementation of Random Latin Hypercube method method using R"; + } + +}