bootstrap

prior-art/Code/10ops.txt

@@ -0,0 +1 @@
+0~*DJORPEL

prior-art/Code/14ops.txt

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

prior-art/Code/19ops.txt

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

prior-art/Code/7ops.txt

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

prior-art/Code/RPN_to_eq.py

@@ -0,0 +1,75 @@
+# 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","O","J"]
+    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)
+            if i=="O":
+                a="(2*("+a+"))"
+                stack = np.append(stack,a)
+            if i=="J":
+                a="(2*("+a+")+1)"
+                stack = np.append(stack,a)
+    return(stack[0])

prior-art/Code/RPN_to_pytorch.py

@@ -0,0 +1,140 @@
+# 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
+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_snapped import get_number_DL_snapped
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+
+# parameters: path to data, RPN expression (obtained from bf)
+def RPN_to_pytorch(data, 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
+
+    # 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_()
+        if torch.isnan(loss):
+            break
+                
+    for nan_i in range(len(trainable_parameters)):
+        if torch.isnan(trainable_parameters[nan_i])==True or abs(trainable_parameters[nan_i])>1e7:
+            return 1000000, 10000000, "1"
+
+    ii = -1
+    for parm in unsnapped_param_dict:
+        if ii == -1:
+            ii = ii + 1
+        else:
+            eq = eq.subs(parm, trainable_parameters[ii])
+            ii = ii + 1
+
+    complexity = 0
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    numbers_expr = [subexpression for subexpression in preorder_traversal(eq) if is_atomic_number(subexpression)]
+    complexity = 0
+    for j in numbers_expr:
+        try:
+            complexity = complexity + get_number_DL_snapped(float(j))
+        except:
+            complexity = complexity + 1000000
+    n_variables = len(eq.free_symbols)
+    n_operations = len(count_ops(eq,visual=True).free_symbols)
+    if n_operations!=0 or n_variables!=0:
+        complexity = complexity + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+
+    error = get_symbolic_expr_error(data,str(eq))
+    return error, complexity, eq
+
+
+
+

prior-art/Code/S_NN_eval.py

@@ -0,0 +1,122 @@
+from __future__ import print_function
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+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)
+
+
+
+

prior-art/Code/S_NN_train.py

@@ -0,0 +1,161 @@
+from __future__ import print_function
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+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=1000, lrs=1e-2, N_red_lr=4, pretrained_path=""):
+    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//N_red_lr
+        epochs = int(epochs)
+
+        if n_variables==0 or n_variables==1:
+            print("Solved!")#, variables[0])
+            return 0
+
+        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)
+
+        if pretrained_path!="":
+            model_feynman.load_state_dict(torch.load(pretrained_path))
+
+        check_es_loss = 10000
+
+        for i_i in range(N_red_lr):
+            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)
+                    loss.backward()
+                    optimizer_feynman.step()
+
+                # Early stopping
+                if epoch%20==0 and epoch>0:
+                    if check_es_loss < loss:
+                        break
+                    else:
+                        torch.save(model_feynman.state_dict(), "results/NN_trained_models/models/" + filename + ".h5")
+                        check_es_loss = loss
+                if epoch==0:
+                    if check_es_loss < loss:
+                        torch.save(model_feynman.state_dict(), "results/NN_trained_models/models/" + filename + ".h5")
+                        check_es_loss = loss
+        
+            print(loss)
+            lrs = lrs/10
+
+        return 1
+
+    except NameError:
+        print("Error in file: %s" %filename)
+        raise
+
+
+
+

prior-art/Code/S_add_bf_on_numbers_on_pareto.py

@@ -0,0 +1,124 @@
+# 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
+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 S_brute_force_number import brute_force_number
+
+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
+from RPN_to_eq import RPN_to_eq
+
+from S_get_number_DL_snapped import get_number_DL_snapped
+
+# parameters: path to data, math (not RPN) expression 
+def add_bf_on_numbers_on_pareto(pathdir, filename, PA, math_expr):
+    input_data = np.loadtxt(pathdir+filename)
+    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 bf on one parameter at a time
+    bf_on_numbers_expr = []
+    for w in range(len(eq_numbers)):
+        try:
+            param_dict = {}
+            unsnapped_param_dict = {'p':1}
+            eq_ = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+            eq = eq_
+            
+            np.savetxt(pathdir+"number_for_bf_%s.txt" %w, [eq_numbers[w]])
+            brute_force_number(pathdir,"number_for_bf_%s.txt" %w)
+            # Load the predictions made by the bf code
+            bf_numbers = np.loadtxt("results.dat",usecols=(1,),dtype="str")
+            new_numbers = copy.deepcopy(eq_numbers)
+            
+            # replace the number under consideration by all the proposed bf numbers
+            for kk in range(len(bf_numbers)):
+                eq = eq_
+                new_numbers[w] = parse_expr(RPN_to_eq(bf_numbers[kk]))
+                
+                jj = 0
+                for parm in unsnapped_param_dict:
+                    if parm!="p":
+                        eq = eq.subs(parm, new_numbers[jj])
+                        jj = jj + 1
+
+                bf_on_numbers_expr = bf_on_numbers_expr + [eq]
+        except:
+            continue
+                
+    for i in range(len(bf_on_numbers_expr)):
+        try:
+            # Calculate the error of the new, snapped expression
+            snapped_error = get_symbolic_expr_error(input_data,str(bf_on_numbers_expr[i]))
+            # Calculate the complexity of the new, snapped expression
+            expr = simplify(powsimp(bf_on_numbers_expr[i]))
+            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)]
+
+            snapped_complexity = 0
+            for j in numbers_expr:
+                snapped_complexity = snapped_complexity + get_number_DL_snapped(float(j))
+            # 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:
+                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)
+        

prior-art/Code/S_add_snap_expr_on_pareto.py

@@ -0,0 +1,203 @@
+# 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
+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, S,Float
+import re
+
+from S_get_number_DL_snapped import get_number_DL_snapped
+
+def intify(expr):
+    floats = S(expr).atoms(Float)
+    ints = [i for i in floats if int(i) == i]
+    return expr.xreplace(dict(zip(ints, [int(i) for i in ints])))
+
+# parameters: path to data, math (not RPN) expression 
+def add_snap_expr_on_pareto(pathdir, filename, math_expr, PA, DR_file=""):
+    input_data = np.loadtxt(pathdir+filename)
+    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, "pp", 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 zero snap one parameter at a time
+#    zero_snapped_expr = []
+#    for w in range(len(eq_numbers)):
+#        try:
+#            param_dict = {}
+#            unsnapped_param_dict = {'pp':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!="pp":
+#                    eq = eq.subs(parm, eq_numbers[jj])
+#                    jj = jj + 1
+#            zero_snapped_expr = zero_snapped_expr + [eq]
+#        except:
+#            continue
+
+
+    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)):
+        try:
+            param_dict = {}
+            unsnapped_param_dict = {'pp':1}
+            eq = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+            del unsnapped_param_dict["pp"]
+            temp_unsnapped_param_dict = copy.deepcopy(unsnapped_param_dict)
+            new_numbers = integerSnap(eq_numbers,w+1)
+            new_numbers = {"pp"+str(k): v for k, v in new_numbers.items()}
+            temp_unsnapped_param_dict.update(new_numbers)
+            #for kk in range(len(new_numbers)):                                                                                                                                   
+            #    eq_numbers[new_numbers[kk][0]] = new_numbers[kk][1]                                                                                                              
+            new_eq = re.sub(r"(pp\d*)",r"{\1}",str(eq))
+            new_eq = new_eq.format_map(temp_unsnapped_param_dict)
+            integer_snapped_expr = integer_snapped_expr + [parse_expr(new_eq)]
+        except:
+            continue
+
+
+
+    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)):
+        try:
+            param_dict = {}
+            unsnapped_param_dict = {'pp':1}
+            eq = unsnap_recur(expr,param_dict,unsnapped_param_dict)
+            del unsnapped_param_dict["pp"]
+            temp_unsnapped_param_dict = copy.deepcopy(unsnapped_param_dict)
+            new_numbers = rationalSnap(eq_numbers,w+1)
+            new_numbers = {"pp"+str(k): v for k, v in new_numbers.items()}
+            temp_unsnapped_param_dict.update(new_numbers)
+            #for kk in range(len(new_numbers)):                                                                                                                                   
+            #    eq_numbers_snap[new_numbers[kk][0]] = new_numbers[kk][1][1:3]                                                                                                    
+            new_eq = re.sub(r"(pp\d*)",r"{\1}",str(eq))
+            new_eq = new_eq.format_map(temp_unsnapped_param_dict)
+            rational_snapped_expr = rational_snapped_expr + [parse_expr(new_eq)]
+        except:
+            continue
+
+    snapped_expr = np.append(integer_snapped_expr,rational_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(input_data,str(snapped_expr[i]))
+            # Calculate the complexity of the new, snapped expression
+            #expr = simplify(powsimp(snapped_expr[i]))
+            expr = snapped_expr[i]
+            for s in (expr.free_symbols):
+                s = symbols(str(s), real = True)
+            expr =  parse_expr(str(snapped_expr[i]),locals())
+            expr = intify(expr)
+            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 da 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()))
+                expr =  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)
+        
+        
+        
+            

