Add files via upload

Code/14ops.txt

@@ -0,0 +1 @@
++*-D><~IRPSCLE

Code/19ops.txt

@@ -0,0 +1 @@
++*-D><~IRPLESCANT01

Code/7ops.txt

@@ -0,0 +1 @@
++*D>~R0

Code/README

@@ -0,0 +1,16 @@
+Run compile.sh to compile the fortran files.
+
+ai_feynman_example.py contains an example of running a code on some examples (found in the example_data directory). The function has the following parameters
+
+pathdir - path to the directory containing the data file
+filename - the name of the file containing the data
+BF_try_time - time limit for each brute force call (set by default to 60 seconds)
+BF_ops_file_type - file containing the symbols to be used in the brute force code (set by default to "14ops.txt")
+polyfit_deg - maximum degree of the polynomial tried by the polynomial fit routine (set be default to 4)
+NN_epochs - number of epochs for the training (set by default to 4000)
+
+The solution file will be saved in the directory called results under the name solution_{filename}.
+
+ai_feynman_terminal_example.py allows calling the aiFeynman function from the command line. 
+e.g. python ai_feynman_terminal_example.py --pathdir=../example_data/ --filename=example1.txt
+python ai_feynman_terminal_example.py --help displays all the available parameters that can be passed to the function.

Code/RPN_to_eq.py

@@ -0,0 +1,69 @@
+# Turns an RPN expression to normal mathematical notation
+
+import numpy as np
+
+def RPN_to_eq(expr):
+
+    variables = ["0","1","a","b","c","d","e","f","g","h","i","j","k","l","m","n","P"]
+    operations_1 = [">","<","~","\\","L","E","S","C","A","N","T","R"]
+    operations_2 = ["+","*","-","/"]
+
+    stack = np.array([])
+
+    for i in (expr):
+        if i in variables:
+            if i == "P":
+                stack = np.append(stack,"pi")
+            elif i == "0":
+                stack = np.append(stack,"0")
+            elif i == "1":
+                stack = np.append(stack,"1")  
+            else:
+                stack = np.append(stack,"x" + str(ord(i)-97))
+        elif i in operations_2:
+            a1 = stack[-1]
+            a2 = stack[-2]
+            stack = np.delete(stack,-1)
+            stack = np.delete(stack,-1)
+            a = "("+a2+i+a1+")"
+            stack = np.append(stack,a)
+        elif i in operations_1:
+            a = stack[-1]
+            stack = np.delete(stack,-1)
+            if i==">":
+                a="("+a+"+1)"
+                stack = np.append(stack,a)
+            if i=="<":
+                a="("+a+"-1)"
+                stack = np.append(stack,a)
+            if i=="~":
+                a="(-"+a+")"
+                stack = np.append(stack,a)
+            if i=="\\":
+                a="("+a+")**(-1)"
+                stack = np.append(stack,a)
+            if i=="L":
+                a="log("+a+")"
+                stack = np.append(stack,a)
+            if i=="E":
+                a="exp("+a+")"
+                stack = np.append(stack,a)
+            if i=="S":
+                a="sin("+a+")"
+                stack = np.append(stack,a)
+            if i=="C":
+                a="cos("+a+")"
+                stack = np.append(stack,a)
+            if i=="A":
+                a="abs("+a+")"
+                stack = np.append(stack,a)
+            if i=="N":
+                a="asin("+a+")"
+                stack = np.append(stack,a)
+            if i=="T":
+                a="atan("+a+")"
+                stack = np.append(stack,a)
+            if i=="R":
+                a="sqrt("+a+")"
+                stack = np.append(stack,a)
+    return(stack[0])

Code/RPN_to_pytorch.py

@@ -0,0 +1,130 @@
+# Turns a mathematical expression (already RPN turned) to pytorch expression, trains the parameters, and returns the new error, complexity and the new symbolic expression
+
+import numpy as np 
+import matplotlib.pyplot as plt
+import pandas as pd
+import torch 
+import torch.nn as nn 
+import torch.nn.functional as F 
+import torch.optim as optim 
+import torch.utils.data as utils
+from torch.autograd import Variable
+from sklearn.metrics import roc_curve, auc
+from sklearn.preprocessing import label_binarize
+from sklearn.manifold import TSNE
+import seaborn as sns
+import warnings
+warnings.filterwarnings("ignore")
+import sympy
+
+from sympy import *
+from sympy.abc import x,y
+from sympy.parsing.sympy_parser import parse_expr
+from sympy import Symbol, lambdify, N
+
+from S_get_number_DL import get_number_DL
+
+# parameters: path to data, RPN expression (obtained from bf)
+def RPN_to_pytorch(data_file, math_expr, lr = 1e-2, N_epochs = 500):
+    param_dict = {}
+    unsnapped_param_dict = {'p':1}
+
+    def unsnap_recur(expr, param_dict, unsnapped_param_dict):
+        """Recursively transform each numerical value into a learnable parameter."""
+        import sympy
+        from sympy import Symbol
+        if isinstance(expr, sympy.numbers.Float) or isinstance(expr, sympy.numbers.Integer) or isinstance(expr, sympy.numbers.Rational) or isinstance(expr, sympy.numbers.Pi):
+            used_param_names = list(param_dict.keys()) + list(unsnapped_param_dict)
+            unsnapped_param_name = get_next_available_key(used_param_names, "p", is_underscore=False)
+            unsnapped_param_dict[unsnapped_param_name] = float(expr)
+            unsnapped_expr = Symbol(unsnapped_param_name)
+            return unsnapped_expr
+        elif isinstance(expr, sympy.symbol.Symbol):
+            return expr
+        else:
+            unsnapped_sub_expr_list = []
+            for sub_expr in expr.args:
+                unsnapped_sub_expr = unsnap_recur(sub_expr, param_dict, unsnapped_param_dict)
+                unsnapped_sub_expr_list.append(unsnapped_sub_expr)
+            return expr.func(*unsnapped_sub_expr_list)
+
+
+    def get_next_available_key(iterable, key, midfix="", suffix="", is_underscore=True):
+        """Get the next available key that does not collide with the keys in the dictionary."""
+        if key + suffix not in iterable:
+            return key + suffix
+        else:
+            i = 0
+            underscore = "_" if is_underscore else ""
+            while "{}{}{}{}{}".format(key, underscore, midfix, i, suffix) in iterable:
+                i += 1
+            new_key = "{}{}{}{}{}".format(key, underscore, midfix, i, suffix)
+            return new_key
+
+    # Load the actual data
+
+    data = np.loadtxt(data_file)
+
+    # Turn BF expression to pytorch expression
+    eq = parse_expr(math_expr)
+    eq = unsnap_recur(eq,param_dict,unsnapped_param_dict)
+
+    N_vars = len(data[0])-1
+    N_params = len(unsnapped_param_dict)
+
+    possible_vars = ["x%s" %i for i in np.arange(0,30,1)]
+    variables = []
+    params = []
+    for i in range(N_vars):
+        variables = variables + [possible_vars[i]]
+    for i in range(N_params-1):
+        params = params + ["p%s" %i]
+
+    symbols = params + variables
+
+    f = lambdify(symbols, N(eq), torch)
+
+    # Set the trainable parameters in the expression
+
+    trainable_parameters = []
+    for i in unsnapped_param_dict:
+        if i!="p":
+            vars()[i] = torch.tensor(unsnapped_param_dict[i])
+            vars()[i].requires_grad=True
+            trainable_parameters = trainable_parameters + [vars()[i]]
+
+    # Prepare the loaded data
+
+    real_variables = []
+    for i in range(len(data[0])-1):
+        real_variables = real_variables + [torch.from_numpy(data[:,i]).float()]
+
+    input = trainable_parameters + real_variables
+
+    y = torch.from_numpy(data[:,-1]).float()
+
+    for i in range(N_epochs):
+        # this order is fixed i.e. first parameters
+        yy = f(*input)
+        loss = torch.mean((yy-y)**2)
+        loss.backward()
+        with torch.no_grad():
+            for j in range(N_params-1):
+                trainable_parameters[j] -= lr * trainable_parameters[j].grad
+                trainable_parameters[j].grad.zero_()
+
+    # get the updated symbolic regression
+    ii = -1
+    complexity = 0
+    for parm in unsnapped_param_dict:
+        if ii == -1:
+            ii = ii + 1
+        else:
+            eq = eq.subs(parm, trainable_parameters[ii])
+            complexity = complexity + get_number_DL(trainable_parameters[ii].detach().numpy())
+            ii = ii+1
+            
+        
+    error = torch.mean((f(*input)-y)**2).data.numpy()*1
+    return error, complexity, eq
+

Code/S_NN_eval.py

@@ -0,0 +1,123 @@
+from __future__ import print_function
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+from torchvision import datasets, transforms
+import pandas as pd
+import numpy as np
+import torch
+from torch.utils import data
+import pickle
+from torch.optim.lr_scheduler import CosineAnnealingLR
+from matplotlib import pyplot as plt
+import time
+
+is_cuda = torch.cuda.is_available()
+
+bs = 2048
+
+class MultDataset(data.Dataset):
+    def __init__(self, factors, product):
+        'Initialization'
+        self.factors = factors
+        self.product = product
+
+    def __len__(self):
+        'Denotes the total number of samples'
+        return len(self.product)
+
+    def __getitem__(self, index):
+        # Load data and get label
+        x = self.factors[index]
+        y = self.product[index]
+
+        return x, y
+
+def rmse_loss(pred, targ):
+    denom = targ**2
+    denom = torch.sqrt(denom.sum()/len(denom))
+
+    return torch.sqrt(F.mse_loss(pred, targ))/denom
+
+
+def NN_eval(pathdir,filename):
+    try:
+        n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+filename, usecols=(0,))
+
+        if n_variables==0:
+            return 0
+        elif n_variables==1:
+            variables = np.reshape(variables,(len(variables),1))
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+
+        f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables[0:int(5*len(variables)/6)])
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+        product = torch.from_numpy(f_dependent[0:int(5*len(f_dependent)/6)])
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        factors_val = torch.from_numpy(variables[int(5*len(variables)/6):int(len(variables))])
+        if is_cuda:
+            factors_val = factors_val.cuda()
+        else:
+            factors_val = factors_val
+        factors_val = factors_val.float()
+        product_val = torch.from_numpy(f_dependent[int(5*len(variables)/6):int(len(variables))])      
+        if is_cuda:
+            product_val = product_val.cuda()
+        else:
+            product_val = product_val
+        product_val = product_val.float()
+
+        class SimpleNet(nn.Module):
+            def __init__(self, ni):
+                super().__init__()
+                self.linear1 = nn.Linear(ni, 128)
+                self.bn1 = nn.BatchNorm1d(128)
+                self.linear2 = nn.Linear(128, 128)
+                self.bn2 = nn.BatchNorm1d(128)
+                self.linear3 = nn.Linear(128, 64)
+                self.bn3 = nn.BatchNorm1d(64)
+                self.linear4 = nn.Linear(64,64)
+                self.bn4 = nn.BatchNorm1d(64)
+                self.linear5 = nn.Linear(64,1)
+            
+            def forward(self, x):
+                x = F.tanh(self.bn1(self.linear1(x)))
+                x = F.tanh(self.bn2(self.linear2(x)))
+                x = F.tanh(self.bn3(self.linear3(x)))
+                x = F.tanh(self.bn4(self.linear4(x)))
+                x = self.linear5(x)
+                return x  
+
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+                    
+        model.load_state_dict(torch.load("results/NN_trained_models/models/"+filename+".h5"))
+        model.eval()
+        return(rmse_loss(model(factors_val),product_val))
+
+    except Exception as e:
+        print(e)
+        return (100)
+
+
+
+

Code/S_NN_train.py

@@ -0,0 +1,151 @@
+from __future__ import print_function
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+from torchvision import datasets, transforms
+import pandas as pd
+import numpy as np
+import torch
+from torch.utils import data
+import pickle
+from matplotlib import pyplot as plt
+import torch.utils.data as utils
+import time
+import os
+
+bs = 2048
+wd = 1e-2
+
+is_cuda = torch.cuda.is_available()
+
+class MultDataset(data.Dataset):
+    def __init__(self, factors, product):
+        'Initialization'
+        self.factors = factors
+        self.product = product
+
+    def __len__(self):
+        'Denotes the total number of samples'
+        return len(self.product)
+
+    def __getitem__(self, index):
+        # Load data and get label
+        x = self.factors[index]
+        y = self.product[index]
+
+        return x, y
+
+def rmse_loss(pred, targ):
+    denom = targ**2
+    denom = torch.sqrt(denom.sum()/len(denom))
+    return torch.sqrt(F.mse_loss(pred, targ))/denom
+
+def NN_train(pathdir, filename, epochs=-1):
+    try:
+        os.mkdir("results/NN_trained_models/")
+    except:
+        pass
+
+    try:
+        os.mkdir("results/NN_trained_models/models/")
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+
+        epochs = epochs//4
+
+        if n_variables==0:
+            print("Solved! ", variables[0])
+            return 0
+        elif n_variables==1:
+            variables = np.reshape(variables,(len(variables),1))
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables)
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        class SimpleNet(nn.Module):
+            def __init__(self, ni):
+                super().__init__()
+                self.linear1 = nn.Linear(ni, 128)
+                self.bn1 = nn.BatchNorm1d(128)
+                self.linear2 = nn.Linear(128, 128)
+                self.bn2 = nn.BatchNorm1d(128)
+                self.linear3 = nn.Linear(128, 64)
+                self.bn3 = nn.BatchNorm1d(64)
+                self.linear4 = nn.Linear(64,64)
+                self.bn4 = nn.BatchNorm1d(64)
+                self.linear5 = nn.Linear(64,1)
+            
+            def forward(self, x):
+                x = F.tanh(self.bn1(self.linear1(x)))
+                x = F.tanh(self.bn2(self.linear2(x)))
+                x = F.tanh(self.bn3(self.linear3(x)))
+                x = F.tanh(self.bn4(self.linear4(x)))
+                x = self.linear5(x)
+                return x
+
+        my_dataset = utils.TensorDataset(factors,product) # create your datset
+        my_dataloader = utils.DataLoader(my_dataset, batch_size=bs, shuffle=True) # create your dataloader
+
+        if is_cuda:
+            model_feynman = SimpleNet(n_variables).cuda()
+        else:
+            model_feynman = SimpleNet(n_variables)
+
+        max_loss = 10000
+
+        lrs = 1e-2
+        for i_i in range(4):
+            optimizer_feynman = optim.Adam(model_feynman.parameters(), lr = lrs)
+            for epoch in range(epochs):
+                model_feynman.train()
+                for i, data in enumerate(my_dataloader):
+                    optimizer_feynman.zero_grad()
+                
+                    if is_cuda:
+                        fct = data[0].float().cuda()
+                        prd = data[1].float().cuda()
+                    else:
+                        fct = data[0].float()
+                        prd = data[1].float()
+                    
+                    loss = rmse_loss(model_feynman(fct),prd)
+                    if loss < max_loss:
+                        torch.save(model_feynman.state_dict(), "results/NN_trained_models/models/" + filename + ".h5")
+                        max_loss = loss
+                    loss = rmse_loss(model_feynman(fct),prd)
+                    loss.backward()
+                    optimizer_feynman.step()
+            print(loss)
+            lrs = lrs/10
+
+        return 1
+
+    except NameError:
+        print("Error in file: %s" %filename)
+        raise
+
+
+
+

