Merge pull request #47 from bpuchala/0.2.X-backlog

0.2.X backlog

Author: bpuchala

INSTALL.md

@@ -5,7 +5,7 @@
 
 **C++11**
 
-CASM must be compiled with a compiler that supports the C++11 standard.
+CASM must be compiled with a compiler that supports the C++11 standard. Testing is done with gcc-4.8.5 and clang-800.0.42.1.
 
 **If using Mac OS X - Xcode command-line tools**
 

python/casm/casm/learn/evolve.py

@@ -8,6 +8,7 @@
 import deap.tools
 import deap.algorithms
 from operator import attrgetter
+from deap.tools import HallOfFame
 
 
 def initNRandomOn(container, n_features, n_features_init):
@@ -204,6 +205,7 @@ def initialize_population(n_population, toolbox, filename=None, verbose=True):
   """
   if filename is not None and os.path.exists(filename):
     load initial population
+    # may be List[Individual] or HallOfFame, which is converted to List[Individual]
   else:
     create random initial population of size n_population via toolbox.population
   """
@@ -214,6 +216,9 @@ def initialize_population(n_population, toolbox, filename=None, verbose=True):
       print "Loading initial population:", filename
     with open(filename, 'rb') as f:
       pop = pickle.load(f)
+    if isinstance(pop, HallOfFame):
+      # convert to List
+      pop = [indiv for indiv in pop]
   else:
     if verbose:
       print "Constructing initial population"

python/casm/casm/learn/fit.py

@@ -474,12 +474,30 @@ def print_input_help():
   # The "problem_specs"/"weight" options specify the method to use for weighting 
   # training data. 
   #
-  #   If weights are included, then the linear model is changed from
-  #     X*b = y  ->  L*X*b = L*y, 
+  #   Ordinary least squares minimizes
+  #     (y-X*b).transpose() * (y-X*b)
+  #  
+  #   where 'X' is the correlation matrix of shape (Nvalue, Nbfunc), and 'y' 
+  #   is a vector of Nvalue calculated properties, and 'b' are the fitting 
+  #   coefficients (ECI).
   #
-  #   where 'X' is the correlation matrix of shape (Nvalue, Nbfunc),
-  #   and 'property' is a vector of Nvalue calculated properties, and 
-  #   W = L*L.transpose() is the weight matrix.
+  #   Weighted least squares minimizes
+  #     (y-X*b).transpose() * W * (y-X*b)
+  #  
+  #   Using the SVD, and given that W is Hermitian:
+  #     U * S * U.transpose() == W
+  #  
+  #   Define L such that:
+  #     L.transpose() = U * sqrt(S)
+  #  
+  #   Then we can write the weighted least squares problem using:
+  #     (y-X*b).transpose() * L.transpose() * L * (y-X*b)
+  #  
+  #   Or:
+  #     (L*y-L*X*b).transpose() * (L*y-L*X*b)
+  #    
+  #   So, if weights are included, then the linear model is changed from
+  #     X*b = y  ->  L*X*b = L*y
   #
   #   By default, W = np.matlib.eye(Nvalue) (unweighted).
   #
@@ -775,8 +793,8 @@ def print_input_help():
   #   Options for "evolve_params_kwargs":
   #
   #     "n_population": int, optional, default=100
-  #        Population size. This many random initial starting individuals are 
-  #        created.
+  #        Initial population size. This many random initial starting individuals 
+  #         are created if no "pop_begin_filename" file exists.
   #     
   #     "n_halloffame": int, optional, default=25
   #        Maxsize of the hall of fame which holds the best individuals 
@@ -797,19 +815,32 @@ def print_input_help():
   #        with.
   #
   #     "pop_begin_filename": string, optional, default="population_begin.pkl"
-  #        Filename where the initial population is read from, if it exists.
-  #
+  #        Filename suffix where the initial population is read from, if it 
+  #        exists. For example, if "filename_prefix" is "Ef_kfold10" and 
+  #        "pop_begin_filename" is "population_begin.pkl", then the initial
+  #        population is read from the file "Ef_kfold10_population_begin.pkl".
+  #
+  #        The population file may contain either a  list of individual,
+  #        as written to the "population_end.pkl" file, or a HallOfFame
+  #        instance, as written to either an "evolve_halloffame.pkl" file or
+  #        overall casm-learn "halloffame.pkl" file.
+  #        
   #     "pop_end_filename": string, optional, default="population_end.pkl"