prior-art/Code/S_add_sym_on_pareto.py

@@ -0,0 +1,26 @@
+# Combines 2 pareto fromtier obtained from the separability test into a new one.
+
+from get_pareto import Point, ParetoSet
+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
+
+
+
+

prior-art/Code/S_brute_force.py

@@ -0,0 +1,42 @@
+# 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 pathlib import Path
+
+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=="*":
+        # 'check=False' because it will return exit status 124 when it time out
+        subprocess.run([Path("./brute_force_oneFile_v2.scr").resolve(),
+                        file_type, "%s" %try_time,
+                        Path(pathdir+filename).resolve()],
+                       shell=False, check=False)
+        #subprocess.call(["./brute_force_oneFile_mdl_v3.scr", file_type, "%s" %try_time, pathdir+filename, "10", "0"])
+    if sep_type=="+":
+        # 'check=False' because it will return exit status 124 when it time out
+        subprocess.run([Path("./brute_force_oneFile_v3.scr").resolve(),
+                        file_type, "%s" %try_time,
+                        Path(pathdir+filename).resolve()],
+                       shell=False, check=False) 
+        #subprocess.call(["./brute_force_oneFile_mdl_v2.scr", file_type, "%s" %try_time, pathdir+filename, "10", "0"])
+    return 1
+

prior-art/Code/S_brute_force_number.py

@@ -0,0 +1,27 @@
+# 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
+
+def brute_force_number(pathdir,filename):
+    try_time = 2
+    file_type = "10ops.txt"
+
+    try:
+        os.remove("results.dat")
+    except:
+        pass
+    
+    subprocess.call(["./brute_force_oneFile_v1.scr", file_type, "%s" %try_time, pathdir+filename])
+
+    return 1
+

prior-art/Code/S_change_output.py