Code/S_add_snap_expr_on_pareto.py

@@ -0,0 +1,192 @@
+# Adds on the pareto all the snapped versions of a given expression (all paramters are snapped in the end)
+
+import numpy as np 
+import matplotlib.pyplot as plt
+import pandas as pd
+import torch 
+import torch.nn as nn 
+import torch.nn.functional as F 
+import torch.optim as optim 
+import torch.utils.data as utils
+from torch.autograd import Variable
+from sklearn.metrics import roc_curve, auc
+from sklearn.preprocessing import label_binarize
+from sklearn.manifold import TSNE
+import seaborn as sns
+import copy
+import warnings
+warnings.filterwarnings("ignore")
+import sympy
+from S_snap import integerSnap
+from S_snap import zeroSnap
+from S_snap import rationalSnap
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+from get_pareto import Point, ParetoSet
+
+from sympy import preorder_traversal, count_ops
+from sympy.abc import x,y
+from sympy.parsing.sympy_parser import parse_expr
+from sympy import Symbol, lambdify, N, simplify, powsimp, Rational, symbols
+
+from S_get_number_DL_snapped import get_number_DL_snapped
+
+# parameters: path to data, math (not RPN) expression 
+def add_snap_expr_on_pareto(pathdir, filename, math_expr, PA, DR_file=""):
+    def unsnap_recur(expr, param_dict, unsnapped_param_dict):
+        """Recursively transform each numerical value into a learnable parameter."""
+        import sympy
+        from sympy import Symbol
+        if isinstance(expr, sympy.numbers.Float) or isinstance(expr, sympy.numbers.Integer) or isinstance(expr, sympy.numbers.Rational) or isinstance(expr, sympy.numbers.Pi):
+            used_param_names = list(param_dict.keys()) + list(unsnapped_param_dict)
+            unsnapped_param_name = get_next_available_key(used_param_names, "p", is_underscore=False)
+            unsnapped_param_dict[unsnapped_param_name] = float(expr)
+            unsnapped_expr = Symbol(unsnapped_param_name)
+            return unsnapped_expr
+        elif isinstance(expr, sympy.symbol.Symbol):
+            return expr
+        else:
+            unsnapped_sub_expr_list = []
+            for sub_expr in expr.args:
+                unsnapped_sub_expr = unsnap_recur(sub_expr, param_dict, unsnapped_param_dict)
+                unsnapped_sub_expr_list.append(unsnapped_sub_expr)
+            return expr.func(*unsnapped_sub_expr_list)
+
+
+    def get_next_available_key(iterable, key, midfix="", suffix="", is_underscore=True):
+        """Get the next available key that does not collide with the keys in the dictionary."""
+        if key + suffix not in iterable:
+            return key + suffix
+        else:
+            i = 0
+            underscore = "_" if is_underscore else ""
+            while "{}{}{}{}{}".format(key, underscore, midfix, i, suffix) in iterable:
+                i += 1
+            new_key = "{}{}{}{}{}".format(key, underscore, midfix, i, suffix)
+            return new_key
+
+    eq = parse_expr(str(math_expr))
+    expr = eq
+    
+    # Get the numbers appearing in the expression
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    eq_numbers = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+       
+    # Do integer snap one parameter at a time
+    integer_snapped_expr = []
+    for w in range(len(eq_numbers)):
+        param_dict = {}
+        unsnapped_param_dict = {'p':1}
+        eq = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+        new_numbers = integerSnap(eq_numbers,w+1)
+        for kk in range(len(new_numbers)):
+            eq_numbers[new_numbers[kk][0]] = new_numbers[kk][1]
+        jj = 0
+        for parm in unsnapped_param_dict:
+            if parm!="p":
+                eq = eq.subs(parm, eq_numbers[jj])
+                jj = jj + 1
+        integer_snapped_expr = integer_snapped_expr + [eq]
+    
+    # Get the numbers appearing in the expression
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    eq_numbers = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+
+    # Do zero snap one parameter at a time
+    zero_snapped_expr = []
+    for w in range(len(eq_numbers)):
+        param_dict = {}
+        unsnapped_param_dict = {'p':1} 
+        eq = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+        new_numbers = zeroSnap(eq_numbers,w+1)
+        for kk in range(len(new_numbers)):
+            eq_numbers[new_numbers[kk][0]] = new_numbers[kk][1]
+        jj = 0
+        for parm in unsnapped_param_dict:
+            if parm!="p":
+                eq = eq.subs(parm, eq_numbers[jj])
+                jj = jj + 1
+        zero_snapped_expr = zero_snapped_expr + [eq]
+
+    # Get the numbers appearing in the expression
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    eq_numbers = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+
+    # Do rational snap one parameter at a time
+    rational_snapped_expr = []
+    for w in range(len(eq_numbers)):
+        eq_numbers_snap = copy.deepcopy(eq_numbers)
+        param_dict = {}
+        unsnapped_param_dict = {'p':1}
+        eq = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+        new_numbers = rationalSnap(eq_numbers,w+1)
+        for kk in range(len(new_numbers)):
+            eq_numbers_snap[new_numbers[kk][0]] = new_numbers[kk][1][1:3]
+        jj = 0
+        for parm in unsnapped_param_dict:
+            if parm!="p":
+                
+                try:
+                    eq = eq.subs(parm, Rational(eq_numbers_snap[jj][0],eq_numbers_snap[jj][1])) 
+                except:
+                    eq = eq.subs(parm, eq_numbers_snap[jj])
+                jj = jj + 1
+        rational_snapped_expr = rational_snapped_expr + [eq]
+        
+    snapped_expr = np.append(integer_snapped_expr,zero_snapped_expr)
+    snapped_expr = np.append(snapped_expr,rational_snapped_expr)
+    
+    
+    for i in range(len(snapped_expr)):
+        try:
+            # Calculate the error of the new, snapped expression
+            snapped_error = get_symbolic_expr_error(pathdir,filename,str(snapped_expr[i]))
+            # Calculate the complexity of the new, snapped expression
+            expr = simplify(powsimp(snapped_expr[i]))
+            for s in (expr.free_symbols):
+                s = symbols(str(s), real = True)
+            expr =  simplify(parse_expr(str(snapped_expr[i]),locals()))
+            is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+            numbers_expr = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+            
+            if DR_file=="":
+                snapped_complexity = 0
+                for j in numbers_expr:
+                    snapped_complexity = snapped_complexity + get_number_DL_snapped(float(j))
+
+                n_variables = len(expr.free_symbols)
+                n_operations = len(count_ops(expr,visual=True).free_symbols)
+                if n_operations!=0 or n_variables!=0:
+                    snapped_complexity = snapped_complexity + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+            
+            # If a bf file is provided, replace the variables with the actual ones before calculating the complexity
+            else:
+                dr_data = np.loadtxt(DR_file,dtype="str",delimiter=",")
+
+                expr = str(expr)
+                old_vars = ["x%s" %k for k in range(len(dr_data)-3)]
+                for i_dr in range(len(old_vars)):
+                    expr = expr.replace(old_vars[i_dr],"("+dr_data[i_dr+2]+")")
+                expr = "("+dr_data[1]+")*(" + expr +")"
+
+                expr = parse_expr(expr)
+                for s in (expr.free_symbols):
+                    s = symbols(str(s), real = True)
+                expr =  simplify(parse_expr(str(expr),locals()))
+
+                snapped_complexity = 0
+                for j in numbers_expr:
+                    snapped_complexity = snapped_complexity + get_number_DL_snapped(float(j))
+                
+                n_variables = len(expr.free_symbols)
+                n_operations = len(count_ops(expr,visual=True).free_symbols)
+                if n_operations!=0 or n_variables!=0:
+                    snapped_complexity = snapped_complexity + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+
+            PA.add(Point(x=snapped_complexity, y=snapped_error, data=str(expr)))
+        except:
+            continue
+    return(PA)
+        
+        
+        
+            

Code/S_add_sym_on_pareto.py

@@ -0,0 +1,29 @@
+# Combines 2 pareto fromtier obtained from the separability test into a new one.
+
+from get_pareto import Point, ParetoSet
+from RPN_to_pytorch import RPN_to_pytorch
+from RPN_to_eq import RPN_to_eq
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+from sympy.parsing.sympy_parser import parse_expr
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from os import path
+from sympy import Symbol, lambdify, N
+from get_pareto import Point, ParetoSet
+
+def add_sym_on_pareto(pathdir,filename,PA1,idx1,idx2,PA,sym_typ):
+    possible_vars = ["x%s" %i for i in np.arange(0,30,1)]
+    PA1 = np.array(PA1.get_pareto_points()).astype('str')
+    for i in range(len(PA1)):
+        exp1 = PA1[i][2]
+        for j in range(len(possible_vars)-2,idx2-1,-1):
+            exp1 = exp1.replace(possible_vars[j],possible_vars[j+1])
+        exp1 = exp1.replace(possible_vars[idx1],"(" + possible_vars[idx1] + sym_typ + possible_vars[idx2] + ")")
+        PA.add(Point(x=float(PA1[i][0]),y=float(PA1[i][1]),data=str(exp1)))
+
+    return PA
+
+
+
+

Code/S_brute_force.py

@@ -0,0 +1,30 @@
+# runs BF on data and saves the best RPN expressions in results.dat
+# all the .dat files are created after I run this script
+# the .scr are needed to run the fortran code
+
+import numpy as np
+import os
+import shutil
+import subprocess
+from subprocess import call
+import sys
+import csv
+import sympy as sp
+from sympy.parsing.sympy_parser import parse_expr
+
+# sep_type = 3 for add and 2 for mult and 1 for normal
+def brute_force(pathdir,filename,BF_try_time,BF_ops_file_type,sep_type="*"):
+    try_time = BF_try_time
+    try_time_prefactor = BF_try_time
+    file_type = BF_ops_file_type
+    try:
+        os.remove("results.dat")
+    except:
+        pass
+    if sep_type=="*":
+        subprocess.call(["./brute_force_oneFile_v2.scr", file_type, "%s" %try_time, pathdir+filename])
+    if sep_type=="+":
+        subprocess.call(["./brute_force_oneFile_v3.scr", file_type, "%s" %try_time, pathdir+filename])
+
+    return 1
+

Code/S_change_output.py

