Initial commit

README.md

@@ -0,0 +1,94 @@
+## PIES
+
+This repository contains code for the simulations in [_Partial Identification by Extending
+Subdistributions_](https://docs.google.com/viewer?a=v&pid=sites&srcid=ZGVmYXVsdGRvbWFpbnxhdG9yZ292aXRza3l8Z3g6MTIwMTM3MDlhYzAzNzlh) (PIES) by
+[Alexander Torgovitsky](https://sites.google.com/site/atorgovitsky/).  
+
+### Software Prerequisites
+
+* [MATLAB](https://www.mathworks.com/products/matlab.html).
+
+  No special toolboxes are required.
+
+* [A Mathematical Programming Language (AMPL)](http://ampl.com/).
+
+ The student version of AMPL is size restricted.  
+ Some simulations will probably run with the student version, but most will require a full license.
+
+* A linear programming solver for AMPL.
+
+ The default is [Gurobi](http://www.gurobi.com/).  
+ It can easily be changed by passing e.g.  `Settings.IdentifiedSet = 'cplex'`
+ when calling `./src/IdentifiedSet.m`.  
+ See the discussion on usage below.
+
+* The [AMPL-MATLAB API](http://ampl.com/api/latest/matlab/getting-started.html).
+
+* Linux (or perhaps OSX).
+
+ I coded these simulations for Linux and made no attempt to be platform-independent.  
+ However, the code is primarily in MATLAB, so should be mostly
+ platform-independent.  
+ Some file operations are used for recording the results of simulations.  
+ These would be likely sources of issues for other operating systems, but should
+ be easy enough to fix.
+
+### Setup and Usage
+
+* **Important first step**
+
+ Open `./cfg/Config.m`. 
+
+ Change `SAVEDIRECTORY` to a directory where you want simulations to be saved.
+
+ Change `AMPLAPISETUPPATH` to the correct location for your installation of the
+ AMPL-MATLAB API.
+
+* The primary code is contained is `./src/IdentifiedSet.m` and the routines
+  called from within.
+
+ It contains many options that can be set, which are given default values in the
+ structure called `Settings` that is defined at the top of that file.
+
+ Another structure called `Assumptions` serves a similar purpose but only
+ contains options that pertain to what assumptions are imposed when constructing
+ the identified set.
+
+* The directory `./run/` contains a file called `./Run.m` that can be used to
+  run `IdentifiedSet.m` under limited pre-set options.
+
+  For example, the command `Run(2,3,5)` would run specification number 2 from
+  the paper with (d\_{1}, d\_{2}) = (3,5) in the notation there.
+
+* Sequentially running all of the simulations in the paper will take a long
+  time.
+ The process can be sped up by using the file `./run/MultiBatchRun.m`, which opens
+ multiple MATLAB threads.
+ (Unfortunately, the AMPL-MATLAB API is not easy to parallelize.)
+ The command to reproduce the results in the paper is
+
+ `MultiBatchRun(1:1:14, 'your-save-dir')`
+
+ Note that this will open 14 MATLAB and AMPL instances at one time, which will strain a typical system.
+ To open fewer threads, pass a smaller array of numbers in the first argument, wait a while, then pass the remaining ones, using the same directory name.
+
+### My Software Versions
+
+The numerical results in the published paper were run with
+
+* MATLAB version 8.6.0.267246 (R2015b)
+* AMPL version 20170711
+* Gurobi version 7.5.0
+* AMPL-MATLAB API version 1.2.2.
+
+### Software Acknowledgment
+
+The code uses the MATLAB function
+[MergeBrackets.m](https://www.mathworks.com/matlabcentral/fileexchange/24254-interval-merging) developed by Bruno Luong.  
+A copy of this code is included in `./src/MergeBrackets.m`.
+
+### Problems or Bugs?
+
+Please use [GitHub] to open an issue and I will be happy to look into it.
+
+[GitHub]: http://www.github.com/a-torgovitsky/pies

cfg/Config.m

@@ -0,0 +1,2 @@
+SAVEDIRECTORY = '~'
+AMPLAPISETUPPATH = '/usr/local/amplapi/matlab/setUp.m'

run/BatchRun.m

@@ -0,0 +1,57 @@
+%*******************************************************************************
+% BatchRun.m
+%
+% Run all of the simulations corresponding to a given batch number, as
+% defined in this script.
+% This is used to simplify the process of running all specifications that
+% are reported in the paper.
+%*******************************************************************************
+
+function [] = BatchRun(BatchNumber, SaveDir)
+
+if ~exist('SaveDir', 'var')
+    SaveDir = [];
+end
+
+%*******************************************************************************
+% Define batches
+%*******************************************************************************
+NList = [3 3; 3 5; 5 3; 5 5; 7 3];
+for b = 1:1:3
+    Batch{b} = [b*ones(size(NList, 1), 1), NList, zeros(size(NList, 1), 1)];
+end
+
+NList = [3 3; 3 5; 5 3; 5 5];
+for b = 4:1:7
+    Batch{b} = [b*ones(size(NList, 1), 1), NList, zeros(size(NList, 1), 1)];
+end
+
+NList1 = [3 3; 5 5];
+NList2 = [3 5; 5 3];
+Batch{8} = [8*ones(size(NList1, 1), 1), NList1, zeros(size(NList1, 1), 1)];
+Batch{9} = [8*ones(size(NList2, 1), 1), NList2, zeros(size(NList2, 1), 1)];
+
+Batch{10} = [9*ones(size(NList1, 1), 1), NList1, zeros(size(NList1, 1), 1)];
+Batch{11} = [9*ones(size(NList2, 1), 1), NList2, zeros(size(NList2, 1), 1)];
+
+S = 5;
+NList  = [3 4; 3 9; 3 14];
+Batch{12} = [S*ones(size(NList, 1), 1), NList, ones(size(NList, 1), 1)];
+
+NList = [3 19];
+Batch{13} = [S*ones(size(NList, 1), 1), NList, ones(size(NList, 1), 1)];
+
+NList = [3 24];
+Batch{14} = [S*ones(size(NList, 1), 1), NList, ones(size(NList, 1), 1)];
+
+assert(BatchNumber <= length(Batch));
+
+for s = 1:1:size(Batch{BatchNumber}, 1)
+    SaveDir = Run(  Batch{BatchNumber}(s,1),...
+                    Batch{BatchNumber}(s,2),...
+                    Batch{BatchNumber}(s,3),...
+                    Batch{BatchNumber}(s,4),...
+                    SaveDir);
+end
+
+end

run/MultiBatchRun.m

@@ -0,0 +1,28 @@
+%*******************************************************************************
+% MultiBatchRun.m
+%
+% Run all batches in BatchList, each one on a separate instance of Matlab.
+% Save them all in the same directory "SaveDir" which is required in
+% this routine.
+%*******************************************************************************
+function [] = MultiBatchRun(BatchList, SaveDir)
+
+    errstr = 'Need to pass nonempty SaveDir for this routine.';
+    if ~exist('SaveDir', 'var')
+        error(errstr);
+    else
+        if isempty(SaveDir)
+            error(errstr);
+        end
+    end
+
+
+    sbase = ['!matlab -nodesktop -nosplash -singleCompThread'...
+         ' -r "BatchRun(%d, ''%s'')" &'];
+
+    for b = 1:1:length(BatchList)
+        s = sprintf(sbase, BatchList(b), SaveDir)
+        eval(s)
+        pause(10);
+    end
+end

run/Run.m

@@ -0,0 +1,85 @@
+%*******************************************************************************
+% Run.m
+%
+% Should generally pass the following (or else defaults will be used):
+%   Spec Number: 1 through 9 as in the paper
+%   N1: Number of support points for X1
+%   N2: Number of points for X2
+%   Continuous: Flag to treat X2 as continuous
+%
+% Optional:
+%   SaveDir: Name of the directory to save the simulation in
+%       If not passed then a name will be generated with a time stamp
+%*******************************************************************************
+function [SaveDir] = Run(SpecNumber, N1, N2, Continuous, SaveDir)
+
+addpath('../src');
+
+if ~exist('SaveDir', 'var');
+    SaveDir = [];
+end
+
+[CleanUpPath, SaveDir] = Setup('pies', SaveDir);
+
+if exist('N1', 'var')
+    Settings.N1 = N1;
+end
+if exist('N2', 'var')
+    Settings.N2 = N2;
+end
+if exist('Continuous', 'var')
+    Settings.Continuous = Continuous;
+end
+if ~exist('SpecNumber', 'var')
+    SpecNumber = 1;
+end
+
+switch SpecNumber
+    case 1
+        Assumptions.U1MedZeroGivenY2X1 = 1;
+
+    case 2
+        Assumptions.U1MedZeroGivenY2X1 = 1;
+        Assumptions.U1SymmetricGivenY2X1 = 1;
+
+    case 3
+        Assumptions.U1MedZeroGivenY2X1 = 1;
+        Assumptions.U1SymmetricGivenY2X1 = 1;
+        Assumptions.U1IndY2X1 = 1;
+
+    case 4
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+
+    case 5
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+        Assumptions.U1IndX1X2 = 1;
+
+    case 6
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+        Assumptions.U1IndX1X2 = 1;
+        Assumptions.U1SymmetricGivenX1X2 = 1;
+
+    case 7
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+        Assumptions.U1U2IndX1X2 = 1;
+
+    case 8
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+        Assumptions.U2MedZeroGivenX1X2 = 1;
+        Assumptions.U1U2IndX1X2 = 1;
+        Settings.ParametricFS = 1;
+
+    case 9
+        Assumptions.U1MedZeroGivenX1X2 = 1;
+        Assumptions.U1SymmetricGivenX1X2 = 1;
+        Assumptions.U2MedZeroGivenX1X2 = 1;
+        Assumptions.U2SymmetricGivenX1X2 = 1;
+        Assumptions.U1U2IndX1X2 = 1;
+        Settings.ParametricFS = 1;
+
+    otherwise
+        error('SpecNumber not recognized.');
+end
+
+Settings.SpecNumber = SpecNumber;
+IdentifiedSet(Settings, Assumptions);