@@ -0,0 +1,182 @@
+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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.arccos(data[:,-1]) 
+        np.savetxt(pathdir_write_to+filename,data)     
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.arcsin(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.arctan(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.cos(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.exp(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = 1/data[:,-1]
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.log(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.sin(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.sqrt(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = data[:,-1]**2
+        np.savetxt(pathdir_write_to+filename,data)
+        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=3):
+    try:
+        os.mkdir(pathdir_write_to)
+    except:
+        pass
+    data = np.loadtxt(pathdir+filename)
+    try:
+        data[:,-1] = np.tan(data[:,-1])
+        np.savetxt(pathdir_write_to+filename,data)
+        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
+
+
+

prior-art/Code/S_combine_pareto.py

@@ -0,0 +1,36 @@
+# Combines 2 pareto fromtier obtained from the separability test into a new one.
+
+from get_pareto import Point, ParetoSet
+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(input_data,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(input_data,new_eq),data=new_eq))
+            except:
+                continue
+    return PA
+
+
+
+

prior-art/Code/S_final_gd.py

@@ -0,0 +1,146 @@
+# 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
+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_snapped import get_number_DL_snapped
+from S_get_symbolic_expr_error import get_symbolic_expr_error
+
+# parameters: path to data, RPN expression (obtained from bf)
+def final_gd(data, math_expr, lr = 1e-2, N_epochs = 5000):
+    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
+
+    # 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_()
+        if torch.isnan(loss):
+            break
+                
+    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/10 * trainable_parameters[j].grad
+                trainable_parameters[j].grad.zero_()  
+        if torch.isnan(loss):
+                break
+
+    for nan_i in range(len(trainable_parameters)):
+        if torch.isnan(trainable_parameters[nan_i])==True or abs(trainable_parameters[nan_i])>1e7:
+            return 1000000, 10000000, "1"
+
+    # get the updated symbolic regression
+    ii = -1
+    for parm in unsnapped_param_dict:
+        if ii == -1:
+            ii = ii + 1
+        else:
+            eq = eq.subs(parm, trainable_parameters[ii])
+            ii = ii + 1
+
+    is_atomic_number = lambda expr: expr.is_Atom and expr.is_number
+    numbers_expr = [subexpression for subexpression in preorder_traversal(eq) if is_atomic_number(subexpression)]
+    complexity = 0
+    for j in numbers_expr:
+        try:
+            complexity = complexity + get_number_DL_snapped(float(j))
+        except:
+            complexity = complexity + 1000000
+    n_variables = len(eq.free_symbols)
+    n_operations = len(count_ops(eq,visual=True).free_symbols)
+    if n_operations!=0 or n_variables!=0:
+        complexity = complexity + (n_variables+n_operations)*np.log2((n_variables+n_operations))
+
+    error = get_symbolic_expr_error(data,str(eq))
+    return error, complexity, eq

prior-art/Code/S_get_number_DL.py

@@ -0,0 +1,18 @@
+# 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))
+    elif np.abs(n - np.pi) < epsilon:
+        return np.log2(1+3)
+    # check if real
+    else:
+        PrecisionFloorLoss = 1e-14
+        return np.log2(1 + (float(n) / PrecisionFloorLoss) ** 2) / 2
+

prior-art/Code/S_get_number_DL_snapped.py

@@ -0,0 +1,23 @@
+# 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))
+    elif np.abs(n - np.pi) < epsilon:
+        return np.log2(1+3)
+    else:
+        PrecisionFloorLoss = 1e-14
+        return np.log2(1 + (float(n) / PrecisionFloorLoss) ** 2) / 2
+    
+    
+    

prior-art/Code/S_get_symbolic_expr_error.py

@@ -0,0 +1,42 @@
+# 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 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(data,expr):
+    try:
+        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]]
+        
+        # Remove accidental nan's
+        good_idx = np.where(np.isnan(f(*real_variables))==False)
+       
+        # 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:
+            try:
+                #return np.sqrt(np.mean((f(*real_variables)[good_idx]-data[good_idx][:,-1])**2))/np.sqrt(np.mean(data[good_idx][:,-1]**2))
+                return np.mean(np.log2(1+abs(f(*real_variables)[good_idx]-data[good_idx][:,-1])*2**30))
+            except:
+                # use this for the case in which the expression is just one number (i.e. not array)
+                #return np.sqrt(np.mean((f(*real_variables)-data[:,-1])**2))/np.sqrt(np.mean(data[:,-1]**2))
+                return np.mean(np.log2(1+abs(f(*real_variables)-data[:,-1])*2**30))
+    except:
+        return 1000000

prior-art/Code/S_polyfit.py

@@ -0,0 +1,89 @@
+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)]
+
+        if np.isnan(C_1_2).any()==False:
+            variables = np.matmul(C_1_2,variables.T).T
+            res = getBest(variables,f_dependent,maxdeg)
+            parameters = res[0]
+            params_error = res[1]
+            deg = res[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:
+            res = getBest(variables,f_dependent,maxdeg)
+            parameters = res[0]
+            params_error = res[1]
+            deg = res[2]
+
+            eq = mk_sympy_function(parameters,n_variables,deg)
+            for i in range(len(x)):
+                eq = eq.subs(z[i],x[i])
+            eq = simplify(eq)
+
+    else:
+        res = getBest(variables,f_dependent,maxdeg)
+        parameters = res[0]
+        params_error = res[1]
+        deg = res[2]
+        eq = mk_sympy_function(parameters,n_variables,deg)
+        try:
+            eq = eq.subs("z0","x0")  
+        except:
+            pass
+
+    return (eq, params_error)
+

prior-art/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)
+

prior-art/Code/S_remove_input_neuron.py

@@ -0,0 +1,33 @@
+# Remove on input neuron from a NN
+
+from __future__ import print_function
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+import torch.optim as optim
+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
+
+is_cuda = torch.cuda.is_available()
+
+def remove_input_neuron(net,n_inp,idx_neuron,ct_median,save_filename):
+    removed_weights = net.linear1.weight[:,idx_neuron]
+    # Remove the weights associated with the removed input neuron                                                                                                               
+    t = torch.transpose(net.linear1.weight,0,1)
+    preserved_ids = torch.LongTensor(np.array(list(set(range(n_inp)) - set([idx_neuron]))))
+    t = nn.Parameter(t[preserved_ids, :])
+    net.linear1.weight = nn.Parameter(torch.transpose(t,0,1))
+    # Adjust the biases
+    if is_cuda:
+        net.linear1.bias = nn.Parameter(net.linear1.bias+torch.tensor(ct_median*removed_weights).float().cuda())
+    else:
+        net.linear1.bias = nn.Parameter(net.linear1.bias+torch.tensor(ct_median*removed_weights).float())
+    torch.save(net.state_dict(), save_filename)
+

prior-art/Code/S_run_aifeynman.py

@@ -0,0 +1,221 @@
+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
+from S_final_gd import final_gd
+from S_add_bf_on_numbers_on_pareto import add_bf_on_numbers_on_pareto
+from dimensionalAnalysis import dimensionalAnalysis
+
+PA = ParetoSet()
+def run_AI_all(pathdir,filename,BF_try_time=60,BF_ops_file_type="14ops", polyfit_deg=3, 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 len(data[0])<3:
+        print("Just one variable!")
+        pass
+    elif 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")
+    elif path.exists("results/NN_trained_models/models/" + filename + "_pretrained.h5"):
+        print("Found pretrained NN \n")
+        NN_train(pathdir,filename,NN_epochs/2,lrs=1e-3,N_red_lr=3,pretrained_path="results/NN_trained_models/models/" + filename + "_pretrained.h5")
+        print("NN loss after training: ", 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_)
+        combine_pareto_data = np.loadtxt(pathdir+filename)
+        PA = combine_pareto(combine_pareto_data,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_)
+        combine_pareto_data = np.loadtxt(pathdir+filename)
+        PA = combine_pareto(combine_pareto_data,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=3, NN_epochs=4000, vars_name=[],test_percentage=0):    
+    # If the variable names are passed, do the dimensional analysis first
+    filename_orig = filename
+    try:
+        if vars_name!=[]:
+            dimensionalAnalysis(pathdir,filename,vars_name)
+            DR_file = filename + "_dim_red_variables.txt"
+            filename = filename + "_dim_red"
+        else:
+            DR_file = ""
+    except:
+        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:(100-test_percentage)*len(input_data)//100]]
+    test_data = input_data[sep_idx[test_percentage*len(input_data)//100:len(input_data)]]
+
+    np.savetxt(pathdir+filename+"_train",train_data)
+    if test_data.size != 0:
+        np.savetxt(pathdir+filename+"_test",test_data)
+
+    PA = ParetoSet()
+    # 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=PA)
+    PA_list = PA.get_pareto_points()
+
+    # Run bf snap on the resulted equations
+    for i in range(len(PA_list)):
+        try:
+            PA = add_bf_on_numbers_on_pareto(pathdir,filename,PA,PA_list[i][-1])
+        except:
+            continue
+    PA_list = PA.get_pareto_points()
+
+    np.savetxt("results/solution_before_snap_%s.txt" %filename,PA_list,fmt="%s")
+
+    # Run zero, integer and rational snap on the resulted equations  
+    for j in range(len(PA_list)):
+        PA = add_snap_expr_on_pareto(pathdir,filename,PA_list[j][-1],PA, "")
+
+    PA_list = PA.get_pareto_points()
+    np.savetxt("results/solution_first_snap_%s.txt" %filename,PA_list,fmt="%s")
+    
+    # Run gradient descent on the data one more time
+    final_gd_data = np.loadtxt(pathdir+filename)
+    for i in range(len(PA_list)):
+        try:
+            gd_update = final_gd(final_gd_data,PA_list[i][-1])
+            PA.add(Point(x=gd_update[1],y=gd_update[0],data=gd_update[2]))
+        except:
+            continue
+  
+    PA_list = PA.get_pareto_points()
+    for j in range(len(PA_list)):
+        PA = add_snap_expr_on_pareto(pathdir,filename,PA_list[j][-1],PA, DR_file)
+
+    list_dt = np.array(PA.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=="" and test_data.size != 0:
+        test_errors = []
+        input_test_data = np.loadtxt(pathdir+filename+"_test")
+        for i in range(len(list_dt)):
+            test_errors = test_errors + [get_symbolic_expr_error(input_test_data,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" %filename_orig,save_data,fmt="%s")
+    return save_data
+    

prior-art/Code/S_run_bf_polyfit.py

@@ -0,0 +1,246 @@
+# 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_get_number_DL_snapped import get_number_DL_snapped
+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
+from S_add_snap_expr_on_pareto import add_snap_expr_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=3, output_type=""):
+    input_data = np.loadtxt(pathdir_transformed+filename)
+#############################################################################################################################
+    if np.isnan(input_data).any()==False:
+        # run BF on the data (+)
+        print("Checking for brute force + \n")
+        brute_force(pathdir_transformed,filename,BF_try_time,BF_ops_file_type,"+")
+        
+        try:
+            # 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 = "asin(" + 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(input_data,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_snapped(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(input_data,eqns[i])
+                    PA.add(Point(x=bf_gd_update[1],y=bf_gd_update[0],data=bf_gd_update[2]))
+                except:
+                    continue
+        except:
+            pass
+
+    #############################################################################################################################
+        # run BF on the data (*)
+        print("Checking for brute force * \n")
+        brute_force(pathdir_transformed,filename,BF_try_time,BF_ops_file_type,"*")
+
+        try:
+            # 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 = "asin(" + 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(input_data,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_snapped(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(input_data,eqns[i])
+                    PA.add(Point(x=bf_gd_update[1],y=bf_gd_update[0],data=bf_gd_update[2]))
+                except:
+                    continue
+        except:
+            pass
+        
+    #############################################################################################################################
+        # run polyfit on the data
+        print("Checking polyfit \n")
+        try:
+            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 = "asin(" + 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(input_data,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_snapped(float(j))
+            try:
+                # 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))
+            except:
+                pass
+
+            #run zero snap on polyfit output
+            PA_poly = ParetoSet()
+            PA_poly.add(Point(x=complexity, y=polyfit_err, data=str(eqn)))
+            PA_poly = add_snap_expr_on_pareto(pathdir, filename, str(eqn), PA_poly)
+
+            for l in range(len(PA_poly.get_pareto_points())):
+                PA.add(Point(PA_poly.get_pareto_points()[l][0],PA_poly.get_pareto_points()[l][1],PA_poly.get_pareto_points()[l][2]))
+        
+        except:
+            pass
+
+        print("Complexity  RMSE  Expression")
+        for pareto_i in range(len(PA.get_pareto_points())):
+            print(PA.get_pareto_points()[pareto_i])
+        
+        return PA
+    else:
+        return PA

prior-art/Code/S_separability.py

@@ -0,0 +1,384 @@
+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
+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,))
+
+        # Pick only data which is close enough to the maximum value (5 times less or higher)                                                                   
+        max_output = np.max(abs(f_dependent))
+        use_idx = np.where(abs(f_dependent)>=max_output/5)
+        f_dependent = f_dependent[use_idx]
+        f_dependent = np.reshape(f_dependent,(len(f_dependent),1))
+        variables = variables[use_idx]
+
+        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)
+
+        

prior-art/Code/S_snap.py

@@ -0,0 +1,85 @@
+# The following are snap functions for finding a best approximated integer or rational number for a real number:
+
+import numpy as np
+from sympy import Rational
+
+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 dict(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 dict(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,10) for x in p))
+    chosen = np.argsort(snaps[:, 3])[:top]    
+    d = dict(list(zip(chosen, snaps[chosen, 1:3])))
+    d = {k:  f"{val[0]}/{val[1]}" for k,val in d.items()}
+    
+    return d
+
+
+

prior-art/Code/S_symmetry.py

@@ -0,0 +1,562 @@
+# 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
+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 S_remove_input_neuron import remove_input_neuron
+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"
+            ct_median = torch.median(torch.from_numpy(variables[:,j]))
+            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)
+            remove_input_neuron(model,n_variables,j,ct_median,"results/NN_trained_models/models/"+filename + "-translated_minus_pretrained.h5")
+            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
+            ct_median =torch.median(torch.from_numpy(variables[:,j]))
+            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)
+            remove_input_neuron(model,n_variables,j,ct_median,"results/NN_trained_models/models/"+filename + "-translated_divide_pretrained.h5")
+            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
+            ct_median =torch.median(torch.from_numpy(variables[:,j]))
+            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)
+            remove_input_neuron(model,n_variables,j,ct_median,"results/NN_trained_models/models/"+filename + "-translated_multiply_pretrained.h5")
+            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
+            ct_median =torch.median(torch.from_numpy(variables[:,j]))
+            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)
+            remove_input_neuron(model,n_variables,j,ct_median,"results/NN_trained_models/models/"+filename + "-translated_plus_pretrained.h5")
+            return ("results/translated_data_plus/", file_name)
+        
+    except Exception as e:
+        print(e)
+        return (-1,-1)