-  #        Filename where the final population is saved.
+  #        Filename where the final population is saved. For example, if 
+  #        "filename_prefix" is "Ef_kfold10" and "pop_end_filename" is 
+  #        "population_end.pkl", then the final population is saved to the 
+  #        file "Ef_kfold10_population_end.pkl".
   #      
   #     "halloffame_filename": string, optional, default="evolve_halloffame.pkl"
   #        Filename where a hall of fame is saved holding the best individuals 
-  #        encountered in any generation.
+  #        encountered in any generation. For example, if "filename_prefix" is 
+  #        "Ef_kfold10" and "halloffame_filename" is "evolve_halloffame.pkl", 
+  #        then it is saved to the file "Ef_kfold10_evolve_halloffame.pkl".
   #      
   #     "filename_prefix": string, optional
   #        Prefix for filenames, default uses input file filename excluding 
   #        extension. For example, if input file is named "Ef_kfold10.json", then
-  #        "Ef_kfold10_population_begin.pkl", "specs2_population_end.pkl", and 
+  #        "Ef_kfold10_population_begin.pkl", "Ef_kfold10_population_end.pkl", and 
   #        "Ef_kfold10_evolve_halloffame.pkl" are used.
   
     "feature_selection" : {
@@ -1414,11 +1445,11 @@ def read_sample_weight(input, tdata, verbose=True):
   # use method to get weights
   if specs["weight"]["method"] == "wCustom":
     if verbose:
-      print "Reading custom weights"
+      print "# Reading custom weights"
     sample_weight = tdata.data["weight"].values
   elif specs["weight"]["method"] == "wCustom2d":
     if verbose:
-      print "Reading custom2d weights"
+      print "# Reading custom2d weights"
     cols = ["weight(" + str(i) + ")" for i in xrange(tdata.n_samples)]
     sample_weight = tdata.data.loc[:,cols].values
   elif specs["weight"]["method"] == "wHullDist":
@@ -1543,7 +1574,9 @@ def check_input(name):
       pickle.dump(fdata, open(fit_data_filename, 'wb'))
 
     if verbose:
-      print "# Writing problem specs to:", fit_data_filename, "\n"
+      print "# Writing problem specs to:", fit_data_filename
+      print "# To inspect or customize the problem specs further, use the '--checkspecs' method\n"
+    
 
   # during runtime only, if LinearRegression and LeaveOneOut, update fdata.cv and fdata.scoring
   # to use optimized LOOCV score method

python/casm/casm/learn/tools.py

@@ -76,12 +76,32 @@ def set_sample_weight(sample_weight, y=None, X=None):
   """ 
   Calculate weighted data and weighted target values.
   
-  Uses sample weights to calculate
+  Ordinary least squares minimizes
+    (y-X*b).transpose() * (y-X*b)
   
-    L*X * b = L*y 
+  where 'X' is the correlation matrix of shape (Nvalue, Nbfunc), and 'y' 
+  is a vector of Nvalue calculated properties, and 'b' are the fitting 
+  coefficients (ECI).
+
+  Weighted least squares minimizes
+    (y-X*b).transpose() * W * (y-X*b)
+  
+  Using the SVD, and given that W is Hermitian:
+    U * S * U.transpose() == W
+  
+  Define L such that:
+    L.transpose() = U * sqrt(S)
   
-  a weighted linear model where the weights are given by W = L * L.transpose().
+  Then we can write the weighted least squares problem using:
+    (y-X*b).transpose() * L.transpose() * L * (y-X*b)
   
+  Or:
+    (L*y-L*X*b).transpose() * (L*y-L*X*b)
+    
+  So, if weights are included, then the linear model is changed from
+    X*b = y  ->  L*X*b = L*y
+    
+    
   Arguments
   ---------
     
@@ -142,7 +162,8 @@ def set_sample_weight(sample_weight, y=None, X=None):
     raise Exception("Error in set_sample_weight: sample_weight dimension > 2")
 
   # weighted data
-  L = np.linalg.cholesky(W)
+  U, S, V = np.linalg.svd(W)
+  L = U.dot(np.diag(np.sqrt(S))).transpose()
 
   if X is not None:
     weighted_X = np.dot(L, X)

src/casm/app/enum.cc

@@ -144,9 +144,9 @@ namespace CASM {
                   "    volume from MINV to MAXV (units: number of primitive \n"
                   "    cells).                                              \n\n"
 
-                  "  casm enum --configs --scellname NAME                   \n"
-                  "  - To enumerate configurations for a particular         \n"
-                  "    supercell.                                           \n\n"
+                  "  casm enum --configs --scelnames NAME1 NAME2            \n"
+                  "  - To enumerate configurations for one or more particular\n"
+                  "    supercells.                                          \n\n"
 
                   "  casm enum --configs [...] --filter '... casm query commands...' \n"
                   "  - To perform restricted enumeration of configurations  \n"
@@ -199,9 +199,9 @@ namespace CASM {
         std::cerr << "Error in 'casm enum'. If --supercells is given, --max must be given." << std::endl;
         return ERR_INVALID_ARG;
       }
-      if(vm.count("configs") && (vm.count("max") + vm.count("all") != 1)) {
+      if(vm.count("configs") && (vm.count("max") + vm.count("scelnames") + vm.count("all") != 1)) {
         std::cerr << "\n" << enum_opt.desc() << "\n" << std::endl;
-        std::cerr << "Error in 'casm enum'. If --configs is given, exactly one of either --max or --all must be given." << std::endl;
+        std::cerr << "Error in 'casm enum'. If --configs is given, exactly one of either --max, --scelnames, or --all must be given." << std::endl;
         return ERR_INVALID_ARG;
       }
     }
@@ -309,7 +309,7 @@ namespace CASM {
             }
           }
         }
-        if(vm.count("scellname")) {
+        if(vm.count("scelnames")) {
           Index j;
           std::cout << "Enumerate configurations for named supercells" << std::endl << std::endl;
           for(int i = 0; i < scellname_list.size(); i++) {

src/casm/app/format.cc

@@ -732,7 +732,7 @@ LCHARG = .FALSE.\n";
                 << "'chebychev'. Otherwise, specifies a JSON object containing a composition vector or\n"
                 << "a JSON array containing multiple composition vectors. A single composition vector\n"
                 << "is formatted as, e.g.\n"
-                << "   \"composition\" : [\"Au\" : 0.25, \"Cu\" : 0.75] \n"
+                << "   \"composition\" : {\"Au\" : 0.25, \"Cu\" : 0.75} \n"
                 << "The site basis functions will then be constructed as to be optimized for that composition.\n\n"
 
                 << "To specify different compositions on multiple sublattices, an array can be used. \n"
@@ -741,11 +741,11 @@ LCHARG = .FALSE.\n";
 
                 << "   \"site_basis_functions\" : [\n"
                 << "                                {\n"
-                << "                                  \"composition\" : [\"Ga\" : 0.3, \"In\" : 0.7],\n"
+                << "                                  \"composition\" : {\"Ga\" : 0.3, \"In\" : 0.7},\n"
                 << "                                  \"sublat_indices\" : [0]\n"
                 << "                                },\n"
                 << "                                {\n"
-                << "                                  \"composition\" : [\"Ga\" : 1.0, \"In\" : 0.0],\n"
+                << "                                  \"composition\" : {\"Ga\" : 1.0, \"In\" : 0.0},\n"
                 << "                                  \"sublat_indices\" : [1,2]\n"
                 << "                                }\n"
                 << "                             ]\n\n"

src/casm/monte_carlo/grand_canonical/GrandCanonical.cc

@@ -423,8 +423,7 @@ namespace CASM {
     auto corr = correlations(config, _clexulator());
     double formation_energy = _eci() * corr.data();
     auto comp_x = primclex().composition_axes().param_composition(CASM::comp_n(config));
-    double Ep = formation_energy - comp_x.dot(m_condition.param_chem_pot());
-    return Ep / supercell().volume();
+    return formation_energy - comp_x.dot(m_condition.param_chem_pot());
   }
 
   /// \brief Calculate delta correlations for an event

src/casm/monte_carlo/grand_canonical/GrandCanonicalIO.cc

@@ -286,7 +286,8 @@ namespace CASM {
 
     DataFormatter<ConstMonteCarloPtr> formatter;
 
-    formatter.push_back(ConstantValueFormatter<std::string, ConstMonteCarloPtr>("configname", configname));
+    bool print_json = true;
+    formatter.push_back(ConstantValueFormatter<std::string, ConstMonteCarloPtr>("configname", configname, print_json));
     formatter.push_back(MonteCarloTFormatter<GrandCanonical>());
     formatter.push_back(GrandCanonicalLTEFormatter(phi_LTE1));
     std::set<std::string> exclude;