Author: jcouteau Date: 2009-04-20 16:25:44 +0000 (Mon, 20 Apr 2009) New Revision: 116 Added: trunk/sensitivity/SensitivityCalculatorRegularFractions.java trunk/sensitivity/regularfractions.R Log: Regular Fractions method working, conform to the new CDC. Added: trunk/sensitivity/SensitivityCalculatorRegularFractions.java =================================================================== --- trunk/sensitivity/SensitivityCalculatorRegularFractions.java (rev 0) +++ trunk/sensitivity/SensitivityCalculatorRegularFractions.java 2009-04-20 16:25:44 UTC (rev 116) @@ -0,0 +1,614 @@ +package sensitivity; + +import java.io.File; +import java.util.List; + +import javax.swing.Box; +import javax.swing.JLabel; +import javax.swing.JOptionPane; +import javax.swing.JScrollPane; +import javax.swing.JTextPane; + +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.AbstractSensitivityCalculator; +import fr.ifremer.isisfish.simulator.sensitivity.DesignPlan; +import fr.ifremer.isisfish.simulator.sensitivity.Domain; +import fr.ifremer.isisfish.simulator.sensitivity.Factor; +import fr.ifremer.isisfish.simulator.sensitivity.Scenario; +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.DiscreteDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.EquationContinuousDomain; +import fr.ifremer.isisfish.simulator.sensitivity.domain.MatrixContinuousDomain; +import fr.ifremer.isisfish.util.Doc; + +public class SensitivityCalculatorRegularFractions extends + AbstractSensitivityCalculator { + + /** to use log facility, just put in your code: log.info("..."); */ + private static Log log = LogFactory + .getLog(SensitivityCalculatorRFrF2.class); + + @Doc("the path of the directory where the R function is stored") + public String param_pathToFunction = ""; + + @Doc("unique prime number of levels of all input and unit factors") + public int param_p = 2; + + @Doc("number of unit factors (so that there are N=p^r units)") + public int param_r = 2; + + @Doc("resolution of the fraction") + public int param_resolution = 2; + + @Doc("True to be able to modify the code sent to R") + public boolean param_modifR = false; + + /** + * 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) { + + setIsisFactorsR(plan, outputdirectory); + + double[] dataframeplan = new double[0]; + int nbExperiments = 0; + int factorNumber = plan.getFactors().size(); + List<Factor> factors = plan.getFactors(); + SensitivityScenarios thisExperiment = new SensitivityScenarios(); + List<Scenario> thisExperimentScenarios = thisExperiment.getScenarios(); + + String factorNames = ""; + + //Test all factors, if one is discrete, return null + //Create a string with all factors names + for (int i = 0; i < factorNumber; i++) { + Domain domain = factors.get(i).getDomain(); + if (i != 0) { + factorNames += ","; + } + + factorNames += "\"" + factors.get(i).getName() + "\""; + + if (domain instanceof DiscreteDomain) { + JOptionPane + .showMessageDialog( + null, + "Error", + factors.get(i).getName() + + " has a discrete domain, this is not acceptable for this method.", + JOptionPane.ERROR_MESSAGE); + return null; + } + } + + REngine engine = new RProxy(); + + try { + + //Set working directory to get Isis R session + log.info("setwd(\"" + outputdirectory.getParent() + "\")"); + engine.voidEval("setwd(\"" + outputdirectory.getParent() + "\")"); + + //Get Isis R session + log.info("load(\".RData\")"); + engine.voidEval("load(\".RData\")"); + + //Set the working directory (to import the R function) + engine.voidEval("setwd(\"" + param_pathToFunction + "\")"); + log.info("Message sent to R : " + "setwd(\"" + param_pathToFunction + + "\")"); + + //Import the function + log.info("source(\"regularfractions.R\")"); + engine.voidEval("source(\"regularfractions.R\")"); + + //Create the instruction + String rInstruction = "x<-regular.fraction(%s,%s,%s,%s)"; + String rCall = String.format(rInstruction, factors.size(), param_p, + param_r, param_resolution); + + if (param_modifR) { + JLabel label = new JLabel( + "Modifier le code R envoyé si vous le souhaitez"); + JTextPane text = new JTextPane(); + text.setText(rCall); + text.setSize(400, 400); + text.setPreferredSize(text.getSize()); + + Box box = Box.createVerticalBox(); + box.add(label); + box.add(new JScrollPane(text)); + + JOptionPane.showMessageDialog(null, box, "R modif", + JOptionPane.QUESTION_MESSAGE); + rCall = text.getText(); + } + + // Run function + engine.voidEval(rCall); + log.info("Message sent to R : " + rCall); + + // Run function + engine.voidEval("call<-\"" + rCall + "\""); + log.info("Message sent to R : " + "call<-\"" + rCall + "\""); + + // Creating the factors vector. + rInstruction = "factornames<-c(%s)"; + rCall = String.format(rInstruction, factorNames); + + engine.voidEval(rCall); + log.info("Message sent to R : " + rCall); + + // Set output directory + 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("dput(x,file=\".regularfractions\")"); + log.info("Message sent to R : " + + "dput(x,file=\".regularfractions\")"); + + // Export the factornames object for the second run in a .factornames file + engine.voidEval("dput(factornames,file=\".factornames\")"); + log.info("Message sent to R : " + + "dput(factornames,file=\".factornames\")"); + + // Get back experiment plan + dataframeplan = (double[]) engine.eval("x$plan"); + log.info("Message sent to R : " + "x$plan"); + + //Get back the simulation number + log.info("length(x[,1])"); + int simulationNumber = (Integer) engine.eval("length(x$plan[,1])"); + + // Transform the result from R in a matrix + MatrixND morris = MatrixFactory.getInstance().create(dataframeplan, + new int[] { factors.size(), simulationNumber }); + + // Setting up the scenarios. + for (int j = 0; j < simulationNumber; j++) { + Scenario experimentScenario = new Scenario(); + for (int i = 0; i < factors.size(); i++) { + Factor factor = plan.getFactors().get(i); + if ((factor.getDomain() instanceof MatrixContinuousDomain) + || (factor.getDomain() instanceof EquationContinuousDomain)) { + factor.setValueForIdentifier(morris.getValue(new int[] { + i, j })); + } else { + Double value = (Double) ((ContinuousDomain) factor + .getDomain()).getMinBound() + + ((Double) ((ContinuousDomain) factor + .getDomain()).getMaxBound() - (Double) ((ContinuousDomain) factor + .getDomain()).getMinBound()) + * (morris.getValue(new int[] { i, j }) / (param_p - 1)); + factor.setValueForIdentifier(value); + } + experimentScenario.addFactor(factor); + } + thisExperimentScenarios.add(experimentScenario); + thisExperiment.setScenarios(thisExperimentScenarios); + } + + String dataframe = "data<-data.frame("; + + //Create the factors vectors and the dataframe instruction + for (int j = 0; j < factorNumber; j++) { + Factor factor = thisExperimentScenarios.get(0).getFactors() + .get(j); + if (factor.getDomain() instanceof EquationContinuousDomain) { + + String vector = factor.getName().replaceAll(" ", "") + + "<-c("; + for (int i = 0; i < simulationNumber; i++) { + if (i < (simulationNumber - 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.getName().replaceAll(" ", "") + + "<-c("; + for (int i = 0; i < simulationNumber; i++) { + if (i < (simulationNumber - 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.getName().replaceAll(" ", "") + + "<-c("; + for (int i = 0; i < simulationNumber; i++) { + if (i < (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); + } + + if (j < factorNumber - 1) { + dataframe = dataframe + + factor.getName().replaceAll(" ", "") + "=factor(" + + factor.getName().replaceAll(" ", "") + "),"; + } else { + dataframe += factor.getName().replaceAll(" ", "") + + "=factor(" + factor.getName().replaceAll(" ", "") + + "))"; + } + + } + engine.voidEval(dataframe); + log.info("Message sent to R : " + dataframe); + + log.info("Message sent to R : " + + "isis.factor.distribution<-c(0.0)"); + engine.voidEval("isis.factor.distribution<-c(0.0)"); + + log + .info("Message sent to R : " + + "isis.MethodExp<-list(\"isis.factors\"=isis.factors,\"isis.factor.distribution\"=isis.factor.distribution,\"call\"=call)"); + engine + .voidEval("isis.MethodExp<-list(\"isis.factors\"=isis.factors,\"isis.factor.distribution\"=isis.factor.distribution,\"call\"=call)"); + + log + .info("Message sent to R : " + + "attributes(isis.MethodExp)<-list(nomModel=\"isis-fish-externe-R\")"); + engine + .voidEval("attributes(isis.MethodExp)<-list(nomModel=\"isis-fish-externe-R\")"); + + //Create isis.Simule + log.info("Message sent to R : " + "isis.simule<-data.frame(data)"); + engine.voidEval("isis.simule<-data.frame(data)"); + + log + .info("Message sent to R : " + + "attributes(isis.simule)<-list(nomModel=\"isis-fish-externe-R\")"); + engine + .voidEval("attributes(isis.simule)<-list(nomModel=\"isis-fish-externe-R\")"); + + log.info("Message sent to R : " + + "names(isis.simule)<-isis.factors[[1]]"); + engine.voidEval("names(isis.simule)<-isis.factors[[1]]"); + + // Export the data.frame object for the second run in a .expandgrid file + engine.voidEval("write.csv(data,file=\".data.csv\")"); + log + .info("Message sent to R : write.csv(data,file=\".data.csv\")"); + + + //Set working directory to save Isis R session + log.info("setwd(\"" + outputdirectory.getParent() + "\")"); + engine.voidEval("setwd(\"" + outputdirectory.getParent() + "\")"); + + // Save Isis R session + log.info("save.image()"); + engine.voidEval("save.image()"); + + engine.terminate(); + + } catch (RException eee) { + eee.printStackTrace(); + throw new RuntimeException("R evaluation failed", eee); + } + + return thisExperiment; + + } + + public void analyzeResult(List<SimulationStorage> simulationStorages, + File outputdirectory) throws SensitivityException { + + REngine engine = new RProxy(); + try { + + //Set working directory to get Isis R session + log.info("setwd(\"" + outputdirectory.getParent() + "\")"); + engine.voidEval("setwd(\"" + outputdirectory.getParent() + "\")"); + + //Get Isis R session + log.info("load(\".RData\")"); + engine.voidEval("load(\".RData\")"); + + // Set output directory + engine.voidEval("setwd(\"" + outputdirectory.getAbsolutePath() + + "\")"); + log.info("Message sent to R : setwd(\"" + + outputdirectory.getAbsolutePath() + "\")"); + + // Get the regularfractions object from the .regularfractions file + engine.voidEval("x<-dget(\".regularfractions\")"); + log.info("Message sent to R : " + "x<-dget(\".regularfractions\")"); + + // Get the factornames object from the .factornames file + engine.voidEval("factornames<-dget(\".factornames\")"); + log.info("Message sent to R : " + + "factornames<-dget(\".factornames\")"); + + //Get back the scenarios + engine.voidEval("factors<-read.csv(\".data.csv\")"); + log + .info("Message sent to R : factors<-read.csv(\".data.csv\")"); + + //Get back the factors number + int factorNumber = ((Double) engine.eval("length(factors)-1")) + .intValue(); + + int sensitivityNumber = simulationStorages.get(0).getParameter() + .getSensitivityExport().size(); + + for (int k = 0; k < sensitivityNumber; k++) { + + String name = simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k).getExportFilename(); + + //Create the results vectors + String result = name + "<-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); + + //Put results in isis.simule + engine.voidEval("isis.simule<-data.frame(isis.simule," + name + + ")"); + log.info("Message sent to R : " + + "isis.simule<-data.frame(isis.simule," + name + ")"); + } + + //adding attribute to isis.Simule + log + .info("Message sent to R : " + + "attributes(isis.simule)<-list(nomModel=\"isis-fish-externe-R\")"); + engine + .voidEval("attributes(isis.simule)<-list(nomModel=\"isis-fish-externe-R\")"); + + for (int k = 0; k < sensitivityNumber; k++) { + + // Creates the R expression to import results in R + String name = simulationStorages.get(0).getParameter() + .getSensitivityExport().get(k).getExportFilename(); + + //Create the dataforaov data.frame + String dataframe = "dataforaov<-data.frame(factors," + name + + "=" + name + ")"; + engine.voidEval(dataframe); + log.info("Message sent to R : " + dataframe); + + //Call aov() + String aovCall = "aovresult<-aov(" + name + "~"; + for (int j = 0; j < factorNumber; j++) { + log.info("Message sent to R : " + "names(factors)[" + + (j + 2) + "]"); + + if (j < (factorNumber - 1)) { + aovCall = aovCall + + engine + .eval("names(factors)[" + (j + 2) + "]") + + "+"; + } else { + aovCall = aovCall + + engine + .eval("names(factors)[" + (j + 2) + "]") + + ",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\")"); + engine.voidEval("row.names(exportsensitivity)<-factornames"); + log.info("Message sent to R : " + + "row.names(exportsensitivity)<-factornames"); + + //Set dataforaov names + engine + .voidEval("resultsnames<-c(\"Simulation\",factornames,\"Result\")"); + log + .info("Message sent to R : " + + "resultsnames<-c(\"Simulation\",factornames,\"Result\")"); + engine.voidEval("names(dataforaov)<-resultsnames"); + log.info("Message sent to R : " + + "names(dataforaov)<-resultsnames"); + + /*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). + + //creating isis.methodAnalyse + log + .info("Message sent to R : " + + "isis.methodAnalyse<-list(\"isis.factors\"=isis.factors,\"isis.factor.distribution\"=isis.factor.distribution,\"isis.simule\"=isis.simule,call_method=\"" + + aovCall + "\"" + ",aovresult)"); + engine + .voidEval("isis.methodAnalyse<-list(\"isis.factors\"=isis.factors,\"isis.factor.distribution\"=isis.factor.distribution,\"isis.simule\"=isis.simule,call_method=\"" + + aovCall + "\"" + ",aovresult)"); + + String renameIsisMethodAnalyse = "%s.isis.methodAnalyse<-isis.methodAnalyse"; + String simulationName = simulationStorages.get(0).getName() + .replaceAll("-", ""); + log.info("Message sent to R : " + + String.format(renameIsisMethodAnalyse, simulationName + + "." + name)); + engine.voidEval(String.format(renameIsisMethodAnalyse, + simulationName + "." + name)); + + } + + //Rename R objects for saving purpose + + String renameIsisSimule = "%s.isis.simule<-isis.simule"; + String renameIsisFactorDistribution = "%s.isis.factor.distribution<-isis.factor.distribution"; + String renameIsisFactor = "%s.isis.factor<-isis.factors"; + String renameIsisMethodExp = "%s.isis.methodExp<-isis.MethodExp"; + + String simulationName = simulationStorages.get(0).getName() + .replaceAll("-", ""); + + log.info("Message sent to R : " + + String.format(renameIsisSimule, simulationName)); + engine.voidEval(String.format(renameIsisSimule, simulationName)); + + log.info("Message sent to R : " + + String.format(renameIsisFactorDistribution, + simulationName)); + engine.voidEval(String.format(renameIsisFactorDistribution, + simulationName)); + + log.info("Message sent to R : " + + String.format(renameIsisFactor, simulationName)); + engine.voidEval(String.format(renameIsisFactor, simulationName)); + + log.info("Message sent to R : " + + String.format(renameIsisMethodExp, simulationName)); + engine.voidEval(String.format(renameIsisMethodExp, simulationName)); + + //Set working directory to save Isis R session + log.info("setwd(\"" + outputdirectory.getParent() + "\")"); + engine.voidEval("setwd(\"" + outputdirectory.getParent() + "\")"); + + // Save Isis R session + log.info("save.image()"); + engine.voidEval("save.image()"); + + engine.terminate(); + + } catch (Exception eee) { + eee.printStackTrace(); + throw new RuntimeException("R evaluation failed", eee); + } + + } + + public String getDescription() { + return "Implementation of Regular fractions method using R"; + } + +} Added: trunk/sensitivity/regularfractions.R =================================================================== --- trunk/sensitivity/regularfractions.R (rev 0) +++ trunk/sensitivity/regularfractions.R 2009-04-20 16:25:44 UTC (rev 116) @@ -0,0 +1,211 @@ +#================================================================================ +# FONCTIONS DE CALCUL D'UNE FRACTION REGULIERE DE RESOLUTION DONNEE +# Auteur: H. Monod, INRA Jouy en Josas +# Copyright INRA 2009 +# Ce programme est une version preliminaire et simplifiee d'une librairie R +# en preparation par A. Kobilinsky, H. Monod, A. Bouvier +#================================================================================ +regular.fraction <- function(s,p,r,resolution){ + # DESCRIPTION + # generates a regular fractional factorial design of given resolution, + # for s factors at p levels in p^r units + # ARGUMENTS + # s : number of input factors + # p : unique prime number of levels of all input and unit factors + # r : number of unit factors (so that there are N=p^r units) + # resolution : resolution of the fraction + # max.sol : maximum number of solutions + # DETAILS + # This is a simplified version of a more general library in preparation. + # In this version, all factors must have the same prime number of levels + # and only fractions with a given resolution can be constructed. The first + # q factors are used as basic factors. The first solution is kept although + # it may not be the most interesting one (no control of aberration). This + # function is programmed entirely in R and so it is not efficient with respect + # to computer time. There is no explicit check on the arguments and so it + # is up to the user to restrict p to a prime number such as 2, 3, 5 or 7. + # OUTPUT: + # a list with two components: plan (the design in base p) and matrice.cle + # (the design key). The design has N=p^r rows (units) and s columns (factors). + # All its elements are integers modulo p that represent the factor levels. + + # ensemble ineligible + cat("Determination des termes ineligibles: ") + ineligible <- diag(s) + for(reso in 2:(resolution-1)){ + combis <- combn(s,reso) + ncombi <- ncol(combis) + select <- cbind( c(combis), rep(seq(ncombi),rep(reso,ncombi)) ) + ineli <- matrix(0,s,ncombi) + ineli[select] <- 1 + ineligible <- cbind(ineligible,ineli) + } + cat(ncol(ineligible)," termes ineligibles.\n") + if( (p!=2) ){ + ineligible <- representative.basep(ineligible,p) + } + # Identification of the last non-zero coefficients in each ineligible trt character + ineligible.lnz <- apply(ineligible, 2, function(x){max(seq(along=x)[x!=0])}) + # initialisation of PhiStar by using the first q factors as basic factors + PhiStar <- diag(r) + # + f <- ncol(PhiStar) + if(s == f){ + check <- !any(apply(((PhiStar %*% ineligible)%%p)==0, 2, all)) + if(check) return(list(PhiStar)) + } + # Calculation of the set of initially admissible elements of U* + admissible <- t(convertinto.basep(seq((p^r)-1),p)) + nb.admissible <- ncol(admissible) + # Backtrack search - preliminaries + eeU <- list(length=s-f) + leeU <- rep(NA,s-f) + neeU <- rep(0,s-f) + # Backtrack search + cat("Recherche d'une solution (algorithme backtrack).\n") + jprev <- 0 ; j <- 1 + solved <- FALSE + while((j > 0)&(!solved)){ + PhiStar <- PhiStar[,seq(f+j-1), drop=FALSE] + if(jprev < j){ + ineligible.j <- ineligible[ seq(f+j-1), ineligible.lnz==(f+j), drop=FALSE ] + admissible.keep <- planor.kernelcheck.basep(PhiStar, admissible, ineligible.j, p) + eeU[[j]] <- seq(nb.admissible)[admissible.keep] + leeU[j] <- length(eeU[[j]]) + neeU[j] <- 0 + } + if(neeU[j] < leeU[j]){ + neeU[j] <- neeU[j]+1 + newcolj <- (eeU[[j]])[neeU[j]] + PhiStar <- cbind(PhiStar,admissible[,newcolj]) + if(j == (s-f)){ + cat("Solution obtenue. ") + solved <- TRUE + jprev <- j ; j <- j + } + else{ + jprev <- j ; j <- j+1 + } + } + else{ + jprev <- j ; j <- j-1 + } + } + if(solved){ + # Construction du plan + plan <- crossing(rep(p,r),start=0) %*% PhiStar %%p + # Sortie + out <- list(plan=plan, matrice.cle=PhiStar, p=p) + } + else{ + cat("Pas de solution. ") + out <- NULL + } + cat("Recherche terminee.\n") + return(out) +} +#--------------------------------------------------------------------------- +planor.kernelcheck.basep <- function(PhiStar, admissible, IneligibleSet, p){ + ImagesIS <- (- PhiStar %*% IneligibleSet)%%p + avoid <- convertfrom.basep( t(ImagesIS), p) + candidate <- convertfrom.basep( t(admissible), p) + test <- !(candidate %in% avoid) + return(test) +} +#--------------------------------------------------------------------------- +convertinto.basep <- function (x, p) { + # Conversion of an integer or integer vector x into base p + # The coefficients are ordered by increasing powers of p + if (!is.numeric(x)) + stop("cannot decompose non-numeric arguments") + if (length(x) > 1) { + l <- matrix(0, length(x), length(Recall(max(x),p))) + for(i in seq(along = x)){ + dec.i <- Recall(x[i],p) + l[i, seq(along=dec.i) ] <- dec.i + } + return(l) + } + if (x != round(x) || x < 0) + return(x) + val <- x%%p + while ( (x <- x%/%p) > 0 ) { + newval <- x%%p + val <- c(val,newval) + } + return(val) +} +#--------------------------------------------------------------------------- +convertfrom.basep <- function (x, p) { + # Conversion of integers x coded as vectors of coefficients in base p + # to classical integers in base 10 + + if (!is.numeric(x)) + stop("cannot recompose non-numeric arguments") + if( (max(x)>p) || (min(x)<0) ) + stop("x must be reduced modulo p") + if (is.matrix(x)) { + l <- rep(NA, nrow(x)) + for(i in seq(along = l)){ + l[i] <- Recall(x[i,],p) + } + return(l) + } + val <- sum( x * p^(seq(along=x)-1) ) + return(val) +} +#--------------------------------------------------------------------------- +inverses.basep <- function(p){ + # Raw calculation of the inverses modulo p + + if(p==2) return(1) + else if(p==3) return(c(1,2)) + products <- outer(seq(2,p-2), seq(2,p-2), "*")%%p + inverses <- 1 + apply(products, 1, function(x){ seq(along=x)[x==1] }) + return( c(1,inverses,p-1) ) +} +#--------------------------------------------------------------------------- +representative.basep <- function(mat,p){ + # generates the minimal set of representatives in base p + # of the columns x of matrix mat + + mat <- as.matrix(mat) + # + if(p==2) return(mat %%2) + # + representative <- NULL + for(j in seq(ncol(mat))){ + x <- mat[,j] + select <- seq(x)[x != 0] + nbtocross <- length(select)-1 + if( nbtocross <= 0 ) mat.j <- x + else{ + select <- select[seq(nbtocross)] + N <- (p-1)^nbtocross + mat.j <- matrix(x, nrow(mat), N) + mat.j[select,] <- t( crossing(rep(p-1,nbtocross),start=1) ) + } + representative <- cbind(representative, mat.j) + } + return(representative %%p) +} +#--------------------------------------------------------------------------- +crossing <- function(n,start=1){ + # Generates all n1 x n2 x ... x ns combinations of size s with n1,...,ns integers + + N <- prod(n) + s <- length(n) + n <- c(n,1) + crosses <- matrix(NA, N, s) + for(i in seq(s)) + { + motif <- start + seq(n[s+1-i])-1 + repet1 <- rep( prod(n[s+1-i+seq(i)]), n[s+1-i] ) + if(i==s){ repet2 <- 1 } + else{ repet2 <- prod(n[seq(s-i)]) } + crosses[,s-i+1] <- rep( rep( motif, repet1 ), repet2 ) + } + return(crosses) +} +#--------------------------------------------------------------------------- +