prior-art/Code/ai_feynman_example.py

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

prior-art/Code/ai_feynman_terminal_example.py

@@ -0,0 +1,19 @@
+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=3, 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")
+parser.add_argument("--vars_name", type=list, default=[], help="List with the names of the variables")
+parser.add_argument("--test_percentage", type=float, default=0, help="Percentage of the input data to be kept as the test set")
+
+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, vars_name=opts.vars_name, test_percentage=opts.test_percentage)
+          

prior-art/Code/brute_force_oneFile_mdl_v2.scr

@@ -0,0 +1,23 @@
+#!/bin/bash
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+opsfile=$1
+maxtime=$2
+f=$3
+sigma=$4
+band=$5
+
+outfile=brute_solutions.dat
+outfile2=brute_constant.dat
+outfile3=brute_formulas.dat
+if [ -f $outfile ]; then /bin/rm $outfile; fi
+if [ -f $outfile2 ]; then /bin/rm $outfile2; fi
+if [ -f $outfile3 ]; then /bin/rm $outfile3; fi
+
+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 $sigma $band >args.dat
+timeout $maxtime ./symbolic_regress_mdl2.x

prior-art/Code/brute_force_oneFile_mdl_v3.scr

@@ -0,0 +1,23 @@
+#!/bin/bash
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+opsfile=$1
+maxtime=$2
+f=$3
+sigma=$4
+band=$5
+
+outfile=brute_solutions.dat
+outfile2=brute_constant.dat
+outfile3=brute_formulas.dat
+if [ -f $outfile ]; then /bin/rm $outfile; fi
+if [ -f $outfile2 ]; then /bin/rm $outfile2; fi
+if [ -f $outfile3 ]; then /bin/rm $outfile3; fi
+
+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 "$sigma" "$band" >args.dat
+timeout $maxtime ./symbolic_regress_mdl3.x;

prior-art/Code/brute_force_oneFile_v1.scr

@@ -0,0 +1,20 @@
+#!/bin/bash
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+opsfile=$1
+maxtime=$2
+f=$3
+
+outfile=brute_solutions.dat
+outfile2=brute_constant.dat
+
+if [ -f $outfile ]; then /bin/rm $outfile; fi
+if [ -f $outfile2 ]; then /bin/rm $outfile2; fi
+
+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 "$sigma" "$band" >args.dat
+timeout $maxtime ./symbolic_regress1.x;

prior-art/Code/brute_force_oneFile_v2.scr

@@ -0,0 +1,22 @@
+#!/bin/bash
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+opsfile=$1
+maxtime=$2
+f=$3
+
+outfile=brute_solutions.dat
+outfile2=brute_constant.dat
+outfile3=brute_formulas.dat
+if [ -f $outfile ]; then /bin/rm $outfile; fi
+if [ -f $outfile2 ]; then /bin/rm $outfile2; fi
+if [ -f $outfile3 ]; then /bin/rm $outfile3; fi
+
+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 ./symbolic_regress2.x;
+

prior-art/Code/brute_force_oneFile_v3.scr

@@ -0,0 +1,21 @@
+#!/bin/bash
+# USAGE EXAMPLE: solve_mysteries.scr ops6.txt 2
+# USAGE EXAMPLE: solve_mysteries.scr allops.txt 1800
+opsfile=$1
+maxtime=$2
+f=$3
+
+outfile=brute_solutions.dat
+outfile2=brute_constant.dat
+outfile3=brute_formulas.dat
+if [ -f $outfile ]; then /bin/rm $outfile; fi
+if [ -f $outfile2 ]; then /bin/rm $outfile2; fi
+if [ -f $outfile3 ]; then /bin/rm $outfile3; fi
+
+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 ./symbolic_regress3.x;

prior-art/Code/compile.sh

@@ -0,0 +1,12 @@
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress1.x symbolic_regress1.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
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress_mdl2.x symbolic_regress_mdl2.f
+gfortran -ffixed-line-length-none -O3 -o symbolic_regress_mdl3.x symbolic_regress_mdl3.f
+
+chmod 555 brute_force_oneFile_v1.scr
+chmod 555 brute_force_oneFile_v2.scr
+chmod 555 brute_force_oneFile_v3.scr
+chmod 555 brute_force_oneFile_mdl_v2.scr
+chmod 555 brute_force_oneFile_mdl_v3.scr
+

prior-art/Code/dimensionalAnalysis.py

