From f9dfbbe35573d66566601f37063fb7289dbe1737 Mon Sep 17 00:00:00 2001 From: Alexander Torgovitsky Date: Tue, 21 Nov 2017 19:36:47 -0600 Subject: [PATCH] Initial commit --- README.md | 94 +++++++ cfg/Config.m | 2 + run/BatchRun.m | 57 +++++ run/MultiBatchRun.m | 28 ++ run/Run.m | 85 +++++++ src/BivBinResp.mod | 341 +++++++++++++++++++++++++ src/CreateNewResultsDir.m | 33 +++ src/GenerateDGP.m | 116 +++++++++ src/IdentifiedSet.m | 520 ++++++++++++++++++++++++++++++++++++++ src/InitializeAMPL.m | 21 ++ src/MergeBrackets.m | 94 +++++++ src/PlotID.m | 83 ++++++ src/PrintStructure.m | 35 +++ src/Setup.m | 106 ++++++++ src/UpdateStruct.m | 45 ++++ 15 files changed, 1660 insertions(+) create mode 100644 README.md create mode 100644 cfg/Config.m create mode 100644 run/BatchRun.m create mode 100644 run/MultiBatchRun.m create mode 100644 run/Run.m create mode 100644 src/BivBinResp.mod create mode 100644 src/CreateNewResultsDir.m create mode 100644 src/GenerateDGP.m create mode 100644 src/IdentifiedSet.m create mode 100644 src/InitializeAMPL.m create mode 100644 src/MergeBrackets.m create mode 100644 src/PlotID.m create mode 100644 src/PrintStructure.m create mode 100644 src/Setup.m create mode 100644 src/UpdateStruct.m diff --git a/README.md b/README.md new file mode 100644 index 0000000..ed51f44 --- /dev/null +++ b/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 diff --git a/cfg/Config.m b/cfg/Config.m new file mode 100644 index 0000000..2f23c1d --- /dev/null +++ b/cfg/Config.m @@ -0,0 +1,2 @@ +SAVEDIRECTORY = '~' +AMPLAPISETUPPATH = '/usr/local/amplapi/matlab/setUp.m' diff --git a/run/BatchRun.m b/run/BatchRun.m new file mode 100644 index 0000000..890dd79 --- /dev/null +++ b/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 diff --git a/run/MultiBatchRun.m b/run/MultiBatchRun.m new file mode 100644 index 0000000..dbc0814 --- /dev/null +++ b/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 diff --git a/run/Run.m b/run/Run.m new file mode 100644 index 0000000..cda64d9 --- /dev/null +++ b/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); diff --git a/src/BivBinResp.mod b/src/BivBinResp.mod new file mode 100644 index 0000000..6a0bcdc --- /dev/null +++ b/src/BivBinResp.mod @@ -0,0 +1,341 @@ +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# BivBinResp.mod +# +# The AMPL model file. +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +#******************************************************************************* +# Data +#******************************************************************************* +param N1, default 3; +param N2, default 3; +set K1 = 1..N1 ordered by Integers; +set K2 = 1..N2 ordered by Integers; +param X1 {k1 in K1} default 0; +param X2 {k2 in K2} default 0; + +# This is redundant but helps keep the indexing clear +# J1, J2 follows the indexing in the paper, while Y1 and Y2 +# are the values the internal variables take in bivariate binary response +set J1 = 1..2; +param Y1 {j1 in J1} = j1 - 1; +set J2 = 1..2; +param Y2 {j2 in J2} = j2 - 1; + +param FYGX {j1 in J1, j2 in J2, k1 in K1, k2 in K2} in [0,1] default .5; +param PX {k1 in K1, k2 in K2} in [0,1] default .5; +param PX2GX1 {k2 in K2, k1 in K1} default .5; + +#******************************************************************************* +# Derived sets/parameters +#******************************************************************************* +param beta0 default .5; +param beta1 default -.75; +param beta2 default 1; + +param g1 {j1 in ({0} union J1), j2 in J2, k1 in K1} = + if (j1 == 1) then (beta0 + beta1*Y2[j2] + beta2*X1[k1]) + else if (j1 == 0) then (-Infinity) + else if (j1 == 2) then (+Infinity); + +#******************************************************************************* +# With no first stage equation just set g2(x) = P[Y2 = 0 | X] +# Observational equivalence will then force +# F[+Infinity,u2,x] = P[Y2 = 0 | X] +# which forces the marginal of U2 to be uniform, independently of X +#******************************************************************************* +param pi0 default .3; +param pi1 default 1; +param pi2 default .2; + +param ParametricFS binary default 0; +param g2 {j2 in ({0} union J2), k1 in K1, k2 in K2} = + if (j2 == 1) then + (if (ParametricFS) then + (pi0 + pi1*X1[k1] + pi2*X2[k2]) + else + (FYGX[2,1,k1,k2]) + ) + else if (j2 == 0) then (-Infinity) + else if (j2 == 2) then (+Infinity); + +set Ug1 = + (setof {j1 in ({0} union J1), j2 in J2, k1 in K1} g1[j1,j2,k1]) + union (setof {j1 in ({0} union J1), j2 in J2, k1 in K1} -1*g1[j1,j2,k1]) + union {-Infinity, +Infinity, 0} + ordered by Reals; + +set Ug2 = + (setof {j2 in ({0} union J2), k1 in K1, k2 in K2} g2[j2,k1,k2]) + union (setof {j2 in ({0} union J2), k1 in K1, k2 in K2} -1*g2[j2,k1,k2]) + union {-Infinity, +Infinity,0} + ordered by Reals; + +#******************************************************************************* +# Variables +# +# F is the key variable of optimization +# The other variables below are derived from it +#******************************************************************************* +var F {u1 in Ug1, u2 in Ug2, k1 in K1, k2 in K2} in [0,1] default 0; + +var ASF0 = 1 - sum {k1 in K1, k2 in K2} ( + F[g1[1,1,k1], +Infinity, k1, k2]*PX[k1, k2]); +var ASF1 = 1 - sum {k1 in K1, k2 in K2} ( + F[g1[1,2,k1], +Infinity, k1, k2]*PX[k1, k2]); +var ATE = ASF1 - ASF0; + +param k1fix, default 1 in K1; +param k2fix, default 1 in K2; +var ASF0Fixed = 1 - F[g1[1,1,k1fix], +Infinity, k1fix, k2fix]; +var ASF1Fixed = 1 - F[g1[1,2,k1fix], +Infinity, k1fix, k2fix]; +var ATEFixed = ASF1Fixed - ASF0Fixed; + +#******************************************************************************* +# Objective functions +#******************************************************************************* +minimize Constant: 1; # For checking feasibility +minimize ObjMinATE: ATE; +maximize ObjMaxATE: ATE; +minimize ObjMinATEFixed: ATEFixed; +maximize ObjMaxATEFixed: ATEFixed; + +#******************************************************************************* +# Observational equivalence +#******************************************************************************* +var ObsEq_Gap {j1 in J1, j2 in J2, k1 in K1, k2 in K2} = + sum {jj2 in J2 : jj2 <= j2} ( + F[g1[j1,jj2,k1], g2[jj2,k1,k2], k1, k2] + - + F[g1[j1,jj2,k1], g2[jj2 - 1,k1,k2], k1, k2] + ) - FYGX[j1,j2,k1,k2]; +subject to ObsEq {j1 in J1, j2 in J2, k1 in K1, k2 in K2}: + ObsEq_Gap[j1,j2,k1,k2] = 0; + +#******************************************************************************* +# Subdistribution shape properties +#******************************************************************************* +subject to Grounded1 {u2 in Ug2, k1 in K1, k2 in K2}: + F[-Infinity,u2,k1,k2] = 0; +subject to Grounded2 {u1 in Ug1, k1 in K1, k2 in K2}: + F[u1,-Infinity,k1,k2] = 0; +subject to HasMargins {k1 in K1, k2 in K2}: + F[+Infinity,+Infinity,k1,k2] = 1; +subject to Increasing {t1 in 2..card(Ug1), t2 in 2..card(Ug2), + k1 in K1, k2 in K2}: + F[member(t1,Ug1),member(t2,Ug2),k1,k2] + - F[member(t1-1,Ug1),member(t2,Ug2),k1,k2] + - F[member(t1,Ug1),member(t2-1,Ug2),k1,k2] + + F[member(t1-1,Ug1),member(t2-1,Ug2),k1,k2] >= 0; + +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +# SINGLE EQUATION RESTRICTIONS +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* + +#******************************************************************************* +# U1 | Y2, X1 has median 0. Not using X2 at all here. +#******************************************************************************* +subject to U1MedZeroGivenY2X1 {j2 in J2, k1 in K1}: + (sum {k2 in K2} ( + (F[0,g2[j2,k1,k2],k1,k2] - F[0,g2[j2-1,k1,k2],k1,k2]) + *PX2GX1[k2,k1] + )) + / + (if (j2 == 2) + then (sum {k2 in K2} + ((1 - FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + else (sum {k2 in K2} + ((FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + ) + = .5; + +#******************************************************************************* +# U1 independent of (Y2, X1)---not using X2 here +# base group is Y2 = 0, X1 = min(X1) +#******************************************************************************* +subject to U1IndY2X1 {u1 in Ug1, j2 in J2, k1 in K1}: + (sum {k2 in K2} ( + (F[u1,g2[j2,k1,k2],k1,k2] - F[u1,g2[j2-1,k1,k2],k1,k2]) + *PX2GX1[k2,k1] + )) + / + (if (j2 == 2) + then (sum {k2 in K2} + ((1 - FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + else (sum {k2 in K2} + ((FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + ) + = + (sum {k2 in K2} ( + F[u1,g2[1,member(1,K1),k2],member(1,K1),k2] + *PX2GX1[k2,member(1,K1)] + )) + / + (sum {k2 in K2} ( + FYGX[2,1,member(1,K1),k2]*PX2GX1[k2,member(1,K1)] + )); + +#******************************************************************************* +# Imposes Symmetry around 0 for U1 given Y2, X1 +# +# P[U1 <= u1 | Y2, X1] = 1 - P[U1 <= -u1 | Y2, X1] +#******************************************************************************* +subject to U1SymmetricGivenY2X1 {u1 in Ug1, j2 in J2, k1 in K1}: + (sum {k2 in K2} ( + (F[u1,g2[j2,k1,k2],k1,k2] - F[u1,g2[j2-1,k1,k2],k1,k2]) + *PX2GX1[k2,k1] + )) + / + (if (j2 == 2) + then (sum {k2 in K2} + ((1 - FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + else (sum {k2 in K2} + ((FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + ) + = + 1 - + (sum {k2 in K2} ( + (F[-u1,g2[j2,k1,k2],k1,k2] - F[-u1,g2[j2-1,k1,k2],k1,k2]) + *PX2GX1[k2,k1] + )) + / + (if (j2 == 2) + then (sum {k2 in K2} + ((1 - FYGX[2,1,k1,k2])*PX2GX1[k2,k1]) + ) + else (sum {k2 in K2} + (FYGX[2,1,k1,k2]*PX2GX1[k2,k1]) + ) + ); + +#******************************************************************************* +# This imposes U1 | Y2, X1, X2 has median 0 +# +# This is +# P[U1 <= 0 | Y2 = j2, X1 = k1, X2 = k2] +# = +# (P[U1<=0, U2<=g(j2,k1,k2) | k1,k2] - P[U1<=0, U2<=g(j2-1,k1,k2)|k1,k2]) +# / +# (P[Y2 <= j2 | X1 = k1, X2 = k2] - P[Y2 <= j2 - 1 | X1 = k1, X2 = k2]) +# = +# 1/2 +#******************************************************************************* +subject to U1MedZeroGivenY2X1X2 {j2 in J2, k1 in K1, k2 in K2}: + (F[0,g2[j2,k1,k2],k1,k2] - F[0,g2[j2-1,k1,k2],k1,k2]) + / + (if (j2 == 2) then (1 - FYGX[2,1,k1,k2]) + else (FYGX[2,1,k1,k2])) + = .5; + +#******************************************************************************* +# This imposes U1 is independent of Y2, X1, X2 +# +# The base group here is Y2 = 0, X1 = min(X1), X2 = min(X2); +#******************************************************************************* +subject to U1IndY2X1X2 {u1 in Ug1, j2 in J2, k1 in K1, k2 in K2}: + (F[u1,g2[j2,k1,k2],k1,k2] - F[u1,g2[j2-1,k1,k2],k1,k2]) + / + (if (j2 == 2) then (1 - FYGX[2,1,k1,k2]) + else (FYGX[2,1,k1,k2])) + = + F[u1,g2[1,member(1,K1),member(1,K2)],member(1,K1),member(1,K2)] + / + FYGX[2,1,member(1,K1),member(1,K2)]; + +#******************************************************************************* +# This imposes U1 is independent of Y2 given X1,X2 +# +# Base group here is Y2 = 0 +#******************************************************************************* +subject to U1IndY2GivenX1X2 {u1 in Ug1, j2 in J2, k1 in K1, k2 in K2}: + (F[u1,g2[j2,k1,k2],k1,k2] - F[u1,g2[j2-1,k1,k2],k1,k2]) + / + (if (j2 == 2) then (1 - FYGX[2,1,k1,k2]) + else (FYGX[2,1,k1,k2])) + = + F[u1,g2[1,k1,k2],k1,k2]/FYGX[2,1,k1,k2]; + +#******************************************************************************* +# Imposes Symmetry around 0 for U1 given Y2, X1, X2 +# +# P[U1 <= u1 | Y2, X1, X2] = 1 - P[U1 <= -u1 | Y2, X1, X2] +#******************************************************************************* +subject to U1SymmetricGivenY2X1X2 {u1 in Ug1, j2 in J2, k1 in K1, k2 in K2}: + (F[u1,g2[j2,k1,k2],k1,k2] - F[u1,g2[j2-1,k1,k2],k1,k2]) + / + (if (j2 == 2) then (1 - FYGX[2,1,k1,k2]) + else (FYGX[2,1,k1,k2])) + = + 1 - + (F[-u1,g2[j2,k1,k2],k1,k2] - F[-u1,g2[j2-1,k1,k2],k1,k2]) + / + (if (j2 == 2) then (1 - FYGX[2,1,k1,k2]) + else (FYGX[2,1,k1,k2])); + +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +# SINGLE EQUATION INSTRUMENT RESTRICTIONS +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* + +#******************************************************************************* +# This is just P[U1 <= 0 | X1, X2] = .5 +# But not conditioning on Y2 as above +#******************************************************************************* +subject to U1MedZeroGivenX1X2 {k1 in K1, k2 in K2}: + F[0,+Infinity,k1,k2] = .5; + +#******************************************************************************* +# U1 independent of (X1,X2) +#******************************************************************************* +subject to U1IndX1X2 {u1 in Ug1, k1 in K1, k2 in K2}: + F[u1,+Infinity,k1,k2] = F[u1,+Infinity,member(1,K1),member(1,K2)]; + +#******************************************************************************* +# U1 given (X1,X2) is symmetric around 0 +#******************************************************************************* +subject to U1SymmetricGivenX1X2 {u1 in Ug1, k1 in K1, k2 in K2}: + F[u1,+Infinity,k1,k2] = 1 - F[-u1,+Infinity,k1,k2]; + +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +# TRIANGULAR MODEL RESTRICTIONS +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* +#******************************************************************************* + +#******************************************************************************* +# (U1,U2) independent of (X1,X2) +#******************************************************************************* +subject to U1U2IndX1X2 {u1 in Ug1, u2 in Ug2, k1 in K1, k2 in K2}: + F[u1,u2,k1,k2] = F[u1,u2,member(1,K1),member(1,K2)]; + +#******************************************************************************* +# U2 median 0 conditional on (X1,X2) +#******************************************************************************* +subject to U2MedZeroGivenX1X2 {k1 in K1, k2 in K2}: + F[+Infinity,0,k1,k2] = .5; + +#******************************************************************************* +# U2 symmetric around 0 conditional on (X1,X2) +#******************************************************************************* +subject to U2SymmetricGivenX1X2 {u2 in Ug2, k1 in K1, k2 in K2}: + F[+Infinity, u2, k1, k2] = 1 - F[+Infinity, -u2, k1, k2]; diff --git a/src/CreateNewResultsDir.m b/src/CreateNewResultsDir.m new file mode 100644 index 0000000..a4d5348 --- /dev/null +++ b/src/CreateNewResultsDir.m @@ -0,0 +1,33 @@ +%******************************************************************************* +% CreateNewResultsDir +% +% Create a new results directory (numbered sequentially) and cd to it. +%******************************************************************************* +function [DirName] = CreateNewResultsDir(DirName) + CurrentDirectory = pwd; + [path namestr] = fileparts(CurrentDirectory); + + % Make sure we are either in results or some directory directly below it + if ~strcmp(namestr, 'results') + cd('..'); + [path namestr] = fileparts(pwd); + if ~strcmp(namestr, 'results') + error('Something is wrong: not results or a subdirectory?'); + end + end + + DirNumber = []; + if ~exist('DirName', 'var') + for i = 0:1:999 + DirName = sprintf('%03d', i + 1); + if ~exist(DirName, 'dir') + DirNumber = i; + break; + end + end + end + + mkdir(DirName); + cd(DirName); + addpath(genpath('../../')); +end diff --git a/src/GenerateDGP.m b/src/GenerateDGP.m new file mode 100644 index 0000000..0113ec0 --- /dev/null +++ b/src/GenerateDGP.m @@ -0,0 +1,116 @@ +%******************************************************************************* +% GenerateDGP.m +% +% Generate DGP quantities from Settings +%******************************************************************************* +function DGP = GenerateDGP(Settings) + %*************************************************************************** + % Distribution of X1 is always uniform discrete + % Distribution of X2 can either be uniform discrete or continuous + % If continuous, then uniform over [-2, 2] + % In either case they are taken to be independent + %*************************************************************************** + DGP.X1 = ... + ((1:1:Settings.N1)' + floor(Settings.N1/2) - Settings.N1)... + /((Settings.N1 - 1)/2); + + if Settings.Continuous + QSpacing = 1/(Settings.N2+1); + QList = QSpacing:QSpacing:(1 - QSpacing); + DGP.X2 = unifinv(QList, -2, 2); + + % Make sure (0,0) is in the support + DGP.X2 = [DGP.X2(:); 0]; + DGP.X2 = unique(DGP.X2); + DGP.X1 = [DGP.X1(:); 0]; + DGP.X1 = unique(DGP.X1); + else + DGP.X2 = (1:1:Settings.N2)' + floor(Settings.N2/2) - Settings.N2; + end + % Joint probability of (X1, X2) + DGP.PX = (1/length(DGP.X1))*(1/length(DGP.X2))... + *ones(length(DGP.X1), length(DGP.X2)); + + % Probability of X2 given X1 + DGP.PX2GX1 = nan(length(DGP.X2), length(DGP.X1)); + for k2 = 1:1:length(DGP.X2) + for k1 = 1:1:length(DGP.X1) + DGP.PX2GX1(k2,k1) = DGP.PX(k1,k2)/sum(DGP.PX(k1,:)); + end + end + + %*************************************************************************** + % Generate conditional probabilities of Y observed in the data + %*************************************************************************** + DGP.g1 = @(y2, x1) Settings.Beta0 + Settings.Beta1*y2 + Settings.Beta2*x1; + for j2 = 1:1:2 + for k1 = 1:1:length(DGP.X1) + DGP.PrU1LTg1(j2, k1) = ... + normcdf(DGP.g1(j2-1,DGP.X1(k1)), 0, Settings.Sigma1); + end + end + + DGP.g2 = @(x1, x2) Settings.Pi0 + Settings.Pi1*x1 + Settings.Pi2*x2; + for k1 = 1:1:length(DGP.X1) + for k2 = 1:1:length(DGP.X2) + DGP.PrU2LTg2(k1,k2) = ... + normcdf(DGP.g2(DGP.X1(k1), DGP.X2(k2)), 0, ... + Settings.Sigma2); + end + end + + DGP.PrYGX = nan(2,2,length(DGP.X1), length(DGP.X2)); + DGP.FYGX = nan(2,2,length(DGP.X1), length(DGP.X2)); + for k1 = 1:1:length(DGP.X1) + for k2 = 1:1:length(DGP.X2) + % Clayton copula with parameter theta + % Taking theta = 1 corresponds to independence + DGP.PrYGX(1,1,k1,k2) = ... + exp(-1*(... + (-log(DGP.PrU1LTg1(1,k1)))^Settings.Theta ... + + ... + (-log(DGP.PrU2LTg2(k1,k2)))^Settings.Theta ... + )^(1/Settings.Theta)); + + DGP.PrYGX(1,2,k1,k2) = DGP.PrU1LTg1(2, k1) - ... + exp(-1*(... + (-log(DGP.PrU1LTg1(2,k1)))^Settings.Theta ... + + ... + (-log(DGP.PrU2LTg2(k1,k2)))^Settings.Theta ... + )^(1/Settings.Theta)); + + DGP.PrYGX(2,1,k1,k2) = DGP.PrU2LTg2(k1, k2) - ... + exp(-1*(... + (-log(DGP.PrU1LTg1(1,k1)))^Settings.Theta ... + + ... + (-log(DGP.PrU2LTg2(k1,k2)))^Settings.Theta ... + )^(1/Settings.Theta)); + + DGP.PrYGX(2,2,k1,k2) = 1 - DGP.PrYGX(1,1,k1,k2) ... + - DGP.PrYGX(1,2,k1,k2) - DGP.PrYGX(2,1,k1,k2); + + DGP.FYGX(1,1,k1,k2) = DGP.PrYGX(1,1,k1,k2); + DGP.FYGX(1,2,k1,k2) = DGP.PrYGX(1,1,k1,k2) + DGP.PrYGX(1,2,k1,k2); + DGP.FYGX(2,1,k1,k2) = DGP.PrYGX(1,1,k1,k2) + DGP.PrYGX(2,1,k1,k2); + DGP.FYGX(2,2,k1,k2) = 1; + end + end + + % Value of the ASF and ATE in the DGP + DGP.ASF = ones(2,1); + if ~Settings.Continuous + for k1 = 1:1:length(DGP.X1) + for k2 = 1:1:length(DGP.X2) + for j2 = 1:1:2 + DGP.ASF(j2) = DGP.ASF(j2) - ... + normcdf(DGP.g1(j2-1, DGP.X1(k1)),... + 0, Settings.Sigma1)*DGP.PX(k1,k2); + end + end + end + else + DGP.ASF(1) = 1 - normcdf(DGP.g1(0, 0), 0, Settings.Sigma1); + DGP.ASF(2) = 1 - normcdf(DGP.g1(1, 0), 0, Settings.Sigma1); + end + DGP.ATE = DGP.ASF(2) - DGP.ASF(1); +end diff --git a/src/IdentifiedSet.m b/src/IdentifiedSet.m new file mode 100644 index 0000000..fce7dc0 --- /dev/null +++ b/src/IdentifiedSet.m @@ -0,0 +1,520 @@ +%******************************************************************************* +% IdentifiedSet.m +% +% The main routine. +%******************************************************************************* +function [] = IdentifiedSet(SettingsIn, AssumptionsIn) + +%******************************************************************************* +% Define default settings +%******************************************************************************* +Settings.N1 = 3; +Settings.N2 = 3; +Settings.Continuous = 0; +Settings.Beta0 = .5; +Settings.Beta1 = -.75; +Settings.Beta2 = 1; +Settings.Pi0 = .3; +Settings.Pi1 = 1; +Settings.Pi2 = .2; +Settings.Sigma1 = 1; +Settings.Sigma2 = 1; +Settings.Theta = 1; % Parameter of the Gumbel Copula (1 = independence) +Settings.ParametricFS = 0; + +Settings.Beta0GridLB = -1.5; +Settings.Beta0GridUB = 3; +Settings.Beta1GridLB = -3; +Settings.Beta1GridUB = 2; +Settings.Pi0GridLB = -.5; +Settings.Pi0GridUB = 1.5; +Settings.Pi2GridLB = -1.5; +Settings.Pi2GridUB = 1.5; + +%Settings.Beta0GridStep = .5; % Testing +%Settings.Beta1GridStep = .5; +%Settings.Pi0GridStep = .2; +%Settings.Pi2GridStep = .2; + +%Settings.Beta0GridStep = .2; % Coarse Beta +%Settings.Beta1GridStep = .2; +%Settings.Pi0GridStep = .2; +%Settings.Pi2GridStep = .2; + +Settings.Beta0GridStep = .05; % Production +Settings.Beta1GridStep = .05; +Settings.Pi0GridStep = .05; +Settings.Pi2GridStep = .05; + +Settings.SpecNumber = []; +Settings.LPSolver = 'gurobi'; +Settings.Noisy = 0; +Settings.VisualDisplayEachIter = 1; +Settings.RecordEachColumn = 1; +Settings.VisualDisplayEnd = 1; + +Settings.Grayscale = .7; +Settings.FontsizeTitle = 20; +Settings.FontsizeLabel = 20; +Settings.FontsizeAxis = 20; +Settings.FontSizeTickLabels = 6; +Settings.MarkersizeBig = 7; +Settings.MarkersizeSmall = 2; +Settings.FigHeight = 7; +Settings.FigWidth = 6; + +Assumptions.U1MedZeroGivenY2X1 = 0; +Assumptions.U1IndY2X1 = 0; +Assumptions.U1SymmetricGivenY2X1 = 0; +Assumptions.U1MedZeroGivenY2X1X2 = 0; +Assumptions.U1IndY2X1X2 = 0; +Assumptions.U1IndY2GivenX1X2 = 0; +Assumptions.U1SymmetricGivenY2X1X2 = 0; +Assumptions.U1MedZeroGivenX1X2 = 0; +Assumptions.U1IndX1X2 = 0; +Assumptions.U1SymmetricGivenX1X2 = 0; +Assumptions.U1U2IndX1X2 = 0; +Assumptions.U2MedZeroGivenX1X2 = 0; +Assumptions.U2SymmetricGivenX1X2 = 0; + +if isempty(SettingsIn) + warning('Empty structure for input. Using defaults.'); + SettingsIn = Settings; +end +if ~isstruct(SettingsIn) + error('Expected structure for input. Quitting.'); +else + Settings = UpdateStruct(Settings, SettingsIn, 1); +end + +Assumptions = UpdateStruct(Assumptions, AssumptionsIn, 1); + + +%******************************************************************************* +% Different grid settings for the continuous case +%******************************************************************************* +if Settings.Continuous + Settings.Beta0GridLB = -1.5; + Settings.Beta0GridUB = 2; + Settings.Beta1GridLB = -2; + Settings.Beta1GridUB = 1.5; + + Settings.Beta0GridStep = .1; + Settings.Beta1GridStep = Settings.Beta0GridStep; +end + +%******************************************************************************* +% Prepare to record +%******************************************************************************* +if ~isempty(Settings.SpecNumber) + Settings.DirName = [ sprintf('Spec-%03d', Settings.SpecNumber) ... + '-X1-' sprintf('%d', Settings.N1) ... + '-X2-' sprintf('%d', Settings.N2)]; + if Settings.Continuous + Settings.DirName = [Settings.DirName '-continuous']; + end + CreateNewResultsDir(Settings.DirName); +else + Settings.DirName = CreateNewResultsDir; +end +fid = fopen('Settings.out', 'w'); +PrintStructure(Settings, fid); +fclose(fid); +fid = fopen('Assumptions.out', 'w'); +PrintStructure(Assumptions, fid); +fclose(fid); +diary('Log.out'); + +%******************************************************************************* +% Quit if this simulation has already been done +% (useful for batching) +%******************************************************************************* +DONEFN = './Done.out'; +if exist(DONEFN, 'file') == 2 + disp(sprintf('Found %s, so quitting.', DONEFN)); + return; +end + +STATUSFN = './Status.out'; +if exist(STATUSFN, 'file') == 2 + disp(sprintf(['Found %s, so quitting. ' ... + 'To overwrite, delete this file and retry.'], STATUSFN)); + return; +end + +%******************************************************************************* +% Initialize AMPL instance +% Set some solver options +%******************************************************************************* +[ampl, CleanUpAMPL] = InitializeAMPL({'BivBinResp.mod'}); + +ampl.setOption('TMPDIR', './'); +ampl.setOption('presolve_eps', '1e-5'); +ampl.setOption('presolve_warnings', '-1'); +ampl.setOption('solver_msg', '0'); +ampl.setOption('print_precision', '6'); +ampl.setOption('show_stats', '0'); +ampl.setOption('show_boundtol', '0'); +ampl.setOption('solver', Settings.LPSolver); + +%******************************************************************************* +% Send data parameters to AMPL +%******************************************************************************* +DGP = GenerateDGP(Settings); + +IdxX1X2 = []; +ValPX = []; +IdxX2X1 = []; +ValPX2GX1 = []; +IdxYX = []; +ValFYGX = []; +for k1 = 1:1:length(DGP.X1) + for k2 = 1:1:length(DGP.X2) + IdxX1X2 = [IdxX1X2; k1 k2]; + ValPX = [ValPX; DGP.PX(k1,k2)]; + IdxX2X1 = [IdxX2X1; k2 k1]; + ValPX2GX1 = [ValPX2GX1; DGP.PX2GX1(k2,k1)]; + + for j1 = 1:1:2 + for j2 = 1:1:2 + IdxYX = [IdxYX; j1 j2 k1 k2]; + ValFYGX = [ValFYGX; DGP.FYGX(j1,j2,k1,k2)]; + end + end + end +end + +ampl.getParameter('N1').setValues(length(DGP.X1)); +ampl.getParameter('N2').setValues(length(DGP.X2)); +ampl.getParameter('X1').setValues((1:1:length(DGP.X1))', DGP.X1(:)); +ampl.getParameter('X2').setValues((1:1:length(DGP.X2))', DGP.X2(:)); +ampl.getParameter('PX').setValues(IdxX1X2, ValPX); +ampl.getParameter('PX2GX1').setValues(IdxX2X1, ValPX2GX1); +ampl.getParameter('FYGX').setValues(IdxYX, ValFYGX); + +[m1 idx1] = min(abs(DGP.X1)); +ampl.getParameter('k1fix').setValues(idx1); +[m2 idx2] = min(abs(DGP.X2)); +ampl.getParameter('k2fix').setValues(idx2); +if Settings.Continuous + if (m1 > 0) | (m2 > 0) + error('Something is wrong -- should be estimated at (0,0).') + end +end + +fidATEDGP = fopen('ATEDGP.out', 'w'); +fprintf(fidATEDGP, '%6.4f', DGP.ATE); +fclose(fidATEDGP); + +%******************************************************************************* +% Turn on the correct assumptions in AMPL +%******************************************************************************* +AssumptionNames = fieldnames(Assumptions); +for i = 1:1:length(AssumptionNames) + if Assumptions.(AssumptionNames{i}) + ampl.getConstraint(AssumptionNames{i}).restore; + else + ampl.getConstraint(AssumptionNames{i}).drop; + end +end + +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* +% Construct identified set +% +% Note that the matrices I am storing results in are set up in a graphical way +% So b0/p0 increases with the column, and b1/p2 decreases with the row +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* + +%******************************************************************************* +% If using a parametric model for the first stage, then create an +% identified set for its parameters. +% One can then build the identified set for the outcome equation at +% each fixed value of these parameters. +% Should be a bit faster than doing the whole 4 dimensional grid +%******************************************************************************* +disp('Suspending diary for identified set construction.') +diary off; +if Settings.ParametricFS + ampl.getParameter('ParametricFS').setValues(1); + + Settings.Pi0Grid = ... + Settings.Pi0GridLB:Settings.Pi0GridStep:Settings.Pi0GridUB; + Settings.Pi2Grid = ... + Settings.Pi2GridUB:(-1*Settings.Pi2GridStep):Settings.Pi2GridLB; + IDSetPi = 2*ones(length(Settings.Pi2Grid), length(Settings.Pi0Grid)); + Str = {'Identified set for (pi0, pi2):'}; + + % Set (beta0, beta1) to their true values so that the outcome equation + % is not the cause of being inside or outside of the identified set. + % This is just for ease in programming. An alternative would be to + % just solve a single equation model with only the first stage. + ampl.getParameter('beta0').setValues(Settings.Beta0); + ampl.getParameter('beta1').setValues(Settings.Beta1); + + for p0 = 1:1:length(Settings.Pi0Grid) + for p2 = 1:1:length(Settings.Pi2Grid) + ampl.getParameter('pi0').setValues(Settings.Pi0Grid(p0)); + ampl.getParameter('pi2').setValues(Settings.Pi2Grid(p2)); + + Obj = SolveProgram(ampl, 'Constant', Settings.Noisy); + + idstr = sprintf('Pi0 = %4.3f, Pi2 = %4.3f',... + Settings.Pi0Grid(p0), Settings.Pi2Grid(p2)); + IDSetPi(p2,p0) = CheckFeasibility(Obj, idstr); + + if Settings.VisualDisplayEachIter + VisualDisplay(IDSetPi, 1, [], Str); + end + end + end + + [IDSetPi2Idx IDSetPi0Idx] = find(IDSetPi == 1); + IDListPi = zeros(size(IDSetPi2Idx, 1), 2); + for p = 1:1:size(IDSetPi2Idx, 1) + IDListPi(p,:) = [Settings.Pi0Grid(IDSetPi0Idx(p)), ... + Settings.Pi2Grid(IDSetPi2Idx(p))]; + end +else + % Also for coding convenience --- knowing these is not necessary + IDListPi = [Settings.Pi0, Settings.Pi2]; +end + +diary on; +if size(IDListPi, 1) == 0 + disp(['Identified set for Pi is empty. '... + 'Maybe the grid is too coarse? '... + 'Nothing to do, so quitting.']); + WriteDone(DONEFN); + return; +end + +Settings.Beta0Grid = ... + Settings.Beta0GridLB:Settings.Beta0GridStep:Settings.Beta0GridUB; +Settings.Beta1Grid = ... + Settings.Beta1GridUB:(-1*Settings.Beta1GridStep):Settings.Beta1GridLB; + +IDSetBeta = 2*ones( length(Settings.Beta1Grid), length(Settings.Beta0Grid),... + size(IDListPi, 1)); +MinATE = inf(size(IDSetBeta)); +MaxATE = -1*MinATE; + +fidTime = fopen('SolveTime.out', 'w'); + +Str{1} = 'Identified set for (beta0, beta1):'; +diary off; +fidStatus = fopen(STATUSFN, 'w'); +StatusStr = ''; +for p = 1:1:size(IDListPi, 1) + ampl.getParameter('pi0').setValues(IDListPi(p,1)); + ampl.getParameter('pi2').setValues(IDListPi(p,2)); + if Settings.ParametricFS + Str{2} = sprintf('Fixing Pi0 = %4.3f, Pi2 = %4.3f',... + IDListPi(p,1), IDListPi(p,2)); + Str{2} = [Str{2} sprintf(' -- point %d of %d.', p, size(IDSetBeta,3))]; + else + Str{2} = []; + end + + for b0 = 1:1:length(Settings.Beta0Grid) + for b1 = 1:1:length(Settings.Beta1Grid) + frewind(fidStatus); + StatusStr = ... + sprintf('%s: Beta point %d of %d, for Pi point %d of %d.\n',... + datestr(datetime('now')),... + (b0-1)*length(Settings.Beta1Grid) + b1,... + length(Settings.Beta0Grid)*length(Settings.Beta1Grid),... + p, size(IDListPi, 1)); + fprintf(fidStatus, StatusStr); + + SolveTime = nan(3,1); + + ampl.getParameter('beta0').setValues(Settings.Beta0Grid(b0)); + ampl.getParameter('beta1').setValues(Settings.Beta1Grid(b1)); + + % Feasibility first + [Obj, SolveTime(1)] = SolveProgram(ampl, 'Constant', Settings.Noisy); + + idstr = sprintf('Beta0 = %4.3f, Beta1 = %4.3f',... + Settings.Beta0Grid(b0), Settings.Beta1Grid(b1)); + IDSetBeta(b1,b0,p) = CheckFeasibility(Obj, idstr); + + % If feasible then do the minimization and maximization problems + if IDSetBeta(b1,b0,p) + + if Settings.Continuous + ObjName = 'ObjMinATEFixed'; + else + ObjName = 'ObjMinATE'; + end + [Obj, SolveTime(2)] = SolveProgram(ampl, ObjName, Settings.Noisy); + MinATE(b1,b0,p) = Obj.value; + + if Settings.Continuous + ObjName = 'ObjMaxATEFixed'; + else + ObjName = 'ObjMaxATE'; + end + [Obj, SolveTime(3)] = SolveProgram(ampl, ObjName, Settings.Noisy); + MaxATE(b1,b0,p) = Obj.value; + end + + if Settings.VisualDisplayEachIter + VisualDisplay(IDSetBeta(:,:,p), 1, SolveTime, Str); + end + + fprintf(fidTime, '%5.3f %5.3f %5.3f\n', SolveTime(:)); + end + end +end +diary on; +fclose(fidStatus); +fclose(fidTime); + +Str = {'Final identified set for (beta0, beta1):'}; +if (Settings.VisualDisplayEnd & ~Settings.VisualDisplayEachIter)... + | (Settings.VisualDisplayEnd & Settings.ParametricFS) + VisualDisplay(max(IDSetBeta, [], 3), 0, [], Str); +end + +%******************************************************************************* +% Record in list format +%******************************************************************************* +fid = fopen('BoundsList.out', 'w'); +ATEList = []; +for p = 1:1:size(IDListPi, 1) + [IDSetBeta1Idx IDSetBeta0Idx] = find(IDSetBeta(:,:,p) == 1); + + if Settings.ParametricFS + fprintf(fid, 'pi0,pi2,beta0,beta1,atemin,atemax\n'); + else + fprintf(fid, 'beta0,beta1,atemin,atemax\n'); + end + + BoundsList = zeros(size(IDSetBeta1Idx, 1), 4); + + for b = 1:1:size(IDSetBeta1Idx, 1) + BoundsList(b,:) = [ Settings.Beta0Grid(IDSetBeta0Idx(b)),... + Settings.Beta1Grid(IDSetBeta1Idx(b)),... + MinATE(IDSetBeta1Idx(b), IDSetBeta0Idx(b)),... + MaxATE(IDSetBeta1Idx(b), IDSetBeta0Idx(b))]; + + if Settings.ParametricFS + fprintf(fid, '%4.3f,%4.3f,', IDListPi(p,:)); + end + fprintf(fid, '%4.3f,%4.3f,%4.3f,%4.3f\n', BoundsList(b,:)); + + ATEList = [ATEList; BoundsList(b,3:4)]; + end +end +fclose(fid); + +%******************************************************************************* +% Determine the identified set for the ATE +% being careful to look for disconnected sets +%******************************************************************************* +if ~isempty(ATEList) + fid = fopen('IDSetATE.out', 'w'); + ATEList = unique(round(ATEList, 6), 'rows'); + [MergedLB MergedUB] = MergeBrackets(ATEList(:,1), ATEList(:,2)); + for m = 1:1:length(MergedLB) + fprintf(fid, '%8.6f %8.6f\n', MergedLB(m), MergedUB(m)); + end + fclose(fid); +end + +%******************************************************************************* +% Make plots for this simulation +%******************************************************************************* +PlotID(max(IDSetBeta, [], 3), MinATE, MaxATE, Settings, 1); +PlotID(max(IDSetBeta, [], 3), MinATE, MaxATE, Settings, 0); + +%******************************************************************************* +% Record this simulation as done +%******************************************************************************* +WriteDone(DONEFN); +diary off; +end + +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* +% Change objective and solve quietly +function [Obj, SolveTime] = SolveProgram(ampl, ObjectiveName, Noisy) + + ampl.eval(['objective ' ObjectiveName ';']); + OptStart = tic; + + if Noisy + eval('ampl.eval(''solve;'')'); + else + evalc('ampl.eval(''solve;'')'); + end + + SolveTime = toc(OptStart); + Obj = ampl.getObjective(ObjectiveName); +end + +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* +% Recover exit code of objective and check if feasible or solved +function Feasible = CheckFeasibility(Obj, idstr) + + SolveResult = char(Obj.result); + + if ismember(SolveResult, {'solved', 'solved?'}) + Feasible = 1; + elseif ismember(SolveResult, {'infeasible', 'infeasible?'}) + Feasible = 0; + else + warning('Problem in %s: SolveResult is %s', idstr, SolveResult); + Feasible = 3; + end + +end + +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* +% Display the current known identified set on screen +function VisualDisplay(IDSet, Clear, SolveTime, Str) + if Clear + clc; + else + disp(sprintf(' ')); % Adds a new line for space + end + disp(pwd); + VisualCode = '01?!'; + IDSetPrint = VisualCode(IDSet + ones(size(IDSet))); + if exist('Str', 'var') + for j = 1:1:length(Str) + disp(Str{j}); + end + end + disp(IDSetPrint); + + if ~exist('SolveTime', 'var') + SolveTime = []; + end + if ~isempty(SolveTime) + str = sprintf('Solve times (in seconds):\n'); + str = [str sprintf(' Feasibility: %5.3f\n', SolveTime(1))]; + str = [str sprintf(' Min: %5.3f\n', SolveTime(2))]; + str = [str sprintf(' Max: %5.3f\n', SolveTime(3))]; + disp(str); + end +end + +%******************************************************************************* +%******************************************************************************* +%******************************************************************************* +function WriteDone(DONEFN) + fiddone = fopen(DONEFN, 'w'); + fprintf(fiddone, '-'); + fclose(fiddone); + disp('Simulation done.'); +end diff --git a/src/InitializeAMPL.m b/src/InitializeAMPL.m new file mode 100644 index 0000000..86b50af --- /dev/null +++ b/src/InitializeAMPL.m @@ -0,0 +1,21 @@ +%******************************************************************************* +% InitializeAMPL.m +% +% Start an AMPL instance and create a destructor that will automatically kill +% the instance when it is cleared. +%******************************************************************************* +function [ampl, CleanUpAMPL] = InitializeAMPL(ModelFiles) + ampl = AMPL; + + CleanUpAMPL = onCleanup(@()ampl.close()); + + for m = 1:1:length(ModelFiles) + if ~exist(ModelFiles{m}, 'file') + error('Model file %d could not be found.', m); + else + % The exist looks through the path--but need to translate for AMPL + fn = which(ModelFiles{m}); + ampl.read(fn); + end + end +end diff --git a/src/MergeBrackets.m b/src/MergeBrackets.m new file mode 100644 index 0000000..6cefc44 --- /dev/null +++ b/src/MergeBrackets.m @@ -0,0 +1,94 @@ +%Copyright (c) 2009, Bruno Luong +%All rights reserved. + +%Redistribution and use in source and binary forms, with or without +%modification, are permitted provided that the following conditions are +%met: + +%* Redistributions of source code must retain the above copyright +%notice, this list of conditions and the following disclaimer. +%* Redistributions in binary form must reproduce the above copyright +%notice, this list of conditions and the following disclaimer in +%the documentation and/or other materials provided with the distribution + +%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +%POSSIBILITY OF SUCH DAMAGE. + + +function [lower upper] = MergeBrackets(left, right) +% function [lower upper] = MergeBrackets(left, right) +% +% Purpose: Interval merging +% +% Given N input closed intervals in braket form: +% Ii := [left(i),right(i)], i = 1,2...,N (mathematical notation) +% The set union{Ii) can be written as a canonical partition by +% intervals Jk; i.e., union{Ii) = union(Jk), where Jk are M intervals +% (with M<=N, so the partition is minimum cardinal), and {Jk} are +% disjoint to each other (their intersections are empty). This function +% returns Jk = [lower(k),upper(k)], k=1,2,...M, in the ascending sorted +% order. +% +% EXAMPLE USAGE: +% >> [lower upper] = MergeBrackets([0 1 2 3 4],[1.5 1.6 3.5 3 5]) +% lower = 0 2 4 +% upper = 1.6 3.5 5 +% +% Algorithm complexity: O(N*log(N)) +% +% Author: Bruno Luong +% Original: 25-May-2009 + +% Detect when right < left (empty Ii), and later remove it (line #29, 30) +notempty = find(right>=left); + +% sort the rest by left bound +[left iorder] = sort(left(notempty)); +right = right(notempty(iorder)); + +% Allocate, as we don't know yet the size, we assume the largest case +lower = zeros(size(left)); +upper = zeros(size(right)); + +% Nothing to do +if isempty(lower) + return +end + +% Initialize +l = left(1); +u = right(1); +k = 0; +% Loop on brakets +for i=1:length(left) + if left(i) > u % new Jk detected + % Stack the old one + k = k+1; + lower(k) = l; + upper(k) = u; + % Reset l and u + l = left(i); + u = right(i); + else + u = max(u, right(i)); + end +end % FOR loop +% Stack the last one +k = k+1; +lower(k) = l; +upper(k) = u; + +% Remove the tails +lower(k+1:end) = []; +upper(k+1:end) = []; + +end % MergeBrackets diff --git a/src/PlotID.m b/src/PlotID.m new file mode 100644 index 0000000..00c7cfa --- /dev/null +++ b/src/PlotID.m @@ -0,0 +1,83 @@ +%******************************************************************************* +% PlotID.m +% +% Create plots of the identified set +%******************************************************************************* +function PlotID(IDSetBeta, MinATE, MaxATE, Settings, Type) + close; + figure('visible', 'off'); + + if (Type == 1) % Grayscale + colormap([1 1 1; Settings.Grayscale*ones(1,3)]); + contourf(Settings.Beta0Grid, Settings.Beta1Grid, IDSetBeta,... + 'LineColor', 'none'); + + Title = 'Identified Set (Shaded Area)'; + Filename = 'IDSetBeta'; + else % With colored partitions + MinATE = round(MinATE,3); + MinATE(abs(MinATE) == Inf) = nan; + + % Recode MinATE to help spread out the numbers + % (so that the colors are more distinct) + U = unique(MinATE); + for k = 1:1:length(U) + I = find(MinATE == U(k)); + MinATE(I) = k; + end + + h = pcolor(Settings.Beta0Grid, Settings.Beta1Grid, MinATE); + set(h, 'EdgeColor', 'none'); + colormap(jet) + + Title = 'Identified Set (Partitioned)'; + Filename = 'IDSetBetaPartitioned'; + end + + xlim([min(Settings.Beta0Grid) max(Settings.Beta0Grid)]); + ylim([min(Settings.Beta1Grid) max(Settings.Beta1Grid)]); + grid on; + grid minor; + hold on; + + if ~(Type == 1) + set(gca, 'layer', 'top'); + end + + xlabel('$\beta_{0}$', 'Interpreter', 'Latex',... + 'FontSize', Settings.FontsizeLabel); + ylabel('$\beta_{1}$', 'Rotation', 0, 'Interpreter', 'Latex',... + 'FontSize', Settings.FontsizeLabel); + plot(Settings.Beta0, Settings.Beta1, 'k+',... + 'MarkerSize', Settings.MarkersizeBig); + + set(gca, 'XTick',... + fix(min(Settings.Beta0Grid)):1:fix(max(Settings.Beta0Grid))); + set(gca, 'YTick',... + fix(min(Settings.Beta1Grid)):1:fix(max(Settings.Beta1Grid))); + set(gca, 'FontSize', Settings.FontsizeAxis) + + daspect([1 Settings.FigWidth/Settings.FigHeight 1]); + + set(gcf, 'InvertHardcopy', 'on'); + set(gcf, 'PaperUnits', 'inches'); + papersize = get(gcf, 'PaperSize'); + left = (papersize(1) - Settings.FigWidth)/2; + bottom = (papersize(2) - Settings.FigHeight)/2; + myfiguresize = [left, bottom, Settings.FigWidth, Settings.FigHeight]; + set(gcf, 'PaperPosition', myfiguresize); + % This gets rid of some extra whitespace + set(gca,'LooseInset',[0 0 0 0]); + + fname = sprintf([Filename '-%s'], Settings.DirName); + savefig(fname); + print(fname, '-dpng', '-r300'); + print(fname,'-depsc','-tiff'); + + title(Title, 'FontSize', Settings.FontsizeTitle); + fname = sprintf([Filename '-Standalone-%s'], Settings.DirName); + savefig(fname); + print(fname, '-dpng'); + print(fname,'-depsc','-tiff'); + +end diff --git a/src/PrintStructure.m b/src/PrintStructure.m new file mode 100644 index 0000000..0708273 --- /dev/null +++ b/src/PrintStructure.m @@ -0,0 +1,35 @@ +%******************************************************************************* +% PrintStructure.m +% +% Recurse through a structure and print its contents to fid +%******************************************************************************* +function PrintStructure(s, fid) + + fields = repmat(fieldnames(s), numel(s), 1); + values = struct2cell(s); + + for i = 1:1:length(fields) + if isstruct(s.(fields{i})) + fprintf(fid, '-----------------\n'); + fprintf(fid, '%s\n', fields{i}); + fprintf(fid, '-----------------\n'); + PrintStructure(s.(fields{i}), fid); % Recursion -- so clever :) + else + if isnumeric(values{i}) + values{i} = num2str(values{i}(:)', '%4.3f'); + end + if isa(values{i}, 'function_handle') + values{i} = func2str(values{i}); + end + stringout = [fields{i} ': ']; + if iscell(values{i}) + stringout = [stringout strjoin(values{i})]; + else + stringout = [stringout values{i}]; + end + + fprintf(fid, '%s\n', stringout); + end + end + fprintf(fid, '-----------------\n'); +end diff --git a/src/Setup.m b/src/Setup.m new file mode 100644 index 0000000..ef9d715 --- /dev/null +++ b/src/Setup.m @@ -0,0 +1,106 @@ +%******************************************************************************* +% Setup.m +% +% Setup actions at the start of Run +%******************************************************************************* +function [CleanUpPath, DirName] = Setup(ScriptTag, DirName) + + OriginalPath = pwd; + CleanUpPath = onCleanup(@()cd(OriginalPath)); + + addpath('../cfg') + evalc('Config'); % Load user-defined paths + + if ~exist(SAVEDIRECTORY, 'dir') + error('SAVEDIRECTORY %s does not exist.', SAVEDIRECTORY); + end + + % Create a directory path name + if ~exist('DirName', 'var') + DirName = []; + end + + if isempty(DirName) + [DirPath DirName] = MakeTargetDirName(SAVEDIRECTORY, ScriptTag, 1); + else + DirPath = [SAVEDIRECTORY DirName]; + end + + % Create a new directory, copy everything there and cd to it + CreateSaveDirectory(DirPath); + + % Create a results subdirectory and move there + if ~exist('./results', 'dir') + mkdir('results'); + end + cd('results'); + + % Add paths + addpath(genpath(DirPath)); + + % Setup AMPL + if ~exist(AMPLAPISETUPPATH, 'file') + error('Path %s for AMPL API SetUp.m file does not exist.',... + AMPLAPISETUPPATH); + else + evalc(['run ' AMPLAPISETUPPATH]); + end +end + +%******************************************************************************* +% MakeTargetDirName +% +% Make a target directory, with option to ask user for a name +%******************************************************************************* +function [DirPath, DirName] ... + = MakeTargetDirName(BASEDIR, STUB, Interactive, tag) +%******************************************************************************* + if ~exist('Interactive') + Interactive = 1; + end + + if Interactive + tag = input('Do you want to add a tag to the directory?\n', 's'); + else + if ~exist('tag') + tag = []; + end + end + if ~isempty(tag) + if ~strcmp(tag(1), '-') + tag = ['-' tag]; + end + end + if ~strcmp(BASEDIR(length(BASEDIR)), '/') + BASEDIR = [BASEDIR '/']; + end + DirName = [STUB tag datestr(now, '-mmddyy-HHMMSS')]; + DirPath = [BASEDIR DirName]; +end + +%******************************************************************************* +% CreateSaveDirectory +% +% Create a directory with a current copy of the code (so that further changes +% don't affect this simulation) and where we can store results. +% Create a subdirectory in that directory to store results. +% And change directories to that subdirectory. +%******************************************************************************* +function [] = CreateSaveDirectory(DirPath) + if isempty(DirPath) + error('DirPath is empty!'); + end + if ~exist(DirPath, 'dir') + disp(sprintf('Saving in new directory %s.', DirPath)); + success = mkdir(DirPath); + if ~success + error('Something is wrong; error creating directory.') + end + copyfile('../*', DirPath); + + elseif (exist(DirPath, 'dir') == 7) + disp(sprintf('Saving in existing directory %s.', DirPath)); + disp(sprintf('\t Directory already exists; not copying code.')); + end + cd(DirPath); +end diff --git a/src/UpdateStruct.m b/src/UpdateStruct.m new file mode 100644 index 0000000..aed7bd2 --- /dev/null +++ b/src/UpdateStruct.m @@ -0,0 +1,45 @@ +%******************************************************************************* +% UpdateStruct.m +% +% Struct is a default structure to be updated (usually default settings) +% StructIn is the new structure (usually user-inputted settings) +% +% Different behavior depending on the flags EnforceSubset and KeepUnknown: +% +% E=0, K=0: +% Don't throw an error if StructIn has fields that Struct does +% not, but don't keep any of these unrecognized fields. +% +% E=0, K=1: +% Don't throw an error if StructIn has fields that Struct does +% not, and add any of the unknown fields to Struct +% +% E=1 +% Throw an error if StructIn has fields that Struct does not +%******************************************************************************* +function Struct = ... + UpdateStruct(Struct, StructIn, EnforceSubset, KeepUnknown) +%******************************************************************************* + if ~isstruct(StructIn) + error('Expected structure for input. Quitting.'); + end + + StructNamesDefault = fieldnames(Struct); + StructNamesIn = fieldnames(StructIn); + + for i = 1:1:length(StructNamesIn) + CurOption = StructNamesIn{i}; + + % Overwrite default with input + if any(strcmp(CurOption, StructNamesDefault)) + Struct.(CurOption) = StructIn.(CurOption); + else + if EnforceSubset + error([ '%s is not a recognized field '... + 'of the structure to be updated.'], CurOption); + elseif KeepUnknown + Struct.(CurOption) = StructIn.(CurOption); + end + end + end +end