bootstrap

Code/generate_claimed_results.py

@@ -1,39 +1,43 @@
 import logging
+import argparse
 import pathlib
-import configparser
+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 Code.S_run_aifeynman import run_aifeynman
 
 
+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 all the whole dataset
+    Run the solver on the whole dataset
     """
 
-    def __init__(self, *, cfg_path: Path):
+    def __init__(self, *,  cfg=_CFG):
-        logging.basicConfig(filename="output.log", level=logging.DEBUG)
+        logging.basicConfig(filename="output_no_units_parallel.log", level=logging.DEBUG)
-        self.config = configparser.ConfigParser()
+        self.cfg = cfg
-        self.config.read(cfg_path)
-        self.cfg = self.config["Default"]
-        self.print_results()
         self.results = {}
 
  
-        self.run_solver()
-
-    def log_results(self):
-        pass
-
     def print_results(self):
-        table = [
+        table = []
-            ["foo", 696000, 1989100000],
+        for file, sol in self.results.items():
-            ["bar", 6371, 5973.6],
+            table.append(sol[-1])
-            ["baz", 1737, 73.5],
-            ["qux", 3390, 641.85],
-        ]
         print(tabulate(
                 table,
                 headers=[
@@ -45,12 +49,23 @@ class RunAll:
             )
         )
     
-    def run_solver(self):
+    def run_solver(self, dirs=None):
+        if not dirs:
             path = Path(self.cfg["dataset_path"])
-        for child in path.iterdir():
+            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="/home/aziz/lambda_lab/AI-Feynman/example_data/",#str(path.resolve()) + "/",
+            pathdir=str(path.resolve()) + "/",
-                filename="example2.txt",#str(child).split("/")[-1],
+            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"]),
@@ -58,12 +73,43 @@ class RunAll:
             vars_name=[],
             test_percentage=int(self.cfg["test_percentage"]),
         )
+
         logging.info(self.results)
-            break
+        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():
+    #cfg_path = pathlib.Path("/home/aziz/lambda_lab/AI-Feynman/configs.cfg")
-        RunAll(cfg_path=cfg_path)
+    #if cfg_path.exists():
-    else:
+    #    RunAll(cfg_path=cfg_path)
-        print(f"No such a file {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)
+

output_no_units_parallel.log

@@ -0,0 +1,154 @@
+INFO:root:{'III.13.18': array([[36.57060766285735, 4.550247640408479, 4550247.640408479, 0.0,
+        23.429392474512944, 8.91505505199440e-11],
+       [0.0003960505387426254, -11.302027839797272, -11302027.839797273,
+        6.321928094887362, 0.00039605053874262544, '4*pi*x0']],
+      dtype=object)}
+INFO:root:{'II.37.1': array([['30.000000001343675', '4.9068905956731355', '4906890.5956731355',
+        '0.0', '30.000000001343675', 'x0'],
+       ['3.289228836945565e-07', '-21.535747281929975',
+        '-21535747.281929974', '3.0', '3.2892288369455636e-07', 'x0 + 1']],
+      dtype='<U32')}
+INFO:root:{'III.19.51': array([['24.081314741250495', '4.589842254381291', '4589842.254381292',
+        '0.0', '24.081314741250495', '0'],
+       ['5.647234375079581e-09', '-27.39980834578584',
+        '-27399808.345785838', '13.584962500721156',
+        '5.6472343750795805e-09', '-0.125/x0**2']], dtype='<U32')}
+INFO:root:{'III.14.14': array([['30.85970357081458', '4.878360910171662', '4878360.910171662',
+        '0.0', '29.412569313163143',
+        'asin(0.000000000000*sqrt(exp(((x0+1))**(-1))))'],
+       ['29.210279221326296', '4.868404243855462', '4868404.243855462',
+        '3.0', '29.21027922132631', '1.50000000000000'],
+       ['28.35236907597377', '4.825397384247454', '4825397.384247454',
+        '4.584962500721156', '28.35236907597379', '2/x0'],
+       ['27.81822835451622', '4.797958637501465', '4797958.637501465',
+        '6.0', '27.818228354516222', '1.5/x0'],
+       ['7.268884829473797e-07', '-20.391762617039728',
+        '-20391762.61703973', '10.0', '7.268884829473792e-07',
+        'exp(1/x0) - 1']], dtype='<U46')}
+INFO:root:{'III.4.33': array([['26.290737189268583', '4.716482690414682', '4716482.690414682',
+        '0.0', '26.290737189268583', 'x0'],
+       ['22.658898217944188', '4.502005807142314', '4502005.8071423145',
+        '11.643856189774725', '22.658898217944188', 'x0 - 0.07'],
+       ['22.658481533490033', '4.501979276544167', '4501979.276544167',
+        '16.965784284662085', '22.658481533490022',
+        '(x0*(x0 - 0.16))**0.5'],
+       ['22.548084074885335', '4.494932946856186', '4494932.946856186',
+        '26.416665599935456', '22.548084074885317',
+        '(x0*(x0 + log(log((2 + 4*pi)**(1/pi)))))**0.5'],
+       ['19.11349123138255', '4.256519417010686', '4256519.417010686',
+        '29.80424344959478', '19.11349123138255',
+        '(x0**2 - 0.16*x0 + 0.01)**0.5'],
+       ['18.032683218363186', '4.172542177060247', '4172542.1770602474',
+        '35.64548429043134', '18.03268321836319',
+        '(x0**2 + x0*log(log((2 + 4*pi)**(1/pi))) + 0.01)**0.5'],
+       ['17.157005285884946', '4.100725850740927', '4100725.8507409273',
+        '66.69395636388344', '17.157005285884953',
+        '(x0**2 - 0.1591549430918953*x0 + 0.01)**0.5']], dtype='<U53')}
+INFO:root:{'II.38.14': array([['27.064393896161725', '4.7583241743201885', '4758324.174320188',
+        '0.0', '27.064393896161725', '0'],
+       ['25.603511150614796', '4.678269763402974', '4678269.763402974',
+        '5.0', '25.603511150614807', '0.5/x0'],
+       ['8.47300466346212e-09', '-26.814479190876447',
+        '-26814479.190876447', '8.754887502163468',
+        '8.47300466346212e-09', '0.500000000000*((x0+1))**(-1)']],
+      dtype='<U32')}
+INFO:root:{'II.6.11': array([[24.91363115902248, 4.63886340443597, 4638863.40443597, 0.0,
+        24.91363115902248, '-8.915072558055073e-11'],
+       [24.91363113734539, 4.638863403180696, 4638863.403180696,
+        10.754887502163468, 24.913631137345384, '-pi*exp(-exp(pi))'],
+       [20.52143209643687, 4.359059508402251, 4359059.508402252,
+        12.584962500721156, 20.521432096436865,
+        'asin(0.0833333333333333*cos(x0))'],
+       [12.27445776778958, 3.617587388836739, 3617587.388836739,
+        17.60964047443681, 12.274457767789572, 'asin(cos(x0)/(4*pi))'],
+       [11.607229189021997, 3.53695171632929, 3536951.71632929,
+        28.236446955124386, 11.607229189021998,
+        'asin(pi*sin(cos(x0)/pi**2)/4)'],
+       [11.577529861964756, 3.5332555730823176, 3533255.5730823176,
+        69.39500893210585, 11.577529861964752,
+        'asin(0.785437643527985*sin(cos(x0)/pi**2))'],
+       [11.577405736751007, 3.533240105552333, 3533240.105552333,
+        100.97220344060344, 11.577405736750997,
+        asin(0.785437643527985*sin(0.1013211780033*cos(x0)))]],
+      dtype=object)}
+INFO:root:{'III.15.27': array([['32.651368992645196', '5.029071576628298', '5029071.5766282985',
+        '0.0', '32.65136899264519', '0'],
+       ['nan', 'nan', 'nan', '5.754887502163468', 'nan', 'acos(-x0)'],
+       ['29.570858703869007', '4.886104233165548', '4886104.233165548',
+        '6.754887502163468', '29.570858703869007', 'acos(1 - x0)'],
+       ['28.17969118815505', '4.8165838966073595', '4816583.89660736',
+        '11.807354922057604', '28.17969118815505', '6*x0/x1'],
+       ['6.643797122366517e-06', '-17.199560550061296',
+        '-17199560.550061297', '12.584962500721156',
+        '6.6437971223665156e-06', '2*pi*x0/x1'],
+       ['6.643797122366517e-06', '-17.19959045660223',
+        '-17199590.45660223', '58.15848945789539',
+        '6.643659400307738e-06', '6.283185307179586*x0/x1']], dtype='<U32')}
+INFO:root:{'III.10.19': array([['29.667577204241404', '4.890815209162709', '4890815.209162709',
+        '0.0', '29.667577204241386', 'x1'],
+       ['29.664692327032483', '4.843528401367819', '4843528.401367819',
+        '5.754887502163468', '28.710934848338415',
+        '1/(0.000000000000+(x0/x1))'],
+       ['27.509606914525047', '4.781863619982525', '4781863.619982526',
+        '9.92481250360578', '27.50960691452505', '(x1**2 + 2)**0.5'],
+       ['26.449356617926124', '4.725160724016205', '4725160.724016205',
+        '12.584962500721156', '26.449356617926107',
+        '(x0 + x1**2 + 1)**0.5'],
+       ['2.6825572107565247e-07', '-21.82988772462039',
+        '-21829887.72462039', '14.169925001442312',
+        '2.6825572107565247e-07', '(x0**2 + x1**2 + 1)**0.5']],
+      dtype='<U32')}
+INFO:root:{'II.38.3': array([['29.20416081803319', '4.86810202440029', '4868102.02440029',
+        '0.0', '29.204160818033195', 'x1'],
+       ['3.132865116530134e-06', '-18.283636592181303',
+        '-18283636.592181303', '3.0', '3.133840984176429e-06',
+        '-0.000000000000+((x0/x1))**(-1)'],
+       ['3.132865116530134e-06', '-18.284085912574707',
+        '-18284085.91257471', '5.754887502163468',
+        '3.132865116530135e-06', 'x1/x0'],
+       ['3.132865116530134e-06', '-18.284085912574707',
+        '-18284085.91257471', '10.0', '3.132865116530134e-06',
+        '1.000000000000*(x1/x0)']], dtype='<U32')}
+INFO:root:{'II.6.15a': array([['28.495761948727154', '4.832675464333084', '4832675.464333084',
+        '0.0', '28.495761948727154', '0'],
+       ['27.380642712447607', '4.775084406610207', '4775084.406610208',
+        '2.0', '27.380642712447603', 'log(log(pi))'],
+       ['27.054294104422542', '4.757785694169851', '4757785.694169851',
+        '6.321928094887362', '27.05429410442252', '4*exp(-3)'],
+       ['25.30155195820018', '4.661153975211263', '4661153.975211264',
+        '7.339850002884624', '25.30155195820018',
+        '0.333333333333333*x0*x2'],
+       ['25.276008725588675', '4.659696763845128', '4659696.763845128',
+        '11.0', '25.276008725588685', 'x0*x2/pi'],
+       ['23.272252622734875', '4.540538957128693', '4540538.957128693',
+        '15.194602975157967', '23.27225262273488',
+        '0.166666666666667*x2*(x0 + x1)'],
+       ['21.63601106924245', '4.435362635591665', '4435362.635591665',
+        '29.0', '21.63601106924245', 'x2*(x0 + x1)*exp(-sqrt(pi))'],
+       ['21.588377090562414', '4.4321828875112645', '4432182.8875112645',
+        '29.558375050011747', '21.588377090562417',
+        '2*x2*(x0 + x1)*log(pi)/(1 + 4*pi)'],
+       ['21.18709874287426', '4.405114140557128', '4405114.140557128',
+        '55.550160087467745', '21.18709874287426',
+        '0.168816319869*(x2*(x1+x0))']], dtype='<U33')}
+INFO:root:{'III.9.52': array([[32.84659416339962, 4.913686532069664, 4913686.532069664, 0.0,
+        30.141650893217463, 'tan(-666.000000000000*sin(pi))'],
+       [32.84182393413422, 4.849625676726862, 4849625.676726862,
+        21.724214059872214, 28.832532911636633, 0.0731742799056452],
+       [28.374965935007765, 4.826546755293563, 4826546.755293563,
+        35.98584193700334, 28.374965935007765,
+        '(12 - 12*cos(x1 - x2))/(x0*(x1 - x2)**2)'],
+       [0.0015229204525354075, -9.35894369787947, -9358943.697879469,
+        38.43621560858933, 0.0015229204525354064,
+        '-4*pi*(cos(x1 - x2) - 1)/(x0*(x1 - x2)**2)']], dtype=object)}
+INFO:root:{'III.17.37': array([['31.225739937707367', '4.964663853808277', '4964663.853808277',
+        '0.0', '31.22573993770736', '0'],
+       ['31.257799461003707', '4.908431509806649', '4908431.509806649',
+        '1.0', '30.032059527985425', '1/(1.000000000000*(x0-(x0+1)))'],
+       ['31.243574989647236', '4.868332410905311', '4868332.410905311',
+        '14.60964047443681', '29.208824854162387',
+        '1/(-1.000000000000*(x0-log(x0)))'],
+       ['1.7431741319824306e-07', '-22.45177997151243',
+        '-22451779.97151243', '16.509775004326936',
+        '1.7431741319824306e-07', '0.000000000000+(x0*((x1*cos(x2))+1))']],
+      dtype='<U36')}

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
+

prior-art/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2020 Silviu
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

prior-art/README.md

@@ -0,0 +1,50 @@
+# AI-Feynman
+
+This code is an improved implementation of AI Feynman: a Physics-Inspired Method for Symbolic Regression, Silviu-Marian Udrescu and Max Tegmark (2019) [[Science Advances](https://advances.sciencemag.org/content/6/16/eaay2631/tab-pdf)] and AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity, Udrescu S.M. et al. (2020) [[arXiv](https://arxiv.org/abs/2006.10782)].
+
+Please check [this Medium article](https://towardsdatascience.com/ai-feynman-2-0-learning-regression-equations-from-data-3232151bd929) for a more detailed eplanation of how to get the code running.
+
+In order to get started, run compile.sh to compile the fortran files used for the brute force code.
+
+ai_feynman_example.py contains an example of running the code on some examples (found in the example_data directory). The examples correspond to the equations I.8.14, I.10.7 and I.50.26 in Table 4 in the paper. More data files on which the code can be tested on can be found in the [Feynman Symbolic Regression Database](https://space.mit.edu/home/tegmark/aifeynman.html). 
+
+The main function of the code, called by the user, has the following parameters:
+
+* pathdir - path to the directory containing the data file
+* filename - the name of the file containing the data
+* BF_try_time - time limit for each brute force call (set by default to 60 seconds)
+* BF_ops_file_type - file containing the symbols to be used in the brute force code (set by default to "14ops.txt")
+* polyfit_deg - maximum degree of the polynomial tried by the polynomial fit routine (set be default to 4)
+* NN_epochs - number of epochs for the training (set by default to 4000)
+* vars_name - name of the variables appearing in the equation (inluding the name ofthe output variable). This should be passed as a list of strings, with the name of the variables appearing in the same order as they are in the file containing the data
+* test_percentage - percentage of the input data to be kept aside and used as the test set
+
+The data file to be analyzed should be a text file with each column containing the numerical values of each (dependent and independent) variable. The solution file will be saved in the directory called "results" under the name solution_{filename}. The solution file will contain several rows (corresponding to each point on the Pareto frontier), each row showing: 
+
+* the mean logarithm in based 2 of the error of the discovered equation applied to the input data (this can be though of as the average error in bits)
+* the cummulative logarithm in based 2 of the error of the discovered equation applied to the input data (this can be though of as the cummulative error in bits)
+* the complexity of the discovered equation (in bits)
+* the error of the discovered equation applied to the input data
+* the symbolic expression of the discovered equation
+
+If test_percentage is different than zero, one more number is added in the beginning of each row, showing the error of the discovered equation on the test set.
+
+ai_feynman_terminal_example.py allows calling the aiFeynman function from the command line. 
+(e.g. python ai_feynman_terminal_example.py --pathdir=../example_data/ --filename=example1.txt). Use python ai_feynman_terminal_example.py --help to display all the available parameters that can be passed to the function.
+
+# Citation
+
+If you compare with, build on, or use aspects of the AI Feynman work, please cite the following:
+
+```
+@article{udrescu2020ai,
+  title={AI Feynman: A physics-inspired method for symbolic regression},
+  author={Udrescu, Silviu-Marian and Tegmark, Max},
+  journal={Science Advances},
+  volume={6},
+  number={16},
+  pages={eaay2631},
+  year={2020},
+  publisher={American Association for the Advancement of Science}
+}
+```

prior-art/configs.cfg

@@ -0,0 +1,7 @@
+[Default]
+dataset_path = /home/aziz/lambda_lab/feynman_dataset/Feynman_without_units/
+operations_file = ./7ops.txt
+polynomial_degree = 3
+number_of_epochs = 500
+bruteforce_time = 10
+test_percentage = 0

prior-art/requirements.txt

@@ -0,0 +1,8 @@
+torch
+numpy
+matplotlib
+sympy
+pandas
+scipy
+sortedcontainers
+tabulate

prior-art/results/solution_I.18.14

@@ -0,0 +1,2 @@
+29.199776515901206 4.86788542219009 4867885.422190091 0.0 29.199776515901203 0
+3.869809801979284e-08 -24.623162098280734 -24623162.098280735 2.0 3.8698098019792816e-08 sin(x0)

prior-art/results/solution_I.39.11

@@ -0,0 +1,3 @@
+28.918403330075854 4.853915993947761 4853915.993947761 1.0 28.91840333007587 1
+27.015303509762383 4.755704985239505 4755704.985239505 3.0 27.01530350976237 1/x0
+3.8490145443977394e-08 -24.630935636380688 -24630935.636380687 6.754887502163468 3.8490145443977394e-08 0.000000000000+((x0-1))**(-1)

prior-art/results/solution_I.41.16

@@ -0,0 +1,9 @@
+26.52341213583394 4.729194479492894 4729194.479492894 0.0 26.52341213583394 0
+22.516263519038066 4.492895532882613 4492895.532882612 6.321928094887363 22.516263519038066 0.1*x0
+22.35215520829344 4.482342038529134 4482342.038529133 11.584962500721156 22.352155208293457 2*x0*exp(-3)
+21.765293633202937 4.442258566416072 4442258.566416072 16.720671786825555 21.739676488086655 asin(0.1*x0 - 0.01)
+21.11880720719005 4.400456448292662 4400456.448292661 22.83845916493269 21.11880720719004 asin(x0/pi**2 - 0.01)
+21.088909897780603 4.398412617913137 4398412.617913137 25.416665599935456 21.088909897780603 asin(0.1*x0 - 1.0*exp(1 - 2*pi))
+21.072971868776285 4.397321883081677 4397321.883081677 25.651484454403228 21.07297186877628 tan(x0*log(log(3))*sin(log(x0 + 1)/2 + 1))
+19.31414652315343 4.27158602231386 4271586.02231386 27.00162810065661 19.314146523153543 asin(0.1*x0 - log(2) + log(log(1 + 2*pi)))
+18.008155737063543 4.170578533475639 4170578.533475639 51.61544621886562 18.008155737063547 asin(0.1*x0 - 0.00673794699908547)

prior-art/results/solution_I.43.43

@@ -0,0 +1,3 @@
+28.91728816127691 4.853860358804314 4853860.358804314 1.0 28.917288161276893 1
+27.01635763399302 4.755761277412207 4755761.277412207 3.0 27.01635763399302 1/x0
+6.040924689221894e-08 -23.98065535814988 -23980655.35814988 6.754887502163468 6.040924689221891e-08 1/(x0 - 1)

prior-art/results/solution_I.6.2a

@@ -0,0 +1,5 @@
+25.5502918088963 4.675267863108713 4675267.863108713 0.0 25.5502918088963 0
+16.991234445196746 4.086718765730164 4086718.765730164 15.491853096329674 16.991234445196753 0.4*exp(-x0**2/2)
+16.991234445196746 4.086718765730163 4086718.765730163 19.101493570766486 16.991234445196746 0.4*exp(-0.5*x0**2)
+7.981114772196693e-08 -23.578834488341045 -23578834.488341045 24.779565475879124 7.981114772196695e-08 sqrt(2)*exp(-0.5*x0**2)/(2*sqrt(pi))
+7.981064652253219e-08 -23.578843548230925 -23578843.548230924 29.67970000576925 7.98106465225322e-08 sqrt(2)*exp(-1.0*x0**2.0)**0.5/(2*sqrt(pi))

prior-art/results/solution_II.10.9

@@ -0,0 +1,3 @@
+28.064858254334787 4.8106928676636205 4810692.86766362 0.0 28.064858254334787 0
+26.604553572874014 4.733601290041219 4733601.2900412185 3.0 26.604553572874018 1/x0
+2.9282883173699072e-08 -25.02536715166485 -25025367.15166485 6.754887502163468 2.928288317369907e-08 1/(x0 + 1)

prior-art/results/solution_II.11.27

@@ -0,0 +1,5 @@
+22.77297713033301 4.509251003419642 4509251.003419642 0.0 22.77297713033301 0.000000000000+(x0*x1)
+22.289843333097 4.437499871419407 4437499.8714194065 2.0 21.668086819815983 tan(0.000000000000+(x0*x1))
+20.993846246888534 4.391894599438973 4391894.599438973 9.0 20.993846246888534 -1.000000000000+exp(sin((x0*x1)))
+9.832766402766665e-08 -23.279571245496566 -23279571.245496567 9.339850002884624 9.820888240763341e-08 -0.000000000000+((x0*x1)/(((x0*x1)/((cos(pi)-1)-1))+1))
+5.129723825496818e-08 -24.54133213519386 -24541332.13519386 16.60964047443681 4.095650542223528e-08 -3.000000000000*((x0*x1)/((((x0*x1)-1)-1)-1))

prior-art/results/solution_II.11.28

@@ -0,0 +1,6 @@
+30.34191981060052 4.923240466352597 4923240.466352597 0.0 30.34191981060053 0
+27.24548632693982 4.767945337873829 4767945.337873829 1.0 27.245486326939798 1
+21.761975506917274 4.443737622278367 4443737.622278367 8.0 21.761975506917274 exp(x0*x1)
+21.438434911857406 4.422127682244653 4422127.682244654 25.97341254929059 21.438434911857414 exp(x0*x1*sqrt(log(log(4*pi))))
+21.345018318999095 4.415827495641937 4415827.4956419375 26.236446955124386 21.34501831899909 exp(x0*x1*log(1 + 2*pi)/2)
+21.312185843725555 4.413606662863952 4413606.662863952 54.4674384247163 21.31218584372553 exp(0.972955074527657*x0*x1)

prior-art/results/solution_II.11.3

@@ -0,0 +1,5 @@
+27.39589624365616 4.775887896382246 4775887.896382246 1.0 27.395896243656154 1
+26.81621912614581 4.745033937928682 4745033.937928681 6.754887502163468 26.816219126145803 exp(x0**3)
+24.541452353081525 4.617148724515219 4617148.724515218 7.339850002884624 24.54145235308151 x0**2 + 1.0
+2.1972568894135705e-07 -22.117793113787034 -22117793.113787033 9.339850002884624 2.1972568894135703e-07 (1 - x0**2)**(-1.0)
+2.1229088077984885e-07 -22.1674542644061 -22167454.2644061 12.0 2.1229088077984885e-07 1.000000000000+(x0/((x0)**(-1)-x0))

prior-art/results/solution_II.13.23

@@ -0,0 +1,4 @@
+26.416902333933454 4.723389399871772 4723389.399871772 0.0 26.416902333933454 x0
+25.3256838054465 4.662529317888433 4662529.317888433 6.754887502163468 25.3256838054465 (x0)*(exp((x1/x2)**3))
+2.8346308778966764e-07 -21.714376478821332 -21714376.478821333 17.194602975157967 2.9061721911169444e-07 0.000000000000+(x0/sqrt((((x1/x2)*(-(x1/x2)))+1)))
+2.8346308778966764e-07 -21.750335782765085 -21750335.782765083 22.67970000576925 2.8346308778966764e-07 x0/sqrt(-x1**2/x2**2 + 1)

prior-art/results/solution_II.13.34

@@ -0,0 +1,4 @@
+30.060513537598077 4.9097977505760255 4909797.750576026 0.0 30.060513537598073 0
+24.955584150305803 4.641290769148962 4641290.769148962 1.0 24.95558415030581 1
+9.497362397801731e-08 -23.25812421267564 -23258124.21267564 13.0 9.967975665930304e-08 (x0/(x0 - 1/x0))**0.5
+8.960019868428877e-08 -23.332283060485445 -23332283.060485445 15.169925001442312 9.468538109174675e-08 (x0**2/(x0**2 - 1))**0.5

prior-art/results/solution_II.15.4

@@ -0,0 +1,4 @@
+28.654323534574495 4.835600391337545 4835600.391337545 0.0 28.55359299998566 asin(-666.000000000000*sin(pi))
+6.099500490531156 4.834770389511966 4834770.389511965 4.754887502163468 28.537170459100448 acos(0.000000000000+cos(cos(x0)))
+23.84319730380234 4.575505804892507 4575505.804892507 9.0 23.84319730380235 -asin(cos(x0))
+5.745117505212346e-08 -24.05308835840271 -24053088.35840271 12.60964047443681 5.745117505212346e-08 -asin(sin(cos(x0)))

prior-art/results/solution_II.2.42

@@ -0,0 +1,3 @@
+27.647022713315675 4.789052220989075 4789052.220989075 0.0 27.647022713315675 (0)*(asin(-1.000000000000+cos((x1*sin(pi)))))
+26.851257452478915 4.746917746523634 4746917.746523634 3.0 26.85125745247892 x0 - 1
+6.772554077576063e-07 -20.493796656088517 -20493796.656088516 10.0 6.772554077576063e-07 0.000000000000+((x0-1)/x1)

prior-art/results/solution_II.21.32

@@ -0,0 +1,5 @@
+26.795605758262347 4.743924525773767 4743924.525773767 0.0 26.79560575826234 0
+25.743683182886485 4.686146571505307 4686146.571505307 7.754887502163468 25.74368318288649 0.222222222222222/x0
+25.725061304488595 4.6851026100465 4685102.6100465 13.868050853745158 25.725061304488584 0.224719101123595/x0
+22.390461830017472 4.484812380693082 4484812.380693082 14.584962500721156 22.390461830017507 0.0833333333333333*x0/(x0 - 1)
+8.244516940322618e-08 -23.531989792812396 -23531989.792812396 16.0 8.24451694032262e-08 x0/(4*pi*(x0 - 1))

prior-art/results/solution_II.34.11

@@ -0,0 +1,3 @@
+30.459303245022873 4.928811035516009 4928811.03551601 0.0 30.45930324502286 0
+28.882093735158097 4.85210342547258 4852103.42547258 2.0 28.882093735158097 log(x0)
+1.9438714084547796e-07 -22.294563879692067 -22294563.879692066 4.0 1.9438714084547778e-07 x0/2

prior-art/results/solution_II.34.29b

@@ -0,0 +1,2 @@
+34.10987896390003 4.694329806427947 4694329.806427947 0.0 25.890121066598777 4.88584612511867e-12
+1.99341380059708e-05 -15.614399252918906 -15614399.252918907 8.339850002884624 1.9934138005970806e-05 2*pi*x0*x1

prior-art/results/solution_II.35.18

@@ -0,0 +1,15 @@
+28.971870848354648 4.856580944027819 4856580.944027819 0.0 28.971870848354648 0
+28.345919874355094 4.82506918276535 4825069.18276535 1.0 28.345919874355086 1
+25.31926820387718 4.662163802422402 4662163.8024224015 12.60964047443681 25.31926820387718 asin(x0*exp(-x1))
+25.284729122683824 4.660194417728512 4660194.417728512 17.509775004326936 25.284729122683824 asin(x0/(exp(x1) + 1))
+24.731971558005636 4.628305346486737 4628305.346486737 25.0 24.731971558005647 asin(sin(x0/(exp(x1) + cos(x1))))
+24.73196283201209 4.605959594935965 4605959.594935965 44.743039599854946 24.35185206993002 asin(0.000000008775+sin((x0/(exp(x1)+cos(x1)))))
+25.691184139593112 4.6037976289129565 4603797.628912956 120.66823055601029 24.31538667221195 asin(sin(x0*(exp(x1) + cos(x1))**(-1.016446352005)) - 0.0207812879234552)
+26.990142594232438 4.557701468518747 4557701.468518747 199.95233207155312 23.550756020717426 ((x0*(exp(x1) + cos(x1))**(-0.994144976139069))**0.304960429668427 + 0.0318155214190483)**1.95718157291412
+23.433957987159474 4.550528740347353 4550528.740347353 247.49648716348534 23.43395798715948 0.19288704556216013*(0.06816739544620876*x0**2 + x0*x1*log(log(pi)) - 0.8872020666474496*x0 - 0.03642442424006058*x1**2 + 0.42211562108625245*x1 - 1)**2
+23.217882186039752 4.537164477984867 4537164.477984867 247.99825724455658 23.21788218603975 0.19288704556216013*(0.06816739544620876*x0**2 + 0.14942362847602014*x0*x1 - 0.8872020666474496*x0 - 0.03642442424006058*x1**2 + pi*x1*log(log(pi)) - 1)**2
+23.201907704524313 4.536171526218943 4536171.526218942 248.926632738039 23.20190770452431 0.19288704556216013*(0.06816739544620876*x0**2 + 0.14942362847602014*x0*x1 - sqrt(pi)*x0/2 - 0.03642442424006058*x1**2 + 0.42211562108625245*x1 - 1)**2
+23.086811364318827 4.528997023887213 4528997.023887212 249.49648716348534 23.086811364318827 0.19288704556216013*(-0.06816739544620876*x0**2 + x0*x1*(1 - log(pi)) + 0.8872020666474496*x0 + 0.03642442424006058*x1**2 - 0.42211562108625245*x1 + 1)**2
+23.050546317717448 4.526729039057007 4526729.039057007 256.06161725630693 23.050546317717437 (0.06816739544620876*x0**2 + 0.14942362847602014*x0*x1 - 0.8872020666474496*x0 - 0.03642442424006058*x1**2 + 0.42211562108625245*x1 - 1)**2*log(log(2) + log(pi))/pi
+22.95314348929238 4.520619842681248 4520619.842681249 261.0059797685666 22.953143489292366 (0.06816739544620876*x0**2 + 0.14942362847602014*x0*x1 - 0.8872020666474496*x0 - 0.03642442424006058*x1**2 + 0.42211562108625245*x1 - 1)**2*exp(-sqrt(pi))*log(pi)
+26.976561099205153 4.513216735201908 4513216.735201908 294.39710052772324 22.83566242154156 0.99098555264132*((x0*(0.983863716837616*exp(x1) + 1)**(-0.97031956911087))**0.293504804372787 + 0.0468885319131299)**1.9545316696167

prior-art/results/solution_II.37.1

@@ -0,0 +1,2 @@
+30.000000001343675 4.9068905956731355 4906890.5956731355 0.0 30.000000001343675 x0
+3.289228836945565e-07 -21.535747281929975 -21535747.281929974 3.0 3.2892288369455636e-07 x0 + 1

prior-art/results/solution_II.38.14

@@ -0,0 +1,3 @@
+27.064393896161725 4.7583241743201885 4758324.174320188 0.0 27.064393896161725 0
+25.603511150614796 4.678269763402974 4678269.763402974 5.0 25.603511150614807 0.5/x0
+8.47300466346212e-09 -26.814479190876447 -26814479.190876447 8.754887502163468 8.47300466346212e-09 0.500000000000*((x0+1))**(-1)

prior-art/results/solution_II.38.3

@@ -0,0 +1,4 @@
+29.20416081803319 4.86810202440029 4868102.02440029 0.0 29.204160818033195 x1
+3.132865116530134e-06 -18.283636592181303 -18283636.592181303 3.0 3.133840984176429e-06 -0.000000000000+((x0/x1))**(-1)
+3.132865116530134e-06 -18.284085912574707 -18284085.91257471 5.754887502163468 3.132865116530135e-06 x1/x0
+3.132865116530134e-06 -18.284085912574707 -18284085.91257471 10.0 3.132865116530134e-06 1.000000000000*(x1/x0)

prior-art/results/solution_II.6.11

@@ -0,0 +1,7 @@
+24.91363115902248 4.63886340443597 4638863.40443597 0.0 24.91363115902248 -8.915072558055073e-11
+24.91363113734539 4.638863403180696 4638863.403180696 10.754887502163468 24.913631137345384 -pi*exp(-exp(pi))
+20.52143209643687 4.359059508402251 4359059.508402252 12.584962500721156 20.521432096436865 asin(0.0833333333333333*cos(x0))
+12.27445776778958 3.617587388836739 3617587.388836739 17.60964047443681 12.274457767789572 asin(cos(x0)/(4*pi))
+11.607229189021997 3.53695171632929 3536951.71632929 28.236446955124386 11.607229189021998 asin(pi*sin(cos(x0)/pi**2)/4)
+11.577529861964756 3.5332555730823176 3533255.5730823176 69.39500893210585 11.577529861964752 asin(0.785437643527985*sin(cos(x0)/pi**2))
+11.577405736751007 3.533240105552333 3533240.105552333 100.97220344060344 11.577405736750997 asin(0.785437643527985*sin(0.1013211780033*cos(x0)))

prior-art/results/solution_II.6.15a

@@ -0,0 +1,9 @@
+28.495761948727154 4.832675464333084 4832675.464333084 0.0 28.495761948727154 0
+27.380642712447607 4.775084406610207 4775084.406610208 2.0 27.380642712447603 log(log(pi))
+27.054294104422542 4.757785694169851 4757785.694169851 6.321928094887362 27.05429410442252 4*exp(-3)
+25.30155195820018 4.661153975211263 4661153.975211264 7.339850002884624 25.30155195820018 0.333333333333333*x0*x2
+25.276008725588675 4.659696763845128 4659696.763845128 11.0 25.276008725588685 x0*x2/pi
+23.272252622734875 4.540538957128693 4540538.957128693 15.194602975157967 23.27225262273488 0.166666666666667*x2*(x0 + x1)
+21.63601106924245 4.435362635591665 4435362.635591665 29.0 21.63601106924245 x2*(x0 + x1)*exp(-sqrt(pi))
+21.588377090562414 4.4321828875112645 4432182.8875112645 29.558375050011747 21.588377090562417 2*x2*(x0 + x1)*log(pi)/(1 + 4*pi)
+21.18709874287426 4.405114140557128 4405114.140557128 55.550160087467745 21.18709874287426 0.168816319869*(x2*(x1+x0))

prior-art/results/solution_III.10.19

@@ -0,0 +1,5 @@
+29.667577204241404 4.890815209162709 4890815.209162709 0.0 29.667577204241386 x1
+29.664692327032483 4.843528401367819 4843528.401367819 5.754887502163468 28.710934848338415 1/(0.000000000000+(x0/x1))
+27.509606914525047 4.781863619982525 4781863.619982526 9.92481250360578 27.50960691452505 (x1**2 + 2)**0.5
+26.449356617926124 4.725160724016205 4725160.724016205 12.584962500721156 26.449356617926107 (x0 + x1**2 + 1)**0.5
+2.6825572107565247e-07 -21.82988772462039 -21829887.72462039 14.169925001442312 2.6825572107565247e-07 (x0**2 + x1**2 + 1)**0.5

prior-art/results/solution_III.12.43

@@ -0,0 +1,6 @@
+28.807896330713724 4.848392407795539 4848392.407795539 0.0 28.807896330713724 0
+27.53380966358232 4.783132334406948 4783132.334406948 4.0 27.533809663582314 0.750000000000000
+24.402752218093024 4.608971963473757 4608971.963473757 5.584962500721156 24.402752218093024 0.166666666666667*x0
+23.233134691280316 4.538111915347569 4538111.915347569 8.339850002884624 23.233134691280316 sin(0.166666666666667*x0)
+21.250729216232127 4.409440442894942 4409440.442894942 8.965784284662087 21.250729216232127 0.16*x0
+1.5306884709397791e-06 -19.31738787705657 -19317387.877056573 9.754887502163468 1.5306884709397796e-06 x0/(2*pi)

prior-art/results/solution_III.13.18

@@ -0,0 +1,2 @@
+36.57060766285735 4.550247640408479 4550247.640408479 0.0 23.429392474512944 8.91505505199440e-11
+0.0003960505387426254 -11.302027839797272 -11302027.839797273 6.321928094887362 0.00039605053874262544 4*pi*x0

prior-art/results/solution_III.14.14

@@ -0,0 +1,5 @@
+30.85970357081458 4.878360910171662 4878360.910171662 0.0 29.412569313163143 asin(0.000000000000*sqrt(exp(((x0+1))**(-1))))
+29.210279221326296 4.868404243855462 4868404.243855462 3.0 29.21027922132631 1.50000000000000
+28.35236907597377 4.825397384247454 4825397.384247454 4.584962500721156 28.35236907597379 2/x0
+27.81822835451622 4.797958637501465 4797958.637501465 6.0 27.818228354516222 1.5/x0
+7.268884829473797e-07 -20.391762617039728 -20391762.61703973 10.0 7.268884829473792e-07 exp(1/x0) - 1

prior-art/results/solution_III.15.27

@@ -0,0 +1,6 @@
+32.651368992645196 5.029071576628298 5029071.5766282985 0.0 32.65136899264519 0
+nan nan nan 5.754887502163468 nan acos(-x0)
+29.570858703869007 4.886104233165548 4886104.233165548 6.754887502163468 29.570858703869007 acos(1 - x0)
+28.17969118815505 4.8165838966073595 4816583.89660736 11.807354922057604 28.17969118815505 6*x0/x1
+6.643797122366517e-06 -17.199560550061296 -17199560.550061297 12.584962500721156 6.6437971223665156e-06 2*pi*x0/x1
+6.643797122366517e-06 -17.19959045660223 -17199590.45660223 58.15848945789539 6.643659400307738e-06 6.283185307179586*x0/x1

prior-art/results/solution_III.17.37

@@ -0,0 +1,4 @@
+31.225739937707367 4.964663853808277 4964663.853808277 0.0 31.22573993770736 0
+31.257799461003707 4.908431509806649 4908431.509806649 1.0 30.032059527985425 1/(1.000000000000*(x0-(x0+1)))
+31.243574989647236 4.868332410905311 4868332.410905311 14.60964047443681 29.208824854162387 1/(-1.000000000000*(x0-log(x0)))
+1.7431741319824306e-07 -22.45177997151243 -22451779.97151243 16.509775004326936 1.7431741319824306e-07 0.000000000000+(x0*((x1*cos(x2))+1))

prior-art/results/solution_III.19.51

@@ -0,0 +1,2 @@
+24.081314741250495 4.589842254381291 4589842.254381292 0.0 24.081314741250495 0
+5.647234375079581e-09 -27.39980834578584 -27399808.345785838 13.584962500721156 5.6472343750795805e-09 -0.125/x0**2

prior-art/results/solution_III.4.33

@@ -0,0 +1,7 @@
+26.290737189268583 4.716482690414682 4716482.690414682 0.0 26.290737189268583 x0
+22.658898217944188 4.502005807142314 4502005.8071423145 11.643856189774725 22.658898217944188 x0 - 0.07
+22.658481533490033 4.501979276544167 4501979.276544167 16.965784284662085 22.658481533490022 (x0*(x0 - 0.16))**0.5
+22.548084074885335 4.494932946856186 4494932.946856186 26.416665599935456 22.548084074885317 (x0*(x0 + log(log((2 + 4*pi)**(1/pi)))))**0.5
+19.11349123138255 4.256519417010686 4256519.417010686 29.80424344959478 19.11349123138255 (x0**2 - 0.16*x0 + 0.01)**0.5
+18.032683218363186 4.172542177060247 4172542.1770602474 35.64548429043134 18.03268321836319 (x0**2 + x0*log(log((2 + 4*pi)**(1/pi))) + 0.01)**0.5
+17.157005285884946 4.100725850740927 4100725.8507409273 66.69395636388344 17.157005285884953 (x0**2 - 0.1591549430918953*x0 + 0.01)**0.5

prior-art/results/solution_III.9.52

@@ -0,0 +1,4 @@
+32.84659416339962 4.913686532069664 4913686.532069664 0.0 30.141650893217463 tan(-666.000000000000*sin(pi))
+32.84182393413422 4.849625676726862 4849625.676726862 21.724214059872214 28.832532911636633 0.0731742799056452
+28.374965935007765 4.826546755293563 4826546.755293563 35.98584193700334 28.374965935007765 (12 - 12*cos(x1 - x2))/(x0*(x1 - x2)**2)
+0.0015229204525354075 -9.35894369787947 -9358943.697879469 38.43621560858933 0.0015229204525354064 -4*pi*(cos(x1 - x2) - 1)/(x0*(x1 - x2)**2)

results/solution_I.18.14

@@ -0,0 +1,2 @@
+29.199776515901206 4.86788542219009 4867885.422190091 0.0 29.199776515901203 0
+3.869809801979284e-08 -24.623162098280734 -24623162.098280735 2.0 3.8698098019792816e-08 sin(x0)

results/solution_I.39.11

@@ -0,0 +1,3 @@
+28.918403330075854 4.853915993947761 4853915.993947761 1.0 28.91840333007587 1
+27.015303509762383 4.755704985239505 4755704.985239505 3.0 27.01530350976237 1/x0
+3.8490145443977394e-08 -24.630935636380688 -24630935.636380687 6.754887502163468 3.8490145443977394e-08 0.000000000000+((x0-1))**(-1)

results/solution_I.41.16

@@ -0,0 +1,9 @@
+26.52341213583394 4.729194479492894 4729194.479492894 0.0 26.52341213583394 0
+22.516263519038066 4.492895532882613 4492895.532882612 6.321928094887363 22.516263519038066 0.1*x0
+22.35215520829344 4.482342038529134 4482342.038529133 11.584962500721156 22.352155208293457 2*x0*exp(-3)
+21.765293633202937 4.442258566416072 4442258.566416072 16.720671786825555 21.739676488086655 asin(0.1*x0 - 0.01)
+21.11880720719005 4.400456448292662 4400456.448292661 22.83845916493269 21.11880720719004 asin(x0/pi**2 - 0.01)
+21.088909897780603 4.398412617913137 4398412.617913137 25.416665599935456 21.088909897780603 asin(0.1*x0 - 1.0*exp(1 - 2*pi))
+21.072971868776285 4.397321883081677 4397321.883081677 25.651484454403228 21.07297186877628 tan(x0*log(log(3))*sin(log(x0 + 1)/2 + 1))
+19.31414652315343 4.27158602231386 4271586.02231386 27.00162810065661 19.314146523153543 asin(0.1*x0 - log(2) + log(log(1 + 2*pi)))
+18.008155737063543 4.170578533475639 4170578.533475639 51.61544621886562 18.008155737063547 asin(0.1*x0 - 0.00673794699908547)

results/solution_I.43.43

@@ -0,0 +1,3 @@
+28.91728816127691 4.853860358804314 4853860.358804314 1.0 28.917288161276893 1
+27.01635763399302 4.755761277412207 4755761.277412207 3.0 27.01635763399302 1/x0
+6.040924689221894e-08 -23.98065535814988 -23980655.35814988 6.754887502163468 6.040924689221891e-08 1/(x0 - 1)

results/solution_I.6.2a

@@ -0,0 +1,5 @@
+25.5502918088963 4.675267863108713 4675267.863108713 0.0 25.5502918088963 0
+16.991234445196746 4.086718765730164 4086718.765730164 15.491853096329674 16.991234445196753 0.4*exp(-x0**2/2)
+16.991234445196746 4.086718765730163 4086718.765730163 19.101493570766486 16.991234445196746 0.4*exp(-0.5*x0**2)
+7.981114772196693e-08 -23.578834488341045 -23578834.488341045 24.779565475879124 7.981114772196695e-08 sqrt(2)*exp(-0.5*x0**2)/(2*sqrt(pi))
+7.981064652253219e-08 -23.578843548230925 -23578843.548230924 29.67970000576925 7.98106465225322e-08 sqrt(2)*exp(-1.0*x0**2.0)**0.5/(2*sqrt(pi))

results/solution_II.10.9

@@ -0,0 +1,3 @@
+28.064858254334787 4.8106928676636205 4810692.86766362 0.0 28.064858254334787 0
+26.604553572874014 4.733601290041219 4733601.2900412185 3.0 26.604553572874018 1/x0
+2.9282883173699072e-08 -25.02536715166485 -25025367.15166485 6.754887502163468 2.928288317369907e-08 1/(x0 + 1)

results/solution_II.11.27

@@ -0,0 +1,5 @@
+22.77297713033301 4.509251003419642 4509251.003419642 0.0 22.77297713033301 0.000000000000+(x0*x1)
+22.289843333097 4.437499871419407 4437499.8714194065 2.0 21.668086819815983 tan(0.000000000000+(x0*x1))
+20.993846246888534 4.391894599438973 4391894.599438973 9.0 20.993846246888534 -1.000000000000+exp(sin((x0*x1)))
+9.832766402766665e-08 -23.279571245496566 -23279571.245496567 9.339850002884624 9.820888240763341e-08 -0.000000000000+((x0*x1)/(((x0*x1)/((cos(pi)-1)-1))+1))
+5.129723825496818e-08 -24.54133213519386 -24541332.13519386 16.60964047443681 4.095650542223528e-08 -3.000000000000*((x0*x1)/((((x0*x1)-1)-1)-1))

results/solution_II.11.28

@@ -0,0 +1,6 @@
+30.34191981060052 4.923240466352597 4923240.466352597 0.0 30.34191981060053 0
+27.24548632693982 4.767945337873829 4767945.337873829 1.0 27.245486326939798 1
+21.761975506917274 4.443737622278367 4443737.622278367 8.0 21.761975506917274 exp(x0*x1)
+21.438434911857406 4.422127682244653 4422127.682244654 25.97341254929059 21.438434911857414 exp(x0*x1*sqrt(log(log(4*pi))))
+21.345018318999095 4.415827495641937 4415827.4956419375 26.236446955124386 21.34501831899909 exp(x0*x1*log(1 + 2*pi)/2)
+21.312185843725555 4.413606662863952 4413606.662863952 54.4674384247163 21.31218584372553 exp(0.972955074527657*x0*x1)

results/solution_II.11.3

@@ -0,0 +1,5 @@
+27.39589624365616 4.775887896382246 4775887.896382246 1.0 27.395896243656154 1
+26.81621912614581 4.745033937928682 4745033.937928681 6.754887502163468 26.816219126145803 exp(x0**3)
+24.541452353081525 4.617148724515219 4617148.724515218 7.339850002884624 24.54145235308151 x0**2 + 1.0
+2.1972568894135705e-07 -22.117793113787034 -22117793.113787033 9.339850002884624 2.1972568894135703e-07 (1 - x0**2)**(-1.0)
+2.1229088077984885e-07 -22.1674542644061 -22167454.2644061 12.0 2.1229088077984885e-07 1.000000000000+(x0/((x0)**(-1)-x0))

results/solution_II.13.23

@@ -0,0 +1,4 @@
+26.416902333933454 4.723389399871772 4723389.399871772 0.0 26.416902333933454 x0
+25.3256838054465 4.662529317888433 4662529.317888433 6.754887502163468 25.3256838054465 (x0)*(exp((x1/x2)**3))
+2.8346308778966764e-07 -21.714376478821332 -21714376.478821333 17.194602975157967 2.9061721911169444e-07 0.000000000000+(x0/sqrt((((x1/x2)*(-(x1/x2)))+1)))
+2.8346308778966764e-07 -21.750335782765085 -21750335.782765083 22.67970000576925 2.8346308778966764e-07 x0/sqrt(-x1**2/x2**2 + 1)

results/solution_II.13.34

@@ -0,0 +1,4 @@
+30.060513537598077 4.9097977505760255 4909797.750576026 0.0 30.060513537598073 0
+24.955584150305803 4.641290769148962 4641290.769148962 1.0 24.95558415030581 1
+9.497362397801731e-08 -23.25812421267564 -23258124.21267564 13.0 9.967975665930304e-08 (x0/(x0 - 1/x0))**0.5
+8.960019868428877e-08 -23.332283060485445 -23332283.060485445 15.169925001442312 9.468538109174675e-08 (x0**2/(x0**2 - 1))**0.5

results/solution_II.15.4

@@ -0,0 +1,4 @@
+28.654323534574495 4.835600391337545 4835600.391337545 0.0 28.55359299998566 asin(-666.000000000000*sin(pi))
+6.099500490531156 4.834770389511966 4834770.389511965 4.754887502163468 28.537170459100448 acos(0.000000000000+cos(cos(x0)))
+23.84319730380234 4.575505804892507 4575505.804892507 9.0 23.84319730380235 -asin(cos(x0))
+5.745117505212346e-08 -24.05308835840271 -24053088.35840271 12.60964047443681 5.745117505212346e-08 -asin(sin(cos(x0)))