@@ -0,0 +1,129 @@
+import numpy as np
+import pandas as pd
+from scipy.sparse.linalg import lsqr
+from scipy.linalg import *
+from sympy import Matrix
+from sympy import symbols, Add, Mul, S
+from getPowers import getPowers
+
+def dimensional_analysis(input,output,units):
+    M = units[input[0]]
+    for i in range(1,len(input)):
+        M = np.c_[M, units[input[i]]]
+    if len(input)==1:
+        M = np.array(M)
+        M = np.reshape(M,(len(M),1))
+    params = getPowers(M,units[output])
+    M = Matrix(M)
+    B = M.nullspace()
+    return (params, B)
+
+# load the data from a file
+def load_data(pathdir, filename):
+    n_variables = np.loadtxt(pathdir+filename, dtype='str').shape[1]-1
+    variables = np.loadtxt(pathdir+filename, usecols=(0,))
+    for i in range(1,n_variables):
+        v = np.loadtxt(pathdir+filename, usecols=(i,))
+        variables = np.column_stack((variables,v))
+    f_dependent = np.loadtxt(pathdir+filename, usecols=(n_variables,))
+    return(variables.T,f_dependent)
+
+def dimensionalAnalysis(pathdir, filename, eq_symbols):
+    file = pd.read_excel("units.xlsx")
+
+    units = {}
+    for i in range(len(file["Variable"])):
+        val = [file["m"][i],file["s"][i],file["kg"][i],file["T"][i],file["V"][i],file["cd"][i]]
+        val = np.array(val)
+        units[file["Variable"][i]] = val
+
+    dependent_var = eq_symbols[-1]
+
+    file_sym = open(filename + "_dim_red_variables.txt" ,"w")
+    file_sym.write(filename)
+    file_sym.write(", ")
+
+    # load the data corresponding to the first line (from mystery_world)
+    varibs = load_data(pathdir,filename)[0]
+    deps = load_data(pathdir,filename)[1]
+
+    # get the data in symbolic form and associate the corresponding values to it
+    input = []
+    for i in range(len(eq_symbols)-1):
+        input = input + [eq_symbols[i]]
+        vars()[eq_symbols[i]] = varibs[i]
+    output = dependent_var
+
+    # Check if all the independent variables are dimensionless
+    ok = 0
+    for j in range(len(input)):
+        if(units[input[j]].any()):
+            ok=1
+
+    if ok==0:
+        dimless_data = load_data(pathdir, filename)[0].T
+        dimless_dep = load_data(pathdir, filename)[1]
+        if dimless_data.ndim==1:
+            dimless_data = np.reshape(dimless_data,(1,len(dimless_data)))
+            dimless_data = dimless_data.T
+        np.savetxt(pathdir + filename + "_dim_red", dimless_data)
+        file_sym.write(", ")
+        for j in range(len(input)):
+            file_sym.write(str(input[j]))
+            file_sym.write(", ")
+        file_sym.write("\n")
+    else:
+        # get the symbolic form of the solved part
+        solved_powers = dimensional_analysis(input,output,units)[0]
+        input_sym = symbols(input)
+        sol = symbols("sol")
+        sol = 1
+        for i in range(len(input_sym)):
+            sol = sol*input_sym[i]**np.round(solved_powers[i],2)
+        file_sym.write(str(sol))
+        file_sym.write(", ")
+
+        # get the symbolic form of the unsolved part
+        unsolved_powers = dimensional_analysis(input,output,units)[1]
+
+        #print(unsolved_powers,unsolved_powers[0])
+        uns = symbols("uns")
+        unsolved = []
+        for i in range(len(unsolved_powers)):
+            uns = 1
+            for j in range(len(unsolved_powers[i])):
+                uns = uns*input_sym[j]**unsolved_powers[i][j]
+            file_sym.write(str(uns))
+            file_sym.write(", ")
+            unsolved = unsolved + [uns]
+        file_sym.write("\n")
+
+        # get the discovered part of the function
+        func = 1
+        for j in range(len(input)):
+            func = func * vars()[input[j]]**dimensional_analysis(input,output,units)[0][j]
+        func = np.array(func)
+
+        # get the new variables needed
+        new_vars = []
+        for i in range(len(dimensional_analysis(input,output,units)[1])):
+            nv = 1
+            for j in range(len(input)):
+                nv = nv*vars()[input[j]]**dimensional_analysis(input,output,units)[1][i][j]
+            new_vars = new_vars + [nv]
+
+        new_vars = np.array(new_vars)
+        new_dependent = deps/func
+
+        if new_vars.size==0:
+            np.savetxt(pathdir + filename + "_dim_red", new_dependent)
+
+        # save this to file
+        all_variables = np.vstack((new_vars, new_dependent)).T
+        np.savetxt(pathdir + filename + "_dim_red", all_variables)
+
+    file_sym.close()
+
+
+#print(dimensionalAnalysis("../_noise_data/", "119_1.24.6", ["m","omega","omega_0","x","E_n"]))
+

prior-art/Code/generate_claimed_results.py

@@ -0,0 +1,115 @@
+import logging
+import argparse
+import pathlib
+import os
+
+from threading import active_count
+from multiprocessing import Pool
+from multiprocessing.pool import ThreadPool
+from random import shuffle
+from tabulate import tabulate
+from pathlib import Path
+from functools import partial
+
+
+from S_run_aifeynman import run_aifeynman
+
+_CFG = {
+  "dataset_path"      : "../Feynman_without_units/",
+  "operations_file"   : "./14ops.txt",
+  "polynomial_degree" : 3,
+  "number_of_epochs"  : 500,
+  "bruteforce_time"   : 60,
+  "test_percentage"   : 0,
+}
+
+class RunAll:
+    """
+    Run the solver on the whole dataset
+    """
+
+    def __init__(self, *,  cfg=_CFG):
+        logging.basicConfig(filename="output_no_units_parallel.log", level=logging.DEBUG)
+        self.cfg = cfg
+        self.results = {}
+
+ 
+    def print_results(self):
+        table = []
+        for file, sol in self.results.items():
+            table.append(sol[-1])
+        print(tabulate(
+                table,
+                headers=[
+                    "Average error",
+                    "Cumulative error",
+                    "Error",
+                    "Symbolic expression",
+                ],
+            )
+        )
+    
+    def run_solver(self, dirs=None):
+        if not dirs:
+            path = Path(self.cfg["dataset_path"])
+            dirs = list(path.iterdir())
+            shuffle(dirs) # Shuffle to sample a different file each time
+            
+        else:
+            path=Path(self.cfg["dataset_path"])
+            child = dirs
+            
+            
+#         for child in dirs:
+#         print(child)
+        print(f"Process PID: {os.getpid()} ---------------- Number of threads: {active_count()}" )
+        self.results[str(child).split("/")[-1]] = run_aifeynman(
+            pathdir=str(path.resolve()) + "/",
+            filename=str(child).split("/")[-1],
+            BF_try_time=int(self.cfg["bruteforce_time"]),
+            BF_ops_file_type=Path(self.cfg["operations_file"]),
+            polyfit_deg=int(self.cfg["polynomial_degree"]),
+            NN_epochs=int(self.cfg["number_of_epochs"]),
+            vars_name=[],
+            test_percentage=int(self.cfg["test_percentage"]),
+        )
+
+        logging.info(self.results)
+        print("@"*120)
+        print("@"*120)
+
+        self.print_results()
+
+
+def get_files(dirs, chunks=5):
+    dirs = list(path.iterdir())
+    dirs = [file for file in dirs if not (str(file).endswith("test") or str(file).endswith("train"))]
+    for i in range(0, len(dirs), chunks):
+        yield dirs[i : i + chunks]
+
+if __name__ == "__main__":
+
+    #cfg_path = pathlib.Path("/home/aziz/lambda_lab/AI-Feynman/configs.cfg")
+    #if cfg_path.exists():
+    #    RunAll(cfg_path=cfg_path)
+    #else:
+    #    print(f"No such a file {cfg_path}")
+
+    solver = RunAll().run_solver
+    path = Path(_CFG["dataset_path"])
+    #dirs = list(path.iterdir())
+    #chunked_dirs = list(get_files(dirs, chunks=24))
+# print(chunked_dirs[0], len(chunked_dirs[0]))
+#    for dd in chunked_dirs:
+ #       pool = Pool(len(dd))
+  #      print(dd, len(dd))
+    #    pool.map(print, dd)
+   #     pool.map(solver, dd)
+   # pool.close()
+   
+    parser = argparse.ArgumentParser(description='Solver')
+    parser.add_argument('--file', help='Enter file path')
+
+    args = parser.parse_args()
+    solver(args.file)
+

prior-art/Code/getPowers.py