@@ -0,0 +1,273 @@
+import numpy as np
+import os
+from S_run_bf_polyfit import run_bf_polyfit
+
+def get_acos(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.arccos(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)     
+                
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "acos")
+
+    except:
+        return PA
+
+    return PA
+
+def get_asin(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.arcsin(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)       
+                
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "asin")
+
+    except:
+        return PA
+
+    return PA
+
+def get_atan(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.arctan(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt) 
+                
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "atan")                
+
+    except:
+        return PA
+
+    return PA
+
+
+def get_cos(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.cos(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)          
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "cos")                
+                
+    except:
+        return PA
+
+    return PA
+
+
+def get_exp(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.exp(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)       
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "exp")                
+                
+    except:
+        return PA
+
+    return PA
+
+
+def get_inverse(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,1/f_dependent))
+        np.savetxt(pathdir_write_to+filename,dt)          
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "inverse")                    
+                    
+    except:
+        return PA
+
+    return PA
+
+
+def get_log(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.log(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)         
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "log")                    
+                    
+    except:
+        return PA
+
+    return PA
+
+
+def get_sin(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.sin(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)         
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "sin")                
+                
+    except:
+        return PA
+
+    return PA
+
+
+def get_sqrt(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.sqrt(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt) 
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "sqrt")                    
+                    
+    except:
+        return PA
+
+    return PA
+
+
+def get_squared(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+   
+        dt = np.column_stack((variables,f_dependent**2))
+        np.savetxt(pathdir_write_to+filename,dt)
+
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "squared")  
+                
+    except:
+        return PA
+
+    return PA
+
+
+def get_tan(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    try:
+        n_variables = np.loadtxt(pathdir+"%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"%s" %filename, usecols=(0,))
+        for j in range(1,n_variables):
+            v = np.loadtxt(pathdir+"%s" %filename, usecols=(j,))
+            variables = np.column_stack((variables,v))
+        f_dependent = np.loadtxt(pathdir+"%s" %filename, usecols=(n_variables,))
+
+        dt = np.column_stack((variables,np.tan(f_dependent)))
+        np.savetxt(pathdir_write_to+filename,dt)        
+        
+        PA = run_bf_polyfit(pathdir,pathdir_write_to,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg, "tan")                
+                
+    except:
+        return PA
+
+    return PA
+
+
+
+
+
+
+
+
+

Code/S_combine_pareto.py

@@ -0,0 +1,38 @@
+# Combines 2 pareto fromtier obtained from the separability test into a new one.
+
+from get_pareto import Point, ParetoSet
+from RPN_to_pytorch import RPN_to_pytorch
+from RPN_to_eq import RPN_to_eq
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+from sympy.parsing.sympy_parser import parse_expr
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from os import path
+from sympy import Symbol, lambdify, N
+from get_pareto import Point, ParetoSet
+
+def combine_pareto(pathdir,filename,PA1,PA2,idx_list_1,idx_list_2,PA,sep_type = "+"):
+    possible_vars = ["x%s" %i for i in np.arange(0,30,1)]
+    PA1 = np.array(PA1.get_pareto_points()).astype('str')    
+    PA2 = np.array(PA2.get_pareto_points()).astype('str')  
+    for i in range(len(PA1)):
+        for j in range(len(PA2)):
+            try:
+                complexity = float(PA1[i][0])+float(PA2[j][0])
+                # replace the variables from the separated parts with the variables reflecting the new combined equation
+                exp1 = PA1[i][2]
+                exp2 = PA2[j][2]
+                for k in range(len(idx_list_1)-1,-1,-1):
+                    exp1 = exp1.replace(possible_vars[k],possible_vars[idx_list_1[k]])
+                for k in range(len(idx_list_2)-1,-1,-1):
+                    exp2 = exp2.replace(possible_vars[k],possible_vars[idx_list_2[k]])
+                new_eq = "(" + exp1 + ")" + sep_type + "(" + exp2 + ")"
+                PA.add(Point(x=complexity,y=get_symbolic_expr_error(pathdir,filename,new_eq),data=new_eq))
+            except:
+                continue
+    return PA
+
+
+
+

Code/S_get_number_DL.py

@@ -0,0 +1,16 @@
+# Calculates the complexity of a number to be used for the Pareto frontier
+
+import numpy as np 
+
+def get_number_DL(n):
+    epsilon = 1e-10
+    # check if integer
+    if np.isnan(n):
+        return 1000000
+    elif np.abs(n - int(n)) < epsilon:
+        return np.log2(1+abs(n))
+    # check if real
+    else:
+        PrecisionFloorLoss = 1e-14
+        return np.log2(1 + (float(n) / PrecisionFloorLoss) ** 2) / 2
+

Code/S_get_number_DL_snapped.py

@@ -0,0 +1,21 @@
+# Calculates the complexity of a number to be used for the Pareto frontier after snapping
+
+import numpy as np 
+from S_snap import bestApproximation
+
+def get_number_DL_snapped(n):
+    epsilon = 1e-10
+    n = float(n)
+    if np.isnan(n):
+        return 1000000
+    elif np.abs(n - int(n)) < epsilon:
+        return np.log2(1 + abs(int(n)))
+    elif np.abs(n - bestApproximation(n,10000)[0]) < epsilon:
+        _, numerator, denominator, _ = bestApproximation(n, 10000)
+        return np.log2((1 + abs(numerator)) * abs(denominator))
+    else:
+        PrecisionFloorLoss = 1e-14
+        return np.log2(1 + (float(n) / PrecisionFloorLoss) ** 2) / 2
+    
+    
+    

Code/S_get_symbolic_expr_error.py

@@ -0,0 +1,33 @@
+# Calculates the error of a given symbolic expression applied to a dataset. The input should be a string of the mathematical expression
+
+from get_pareto import Point, ParetoSet
+from RPN_to_pytorch import RPN_to_pytorch
+from RPN_to_eq import RPN_to_eq
+from sympy.parsing.sympy_parser import parse_expr
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from os import path
+from sympy import Symbol, lambdify, N
+
+def get_symbolic_expr_error(pathdir,filename,expr):
+    data = np.loadtxt(pathdir+filename)
+    N_vars = len(data[0])-1
+    possible_vars = ["x%s" %i for i in np.arange(0,30,1)]
+    variables = []
+    for i in range(N_vars):
+        variables = variables + [possible_vars[i]]
+    eq = parse_expr(expr)
+    f = lambdify(variables, N(eq))
+    real_variables = []
+
+    for i in range(len(data[0])-1):
+        check_var = "x"+str(i)
+        if check_var in np.array(variables).astype('str'):
+            real_variables = real_variables + [data[:,i]]
+
+    # use this to get rid of cases where the loss gets complex because of transformations of the output variable
+    if isinstance(np.mean((f(*real_variables)-data[:,-1])**2), complex):
+        return 1000000
+    else:
+        return np.sqrt(np.mean((f(*real_variables)-data[:,-1])**2))/np.sqrt(np.mean(data[:,-1]**2))

Code/S_polyfit.py

@@ -0,0 +1,74 @@
+import numpy as np
+import os
+from S_polyfit_utils import getBest
+from S_polyfit_utils import basis_vector
+import itertools
+import sys
+import csv
+import sympy
+from sympy import symbols, Add, Mul, S, simplify
+from scipy.linalg import fractional_matrix_power
+
+def mk_sympy_function(coeffs, num_covariates, deg):
+    generators = [basis_vector(num_covariates+1, i) for i in range(num_covariates+1)]
+    powers = map(sum, itertools.combinations_with_replacement(generators, deg))
+    
+    coeffs = np.round(coeffs,2)
+    
+    xs = (S.One,) + symbols('z0:%d'%num_covariates)
+    if len(coeffs)>1:
+        return Add(*[coeff * Mul(*[x**deg for x, deg in zip(xs, power)])
+                     for power, coeff in zip(powers, coeffs)])
+    else:
+        return coeffs[0]
+
+def polyfit(maxdeg, filename):
+    n_variables = np.loadtxt(filename, dtype='str').shape[1]-1
+    variables = np.loadtxt(filename, usecols=(0,))
+    means = [np.mean(variables)]
+
+    for j in range(1,n_variables):
+        v = np.loadtxt(filename, usecols=(j,))
+        means = means + [np.mean(v)]
+        variables = np.column_stack((variables,v))
+       
+    f_dependent = np.loadtxt(filename, usecols=(n_variables,))
+
+    if n_variables>1:
+        C_1_2 = fractional_matrix_power(np.cov(variables.T),-1/2)
+        
+        x = []
+        z = []
+        for ii in range(len(variables[0])):
+            variables[:,ii] = variables[:,ii] - np.mean(variables[:,ii])
+            x = x + ["x"+str(ii)]
+            z = z + ["z"+str(ii)]
+
+        variables = np.matmul(C_1_2,variables.T).T
+
+        parameters = getBest(variables,f_dependent,maxdeg)[0]
+        params_error = getBest(variables,f_dependent,maxdeg)[1]
+        deg = getBest(variables,f_dependent,maxdeg)[2]     
+
+        x = sympy.Matrix(x)
+        M = sympy.Matrix(C_1_2)
+        b = sympy.Matrix(means)
+        M_x = M*(x-b)
+
+        eq = mk_sympy_function(parameters,n_variables,deg)
+        symb = sympy.Matrix(z)
+
+        for i in range(len(symb)):
+            eq = eq.subs(symb[i],M_x[i])
+
+        eq = simplify(eq)
+        
+    else:
+        parameters = getBest(variables,f_dependent,maxdeg)[0]
+        params_error = getBest(variables,f_dependent,maxdeg)[1]
+        deg = getBest(variables,f_dependent,maxdeg)[2]
+        eq = mk_sympy_function(parameters,n_variables,deg)
+        eq = eq.subs("z0","x0")  
+    
+    return (eq, params_error)
+

Code/S_polyfit_utils.py

@@ -0,0 +1,55 @@
+import numpy as np
+from numpy import linalg, zeros, ones, hstack, asarray
+import itertools
+from matplotlib import pyplot as plt
+from scipy.sparse.linalg import lsqr
+import os
+from sympy import symbols, Add, Mul, S
+
+
+def basis_vector(n, i):
+    x = zeros(n, dtype=int)
+    x[i] = 1
+    return x
+
+def as_tall(x):
+    return x.reshape(x.shape + (1,))
+
+
+def multipolyfit(xs, y, deg):
+    
+    y = asarray(y).squeeze()
+    rows = y.shape[0]
+    xs = asarray(xs)
+    try:
+        num_covariates = xs.shape[1]
+    except:
+        num_covariates = 1
+        xs = np.reshape(xs,(len(xs),1))
+
+    xs = hstack((ones((xs.shape[0], 1), dtype=xs.dtype) , xs))
+    
+    generators = [basis_vector(num_covariates+1, i) for i in range(num_covariates+1)]
+        
+    # All combinations of degrees
+    powers = map(sum, itertools.combinations_with_replacement(generators, deg))
+
+    # Raise data to specified degree pattern, stack in order
+    A = hstack(asarray([as_tall((xs**p).prod(1)) for p in powers]))
+    params = lsqr(A, y)[0] # get the best params of the fit
+    rms = lsqr(A, y)[4] # get the rms params of the fit
+    
+    return (params, rms)
+
+
+def getBest(xs,y,max_deg):
+    results = []
+    for i in range(0,max_deg+1):
+        results = results + [multipolyfit(xs,y,i)]
+    results = np.array(results)
+    # get the parameters and error of the fit with the lowest rms error
+    params = results[np.argmin(results[:,1:])][0]
+    error = results[np.argmin(results[:,1:])][1]
+    deg = np.argmin(results[:,1:])
+    return (params, error, deg)
+

Code/S_run_aifeynman.py

@@ -0,0 +1,173 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from os import path
+from get_pareto import Point, ParetoSet
+from RPN_to_pytorch import RPN_to_pytorch
+from RPN_to_eq import RPN_to_eq
+from S_NN_train import NN_train
+from S_NN_eval import NN_eval
+from S_symmetry import *
+from S_separability import *
+from S_change_output import *
+from S_brute_force import brute_force
+from S_combine_pareto import combine_pareto
+from S_get_number_DL import get_number_DL
+from sympy.parsing.sympy_parser import parse_expr
+from sympy import preorder_traversal, count_ops
+from S_polyfit import polyfit
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+from S_add_snap_expr_on_pareto import add_snap_expr_on_pareto
+from S_add_sym_on_pareto import add_sym_on_pareto
+from S_run_bf_polyfit import run_bf_polyfit
+
+
+PA = ParetoSet()
+
+def run_AI_all(pathdir,filename,BF_try_time=60,BF_ops_file_type="14ops", polyfit_deg=4, NN_epochs=4000, PA = PA):
+    try:
+        os.mkdir("results/")
+    except:
+        pass
+    
+    # load the data for different checks
+    data = np.loadtxt(pathdir+filename)
+
+    # Run bf and polyfit
+    PA = run_bf_polyfit(pathdir,pathdir,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+
+
+
+    # Run bf and polyfit on modified output
+    PA = get_acos(pathdir,"results/mystery_world_acos/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_asin(pathdir,"results/mystery_world_asin/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_atan(pathdir,"results/mystery_world_atan/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_cos(pathdir,"results/mystery_world_cos/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_exp(pathdir,"results/mystery_world_exp/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_inverse(pathdir,"results/mystery_world_inverse/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_log(pathdir,"results/mystery_world_log/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_sin(pathdir,"results/mystery_world_sin/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_sqrt(pathdir,"results/mystery_world_sqrt/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_squared(pathdir,"results/mystery_world_squared/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+    PA = get_tan(pathdir,"results/mystery_world_tan/",filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg)
+
+#############################################################################################################################  
+    # check if the NN is trained. If it is not, train it on the data.
+    print("Checking for symmetry \n", filename)
+    if path.exists("results/NN_trained_models/models/" + filename + ".h5") or len(data[0])<3:
+        print("NN already trained \n")
+        print("NN loss: ", NN_eval(pathdir,filename), "\n")
+    else:
+        print("Training a NN on the data... \n")
+        NN_train(pathdir,filename,NN_epochs)
+        print("NN loss: ", NN_eval(pathdir,filename), "\n")
+        
+    # Check which symmetry/separability is the best
+    
+    # Symmetries
+    symmetry_minus_result = check_translational_symmetry_minus(pathdir,filename)
+    symmetry_divide_result = check_translational_symmetry_divide(pathdir,filename)
+    symmetry_multiply_result = check_translational_symmetry_multiply(pathdir,filename)
+    symmetry_plus_result = check_translational_symmetry_plus(pathdir,filename)
+    
+    # Separabilities
+    separability_plus_result = check_separability_plus(pathdir,filename)
+    separability_multiply_result = check_separability_multiply(pathdir,filename)
+    
+    if symmetry_plus_result[0]==-1:
+        idx_min = -1
+    else:
+        idx_min = np.argmin(np.array([symmetry_plus_result[0], symmetry_minus_result[0], symmetry_multiply_result[0], symmetry_divide_result[0], separability_plus_result[0], separability_multiply_result[0]]))
+
+    # Apply the best symmetry/separability and rerun the main function on this new file
+    if idx_min == 0:
+        new_pathdir, new_filename = do_translational_symmetry_plus(pathdir,filename,symmetry_plus_result[1],symmetry_plus_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir,new_filename,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA = add_sym_on_pareto(pathdir,filename,PA1,symmetry_plus_result[1],symmetry_plus_result[2],PA,"+")
+        return PA
+    
+    elif idx_min == 1:
+        new_pathdir, new_filename = do_translational_symmetry_minus(pathdir,filename,symmetry_minus_result[1],symmetry_minus_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir,new_filename,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA = add_sym_on_pareto(pathdir,filename,PA1,symmetry_minus_result[1],symmetry_minus_result[2],PA,"-")
+        return PA
+    
+    elif idx_min == 2:
+        new_pathdir, new_filename = do_translational_symmetry_multiply(pathdir,filename,symmetry_multiply_result[1],symmetry_multiply_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir,new_filename,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA = add_sym_on_pareto(pathdir,filename,PA1,symmetry_multiply_result[1],symmetry_multiply_result[2],PA,"*")
+        return PA
+    
+    elif idx_min == 3:
+        new_pathdir, new_filename = do_translational_symmetry_divide(pathdir,filename,symmetry_divide_result[1],symmetry_divide_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir,new_filename,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA = add_sym_on_pareto(pathdir,filename,PA1,symmetry_divide_result[1],symmetry_divide_result[2],PA,"/")
+        return PA
+    
+    elif idx_min == 4:
+        new_pathdir1, new_filename1, new_pathdir2, new_filename2,  = do_separability_plus(pathdir,filename,separability_plus_result[1],separability_plus_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir1,new_filename1,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA2_ = ParetoSet()
+        PA2 = run_AI_all(new_pathdir2,new_filename2,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA2_)
+        PA = combine_pareto(pathdir,filename,PA1,PA2,separability_plus_result[1],separability_plus_result[2],PA,"+")
+        return PA 
+    
+    elif idx_min == 5:
+        new_pathdir1, new_filename1, new_pathdir2, new_filename2,  = do_separability_multiply(pathdir,filename,separability_multiply_result[1],separability_multiply_result[2])
+        PA1_ = ParetoSet()
+        PA1 = run_AI_all(new_pathdir1,new_filename1,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA1_)
+        PA2_ = ParetoSet()
+        PA2 = run_AI_all(new_pathdir2,new_filename2,BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs, PA2_)
+        PA = combine_pareto(pathdir,filename,PA1,PA2,separability_multiply_result[1],separability_multiply_result[2],PA,"*")
+        return PA 
+    else:
+        return PA
+
+# this runs snap on the output of aifeynman
+def run_aifeynman(pathdir,filename,BF_try_time,BF_ops_file_type, polyfit_deg=4, NN_epochs=4000, DR_file=""):
+    # Split the data into train and test set                                                                                                                                      
+    input_data = np.loadtxt(pathdir+filename)
+    sep_idx = np.random.permutation(len(input_data))
+    train_data = input_data[sep_idx[0:8*len(input_data)//10]]
+    test_data = input_data[sep_idx[8*len(input_data)//10:len(input_data)]]
+
+    np.savetxt(pathdir+filename+"_train",train_data)
+    np.savetxt(pathdir+filename+"_test",test_data)
+
+    # Run the code on the train data 
+    PA = run_AI_all(pathdir,filename+"_train",BF_try_time,BF_ops_file_type, polyfit_deg, NN_epochs)
+    PA_list = PA.get_pareto_points()
+    PA_snapped = ParetoSet()
+    
+    np.savetxt("results/solution_before_snap_%s.txt" %filename,PA_list,fmt="%s")
+    for j in range(len(PA_list)):
+        PA_snapped = add_snap_expr_on_pareto(pathdir,filename,PA_list[j][-1],PA_snapped, DR_file)
+
+    list_dt = np.array(PA_snapped.get_pareto_points())
+    
+    data_file_len = len(np.loadtxt(pathdir+filename))
+    log_err = []
+    log_err_all = []
+    for i in range(len(list_dt)):
+        log_err = log_err + [np.log2(float(list_dt[i][1]))]
+        log_err_all = log_err_all + [data_file_len*np.log2(float(list_dt[i][1]))]
+    log_err = np.array(log_err)
+    log_err_all = np.array(log_err_all)
+
+    # Try the found expressions on the test data                                                                                                                                  
+    if DR_file=="":
+        test_errors = []
+        for i in range(len(list_dt)):
+            test_errors = test_errors + [get_symbolic_expr_error(pathdir,filename+"_test",str(list_dt[i][-1]))]
+        test_errors = np.array(test_errors)
+        # Save all the data to file                                                                                                                                               
+        save_data = np.column_stack((test_errors,log_err,log_err_all,list_dt))
+    else:
+        save_data = np.column_stack((log_err,log_err_all,list_dt))
+
+    np.savetxt("results/solution_%s.txt" %filename,save_data,fmt="%s")

Code/S_run_bf_polyfit.py

@@ -0,0 +1,225 @@
+# add a function to compte complexity
+
+from get_pareto import Point, ParetoSet
+from RPN_to_pytorch import RPN_to_pytorch
+from RPN_to_eq import RPN_to_eq
+import numpy as np
+import matplotlib.pyplot as plt
+from S_brute_force import brute_force
+from S_combine_pareto import combine_pareto
+from S_get_number_DL import get_number_DL
+from sympy.parsing.sympy_parser import parse_expr
+from sympy import preorder_traversal, count_ops
+from S_polyfit import polyfit
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+from S_add_sym_on_pareto import add_sym_on_pareto
+import os
+from os import path
+
+
+def run_bf_polyfit(pathdir,pathdir_transformed,filename,BF_try_time,BF_ops_file_type, PA, polyfit_deg=4, output_type=""):
+    
+#############################################################################################################################
+    
+    # run BF on the data (+)
+    print("Checking for brute force + \n")
+    brute_force(pathdir_transformed,filename,BF_try_time,BF_ops_file_type,"+")
+    
+    # load the BF output data
+    bf_all_output = np.loadtxt("results.dat", dtype="str")
+    express = bf_all_output[:,2]
+    prefactors = bf_all_output[:,1]
+    prefactors = [str(i) for i in prefactors] 
+    
+    # Calculate the complexity of the bf expression the same way as for gradient descent case
+    complexity = []
+    errors = []
+    eqns = []
+    for i in range(len(prefactors)):
+        try:            
+            if output_type=="":
+                eqn = prefactors[i] + "+" + RPN_to_eq(express[i])
+            elif output_type=="acos":
+                eqn = "cos(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="asin":
+                eqn = "sin(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")" 
+            elif output_type=="atan":
+                eqn = "tan(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="cos":
+                eqn = "acos(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="exp":
+                eqn = "log(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="inverse":
+                eqn = "1/(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="log":
+                eqn = "exp(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="sin":
+                eqn = "acos(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="sqrt":
+                eqn = "(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")**2"
+            elif output_type=="squared":
+                eqn = "sqrt(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="tan":
+                eqn = "atan(" + prefactors[i] + "+" + RPN_to_eq(express[i]) + ")"
+                
+            eqns = eqns + [eqn]
+            errors = errors + [get_symbolic_expr_error(pathdir,filename,eqn)]
+            expr = parse_expr(eqn)
+            is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+            numbers_expr = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+            compl = 0
+            for j in numbers_expr:
+                try:
+                    compl = compl + get_number_DL(float(j))
+                except:
+                    compl = compl + 1000000
+
+            # Add the complexity due to symbols
+            n_variables = len(expr.free_symbols)
+            n_operations = len(count_ops(expr,visual=True).free_symbols)
+            if n_operations!=0 or n_variables!=0:
+                compl = compl + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+
+            complexity = complexity + [compl]
+        except:
+            continue
+
+    for i in range(len(complexity)):
+        PA.add(Point(x=complexity[i], y=errors[i], data=eqns[i]))
+    
+    # run gradient descent of BF output parameters and add the results to the Pareto plot
+    for i in range(len(express)):
+        try:
+            bf_gd_update = RPN_to_pytorch(pathdir+filename,eqns[i])
+            PA.add(Point(x=bf_gd_update[1],y=bf_gd_update[0],data=bf_gd_update[2]))
+        except:
+            continue
+
+
+#############################################################################################################################
+    # run BF on the data (*)
+    print("Checking for brute force * \n")
+    brute_force(pathdir_transformed,filename,BF_try_time,BF_ops_file_type,"*")
+    
+    # load the BF output data
+    bf_all_output = np.loadtxt("results.dat", dtype="str")
+    express = bf_all_output[:,2]
+    prefactors = bf_all_output[:,1]
+    prefactors = [str(i) for i in prefactors] 
+                
+    # Calculate the complexity of the bf expression the same way as for gradient descent case
+    complexity = []
+    errors = []
+    eqns = []
+    for i in range(len(prefactors)):
+        try: 
+            if output_type=="":
+                eqn = prefactors[i] + "*" + RPN_to_eq(express[i])
+            elif output_type=="acos":
+                eqn = "cos(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="asin":
+                eqn = "sin(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")" 
+            elif output_type=="atan":
+                eqn = "tan(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="cos":
+                eqn = "acos(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="exp":
+                eqn = "log(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="inverse":
+                eqn = "1/(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="log":
+                eqn = "exp(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="sin":
+                eqn = "acos(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="sqrt":
+                eqn = "(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")**2"
+            elif output_type=="squared":
+                eqn = "sqrt(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            elif output_type=="tan":
+                eqn = "atan(" + prefactors[i] + "*" + RPN_to_eq(express[i]) + ")"
+            
+            eqns = eqns + [eqn]
+            errors = errors + [get_symbolic_expr_error(pathdir,filename,eqn)]
+            expr = parse_expr(eqn)
+            is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+            numbers_expr = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+            compl = 0
+            for j in numbers_expr:
+                try:
+                    compl = compl + get_number_DL(float(j))
+                except:
+                    compl = compl + 1000000
+
+            # Add the complexity due to symbols
+            n_variables = len(expr.free_symbols)
+            n_operations = len(count_ops(expr,visual=True).free_symbols)
+            if n_operations!=0 or n_variables!=0:
+                compl = compl + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+
+            complexity = complexity + [compl]
+        except:
+            continue
+
+    # add the BF output to the Pareto plot
+    for i in range(len(complexity)):
+        PA.add(Point(x=complexity[i], y=errors[i], data=eqns[i]))
+
+    # run gradient descent of BF output parameters and add the results to the Pareto plot
+    for i in range(len(express)):
+        try:
+            bf_gd_update = RPN_to_pytorch(pathdir+filename,eqns[i])
+            PA.add(Point(x=bf_gd_update[1],y=bf_gd_update[0],data=bf_gd_update[2]))
+        except:
+            continue
+
+#############################################################################################################################
+    # run polyfit on the data
+    print("Checking polyfit \n")
+    polyfit_result = polyfit(polyfit_deg, pathdir_transformed+filename)
+    eqn = str(polyfit_result[0])
+    
+    # Calculate the complexity of the polyfit expression the same way as for gradient descent case    
+    if output_type=="":
+        eqn = eqn
+    elif output_type=="acos":
+        eqn = "cos(" + eqn + ")"
+    elif output_type=="asin":
+        eqn = "sin(" + eqn + ")" 
+    elif output_type=="atan":
+        eqn = "tan(" + eqn + ")"
+    elif output_type=="cos":
+        eqn = "acos(" + eqn + ")"
+    elif output_type=="exp":
+        eqn = "log(" + eqn + ")"
+    elif output_type=="inverse":
+        eqn = "1/(" + eqn + ")"
+    elif output_type=="log":
+        eqn = "exp(" + eqn + ")"
+    elif output_type=="sin":
+        eqn = "acos(" + eqn + ")"
+    elif output_type=="sqrt":
+        eqn = "(" + eqn + ")**2"
+    elif output_type=="squared":
+        eqn = "sqrt(" + eqn + ")"
+    elif output_type=="tan":
+        eqn = "atan(" + eqn + ")"
+    
+    polyfit_err = get_symbolic_expr_error(pathdir,filename,eqn)
+    expr = parse_expr(eqn)
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    numbers_expr = [subexpression for subexpression in preorder_traversal(expr) if is_atomic_number(subexpression)]
+    complexity = 0
+    for j in numbers_expr:
+        complexity = complexity + get_number_DL(float(j)) 
+    # Add the complexity due to symbols
+    n_variables = len(polyfit_result[0].free_symbols)
+    n_operations = len(count_ops(polyfit_result[0],visual=True).free_symbols)
+    if n_operations!=0 or n_variables!=0:
+        complexity = complexity + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+    
+    PA.add(Point(x=complexity, y=polyfit_err, data=str(eqn)))
+
+    for pareto_i in range(len(PA.get_pareto_points())):
+        print(PA.get_pareto_points()[pareto_i])
+    
+    return PA

Code/S_separability.py

@@ -0,0 +1,379 @@
+from __future__ import print_function
+import torch
+import os
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+from torchvision import datasets, transforms
+import pandas as pd
+import numpy as np
+import torch
+from torch.utils import data
+import pickle
+from torch.optim.lr_scheduler import CosineAnnealingLR
+from matplotlib import pyplot as plt
+from itertools import combinations
+import time
+
+is_cuda = torch.cuda.is_available()
+
+class SimpleNet(nn.Module):
+    def __init__(self, ni):
+        super().__init__()
+        self.linear1 = nn.Linear(ni, 128)
+        self.bn1 = nn.BatchNorm1d(128)
+        self.linear2 = nn.Linear(128, 128)
+        self.bn2 = nn.BatchNorm1d(128)
+        self.linear3 = nn.Linear(128, 64)
+        self.bn3 = nn.BatchNorm1d(64)
+        self.linear4 = nn.Linear(64,64)
+        self.bn4 = nn.BatchNorm1d(64)
+        self.linear5 = nn.Linear(64,1)
+        
+    def forward(self, x):
+        x = F.tanh(self.bn1(self.linear1(x)))
+        x = F.tanh(self.bn2(self.linear2(x)))
+        x = F.tanh(self.bn3(self.linear3(x)))
+        x = F.tanh(self.bn4(self.linear4(x)))
+        x = self.linear5(x)
+        return x  
+
+def rmse_loss(pred, targ):
+    denom = targ**2
+    denom = torch.sqrt(denom.sum()/len(denom))
+    return torch.sqrt(F.mse_loss(pred, targ))/denom
+
+def check_separability_plus(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        # make some variables at the time equal to the median of factors
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():
+            fact_vary = factors.clone()
+            for k in range(len(factors[0])):
+                fact_vary[:,k] = torch.full((len(factors),),torch.median(factors[:,k]))
+
+            # loop through all indices combinations
+            var_indices_list = np.arange(0,n_variables,1)
+            min_error = 1000
+            best_i = []
+            best_j = []
+            for i in range(1,n_variables):
+                c = combinations(var_indices_list, i)
+                for j in c:
+                    fact_vary_one = factors.clone()
+                    fact_vary_rest = factors.clone()
+                    rest_indx = list(filter(lambda x: x not in j, var_indices_list))
+                    for t1 in rest_indx:
+                        fact_vary_one[:,t1] = torch.full((len(factors),),torch.median(factors[:,t1]))
+                    for t2 in j:
+                        fact_vary_rest[:,t2] = torch.full((len(factors),),torch.median(factors[:,t2]))
+                    # check if the equation is separable
+                    sm = model(fact_vary_one)+model(fact_vary_rest)
+                    #error = torch.sqrt(torch.mean((product-sm+model(fact_vary))**2))/torch.sqrt(torch.mean(product**2))
+                    error = 2*torch.median(abs(product-sm+model(fact_vary)))
+                    if error<min_error:
+                        min_error = error
+                        best_i = j
+                        best_j = rest_indx
+                        
+        return min_error, best_i, best_j
+                    
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)                    
+                    
+                                           
+def do_separability_plus(pathdir, filename, list_i,list_j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        # make some variables at the time equal to the median of factors
+        models_one = []
+        models_rest = []
+        
+        fact_vary = factors.clone()
+        for k in range(len(factors[0])):
+            fact_vary[:,k] = torch.full((len(factors),),torch.median(factors[:,k]))
+        fact_vary_one = factors.clone()
+        fact_vary_rest = factors.clone()
+        for t1 in list_j:
+            fact_vary_one[:,t1] = torch.full((len(factors),),torch.median(factors[:,t1]))
+        for t2 in list_i:
+            fact_vary_rest[:,t2] = torch.full((len(factors),),torch.median(factors[:,t2]))
+
+        with torch.no_grad():
+            str1 = filename+"-add_a"
+            str2 = filename+"-add_b"
+            # save the first half
+            data_sep_1 = variables
+            data_sep_1 = np.delete(data_sep_1,list_j,axis=1)
+            data_sep_1 = np.column_stack((data_sep_1,model(fact_vary_one).cpu()))
+            # save the second half  
+            data_sep_2 = variables
+            data_sep_2 = np.delete(data_sep_2,list_i,axis=1)
+            data_sep_2 = np.column_stack((data_sep_2,model(fact_vary_rest).cpu()-model(fact_vary).cpu()))
+            try:
+                os.mkdir("results/separable_add/")
+            except:
+                pass
+            np.savetxt("results/separable_add/"+str1,data_sep_1)
+            np.savetxt("results/separable_add/"+str2,data_sep_2)
+            # if it is separable, return the 2 new files created and the index of the column with the separable variable
+            return ("results/separable_add/",str1,"results/separable_add/",str2)
+
+    except Exception as e:
+        print(e)
+        return (-1,-1)
+
+    
+
+
+    
+def check_separability_multiply(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables)
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        # make some variables at the time equal to the median of factors
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():
+            fact_vary = factors.clone()
+            for k in range(len(factors[0])):
+                fact_vary[:,k] = torch.full((len(factors),),torch.median(factors[:,k]))
+
+            # loop through all indices combinations
+            var_indices_list = np.arange(0,n_variables,1)
+            min_error = 1000
+            best_i = []
+            best_j = []
+            for i in range(1,n_variables):
+                c = combinations(var_indices_list, i)
+                for j in c:
+                    fact_vary_one = factors.clone()
+                    fact_vary_rest = factors.clone()
+                    rest_indx = list(filter(lambda x: x not in j, var_indices_list))
+                    for t1 in rest_indx:
+                        fact_vary_one[:,t1] = torch.full((len(factors),),torch.median(factors[:,t1]))
+                    for t2 in j:
+                        fact_vary_rest[:,t2] = torch.full((len(factors),),torch.median(factors[:,t2]))
+                    # check if the equation is separable
+                    pd = model(fact_vary_one)*model(fact_vary_rest)
+                    #error = torch.sqrt(torch.mean((product-pd/model(fact_vary))**2))/torch.sqrt(torch.mean(product**2))
+                    error = 2*torch.median(abs(product-pd/model(fact_vary)))
+                    if error<min_error:
+                        min_error = error
+                        best_i = j
+                        best_j = rest_indx
+                        
+        return min_error, best_i, best_j
+                    
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)                         
+
+                    
+                    
+def do_separability_multiply(pathdir, filename, list_i,list_j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables)
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        # make some variables at the time equal to the median of factors
+        models_one = []
+        models_rest = []                    
+
+        fact_vary = factors.clone()
+        for k in range(len(factors[0])):
+            fact_vary[:,k] = torch.full((len(factors),),torch.median(factors[:,k]))
+        fact_vary_one = factors.clone()
+        fact_vary_rest = factors.clone()
+        for t1 in list_j:
+            fact_vary_one[:,t1] = torch.full((len(factors),),torch.median(factors[:,t1]))
+        for t2 in list_i:
+            fact_vary_rest[:,t2] = torch.full((len(factors),),torch.median(factors[:,t2]))        
+        
+        with torch.no_grad():
+            str1 = filename+"-mult_a"
+            str2 = filename+"-mult_b"
+            # save the first half
+            data_sep_1 = variables
+            data_sep_1 = np.delete(data_sep_1,list_j,axis=1)
+            data_sep_1 = np.column_stack((data_sep_1,model(fact_vary_one).cpu()))
+            # save the second half  
+            data_sep_2 = variables
+            data_sep_2 = np.delete(data_sep_2,list_i,axis=1)
+            data_sep_2 = np.column_stack((data_sep_2,model(fact_vary_rest).cpu()/model(fact_vary).cpu()))
+            try:
+                os.mkdir("results/separable_mult/")
+            except:
+                pass
+            np.savetxt("results/separable_mult/"+str1,data_sep_1)
+            np.savetxt("results/separable_mult/"+str2,data_sep_2)
+            # if it is separable, return the 2 new files created and the index of the column with the separable variable
+            return ("results/separable_mult/",str1,"results/separable_mult/",str2)
+
+    except Exception as e:
+        print(e)
+        return (-1,-1)
+
+        

Code/S_snap.py

@@ -0,0 +1,78 @@
+# The following are snap functions for finding a best approximated integer or rational number for a real number:
+
+import numpy as np
+
+def bestApproximation(x,imax):
+    # The input is a numpy parameter vector p.
+    # The output is an integer specifying which parameter to change,
+    # and a float specifying the new value.
+    def float2contfrac(x,nmax):
+        x = float(x)
+        c = [np.floor(x)];
+        y = x - np.floor(x)
+        k = 0
+        while np.abs(y)!=0 and k<nmax:
+            y = 1 / float(y)
+            i = np.floor(y)
+            c.append(i)
+            y = y - i
+            k = k + 1
+        return c
+    
+    def contfrac2frac(seq):
+        ''' Convert the simple continued fraction in `seq`
+            into a fraction, num / den
+            '''
+        num, den = 1, 0
+        for u in reversed(seq):
+            num, den = den + num*u, num
+        return num, den
+    
+    def contFracRationalApproximations(c):
+        return np.array(list(contfrac2frac(c[:i+1]) for i in range(len(c))))
+    
+    def contFracApproximations(c):
+        q = contFracRationalApproximations(c)
+        return q[:,0] / float(q[:,1])
+    
+    def truncateContFrac(q,imax):
+        k = 0
+        while k < len(q) and np.maximum(np.abs(q[k,0]), q[k,1]) <= imax:
+            k = k + 1
+        return q[:k]
+    
+    def pval(p):
+        p = p.astype(float)
+        return 1 - np.exp(-p ** 0.87 / 0.36)
+    
+    xsign = np.sign(x)
+    q = truncateContFrac(contFracRationalApproximations(float2contfrac(abs(x),20)),imax)
+    
+    if len(q) > 0:
+        p = np.abs(q[:,0] / q[:,1] - abs(x)).astype(float) * (1 + np.abs(q[:,0])) * q[:,1]
+        p = pval(p)
+        i = np.argmin(p)
+        return (xsign * q[i,0] / float(q[i,1]), xsign* q[i,0], q[i,1], p[i])
+    else:
+        return (None, 0, 0, 1)
+
+def integerSnap(p, top=1):
+    p = np.array(p)
+    metric = np.abs(p - np.round(p.astype(np.double)))
+    chosen = np.argsort(metric)[:top]
+    return list(zip(chosen, np.round(p.astype(np.double))[chosen]))
+
+
+def zeroSnap(p, top=1):
+    p = np.array(p)
+    metric = np.abs(p)
+    chosen = np.argsort(metric)[:top]
+    return list(zip(chosen, np.zeros(len(chosen))))
+
+
+def rationalSnap(p, top=1):
+    """Snap to nearest rational number using continued fraction."""
+    p = np.array(p)
+    snaps = np.array(list(bestApproximation(x,100) for x in p))
+    chosen = np.argsort(snaps[:, 3])[:top]
+    return list(zip(chosen, snaps[chosen, 0:3]))

Code/S_symmetry.py

@@ -0,0 +1,554 @@
+# checks for symmetries in the data
+
+from __future__ import print_function
+import torch
+import os
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+from torchvision import datasets, transforms
+import pandas as pd
+import numpy as np
+import torch
+from torch.utils import data
+import pickle
+from torch.optim.lr_scheduler import CosineAnnealingLR
+from matplotlib import pyplot as plt
+import time
+
+is_cuda = torch.cuda.is_available()
+
+class SimpleNet(nn.Module):
+    def __init__(self, ni):
+        super().__init__()
+        self.linear1 = nn.Linear(ni, 128)
+        self.bn1 = nn.BatchNorm1d(128)
+        self.linear2 = nn.Linear(128, 128)
+        self.bn2 = nn.BatchNorm1d(128)
+        self.linear3 = nn.Linear(128, 64)
+        self.bn3 = nn.BatchNorm1d(64)
+        self.linear4 = nn.Linear(64,64)
+        self.bn4 = nn.BatchNorm1d(64)
+        self.linear5 = nn.Linear(64,1)
+    
+    def forward(self, x):
+        x = F.tanh(self.bn1(self.linear1(x)))
+        x = F.tanh(self.bn2(self.linear2(x)))
+        x = F.tanh(self.bn3(self.linear3(x)))
+        x = F.tanh(self.bn4(self.linear4(x)))
+        x = self.linear5(x)
+        return x
+
+def rmse_loss(pred, targ):
+    denom = targ**2
+    denom = torch.sqrt(denom.sum()/len(denom))
+    return torch.sqrt(F.mse_loss(pred, targ))/denom
+
+# checks if f(x,y)=f(x-y)
+def check_translational_symmetry_minus(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD \n")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+"/%s" %filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():            
+            # make the shift x->x+a for 2 variables at a time (different variables)
+            min_error = 1000
+            best_i = -1
+            best_j = -1
+            for i in range(0,n_variables,1):
+                for j in range(0,n_variables,1):
+                    if i<j:
+                        fact_translate = factors.clone()
+                        a = 0.5*min(torch.std(fact_translate[:,i]),torch.std(fact_translate[:,j]))
+                        fact_translate[:,i] = fact_translate[:,i] + a
+                        fact_translate[:,j] = fact_translate[:,j] + a
+                        error = torch.median(abs(product-model(fact_translate)))
+                        if error<min_error:
+                            min_error = error
+                            best_i = i
+                            best_j = j
+        return min_error, best_i, best_j
+
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)
+    
+def do_translational_symmetry_minus(pathdir, filename, i,j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        for k in range(1,n_variables):
+            v = np.loadtxt(pathdir+"/%s" %filename, usecols=(k,))
+            variables = np.column_stack((variables,v))
+        
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():   
+            file_name = filename + "-translated_minus"
+            data_translated = variables
+            data_translated[:,i] = variables[:,i]-variables[:,j]
+            data_translated =  np.delete(data_translated, j, axis=1)
+            data_translated = np.column_stack((data_translated,f_dependent))
+            try:
+                os.mkdir("results/translated_data_minus/")
+            except:
+                pass
+            np.savetxt("results/translated_data_minus/"+file_name , data_translated)
+            return ("results/translated_data_minus/",file_name)
+        
+    except Exception as e:
+        print(e)
+        return (-1,-1)
+        
+
+# checks if f(x,y)=f(x/y)
+def check_translational_symmetry_divide(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD \n")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+"/%s" %filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():
+            a = 1.2
+            min_error = 1000
+            best_i = -1
+            best_j = -1
+            # make the shift x->x*a and y->y*a for 2 variables at a time (different variables)
+            for i in range(0,n_variables,1):
+                for j in range(0,n_variables,1):
+                    if i<j:
+                        fact_translate = factors.clone()
+                        fact_translate[:,i] = fact_translate[:,i]*a
+                        fact_translate[:,j] = fact_translate[:,j]*a
+                        error = torch.median(abs(product-model(fact_translate)))
+                        if error<min_error:
+                            min_error = error
+                            best_i = i
+                            best_j = j
+        return min_error, best_i, best_j                        
+
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)
+    
+    
+def do_translational_symmetry_divide(pathdir, filename, i,j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        for k in range(1,n_variables):
+            v = np.loadtxt(pathdir+"/%s" %filename, usecols=(k,))
+            variables = np.column_stack((variables,v))
+        
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():             
+            file_name = filename + "-translated_divide"
+            data_translated = variables
+            data_translated[:,i] = variables[:,i]/variables[:,j]
+            data_translated =  np.delete(data_translated, j, axis=1)
+            data_translated = np.column_stack((data_translated,f_dependent))
+            try:
+                os.mkdir("results/translated_data_divide/")
+            except:
+                pass
+            np.savetxt("results/translated_data_divide/"+file_name , data_translated)
+            return ("results/translated_data_divide/",file_name)
+        
+    except Exception as e:
+        print(e)
+        return (-1,1)
+
+# checks if f(x,y)=f(x*y)
+def check_translational_symmetry_multiply(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD \n")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+"/%s" %filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():
+            a = 1.2
+            min_error = 1000
+            best_i = -1
+            best_j = -1
+            # make the shift x->x*a and y->y/a for 2 variables at a time (different variables)
+            for i in range(0,n_variables,1):
+                for j in range(0,n_variables,1):
+                    if i<j:
+                        fact_translate = factors.clone()
+                        fact_translate[:,i] = fact_translate[:,i]*a
+                        fact_translate[:,j] = fact_translate[:,j]/a
+                        error = torch.median(abs(product-model(fact_translate)))
+                        if error<min_error:
+                            min_error = error
+                            best_i = i
+                            best_j = j
+        return min_error, best_i, best_j                      
+
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)
+
+def do_translational_symmetry_multiply(pathdir, filename, i,j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        for k in range(1,n_variables):
+            v = np.loadtxt(pathdir+"/%s" %filename, usecols=(k,))
+            variables = np.column_stack((variables,v))
+        
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():             
+            file_name = filename + "-translated_multiply"
+            data_translated = variables
+            data_translated[:,i] = variables[:,i]*variables[:,j]
+            data_translated =  np.delete(data_translated, j, axis=1)
+            data_translated = np.column_stack((data_translated,f_dependent))
+            try:
+                os.mkdir("results/translated_data_multiply/")
+            except:
+                pass
+            np.savetxt("results/translated_data_multiply/"+file_name , data_translated)
+            return ("results/translated_data_multiply/",file_name)
+        
+    except Exception as e:
+        print(e)
+        return (-1,1)
+
+# checks if f(x,y)=f(x+y)
+def check_translational_symmetry_plus(pathdir, filename):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        if n_variables==1:
+            print(filename, "just one variable for ADD \n")
+            # if there is just one variable you have nothing to separate
+            return (-1,-1,-1)
+        else:
+            for j in range(1,n_variables):
+                v = np.loadtxt(pathdir+"/%s" %filename, usecols=(j,))
+                variables = np.column_stack((variables,v))
+        
+
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+                product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():
+            min_error = 1000
+            best_i = -1
+            best_j = -1
+            for i in range(0,n_variables,1):
+                for j in range(0,n_variables,1):
+                    if i<j:
+                        fact_translate = factors.clone()
+                        a = 0.5*min(torch.std(fact_translate[:,i]),torch.std(fact_translate[:,j]))
+                        fact_translate[:,i] = fact_translate[:,i] + a
+                        fact_translate[:,j] = fact_translate[:,j] - a
+                        error = torch.median(abs(product-model(fact_translate)))
+                        if error<min_error:
+                            min_error = error
+                            best_i = i
+                            best_j = j
+        return min_error, best_i, best_j                       
+
+    except Exception as e:
+        print(e)
+        return (-1,-1,-1)
+    
+def do_translational_symmetry_plus(pathdir, filename, i,j):
+    try:
+        pathdir_weights = "results/NN_trained_models/models/"
+
+        # load the data
+        n_variables = np.loadtxt(pathdir+"/%s" %filename, dtype='str').shape[1]-1
+        variables = np.loadtxt(pathdir+"/%s" %filename, usecols=(0,))
+
+        for k in range(1,n_variables):
+            v = np.loadtxt(pathdir+"/%s" %filename, usecols=(k,))
+            variables = np.column_stack((variables,v))
+        
+        f_dependent = np.loadtxt(pathdir+"/%s" %filename, usecols=(n_variables,))
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+
+        factors = torch.from_numpy(variables) 
+        if is_cuda:
+            factors = factors.cuda()
+        else:
+            factors = factors
+        factors = factors.float()
+
+        product = torch.from_numpy(f_dependent)
+        if is_cuda:
+            product = product.cuda()
+        else:
+            product = product
+        product = product.float()
+
+        # load the trained model and put it in evaluation mode
+        if is_cuda:
+            model = SimpleNet(n_variables).cuda()
+        else:
+            model = SimpleNet(n_variables)
+        model.load_state_dict(torch.load(pathdir_weights+filename+".h5"))
+        model.eval()
+
+        models_one = []
+        models_rest = []
+
+        with torch.no_grad():             
+            file_name = filename + "-translated_plus"
+            data_translated = variables
+            data_translated[:,i] = variables[:,i]+variables[:,j]
+            data_translated =  np.delete(data_translated, j, axis=1)
+            data_translated = np.column_stack((data_translated,f_dependent))
+            try:
+                os.mkdir("results/translated_data_plus/")
+            except:
+                pass
+            np.savetxt("results/translated_data_plus/"+file_name , data_translated)
+            return ("results/translated_data_plus/", file_name)
+        
+    except Exception as e:
+        print(e)
+        return (-1,-1)

Code/ai_feynman_example.py

@@ -0,0 +1,3 @@
+from S_run_aifeynman import run_aifeynman
+
+run_aifeynman("../example_data/","example3.txt",30,"14ops.txt", polyfit_deg=4, NN_epochs=400)

Code/ai_feynman_terminal_example.py

@@ -0,0 +1,17 @@
+import argparse
+from S_run_aifeynman import run_aifeynman
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument("--pathdir", type=str, help="Path to the directory containing the data file")
+parser.add_argument("--filename", type=str, help="Name of the file containing the data")
+parser.add_argument("--BF_try_time", type=float, default=60, help="Time limit for each brute force code call")
+parser.add_argument("--BF_ops_file_type", type=str, default="14ops.txt", help="File containing the symbols to be used in the brute force code")
+parser.add_argument("--polyfit_deg", type=int, default=4, help="Maximum degree of the polynomial tried by the polynomial fit routine")
+parser.add_argument("--NN_epochs", type=int, default=2000, help="Number of epochs for the training")
+
+opts = parser.parse_args()
+
+run_aifeynman(opts.pathdir, opts.filename, BF_try_time=opts.BF_try_time, BF_ops_file_type=opts.BF_ops_file_type, polyfit_deg=opts.polyfit_deg,
+          NN_epochs=opts.NN_epochs)
+          

Code/brute_force_oneFile.scr

@@ -0,0 +1,19 @@
+#!/bin/csh
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+set opsfile = $1
+set maxtime = $2
+set f = $3
+
+set outfile = brute_solutions.dat
+set outfile2 = brute_formulas.dat
+if -f $outfile /bin/rm $outfile
+if -f $outfile2 /bin/rm $outfile2
+
+echo Trying to solve mysteries with brute force...
+
+echo Trying to solve $f...
+echo /bin/cp -p $f mystery.dat
+/bin/cp -p $f mystery.dat
+echo $opsfile arity2templates.txt mystery.dat results.dat >args.dat
+timeout {$maxtime}s ./symbolic_regress.x  

Code/brute_force_oneFile_v2.scr

@@ -0,0 +1,21 @@
+#!/bin/csh
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+set opsfile = $1
+set maxtime = $2
+set f = $3
+
+set outfile = brute_solutions.dat
+set outfile2 = brute_constant.dat
+set outfile3 = brute_formulas.dat
+if -f $outfile /bin/rm $outfile
+if -f $outfile2 /bin/rm $outfile2
+if -f $outfile3 /bin/rm $outfile3
+
+echo Trying to solve mysteries with brute force...
+
+echo Trying to solve $f...
+echo /bin/cp -p $f mystery.dat
+/bin/cp -p $f mystery.dat
+echo $opsfile arity2templates.txt mystery.dat results.dat >args.dat
+timeout {$maxtime}s ./symbolic_regress2.x 

Code/brute_force_oneFile_v3.scr

@@ -0,0 +1,23 @@
+#!/bin/csh
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+set opsfile = $1
+set maxtime = $2
+set f = $3
+
+set outfile = brute_solutions.dat
+set outfile2 = brute_constant.dat
+set outfile3 = brute_formulas.dat
+if -f $outfile /bin/rm $outfile
+if -f $outfile2 /bin/rm $outfile2
+if -f $outfile3 /bin/rm $outfile3
+
+echo Trying to solve mysteries with brute force...
+
+echo Trying to solve $f...
+echo /bin/cp -p $f mystery.dat
+/bin/cp -p $f mystery.dat
+echo $opsfile arity2templates.txt mystery.dat results.dat >args.dat
+
+timeout {$maxtime}s ./symbolic_regress3.x
+

Code/compile.sh

@@ -0,0 +1,3 @@
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress.x symbolic_regress.f
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress2.x symbolic_regress2.f
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress3.x symbolic_regress3.f

Code/get_pareto.py

@@ -0,0 +1,267 @@
+from collections import namedtuple
+
+import matplotlib.pyplot as plt
+import numpy as np
+from sortedcontainers import SortedKeyList
+
+
+class Point(object):
+    def __init__(self, x, y, data=None, id=None):
+        self.x = x
+        self.y = y
+        self.data = data
+        self.id = id
+
+
+    def __getitem__(self, index):
+        """Indexing: get item according to index."""
+        if index == 0:
+            return self.x
+        elif index == 1:
+            return self.y
+        elif index == 2:
+            return self.data
+        elif index == 3:
+            return self.id
+        else:
+            raise Exception("Index {} is out of range!".format(index))
+
+
+    def __setitem__(self, index, value):
+        """Indexing: set item according to index."""
+        if index == 0:
+            self.x = value
+        elif index == 1:
+            self.y = value
+        elif index == 2:
+            self.data = value
+        elif index == 3:
+            raise Exception("Cannot set Id!")
+        else:
+            raise Exception("Index {} is out of range!".format(index))
+
+
+# In[2]:
+
+
+class ParetoSet(SortedKeyList):
+    """Maintained maximal set with efficient insertion. Note that we use the convention of smaller the better."""
+
+    def __init__(self):
+        super().__init__(key=lambda p: p.x)
+
+
+    def _input_check(self, p):
+        """Check that input is in the correct format.
+
+        Args:
+            p: input
+
+        Returns:
+            Point:
+
+        Raises:
+            TypeError if cannot be converted.
+        """
+
+        if isinstance(p, Point):
+            return p
+        elif isinstance(p, tuple) and len(p) == 2:
+            return Point(x=p[0], y=p[1], data=None)
+        else:
+            raise TypeError("Must be instance of Point or 2-tuple.")
+    
+    
+    def get_id_list(self):
+        id_list = []
+        for point in self:
+            id_list.append(point.id)
+        return id_list
+
+    
+    def add(self, p):
+        """Insert Point into set if minimal in first two indices.
+
+        Args:
+            p (Point): Point to insert
+
+        Returns:
+            bool: True only if point is inserted
+
+        """
+        p = self._input_check(p)
+
+        is_pareto = False
+        # check right for dominated points:
+        right = self.bisect_left(p)
+
+        while len(self) > right and self[right].y >= p.y and not (self[right].x == p.x and self[right].y == p.y):
+            self.pop(right)
+            is_pareto = True
+
+        # check left for dominating points:
+        left = self.bisect_right(p) - 1
+
+        if left == -1 or self[left][1] > p[1]:
+            is_pareto = True
+
+        # if it's the only point it's maximal
+        if len(self) == 0:
+            is_pareto = True
+
+        if is_pareto:
+            super().add(p)
+
+        return is_pareto
+
+
+    def __contains__(self, p):
+        p = self._input_check(p)
+
+        left = self.bisect_left(p)
+
+        while len(self) > left and self[left].x == p.x:
+            if self[left].y == p.y:
+                return True
+
+            left += 1
+
+        return False
+
+
+    def __add__(self, other):
+        """Merge another pareto set into self.
+
+        Args:
+            other (ParetoSet): set to merge into self
+
+        Returns:
+            ParetoSet: self
+
+        """
+
+        for item in other:
+            self.add(item)
+
+        return self
+
+
+    def distance(self, p):
+        """Given a Point, calculate the minimum Euclidean distance to pareto
+        frontier (in first two indices).
+
+        Args:
+            p (Point): point
+
+        Returns:
+            float: minimum Euclidean distance to pareto frontier
+
+        """
+        p = self._input_check(p)
+
+        point = np.array((p.x, p.y))
+        dom = self.dominant_array(p)
+
+        # distance is zero if pareto optimal
+        if dom.shape[0] == 0:
+            return 0.
+
+        # add corners of all adjacent pairs
+        candidates = np.zeros((dom.shape[0] + 1, 2))
+        for i in range(dom.shape[0] - 1):
+            candidates[i, :] = np.max(dom[[i, i+1], :], axis=0)
+
+        # add top and right bounds
+        candidates[-1, :] = (p.x, np.min(dom[:, 1]))
+        candidates[-2, :] = (np.min(dom[:, 0]), p.y)
+
+        return np.min(np.sqrt(np.sum(np.square(candidates - point), axis=1)))
+
+
+    def dominant_array(self, p):
+        """Given a Point, return the set of dominating points in the set (in
+        the first two indices).
+
+        Args:
+            p (Point): point
+
+        Returns:
+            numpy.ndarray: array of dominating points
+
+        """
+        p = self._input_check(p)
+
+        idx = self.bisect_left(p) - 1
+
+        domlist = []
+
+        while idx >= 0 and self[idx][1] < p[1]:
+            domlist.append(self[idx])
+            idx -= 1
+
+        return np.array([x[0:2] for x in domlist])
+
+
+    def to_array(self):
+        """Convert first two indices to numpy.ndarray
+
+        Args:
+            None
+
+        Returns:
+            numpy.ndarray: array of shape (len(self), 2)
+
+        """
+        A = np.zeros((len(self), 2))
+        for i, p in enumerate(self):
+            A[i, :] = p.x, p.y
+
+        return A
+
+    def get_pareto_points(self):
+        """Returns the x, y and data for each point in the pareto frontier
+        
+        """
+        pareto_points = []
+        for i, p in enumerate(self):
+            pareto_points = pareto_points + [[p.x, p.y, p.data]]
+        
+        return pareto_points
+        
+
+    def from_list(self, A):
+        """Convert iterable of Points into ParetoSet.
+
+        Args:
+            A (iterator): iterator of Points
+
+        Returns:
+            None
+
+        """
+        for a in A:
+            self.add(a)
+    
+    
+    def plot(self):
+        """Plotting the Pareto frontier."""
+        array = self.to_array()
+        plt.figure(figsize=(8, 6))
+        plt.plot(array[:, 0], array[:, 1], 'r.')
+        plt.show()
+
+
+if __name__ == "__main__":
+    PA = ParetoSet()
+    A = np.zeros((40, 2))
+    
+    for i in range(40):
+        x = np.random.rand()
+        y = np.random.rand()
+        
+        A[i, 0] = x
+        A[i, 1] = y
+        
+        PA.add(Point(x=x, y=y, data=None))
+    paretoA = PA.to_array()
+

Code/symbolic_regress.f

@@ -0,0 +1,271 @@
+	! Max Tegmark 171119, 190128-31, 190506
+	! Loads templates.csv functions.dat and mystery.dat, returns winner.
+	! scp -P2222 symbolic_regress.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress.x symbolic_regress.f |& more'
+	! SAMPLE USAGE: call symbolic_regress.x 19ops.txt arity2templates.txt mystery16.dat results.dat
+	! functions.dat contains a single line (say "0>+*-/") with the single-character symbols 
+       ! that will be used, drawn from this list: 
+	!
+	! Binary:
+	! +: add
+	! *: multiply
+	! -: subtract
+	! /: divide	(Put "D" instead of "/" in file, since f77 can't load backslash
+	!
+	! Unary:
+	!  >: increment (x -> x+1) 
+	!  <: decrement (x -> x-1)
+	!  ~: negate  	(x-> -x)
+	!  \: invert    (x->1/x) (Put "I" instead of "\" in file, since f77 can't load backslash
+	!  L: logaritm: (x-> ln(x)
+	!  E: exponentiate (x->exp(x))
+	!  S: sin:      (x->sin(x))       
+	!  C: cos:      (x->cos(x))       
+	!  A: abs:      (x->abs(x))
+	!  N: arcsin:   (x->arcsin(x))
+	!  T: arctan:   (x->arctan(x))
+	!  R: sqrt	(x->sqrt(x))
+	!
+	! nonary: 
+	!  0
+	!  1
+	!  P = pi
+	!  a, b, c, ...: input variables for function (need not be listed in functions.dat)
+
+	program symbolic_regress
+	call go
+	end
+	
+	subroutine go
+	implicit none
+	character*60 opsfile, templatefile, mysteryfile, outfile, usedfuncs
+	character*60 comline, functions, ops, formula
+	integer arities(19), nvar, nvarmax, nmax, lnblnk
+	parameter(nvarmax=20, nmax=10000000)
+	real*8 f, newloss, minloss, maxloss, rmsloss, xy(nvarmax+1,nmax), epsilon, DL, DL2, DL3
+	parameter(epsilon=0.000000000000000000001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\LESCANTR01P"/
+      	integer nn(0:2), ii(nmax), kk(nmax), radix(nmax)
+	integer ndata, i, j, n
+	integer*8 nformulas
+	logical done
+	character*60 func(0:2), template
+
+	open(2,file='args.dat',status='old',err=666)
+	read(2,*) opsfile, templatefile, mysteryfile, outfile
+	close(2)
+
+        comline = 'head -1 '//mysteryfile(1:lnblnk(mysteryfile))//' | wc > qaz.dat'
+        if (system(comline).ne.0) stop 'DEATH ERROR counting columns'
+        open(2,file='qaz.dat')
+        read(2,*) i, nvar
+        close(2)
+	nvar = nvar - 1
+	if (nvar.gt.nvarmax) stop 'DEATH ERROR: TOO MANY VARIABLES'
+	write(*,'(1a24,i8)') 'Number of variables.....',nvar
+
+	open(2,file=opsfile,status='old',err=668)
+	read(2,*) usedfuncs
+	close(2)
+	nn(0)=0
+	nn(1)=0
+	nn(2)=0
+	do i=1,lnblnk(usedfuncs) 
+	  if (usedfuncs(i:i).eq.'D') usedfuncs(i:i)='/'
+	  if (usedfuncs(i:i).eq.'I') usedfuncs(i:i)='\'
+	  j = index(functions,usedfuncs(i:i))
+	  if (j.eq.0) then
+            print *,'DEATH ERROR: Unknown function requested: ',usedfuncs(i:i)
+	    stop
+	  else
+	    nn(arities(j)) = nn(arities(j)) + 1
+	    func(arities(j))(nn(arities(j)):nn(arities(j))) = functions(j:j)
+         end if
+	end do
+	! Add nonary ops to retrieve each of the input variables:
+	do i=1,nvar
+          nn(0) = nn(0) + 1
+	  func(0)(nn(0):nn(0)) = char(96+i)
+	end do
+	write(*,'(1a24,1a22)') 'Functions used..........',usedfuncs(1:lnblnk(usedfuncs))
+	do i=0,2
+	  write(*,*) 'Arity ',i,': ',func(i)(1:nn(i))
+	end do
+
+	write(*,'(1a24)') 'Loading mystery data....'
+	call LoadMatrixTranspose(nvarmax+1,nvar+1,nmax,ndata,xy,mysteryfile)	
+	write(*,'(1a24,i8)') 'Number of examples......',ndata
+
+	print *,'Searching for best fit...'	
+	nformulas = 0
+	minloss = 1.e6
+	template = ''
+	ops='===================='
+	open(2,file=templatefile,status='old',err=670)
+	open(3,file=outfile)
+555	read(2,'(1a60)',end=665) template
+	n = lnblnk(template)
+	!print *,"template:",template(1:n),"#####"
+	do i=1,n
+	  ii(i) = ichar(template(i:i))-48
+	  radix(i) = nn(ii(i))
+	  kk(i) = 0
+	  !print *,'ASILOMAR ', i,ii(i),kk(i),radix(i)
+	end do
+	done = .false.
+	do while ((minloss.gt.epsilon).and.(.not.done))
+	  nformulas = nformulas + 1
+	  ! Analyze structure ii:
+	  do i=1,n
+	    ops(i:i) = func(ii(i))(1+kk(i):1+kk(i))
+            !print *,'TEST ',i,ii(i), func(ii(i))
+	  end do
+	  !write(*,'(1f20.12,99i3)') minloss, (ii(i),i=1,n), (kk(i),i=1,n)
+          !write(*,'(1a24)') ops(1:n)
+	  j = 1
+	  maxloss = 0.
+	  do while ((maxloss.lt.minloss).and.(j.le.ndata))
+	    newloss = abs(xy(nvar+1,j) - f(n,ii,ops,xy(1,j)))
+            !!!!!print *,'newloss: ',j,newloss,xy(nvar,j),f(n,ii,ops,xy(1,j))
+	    if (.not.((newloss.ge.0).or.(newloss.le.0))) newloss = 1.e30 ! This was a NaN :-)
+	    if (maxloss.lt.newloss) maxloss = newloss
+	    j = j + 1
+	  end do
+	  if (maxloss.lt.minloss) then ! We have a new best fit
+	    minloss = maxloss
+	    rmsloss = 0.
+	    do j=1,ndata
+	      rmsloss = rmsloss + (xy(nvar+1,j) - f(n,ii,ops,xy(1,j)))**2
+	    end do
+	    rmsloss = sqrt(rmsloss/ndata)
+	    DL  = log(nformulas*max(1.,rmsloss/epsilon))/log(2.)
+	    DL2 = log(nformulas*max(1.,rmsloss/1.e-15))/log(2.)
+	    DL3 = (log(1.*nformulas) + sqrt(1.*ndata)*log(max(1.,rmsloss/1.e-15)))/log(2.)
+	    write(*,'(1f20.12,x,1a22,1i16,4f19.4)') minloss, ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    write(3,'(1f20.12,x,1a22,1i16,4f19.4)') minloss, ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    flush(3)
+	  end if
+	  call multiloop(n,radix,kk,done) 
+	end do
+	goto 555
+665	close(3)
+	close(2)
+	print *,'All done: results in ',outfile
+	return
+666	stop 'DEATH ERROR: missing file args.dat'
+668	print *,'DEATH ERROR: missing file ',opsfile(1:lnblnk(opsfile))
+	stop
+670	print *,'DEATH ERROR: missing file ',templatefile(1:lnblnk(templatefile))
+	stop
+	end
+
+
+
+	real*8 function f(n,arities,ops,x) ! n=number of ops, x=arg vector
+	implicit none
+	integer nmax, n, i, j, arities(n), arity, lnblnk
+	character*60 ops
+	parameter(nmax=100)
+	real*8 x(nmax), y, stack(nmax)
+	character op
+	!write(*,*) 'Evaluating function with ops = ',ops(1:n)
+	!write(*,'(3f10.5,99i3)') (x(i),i=1,3), (arities(i),i=1,n)
+	j = 0 ! Number of numbers on the stack
+	do i=1,n
+	  arity = arities(i)
+	  op = ops(i:i)
+	  if (arity.eq.0) then ! This is a nonary function
+	    if (op.eq."0") then
+              y = 0.
+	    else if (op.eq."1") then
+	      y = 1.
+	    else if (op.eq."P") then
+ 	      y = 4.*atan(1.) ! pi
+           else 
+    	      y = x(ichar(op)-96)
+	    end if
+	  else if (arity.eq.1) then ! This is a unary function
+	    if (op.eq.">") then
+             y = stack(j) + 1
+	    else if (op.eq."<") then
+             y = stack(j) - 1
+	    else if (op.eq."~") then
+              y = -stack(j)
+	    else if (op.eq."\") then
+             y = 1./stack(j)
+	    else if (op.eq."L") then
+             y = log(stack(j))
+	    else if (op.eq."E") then
+             y = exp(stack(j))
+	    else if (op.eq."S") then
+             y = sin(stack(j))
+	    else if (op.eq."C") then
+             y =cos(stack(j))
+	    else if (op.eq."A") then
+             y = abs(stack(j))
+	    else if (op.eq."N") then
+             y = asin(stack(j))
+	    else if (op.eq."T") then
+             y = atan(stack(j))
+	    else
+             y = sqrt(stack(j))
+	    end if
+	  else ! This is a binary function
+	  if (op.eq."+") then
+	      y = stack(j-1)+stack(j)
+	    else if (op.eq."-") then
+	      y = stack(j-1)-stack(j)
+	    else if (op.eq."*") then
+	      y = stack(j-1)*stack(j)
+	    else
+	      y = stack(j-1)/stack(j)
+	    end if
+          end if
+	  j = j + 1 - arity
+          stack(j) = y  
+	  ! write(*,'(9f10.5)') (stack(k),k=1,j)
+	end do
+	if (j.ne.1) stop 'DEATH ERROR: STACK UNBALANCED'
+	f = stack(1)
+        !write(*,'(9f10.5)') 666.,x(1),x(2),x(3),f
+	return
+	end
+
+        subroutine multiloop(n,bases,i,done)
+        ! Handles <n> nested loops with loop variables i(1),...i(n).
+        ! Example: With n=3, bases=2, repeated calls starting with i=(000) will return
+        ! 001, 010, 011, 100, 101, 110, 111, 000 (and done=.true. the last time).
+        ! All it's doing is counting in mixed radix specified by the array <bases>.
+        implicit none
+        integer n, bases(n), i(n), k
+        logical done
+        done = .false.
+        k = 1
+555     i(k) = i(k) + 1
+        if (i(k).lt.bases(k)) return
+        i(k) = 0
+        k = k + 1
+        if (k.le.n) goto 555
+        done = .true.
+        return
+        end
+
+	subroutine LoadMatrixTranspose(nd,n,mmax,m,A,f)
+	! Reads the n x m matrix A from the file named f, stored as its transpose
+	implicit none
+	integer nd,mmax,n,m,j
+	real*8 A(nd,mmax)
+	character*60 f
+	open(2,file=f,status='old')
+	m = 0
+555	m = m + 1
+	if (m.gt.mmax) stop 'DEATH ERROR: m>mmax in LoadVectorTranspose'
+	read(2,*,end=666) (A(j,m),j=1,n)
+	goto 555
+666	close(2)
+	m = m - 1
+	print *,m,' rows read from file ',f
+	return
+	end
+

Code/symbolic_regress2.f

@@ -0,0 +1,292 @@
+	! Max Tegmark 171119, 190128-31, 190218
+	! Same as symbolic_regress2.f except that it fits for the symbolic formula times an arbitrary constant.
+	! Loads templates.csv functions.dat and mystery.dat, returns winner.
+	! scp -P2222 symbolic_regress2.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress2.x symbolic_regress2.f |& more'
+	! SAMPLE USAGE: call symbolic_regress2.x 4ops.txt arity2templates.txt mysteryB3.dat results.dat
+	! functions.dat contains a single line (say "0>+*-/") with the single-character symbols 
+        ! that will be used, drawn from this list: 
+	!
+	! Binary:
+	! +: add
+	! *: multiply
+	! -: subtract
+	! /: divide	(Put "D" instead of "/" in file, since f77 can't load backslash
+	!
+	! Unary:
+	!  >: increment (x -> x+1) 
+	!  <: decrement (x -> x-1)
+	!  ~: negate  	(x-> -x)
+	!  \: invert    (x->1/x) (Put "I" instead of "\" in file, since f77 can't load backslash
+	!  L: logaritm: (x-> ln(x)
+	!  E: exponentiate (x->exp(x))
+	!  S: sin:      (x->sin(x))       
+	!  C: cos:      (x->cos(x))       
+	!  A: abs:      (x->abs(x))
+	!  N: arcsin:   (x->arcsin(x))
+	!  T: arctan:   (x->arctan(x))
+	!  R: sqrt	(x->sqrt(x))
+	!
+	! nonary: 
+	!  0
+	!  1
+	!  a, b, c, ...: input variables for function (need not be listed in functions.dat)
+
+	program symbolic_regress
+	call go
+	end
+	
+	subroutine go
+	implicit none
+	character*60 opsfile, templatefile, mysteryfile, outfile, usedfuncs
+	character*60 comline, functions, ops, formula
+	integer arities(19), nvar, nvarmax, nmax, lnblnk
+	parameter(nvarmax=20, nmax=10000000)
+	real*8 f, newloss, minloss, maxloss, rmsloss, xy(nvarmax+1,nmax), epsilon
+	real*8 ymax, prefactor, DL, DL2, DL3, limit
+	parameter(epsilon=0.000000000000000000001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\LESCANTR01P"/
+      	integer nn(0:2), ii(nmax), kk(nmax), radix(nmax)
+	integer ndata, i, j, n, jmax
+	integer*8 nformulas
+	logical done
+	character*60 func(0:2), template
+
+	open(2,file='args.dat',status='old',err=666)
+	read(2,*) opsfile, templatefile, mysteryfile, outfile
+	close(2)
+
+        comline = 'head -1 '//mysteryfile(1:lnblnk(mysteryfile))//' | wc > qaz.dat'
+        if (system(comline).ne.0) stop 'DEATH ERROR counting columns'
+        open(2,file='qaz.dat')
+        read(2,*) i, nvar
+        close(2)
+	nvar = nvar - 1
+	if (nvar.gt.nvarmax) stop 'DEATH ERROR: TOO MANY VARIABLES'
+	write(*,'(1a24,i8)') 'Number of variables.....',nvar
+
+	open(2,file=opsfile,status='old',err=668)
+	read(2,*) usedfuncs
+	close(2)
+	nn(0)=0
+	nn(1)=0
+	nn(2)=0
+	do i=1,lnblnk(usedfuncs) 
+	  if (usedfuncs(i:i).eq.'D') usedfuncs(i:i)='/'
+	  if (usedfuncs(i:i).eq.'I') usedfuncs(i:i)='\'
+	  j = index(functions,usedfuncs(i:i))
+	  if (j.eq.0) then
+            print *,'DEATH ERROR: Unknown function requested: ',usedfuncs(i:i)
+	    stop
+	  else
+	    nn(arities(j)) = nn(arities(j)) + 1
+	    func(arities(j))(nn(arities(j)):nn(arities(j))) = functions(j:j)
+         end if
+	end do
+	! Add nonary ops to retrieve each of the input variables:
+	do i=1,nvar
+          nn(0) = nn(0) + 1
+	  func(0)(nn(0):nn(0)) = char(96+i)
+	end do
+	write(*,'(1a24,1a22)') 'Functions used..........',usedfuncs(1:lnblnk(usedfuncs))
+	do i=0,2
+	  write(*,*) 'Arity ',i,': ',func(i)(1:nn(i))
+	end do
+
+	write(*,'(1a24)') 'Loading mystery data....'
+	call LoadMatrixTranspose(nvarmax+1,nvar+1,nmax,ndata,xy,mysteryfile)	
+	write(*,'(1a24,i8)') 'Number of examples......',ndata
+	! Find max(abs(y)) to use for normalization estimation (crucial to avoid data point where y~0):
+	jmax=1
+	ymax = abs(xy(1,nvar+1))
+	do j=2,ndata
+	  if (ymax < abs(xy(nvar+1,j))) then 
+             ymax = abs(xy(nvar+1,j))
+	     jmax = j
+	  end if
+	end do
+	print *,'Mystery data has largest magnitude ',ymax,' at j=',jmax
+	print *,'Searching for best fit...'	
+	nformulas = 0
+	minloss = 1.e6
+	template = ''
+	ops='===================='
+	open(2,file=templatefile,status='old',err=670)
+	open(3,file=outfile)
+555	read(2,'(1a60)',end=665) template
+	n = lnblnk(template)
+	!print *,"template:",template(1:n),"#####"
+	do i=1,n
+	  ii(i) = ichar(template(i:i))-48
+	  radix(i) = nn(ii(i))
+	  kk(i) = 0
+	end do
+	done = .false.
+	do while ((minloss.gt.epsilon).and.(.not.done))
+	  nformulas = nformulas + 1
+	  ! Analyze structure ii:
+	  do i=1,n
+	    ops(i:i) = func(ii(i))(1+kk(i):1+kk(i))
+            !print *,'TEST ',i,ii(i), func(ii(i))
+	  end do
+	  !write(*,'(1f20.12,99i3)') minloss, (ii(i),i=1,n), (kk(i),i=1,n)
+          !write(*,'(1a24)') ops(1:n)
+
+	  prefactor = xy(nvar+1,jmax)/f(n,ii,ops,xy(1,jmax))
+	  j = 1
+	  maxloss = 0.
+	  do while ((maxloss.lt.minloss).and.(j.le.ndata))
+	    newloss = abs(xy(nvar+1,j) - prefactor*f(n,ii,ops,xy(1,j)))
+            !!!!!print *,'newloss: ',j,newloss,xy(nvar,j),f(n,ii,ops,xy(1,j))
+	    if (.not.((newloss.ge.0).or.(newloss.le.0))) newloss = 1.e30 ! This was a NaN :-)
+	    if (maxloss.lt.newloss) maxloss = newloss
+	    j = j + 1
+	  end do
+	  if (maxloss.lt.minloss) then ! We have a new best fit
+	    minloss = maxloss
+	    rmsloss = 0.
+	    do j=1,ndata
+	      rmsloss = rmsloss + (xy(nvar+1,j) - prefactor*f(n,ii,ops,xy(1,j)))**2
+	    end do
+	    rmsloss = sqrt(rmsloss/ndata)
+	    DL  = log(nformulas*max(1.,minloss/epsilon))/log(2.)
+	    DL2 = log(nformulas*max(1.,minloss/1.e-15))/log(2.)
+	    DL3 = (log(1.*nformulas) + sqrt(1.*ndata)*log(max(1.,rmsloss/1.e-15)))/log(2.)
+	    write(*,'(2f20.12,x,1a22,1i16,4f19.4)') limit(minloss), limit(prefactor), ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    write(3,'(2f20.12,x,1a22,1i16,4f19.4)') limit(minloss), limit(prefactor), ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    flush(3)
+	  end if
+	  call multiloop(n,radix,kk,done) 
+	end do
+	goto 555
+665	close(3)
+	close(2)
+	print *,'All done: results in ',outfile
+	return
+666	stop 'DEATH ERROR: missing file args.dat'
+668	print *,'DEATH ERROR: missing file ',opsfile(1:lnblnk(opsfile))
+	stop
+670	print *,'DEATH ERROR: missing file ',templatefile(1:lnblnk(templatefile))
+	stop
+	end
+
+	real*8 function limit(x)
+	implicit none
+	real*8 x, xmax
+	parameter(xmax=666.)
+	if (abs(x).lt.xmax) then
+ 	  limit = x
+	else
+	  limit = sign(xmax,x)
+	end if
+	return
+	end
+
+	real*8 function f(n,arities,ops,x) ! n=number of ops, x=arg vector
+	implicit none
+	integer nmax, n, i, j, arities(n), arity, lnblnk
+	character*60 ops
+	parameter(nmax=100)
+	real*8 x(nmax), y, stack(nmax)
+	character op
+	!write(*,*) 'Evaluating function with ops = ',ops(1:n)
+	!write(*,'(3f10.5,99i3)') (x(i),i=1,3), (arities(i),i=1,n)
+	j = 0 ! Number of numbers on the stack
+	do i=1,n
+	  arity = arities(i)
+	  op = ops(i:i)
+	  if (arity.eq.0) then ! This is a nonary function
+	    if (op.eq."0") then
+              y = 0.
+	    else if (op.eq."1") then
+	      y = 1.
+	    else if (op.eq."P") then
+ 	      y = 4.*atan(1.) ! pi
+           else 
+    	      y = x(ichar(op)-96)
+	    end if
+	  else if (arity.eq.1) then ! This is a unary function
+	    if (op.eq.">") then
+             y = stack(j) + 1
+	    else if (op.eq."<") then
+             y = stack(j) - 1
+	    else if (op.eq."~") then
+              y = -stack(j)
+	    else if (op.eq."\") then
+             y = 1./stack(j)
+	    else if (op.eq."L") then
+             y = log(stack(j))
+	    else if (op.eq."E") then
+             y = exp(stack(j))
+	    else if (op.eq."S") then
+             y = sin(stack(j))
+	    else if (op.eq."C") then
+             y =cos(stack(j))
+	    else if (op.eq."A") then
+             y = abs(stack(j))
+	    else if (op.eq."N") then
+             y = asin(stack(j))
+	    else if (op.eq."T") then
+             y = atan(stack(j))
+	    else
+             y = sqrt(stack(j))
+	    end if
+	  else ! This is a binary function
+	  if (op.eq."+") then
+	      y = stack(j-1)+stack(j)
+	    else if (op.eq."-") then
+	      y = stack(j-1)-stack(j)
+	    else if (op.eq."*") then
+	      y = stack(j-1)*stack(j)
+	    else
+	      y = stack(j-1)/stack(j)
+	    end if
+          end if
+	  j = j + 1 - arity
+          stack(j) = y  
+	  ! write(*,'(9f10.5)') (stack(k),k=1,j)
+	end do
+	if (j.ne.1) stop 'DEATH ERROR: STACK UNBALANCED'
+	f = stack(1)
+        !write(*,'(9f10.5)') 666.,x(1),x(2),x(3),f
+	return
+	end
+
+        subroutine multiloop(n,bases,i,done)
+        ! Handles <n> nested loops with loop variables i(1),...i(n).
+        ! Example: With n=3, bases=2, repeated calls starting with i=(000) will return
+        ! 001, 010, 011, 100, 101, 110, 111, 000 (and done=.true. the last time).
+        ! All it's doing is counting in mixed radix specified by the array <bases>.
+        implicit none
+        integer n, bases(n), i(n), k
+        logical done
+        done = .false.
+        k = 1
+555     i(k) = i(k) + 1
+        if (i(k).lt.bases(k)) return
+        i(k) = 0
+        k = k + 1
+        if (k.le.n) goto 555
+        done = .true.
+        return
+        end
+
+	subroutine LoadMatrixTranspose(nd,n,mmax,m,A,f)
+	! Reads the n x m matrix A from the file named f, stored as its transpose
+	implicit none
+	integer nd,mmax,n,m,j
+	real*8 A(nd,mmax)
+	character*60 f
+	open(2,file=f,status='old')
+	m = 0
+555	m = m + 1
+	if (m.gt.mmax) stop 'DEATH ERROR: m>mmax in LoadVectorTranspose'
+	read(2,*,end=666) (A(j,m),j=1,n)
+	goto 555
+666	close(2)
+	m = m - 1
+	print *,m,' rows read from file ',f
+	return
+	end
+

Code/symbolic_regress3.f

@@ -0,0 +1,293 @@
+	! Max Tegmark 171119, 190128-31, 190218, 25
+	! Same as symbolic_regress3.f except that it fits for the symbolic formula plus an arbitrary constant.
+	! Loads templates.csv functions.dat and mystery.dat, returns winner.
+	! scp -P2222 symbolic_regress3.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress3.x symbolic_regress3.f |& more'
+	! SAMPLE USAGE: call symbolic_regress3.x 6ops.txt arity2templates.txt mystery012.dat results.dat
+	! functions.dat contains a single line (say "0>+*-/") with the single-character symbols 
+        ! that will be used, drawn from this list: 
+	!
+	! Binary:
+	! +: add
+	! *: multiply
+	! -: subtract
+	! /: divide	(Put "D" instead of "/" in file, since f77 can't load backslash
+	!
+	! Unary:
+	!  >: increment (x -> x+1) 
+	!  <: decrement (x -> x-1)
+	!  ~: negate  	(x-> -x)
+	!  \: invert    (x->1/x) (Put "I" instead of "\" in file, since f77 can't load backslash
+	!  L: logaritm: (x-> ln(x)
+	!  E: exponentiate (x->exp(x))
+	!  S: sin:      (x->sin(x))       
+	!  C: cos:      (x->cos(x))       
+	!  A: abs:      (x->abs(x))
+	!  N: arcsin:   (x->arcsin(x))
+	!  T: arctan:   (x->arctan(x))
+	!  R: sqrt	(x->sqrt(x))
+	!
+	! nonary: 
+	!  0
+	!  1
+	!  a, b, c, ...: input variables for function (need not be listed in functions.dat)
+
+	program symbolic_regress
+	call go
+	end
+	
+	subroutine go
+	implicit none
+	character*60 opsfile, templatefile, mysteryfile, outfile, usedfuncs
+	character*60 comline, functions, ops, formula
+	integer arities(19), nvar, nvarmax, nmax, lnblnk
+	parameter(nvarmax=20, nmax=10000000)
+	real*8 f, newloss, minloss, maxloss, rmsloss, xy(nvarmax+1,nmax), epsilon
+	real*8 ymin, prefactor, DL, DL2, DL3, limit
+	parameter(epsilon=0.000000000000000000001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\LESCANTR01P"/
+      	integer nn(0:2), ii(nmax), kk(nmax), radix(nmax)
+	integer ndata, i, j, n, jmin
+	integer*8 nformulas
+	logical done
+	character*60 func(0:2), template
+
+	open(2,file='args.dat',status='old',err=666)
+	read(2,*) opsfile, templatefile, mysteryfile, outfile
+	close(2)
+
+        comline = 'head -1 '//mysteryfile(1:lnblnk(mysteryfile))//' | wc > qaz.dat'
+        if (system(comline).ne.0) stop 'DEATH ERROR counting columns'
+        open(2,file='qaz.dat')
+        read(2,*) i, nvar
+        close(2)
+	nvar = nvar - 1
+	if (nvar.gt.nvarmax) stop 'DEATH ERROR: TOO MANY VARIABLES'
+	write(*,'(1a24,i8)') 'Number of variables.....',nvar
+
+	open(2,file=opsfile,status='old',err=668)
+	read(2,*) usedfuncs
+	close(2)
+	nn(0)=0
+	nn(1)=0
+	nn(2)=0
+	do i=1,lnblnk(usedfuncs) 
+	  if (usedfuncs(i:i).eq.'D') usedfuncs(i:i)='/'
+	  if (usedfuncs(i:i).eq.'I') usedfuncs(i:i)='\'
+	  j = index(functions,usedfuncs(i:i))
+	  if (j.eq.0) then
+            print *,'DEATH ERROR: Unknown function requested: ',usedfuncs(i:i)
+	    stop
+	  else
+	    nn(arities(j)) = nn(arities(j)) + 1
+	    func(arities(j))(nn(arities(j)):nn(arities(j))) = functions(j:j)
+         end if
+	end do
+	! Add nonary ops to retrieve each of the input variables:
+	do i=1,nvar
+          nn(0) = nn(0) + 1
+	  func(0)(nn(0):nn(0)) = char(96+i)
+	end do
+	write(*,'(1a24,1a22)') 'Functions used..........',usedfuncs(1:lnblnk(usedfuncs))
+	do i=0,2
+	  write(*,*) 'Arity ',i,': ',func(i)(1:nn(i))
+	end do
+
+	write(*,'(1a24)') 'Loading mystery data....'
+	call LoadMatrixTranspose(nvarmax+1,nvar+1,nmax,ndata,xy,mysteryfile)	
+	write(*,'(1a24,i8)') 'Number of examples......',ndata
+	! Find min(abs(y)) to use for offset estimation:
+	jmin=1
+	ymin = abs(xy(1,nvar+1))
+	do j=2,ndata
+	  if (ymin > abs(xy(nvar+1,j))) then 
+             ymin = abs(xy(nvar+1,j))
+	     jmin = j
+	  end if
+	end do
+	print *,'Mystery data has largest magnitude ',ymin,' at j=',jmin
+	print *,'Searching for best fit...'	
+	nformulas = 0
+	minloss = 1.e6
+	template = ''
+	ops='===================='
+	open(2,file=templatefile,status='old',err=670)
+	open(3,file=outfile)
+555	read(2,'(1a60)',end=665) template
+	n = lnblnk(template)
+	!print *,"template:",template(1:n),"#####"
+	do i=1,n
+	  ii(i) = ichar(template(i:i))-48
+	  radix(i) = nn(ii(i))
+	  kk(i) = 0
+	end do
+	done = .false.
+	do while ((minloss.gt.epsilon).and.(.not.done))
+	  nformulas = nformulas + 1
+	  ! Analyze structure ii:
+	  do i=1,n
+	    ops(i:i) = func(ii(i))(1+kk(i):1+kk(i))
+            !print *,'TEST ',i,ii(i), func(ii(i))
+	  end do
+	  !write(*,'(1f20.12,99i3)') minloss, (ii(i),i=1,n), (kk(i),i=1,n)
+          !write(*,'(1a24)') ops(1:n)
+
+	  prefactor = xy(nvar+1,jmin)-f(n,ii,ops,xy(1,jmin))
+	  j = 1
+	  maxloss = 0.
+	  do while ((maxloss.lt.minloss).and.(j.le.ndata))
+	    newloss = abs(xy(nvar+1,j) - (prefactor+f(n,ii,ops,xy(1,j))))
+            !!!!!print *,'newloss: ',j,newloss,xy(nvar,j),f(n,ii,ops,xy(1,j))
+	    if (.not.((newloss.ge.0).or.(newloss.le.0))) newloss = 1.e30 ! This was a NaN :-)
+	    if (maxloss.lt.newloss) maxloss = newloss
+	    j = j + 1
+	  end do
+	  if (maxloss.lt.minloss) then ! We have a new best fit
+	    minloss = maxloss
+	    rmsloss = 0.
+	    do j=1,ndata
+	      newloss = abs(xy(nvar+1,j) - (prefactor+f(n,ii,ops,xy(1,j))))
+	      rmsloss = rmsloss + newloss**2
+ 	    end do
+	    rmsloss = sqrt(rmsloss/ndata)
+	    DL  = log(nformulas*max(1.,minloss/epsilon))/log(2.)
+	    DL2 = log(nformulas*max(1.,minloss/1.e-15))/log(2.)
+	    DL3 = (log(1.*nformulas) + sqrt(1.*ndata)*log(max(1.,rmsloss/1.e-15)))/log(2.)
+	    write(*,'(2f20.12,x,1a22,1i16,4f19.4)') limit(minloss), limit(prefactor), ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    write(3,'(2f20.12,x,1a22,1i16,4f19.4)') limit(minloss), limit(prefactor), ops(1:n), nformulas, rmsloss, DL, DL2, DL3
+	    flush(3)
+	  end if
+	  call multiloop(n,radix,kk,done) 
+	end do
+	goto 555
+665	close(3)
+	close(2)
+	print *,'All done: results in ',outfile
+	return
+666	stop 'DEATH ERROR: missing file args.dat'
+668	print *,'DEATH ERROR: missing file ',opsfile(1:lnblnk(opsfile))
+	stop
+670	print *,'DEATH ERROR: missing file ',templatefile(1:lnblnk(templatefile))
+	stop
+	end
+
+	real*8 function limit(x)
+	implicit none
+	real*8 x, xmax
+	parameter(xmax=666.)
+	if (abs(x).lt.xmax) then
+ 	  limit = x
+	else
+	  limit = sign(xmax,x)
+	end if
+	return
+	end
+
+	real*8 function f(n,arities,ops,x) ! n=number of ops, x=arg vector
+	implicit none
+	integer nmax, n, i, j, arities(n), arity, lnblnk
+	character*60 ops
+	parameter(nmax=100)
+	real*8 x(nmax), y, stack(nmax)
+	character op
+	!write(*,*) 'Evaluating function with ops = ',ops(1:n)
+	!write(*,'(3f10.5,99i3)') (x(i),i=1,3), (arities(i),i=1,n)
+	j = 0 ! Number of numbers on the stack
+	do i=1,n
+	  arity = arities(i)
+	  op = ops(i:i)
+	  if (arity.eq.0) then ! This is a nonary function
+	    if (op.eq."0") then
+              y = 0.
+	    else if (op.eq."1") then
+	      y = 1.
+	    else if (op.eq."P") then
+ 	      y = 4.*atan(1.) ! pi
+           else 
+    	      y = x(ichar(op)-96)
+	    end if
+	  else if (arity.eq.1) then ! This is a unary function
+	    if (op.eq.">") then
+             y = stack(j) + 1
+	    else if (op.eq."<") then
+             y = stack(j) - 1
+	    else if (op.eq."~") then
+              y = -stack(j)
+	    else if (op.eq."\") then
+             y = 1./stack(j)
+	    else if (op.eq."L") then
+             y = log(stack(j))
+	    else if (op.eq."E") then
+             y = exp(stack(j))
+	    else if (op.eq."S") then
+             y = sin(stack(j))
+	    else if (op.eq."C") then
+             y =cos(stack(j))
+	    else if (op.eq."A") then
+             y = abs(stack(j))
+	    else if (op.eq."N") then
+             y = asin(stack(j))
+	    else if (op.eq."T") then
+             y = atan(stack(j))
+	    else
+             y = sqrt(stack(j))
+	    end if
+	  else ! This is a binary function
+	  if (op.eq."+") then
+	      y = stack(j-1)+stack(j)
+	    else if (op.eq."-") then
+	      y = stack(j-1)-stack(j)
+	    else if (op.eq."*") then
+	      y = stack(j-1)*stack(j)
+	    else
+	      y = stack(j-1)/stack(j)
+	    end if
+          end if
+	  j = j + 1 - arity
+          stack(j) = y  
+	  ! write(*,'(9f10.5)') (stack(k),k=1,j)
+	end do
+	if (j.ne.1) stop 'DEATH ERROR: STACK UNBALANCED'
+	f = stack(1)
+        !write(*,'(9f10.5)') 666.,x(1),x(2),x(3),f
+	return
+	end
+
+        subroutine multiloop(n,bases,i,done)
+        ! Handles <n> nested loops with loop variables i(1),...i(n).
+        ! Example: With n=3, bases=2, repeated calls starting with i=(000) will return
+        ! 001, 010, 011, 100, 101, 110, 111, 000 (and done=.true. the last time).
+        ! All it's doing is counting in mixed radix specified by the array <bases>.
+        implicit none
+        integer n, bases(n), i(n), k
+        logical done
+        done = .false.
+        k = 1
+555     i(k) = i(k) + 1
+        if (i(k).lt.bases(k)) return
+        i(k) = 0
+        k = k + 1
+        if (k.le.n) goto 555
+        done = .true.
+        return
+        end
+
+	subroutine LoadMatrixTranspose(nd,n,mmax,m,A,f)
+	! Reads the n x m matrix A from the file named f, stored as its transpose
+	implicit none
+	integer nd,mmax,n,m,j
+	real*8 A(nd,mmax)
+	character*60 f
+	open(2,file=f,status='old')
+	m = 0
+555	m = m + 1
+	if (m.gt.mmax) stop 'DEATH ERROR: m>mmax in LoadVectorTranspose'
+	read(2,*,end=666) (A(j,m),j=1,n)
+	goto 555
+666	close(2)
+	m = m - 1
+	print *,m,' rows read from file ',f
+	return
+	end
+