@@ -0,0 +1,62 @@
+import numpy as np
+import pandas as pd
+from scipy.sparse.linalg import lsqr
+from scipy.linalg import *
+from sympy import Matrix
+from sympy import symbols, Add, Mul, S
+from numpy.linalg import matrix_rank
+from itertools import combinations
+
+
+N = np.array([[ 0,  1,  1],
+              [ 0, -1, -1],
+              [ 1,  0, 0],
+              [ 0,  0,  0],
+              [ 0,  0,  0],
+              [ 0,  0, 0],])
+
+N = np.array([[ 0,  0,  3,  1,  1,  1,  1,  1,  1],
+              [ 0,  0, -2,  0,  0,  0,  0,  0,  0],
+              [ 1,  1, -1,  0,  0,  0,  0,  0,  0],
+              [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
+              [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
+              [ 0,  0,  0,  0,  0,  0,  0,  0,  0]])
+
+a = np.array([ 1, -2,  1,  0,  0,  0])
+
+def getPowers(N,a):
+    rand_drop_cols = np.arange(0,len(N[0]),1)
+    rand_drop_rows = np.arange(0,len(N),1)
+    rand_drop_rows = np.flip(rand_drop_rows)
+    rank = matrix_rank(N)
+    d_cols = list(combinations(rand_drop_cols,len(N[0])-rank))
+    d_rows = list(combinations(rand_drop_rows,len(N)-rank))
+    for i in d_cols:
+        M = N
+        M = np.delete(M,i,1)
+        M = np.transpose(M)
+        for j in d_rows:
+            P = M
+            P = np.delete(P,j,1)
+            if np.linalg.det(P)!=0:
+                solved_M = np.transpose(P)
+                indices_sol = j
+                indices_powers = i
+                break
+
+    b = np.delete(a,indices_sol)
+    params = np.linalg.solve(solved_M,b)
+
+    sol = []
+    for i in range(len(N[0])):
+        if i in indices_powers:
+            sol = sol + [0]
+        else:
+            sol = sol + [params[0]]
+            params = np.delete(params,0)
+
+    # this is the solution:
+    sol = np.array(sol)
+    return(sol)
+
+

prior-art/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()
+

prior-art/Code/symbolic_regress1.f

@@ -0,0 +1,157 @@
+	! Max Tegmark 171119, 190128-31, 190506
+	! Loads templates.csv functions.dat and mystery.dat, returns winner.
+	! scp -P2222 symbolic_regress1.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress1.x symbolic_regress1.f |& more'
+	! SAMPLE USAGE: call symbolic_regress1.x 10ops.txt arity2templates.txt mystery_constant.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:
+        !  O: double    (x->2*x); note that this is the letter "O", not zero
+	!  J: double+1  (x->2*x+1)
+	!  >: 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(21), 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.00000001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\OJLESCANTR01P"/
+      	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)
+
+        nvar = 0
+ 	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
+
+        include 'tools.f'

prior-art/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(21), 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.00001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\OJLESCANTR01P"/
+      	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
+

prior-art/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(21), 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.00001)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\OJLESCANTR01P"/
+      	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
+

prior-art/Code/symbolic_regress_mdl2.f

@@ -0,0 +1,220 @@
+	! Max Tegmark 171119, 190128-31, 190506, 200427-29
+	! Loads templates.csv functions.dat and mystery.dat, returns winners.
+        ! Rejects Pareto-dominated formulas not based on hard sup-norm cut, but using a 
+        !   hypothesis-testing framework with a z-score z_n = sqrt(n)*(<b_n>-<b_best>)/sigma_best
+	! scp -P2222 symbolic_regress.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress_mdl2.x symbolic_regress_mdl2.f |& more'
+	! SAMPLE USAGE: call symbolic_regress_mdl2.x 7ops.txt arity2templates.txt mystery2.dat results.dat 10 0
+        !               call symbolic_regress_mdl2.x 6ops.txt arity2templates.txt mysteryB3.dat results.dat 10 0 (takes a few minutes)
+	! 		call symbolic_regress_mdl2.x 14ops.txt arity2templates.txt mystery.dat results.dat 10 0
+	! 		call symbolic_regress_mdl2.x 14ops.txt arity2templates.txt mystery.dat results.dat 1000 0 (if skips over correct formula)
+	! 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(21), nvar, nvarmax, nmax, lnblnk
+	parameter(nvarmax=20, nmax=5000000)
+	real*8 f, newloss, minloss, maxloss, rmsloss, limit
+        real*8 xy0(nvarmax+1,nmax), xy(nvarmax+1,nmax), offset(nmax), offst, bestoffset
+        real*8 epsilon, DL, nu, z
+        real*8 lossbits, bitmean, bitsdev, bestbits, bitmargin, sigma, bitexcess, ev
+	parameter(epsilon=1/2.**30)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\OJLESCANTR01P"/
+      	integer nn(0:2), ii(nmax), kk(nmax), radix(nmax), iarr(nmax)
+	integer ndata, i, j, jtest, n
+	integer*8 nformulas, nevals
+	logical done, rejected
+	character*60 func(0:2), template
+
+	nu = 5.
+        bitmargin = 0. ! "Thickness" of pareto frontier; default 0
+ 	open(2,file='args.dat',status='old',err=666)
+        read(2,*) opsfile, templatefile, mysteryfile, outfile, nu, bitmargin
+ 	write(*,'(1a24,f10.3)') 'Rejection threshold.....',nu
+	write(*,'(1a24,f10.3)') 'Bit margin..............',bitmargin
+
+        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,xy0,mysteryfile)	
+	write(*,'(1a24,i8)') 'Number of examples......',ndata
+
+	write(*,'(1a24)') 'Shuffling mystery data....'
+        call permutation(ndata,iarr)
+	do i=1,ndata
+          do j=1,nvar+1
+            xy(j,i) = xy0(j,iarr(i))
+          end do
+        end do    	
+
+	print *,'Searching for best fit...'
+	nformulas = 0
+        nevals = 0
+	bestbits = 1.e6
+        sigma = 1.d40 ! So that 1st function gets accepted
+	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 ((bestbits.gt.0).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))
+	  end do
+	  j = 1
+ 	  jtest = 2  ! Will test after j=2, 3, 5, 9, 17, ... data points
+	  rejected = .false.
+          do while ((.not.rejected).and.(j.le.ndata)) ! Keep going as long as you can't reject this formula
+ 	    nevals = nevals + 1
+            offst = xy(nvar+1,j) - f(n,ii,ops,xy(1,j))
+            rejected = (.not.((offst.ge.0).or.(offst.le.0))) ! This was a NaN, so reject the formula :-)
+            !if (rejected) print *,"NaN!"
+            if (rejected) exit
+            rejected = abs(offst).gt.(1./epsilon) ! Otherwise numerical cancellation can masquerade as successss
+            !if (rejected) print *,"Infinity!"
+            if (rejected) exit     
+	    offset(j) = offst    
+            if (j.ge.jtest) then ! Time for another test                
+                call analyze_offset(j,offset,epsilon,bestoffset,bitmean,bitsdev)
+                bitexcess = bitmean - bestbits - bitmargin
+                z = sqrt(1.*j)*bitexcess/sigma ! This sigma is for previous winner, not for this candidate
+                rejected = (z.gt.nu)
+                jtest = min(2*jtest-1,ndata)
+            end if  
+    	    j = j + 1
+	  end do
+	  if (.not.rejected.and.(bitexcess.lt.0.)) then ! We have a new point on the Pareto frontier
+	    bestbits = min(bitmean,bestbits)
+	    rmsloss = 0.
+            maxloss = 0.
+            sigma = 0.
+	    do j=1,ndata
+              newloss = abs(xy(nvar+1,j) - f(n,ii,ops,xy(1,j)) - bestoffset)
+	      rmsloss = rmsloss + newloss**2
+              if (maxloss.lt.newloss) maxloss = newloss
+	    end do
+	    rmsloss = sqrt(rmsloss/ndata)
+            sigma = bitsdev
+	    DL  = log(1.*nformulas)/log(2.)
+            ev = (1.*nevals)/nformulas
+	    write(*,'(2f20.12,x,1a22,1i16,6f19.4)') bitmean, limit(bestoffset), ops(1:n), nformulas, DL, DL+ndata*bitmean, rmsloss, maxloss, bitsdev, ev
+	    write(3,'(2f20.12,x,1a22,1i16,6f19.4)') bitmean, limit(bestoffset), ops(1:n), nformulas, DL, DL+ndata*bitmean, rmsloss, maxloss, bitsdev, ev
+	    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
+
+	subroutine analyze_offset(n,offset,epsilon,median,bitmean,bitsdev) ! Check how much an array departs from its median
+	implicit none
+        integer n, i
+        real*8 offset(n), epsilon, bitmean, bitsdev
+        real*8 median, mymedian, x, bits, sum1, sum2
+        median = mymedian(n,offset)
+        sum1 = 0.
+        sum2 = 0.
+        do i=1,n
+          x = abs(offset(i)-median)/epsilon        
+          if (x.gt.1) then 
+            bits = 1.44269504089*log(x)  ! = log2(x)
+          else
+            bits = 0.
+          end if
+	  sum1 = sum1 + bits
+          sum2 = sum2 + bits*bits          
+        end do
+        bitmean = sum1/n
+        bitsdev = sqrt(abs(sum2/n-bitmean**2))
+     	return
+        end
+
+	include "tools.f"

prior-art/Code/symbolic_regress_mdl3.f

@@ -0,0 +1,244 @@
+	! Max Tegmark 171119, 190128-31, 190506, 200427-29
+	! Loads templates.csv functions.dat and mystery.dat, returns winners.
+        ! Rejects Pareto-dominated formulas not based on hard sup-norm cut, but using a 
+        !   hypothesis-testing framework with a z-score z_n = sqrt(n)*(<b_n>-<b_best>)/sigma_best
+	! scp -P2222 symbolic_regress.f euler@tor.mit.edu:FEYNMAN
+	! COMPILATION: a f 'f77 -O3 -o symbolic_regress_mdl3.x symbolic_regress_mdl3.f |& more'
+	! SAMPLE USAGE: call symbolic_regress_mdl3.x 7ops.txt arity2templates.txt mystery2.dat results.dat 10 0
+        !               call symbolic_regress_mdl3.x 6ops.txt arity2templates.txt mysteryB3.dat results.dat 10 0 (takes a few minutes)
+	! 		call symbolic_regress_mdl3.x 14ops.txt arity2templates.txt mystery.dat results.dat 10 0
+	! 		call symbolic_regress_mdl3.x 14ops.txt arity2templates.txt mystery.dat results.dat 1000 0 (if skips over correct formula)
+	! 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(21), nvar, nvarmax, nmax, lnblnk
+	parameter(nvarmax=20, nmax=5000000)
+	real*8 f, newloss, minloss, maxloss, rmsloss, limit
+        real*8 xy0(nvarmax+1,nmax), xy(nvarmax+1,nmax), y(nmax), offset(nmax), offst, bestoffset
+        real*8 epsilon, DL, nu, z
+        real*8 lossbits, bitmean, bitsdev, bestbits, bitmargin, sigma, bitexcess, ev
+        real*8 ymin, ymax
+	parameter(epsilon=1/2.**30)
+	data arities /2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0/
+	data functions /"+*-/><~\OJLESCANTR01P"/
+      	integer nn(0:2), ii(nmax), kk(nmax), radix(nmax), iarr(nmax)
+	integer ndata, i, i1, j, jtest, n
+	integer*8 nformulas, nevals
+	logical done, rejected
+	character*60 func(0:2), template
+	nu = 5.
+        bitmargin = 0. ! "Thickness" of pareto frontier; default 0
+ 	open(2,file='args.dat',status='old',err=666)
+        read(2,*) opsfile, templatefile, mysteryfile, outfile, nu, bitmargin
+ 	write(*,'(1a24,f10.3)') 'Rejection threshold.....',nu
+	write(*,'(1a24,f10.3)') 'Bit margin..............',bitmargin
+
+        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,xy0,mysteryfile)	
+	write(*,'(1a24,i8)') 'Number of examples......',ndata
+
+	write(*,'(1a24)') 'Removing problematically small data points....'
+        ymax = 0.
+        do i=1,ndata
+          if (ymax.lt.abs(xy0(nvar+1,i))) ymax=abs(xy0(nvar+1,i))
+        end do
+        ymin = 0.001*ymax ! Require all data to exceed this
+        i1 = 0
+        print *,ymax,ymin
+        do i=1,ndata
+          if (abs(xy0(nvar+1,i)).gt.ymin) then ! Keep this data point
+            i1 = i1 + 1
+            do j=1,nvar+1
+              xy0(j,i1) = xy0(j,i1)
+            end do
+          end if
+        end do
+        write(*,*) ndata-i1," out of ",ndata," data points discarded for being too close to zero"
+        ndata = i1  
+        
+	write(*,'(1a24)') 'Shuffling mystery data....'
+        call permutation(ndata,iarr)
+	do i=1,ndata
+          do j=1,nvar+1
+            xy(j,i) = xy0(j,iarr(i))
+          end do
+            y(i) = xy(nvar+1,i)
+        end do    	
+
+	print *,'Searching for best fit...'
+	nformulas = 0
+        nevals = 0
+	bestbits = 1.e6
+        sigma = 1.d40 ! So that 1st function gets accepted
+	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 ((bestbits.gt.0).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))
+	  end do
+	  j = 1
+ 	  jtest = 2  ! Will test after j=2, 3, 5, 9, 17, ... data points
+	  rejected = .false.
+          do while ((.not.rejected).and.(j.le.ndata)) ! Keep going as long as you can't reject this formula
+ 	    nevals = nevals + 1
+            offst = y(j)/f(n,ii,ops,xy(1,j))
+            rejected = (.not.((offst.ge.0).or.(offst.le.0))) ! This was a NaN, so reject the formula :-)
+            !if (rejected) print *,"NaN!"
+            if (rejected) exit
+            rejected = (abs(offst).lt.epsilon).or.(abs(offst).gt.1./epsilon) ! Otherwise numerical cancellation can masquerade as success
+            !rejected = abs(log(abs(offst))).gt.(1./epsilon) ! Otherwise numerical cancellation can masquerade as successss
+            !if (rejected) print *,"Infinity!"
+            if (rejected) exit     
+	    offset(j) = offst    
+            if (j.ge.jtest) then ! Time for another test                
+                call analyze_offset(j,y,offset,epsilon,bestoffset,bitmean,bitsdev)
+                bitexcess = bitmean - bestbits - bitmargin
+                z = sqrt(1.*j)*bitexcess/sigma ! This sigma is for previous winner, not for this candidate
+                rejected = (z.gt.nu)
+                jtest = min(2*jtest-1,ndata)
+            end if  
+    	    j = j + 1
+	  end do
+	  if (.not.rejected.and.(bitexcess.lt.0.)) then ! We have a new point on the Pareto frontier
+	    bestbits = min(bitmean,bestbits)
+	    rmsloss = 0.
+            maxloss = 0.
+            sigma = 0.
+	    do j=1,ndata
+              newloss = abs(y(j) - f(n,ii,ops,xy(1,j))*bestoffset)
+	      rmsloss = rmsloss + newloss**2
+              if (maxloss.lt.newloss) maxloss = newloss
+	    end do
+	    rmsloss = sqrt(rmsloss/ndata)
+            sigma = bitsdev
+	    DL  = log(1.*nformulas)/log(2.)
+            ev = (1.*nevals)/nformulas
+	    write(*,'(2f20.12,x,1a22,1i16,6f19.4)') bitmean, limit(bestoffset), ops(1:n), nformulas, DL, DL+ndata*bitmean, rmsloss, maxloss, bitsdev, ev
+	    write(3,'(2f20.12,x,1a22,1i16,6f19.4)') bitmean, limit(bestoffset), ops(1:n), nformulas, DL, DL+ndata*bitmean, rmsloss, maxloss, bitsdev, ev
+	    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
+
+	subroutine analyze_offset(n,y,offset,epsilon,median,bitmean,bitsdev) ! Check how much an array departs from its median
+	implicit none
+        integer n, i
+        real*8 y(n), offset(n), epsilon, bitmean, bitsdev
+        real*8 median, mymedian, x, f, bits, sum1, sum2
+        median = mymedian(n,offset)
+        sum1 = 0.
+        sum2 = 0.
+        do i=1,n
+          f = y(i)/offset(i)
+          x = abs(y(i)-median*f)/epsilon        
+          if (x.gt.1) then 
+            bits = 1.44269504089*log(x)  ! = log2(x)
+          else
+            bits = 0.
+          end if
+	  sum1 = sum1 + bits
+          sum2 = sum2 + bits*bits          
+          !print *,i,y(i),offset(i),abs(y(i)*(offset(i)-median)),x
+          !read *
+        end do
+        bitmean = sum1/n
+        bitsdev = sqrt(abs(sum2/n-bitmean**2))
+     	return
+        end
+
+       include "tools.f"

prior-art/Code/tools.f

@@ -0,0 +1,380 @@
+! Max Tegmark 171119, 190128-31, 190506, May 2020
+
+	! 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))
+        !  O: double    (x->2*x); note that this is the letter "O", not zero
+	!  J: double+1  (x->2*x+1)
+	! nonary: 
+	!  0
+	!  1
+	!  P: pi
+	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 if (op.eq."O") then
+             y = 2.*stack(j)
+	    else if (op.eq."J") then
+             y = 1+2.*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
+
+	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
+
+	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
+
+	real*8 function mymedian(n,a)
+        implicit none
+        integer n,nmax, i
+        parameter(nmax=10000000)
+        real*8 a(n), b(nmax)
+        if (n.gt.nmax) stop 'DEATH ERROR: n>nmax in mymedian'
+        do i=1,n
+          b(i) = a(i)
+        end do
+        call sort(n,b)
+        i = nint((n+.5)/2)
+        if (i.eq.0) i=1
+        mymedian = b(i)
+        return
+        end
+
+	subroutine permutation(n,iarr) ! Return a random permutation of the first n integer:
+	integer iarr(n), idum, nmax, i
+        parameter(nmax=10000000)
+        real*8 arr(nmax), brr(nmax), ran1
+        if (n.gt.nmax) stop "PERMUTATION DEATH ERROR: nmax TOO SMALL"
+        idum = -666
+        do i=1,n
+          arr(i) = ran1(idum)
+          brr(i) = i
+        end do
+        call sort2(n,arr,brr)
+        do i=1,n
+          iarr(i) = nint(brr(i))
+        end do
+	return
+        end        
+        
+      SUBROUTINE sort(n,arr) ! Numerical Recipes Quicksort:
+      INTEGER n,M,NSTACK
+      REAL*8 arr(n)
+      PARAMETER (M=7,NSTACK=50)
+      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
+      REAL*8 a,temp
+      jstack=0
+      l=1
+      ir=n
+1     if(ir-l.lt.M)then
+        do 12 j=l+1,ir
+          a=arr(j)
+          do 11 i=j-1,1,-1
+            if(arr(i).le.a)goto 2
+            arr(i+1)=arr(i)
+11        continue
+          i=0
+2         arr(i+1)=a
+12      continue
+        if(jstack.eq.0)return
+        ir=istack(jstack)
+        l=istack(jstack-1)
+        jstack=jstack-2
+      else
+        k=(l+ir)/2
+        temp=arr(k)
+        arr(k)=arr(l+1)
+        arr(l+1)=temp
+        if(arr(l+1).gt.arr(ir))then
+          temp=arr(l+1)
+          arr(l+1)=arr(ir)
+          arr(ir)=temp
+        endif
+        if(arr(l).gt.arr(ir))then
+          temp=arr(l)
+          arr(l)=arr(ir)
+          arr(ir)=temp
+        endif
+        if(arr(l+1).gt.arr(l))then
+          temp=arr(l+1)
+          arr(l+1)=arr(l)
+          arr(l)=temp
+        endif
+        i=l+1
+        j=ir
+        a=arr(l)
+3       continue
+          i=i+1
+        if(arr(i).lt.a)goto 3
+4       continue
+          j=j-1
+        if(arr(j).gt.a)goto 4
+        if(j.lt.i)goto 5
+        temp=arr(i)
+        arr(i)=arr(j)
+        arr(j)=temp
+        goto 3
+5       arr(l)=arr(j)
+        arr(j)=a
+        jstack=jstack+2
+        if(jstack.gt.NSTACK) stop 'NSTACK too small in sort'
+        if(ir-i+1.ge.j-l)then
+          istack(jstack)=ir
+          istack(jstack-1)=i
+          ir=j-1
+        else
+          istack(jstack)=j-1
+          istack(jstack-1)=l
+          l=i
+        endif
+      endif
+      goto 1
+      END
+
+      SUBROUTINE sort2(n,arr,brr)  ! Numerical Recipes Quicksort:
+      INTEGER n,M,NSTACK
+      REAL*8 arr(n),brr(n)
+      PARAMETER (M=7,NSTACK=50)
+      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
+      REAL*8 a,b,temp
+      jstack=0
+      l=1
+      ir=n
+1     if(ir-l.lt.M)then
+        do 12 j=l+1,ir
+          a=arr(j)
+          b=brr(j)
+          do 11 i=j-1,1,-1
+            if(arr(i).le.a)goto 2
+            arr(i+1)=arr(i)
+            brr(i+1)=brr(i)
+11        continue
+          i=0
+2         arr(i+1)=a
+          brr(i+1)=b
+12      continue
+        if(jstack.eq.0)return
+        ir=istack(jstack)
+        l=istack(jstack-1)
+        jstack=jstack-2
+      else
+        k=(l+ir)/2
+        temp=arr(k)
+        arr(k)=arr(l+1)
+        arr(l+1)=temp
+        temp=brr(k)
+        brr(k)=brr(l+1)
+        brr(l+1)=temp
+        if(arr(l+1).gt.arr(ir))then
+          temp=arr(l+1)
+          arr(l+1)=arr(ir)
+          arr(ir)=temp
+          temp=brr(l+1)
+          brr(l+1)=brr(ir)
+          brr(ir)=temp
+        endif
+        if(arr(l).gt.arr(ir))then
+          temp=arr(l)
+          arr(l)=arr(ir)
+          arr(ir)=temp
+          temp=brr(l)
+          brr(l)=brr(ir)
+          brr(ir)=temp
+        endif
+        if(arr(l+1).gt.arr(l))then
+          temp=arr(l+1)
+          arr(l+1)=arr(l)
+          arr(l)=temp
+          temp=brr(l+1)
+          brr(l+1)=brr(l)
+          brr(l)=temp
+        endif
+        i=l+1
+        j=ir
+        a=arr(l)
+        b=brr(l)
+3       continue
+          i=i+1
+        if(arr(i).lt.a)goto 3
+4       continue
+          j=j-1
+        if(arr(j).gt.a)goto 4
+        if(j.lt.i)goto 5
+        temp=arr(i)
+        arr(i)=arr(j)
+        arr(j)=temp
+        temp=brr(i)
+        brr(i)=brr(j)
+        brr(j)=temp
+        goto 3
+5       arr(l)=arr(j)
+        arr(j)=a
+        brr(l)=brr(j)
+        brr(j)=b
+        jstack=jstack+2
+        if(jstack.gt.NSTACK)stop 'NSTACK too small in sort2'
+        if(ir-i+1.ge.j-l)then
+          istack(jstack)=ir
+          istack(jstack-1)=i
+          ir=j-1
+        else
+          istack(jstack)=j-1
+          istack(jstack-1)=l
+          l=i
+        endif
+      endif
+      goto 1
+      END
+
+      ! Numerical Recipes random number generator:
+      FUNCTION ran1(idum)
+      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
+      REAL*8 ran1,AM,EPS,RNMX
+      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
+     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
+      INTEGER j,k,iv(NTAB),iy
+      SAVE iv,iy
+      DATA iv /NTAB*0/, iy /0/
+      if (idum.le.0.or.iy.eq.0) then
+        idum=max(-idum,1)
+        do 11 j=NTAB+8,1,-1
+          k=idum/IQ
+          idum=IA*(idum-k*IQ)-IR*k
+          if (idum.lt.0) idum=idum+IM
+          if (j.le.NTAB) iv(j)=idum
+11      continue
+        iy=iv(1)
+      endif
+      k=idum/IQ
+      idum=IA*(idum-k*IQ)-IR*k
+      if (idum.lt.0) idum=idum+IM
+      j=1+iy/NDIV
+      iy=iv(j)
+      iv(j)=idum
+      ran1=min(AM*iy,RNMX)
+      return
+      END
+