afdee1c/symbolic-regression
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Compare commits

base: afdee1c:7aa7ce80b99a5fbed30bfbaf924a57671e96bc12
Branches Tags
afdee1c:master
..
compare: afdee1c:60805f477fc20bfab4fe05f224bf4f23b9f8a16e
Branches Tags
afdee1c:master

No commits in common. "7aa7ce80b99a5fbed30bfbaf924a57671e96bc12" and "60805f477fc20bfab4fe05f224bf4f23b9f8a16e" have entirely different histories.
7aa7ce80b9 ... 60805f477f

119 changed files with 68 additions and 435772 deletions
Download patch file Download diff file Expand all files Collapse all files

130
Code/generate_claimed_results.py
View file

@ -1,43 +1,39 @@
import logging import logging
import argparse
import pathlib import pathlib
import os import configparser
from threading import active_count
from multiprocessing import Pool
from multiprocessing.pool import ThreadPool
from random import shuffle
from tabulate import tabulate from tabulate import tabulate
from pathlib import Path 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: class RunAll:
""" """
Run the solver on the whole dataset Run the solver on all the whole dataset
""" """
def __init__(self, *, cfg=_CFG): def __init__(self, *, cfg_path: Path):
logging.basicConfig(filename="output_no_units_parallel.log", level=logging.DEBUG) logging.basicConfig(filename="output.log", level=logging.DEBUG)
self.cfg = cfg self.config = configparser.ConfigParser()
self.config.read(cfg_path)
self.cfg = self.config["Default"]
self.print_results()
self.results = {} self.results = {}
self.run_solver()
def log_results(self):
pass
def print_results(self): def print_results(self):
table = [] table = [
for file, sol in self.results.items(): ["foo", 696000, 1989100000],
table.append(sol[-1]) ["bar", 6371, 5973.6],
["baz", 1737, 73.5],
["qux", 3390, 641.85],
]
print(tabulate( print(tabulate(
table, table,
headers=[ headers=[
@ -48,68 +44,26 @@ class RunAll:
], ],
) )
) )
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) def run_solver(self):
print("@"*120) path = Path(self.cfg["dataset_path"])
print("@"*120) for child in path.iterdir():
self.results[str(child).split("/")[-1]] = run_aifeynman(
self.print_results() pathdir="/home/aziz/lambda_lab/AI-Feynman/example_data/",#str(path.resolve()) + "/",
filename="example2.txt",#str(child).split("/")[-1],
BF_try_time=int(self.cfg["bruteforce_time"]),
def get_files(dirs, chunks=5): BF_ops_file_type=Path(self.cfg["operations_file"]),
dirs = list(path.iterdir()) polyfit_deg=int(self.cfg["polynomial_degree"]),
dirs = [file for file in dirs if not (str(file).endswith("test") or str(file).endswith("train"))] NN_epochs=int(self.cfg["number_of_epochs"]),
for i in range(0, len(dirs), chunks): vars_name=[],
yield dirs[i : i + chunks] test_percentage=int(self.cfg["test_percentage"]),
)
logging.info(self.results)
break
if __name__ == "__main__": if __name__ == "__main__":
cfg_path = pathlib.Path("/home/aziz/lambda_lab/AI-Feynman/configs.cfg")
#cfg_path = pathlib.Path("/home/aziz/lambda_lab/AI-Feynman/configs.cfg") if cfg_path.exists():
#if cfg_path.exists(): RunAll(cfg_path=cfg_path)
# RunAll(cfg_path=cfg_path) else:
#else: print(f"No such a file {cfg_path}")
# 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)

1185
notebook_1.ipynb
View file

File diff suppressed because it is too large Load diff

154
output_no_units_parallel.log
View file

@ -1,154 +0,0 @@
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')}

1
prior-art/Code/10ops.txt
View file

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

1
prior-art/Code/14ops.txt
View file

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

1
prior-art/Code/19ops.txt
View file

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

1
prior-art/Code/7ops.txt
View file

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

75
prior-art/Code/RPN_to_eq.py
View file

@ -1,75 +0,0 @@
# 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])

140
prior-art/Code/RPN_to_pytorch.py
View file

@ -1,140 +0,0 @@
# 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

122
prior-art/Code/S_NN_eval.py
View file

@ -1,122 +0,0 @@
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)

161
prior-art/Code/S_NN_train.py
View file

@ -1,161 +0,0 @@
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

124
prior-art/Code/S_add_bf_on_numbers_on_pareto.py
View file

@ -1,124 +0,0 @@
# 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)

203
prior-art/Code/S_add_snap_expr_on_pareto.py
View file

@ -1,203 +0,0 @@
# 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)

26
prior-art/Code/S_add_sym_on_pareto.py
View file

@ -1,26 +0,0 @@
# 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

42
prior-art/Code/S_brute_force.py
View file

@ -1,42 +0,0 @@
# 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

27
prior-art/Code/S_brute_force_number.py
View file

@ -1,27 +0,0 @@
# 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

182
prior-art/Code/S_change_output.py
View file

@ -1,182 +0,0 @@
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

36
prior-art/Code/S_combine_pareto.py
View file

@ -1,36 +0,0 @@
# 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

146
prior-art/Code/S_final_gd.py
View file

@ -1,146 +0,0 @@
# 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

18
prior-art/Code/S_get_number_DL.py
View file

@ -1,18 +0,0 @@
# 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

23
prior-art/Code/S_get_number_DL_snapped.py
View file

@ -1,23 +0,0 @@
# 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

42
prior-art/Code/S_get_symbolic_expr_error.py
View file

@ -1,42 +0,0 @@
# 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

89
prior-art/Code/S_polyfit.py
View file

@ -1,89 +0,0 @@
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)

55
prior-art/Code/S_polyfit_utils.py
View file

@ -1,55 +0,0 @@
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)

33
prior-art/Code/S_remove_input_neuron.py
View file

@ -1,33 +0,0 @@
# 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)

221
prior-art/Code/S_run_aifeynman.py
View file

@ -1,221 +0,0 @@
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

246
prior-art/Code/S_run_bf_polyfit.py
View file

@ -1,246 +0,0 @@
# 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

384
prior-art/Code/S_separability.py
View file

@ -1,384 +0,0 @@
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)

85
prior-art/Code/S_snap.py
View file

@ -1,85 +0,0 @@
# 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

562
prior-art/Code/S_symmetry.py
View file

@ -1,562 +0,0 @@
# 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)

4
prior-art/Code/ai_feynman_example.py
View file

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

19
prior-art/Code/ai_feynman_terminal_example.py
View file

@ -1,19 +0,0 @@
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)

66706
prior-art/Code/arity2templates.txt
View file

File diff suppressed because it is too large Load diff

23
prior-art/Code/brute_force_oneFile_mdl_v2.scr
View file

@ -1,23 +0,0 @@
#!/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

23
prior-art/Code/brute_force_oneFile_mdl_v3.scr
View file

@ -1,23 +0,0 @@
#!/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;

20
prior-art/Code/brute_force_oneFile_v1.scr
View file

@ -1,20 +0,0 @@
#!/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;

22
prior-art/Code/brute_force_oneFile_v2.scr
View file

@ -1,22 +0,0 @@
#!/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;

21
prior-art/Code/brute_force_oneFile_v3.scr
View file

@ -1,21 +0,0 @@
#!/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;

12
prior-art/Code/compile.sh
View file

@ -1,12 +0,0 @@
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

129
prior-art/Code/dimensionalAnalysis.py
View file

@ -1,129 +0,0 @@
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"]))

115
prior-art/Code/generate_claimed_results.py
View file

@ -1,115 +0,0 @@
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)

62
prior-art/Code/getPowers.py
View file

@ -1,62 +0,0 @@
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)

267
prior-art/Code/get_pareto.py
View file

@ -1,267 +0,0 @@
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()

157
prior-art/Code/symbolic_regress1.f
View file

@ -1,157 +0,0 @@
! 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'

292
prior-art/Code/symbolic_regress2.f
View file

@ -1,292 +0,0 @@
! 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

293
prior-art/Code/symbolic_regress3.f
View file

@ -1,293 +0,0 @@
! 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

220
prior-art/Code/symbolic_regress_mdl2.f
View file

@ -1,220 +0,0 @@
! 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"

244
prior-art/Code/symbolic_regress_mdl3.f
View file

@ -1,244 +0,0 @@
! 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"

380
prior-art/Code/tools.f
View file

@ -1,380 +0,0 @@
! 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

21
prior-art/LICENSE
View file

@ -1,21 +0,0 @@
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.

50
prior-art/README.md
View file

@ -1,50 +0,0 @@
# 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}
}
```

7
prior-art/configs.cfg
View file

@ -1,7 +0,0 @@
[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

120000
prior-art/example_data/example1.txt
View file

File diff suppressed because it is too large Load diff

120000
prior-art/example_data/example2.txt
View file

File diff suppressed because it is too large Load diff

120000
prior-art/example_data/example3.txt
View file

File diff suppressed because it is too large Load diff

1936
prior-art/notebook_1.ipynb
View file

File diff suppressed because it is too large Load diff

8
prior-art/requirements.txt
View file

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

2
prior-art/results/solution_I.18.14
View file

@ -1,2 +0,0 @@
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)

3
prior-art/results/solution_I.39.11
View file

@ -1,3 +0,0 @@
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)

9
prior-art/results/solution_I.41.16
View file

@ -1,9 +0,0 @@
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)

3
prior-art/results/solution_I.43.43
View file

@ -1,3 +0,0 @@
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)

5
prior-art/results/solution_I.6.2a
View file

@ -1,5 +0,0 @@
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))

3
prior-art/results/solution_II.10.9
View file

@ -1,3 +0,0 @@
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)

5
prior-art/results/solution_II.11.27
View file

@ -1,5 +0,0 @@
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))

6
prior-art/results/solution_II.11.28
View file

@ -1,6 +0,0 @@
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)

5
prior-art/results/solution_II.11.3
View file

@ -1,5 +0,0 @@
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))

4
prior-art/results/solution_II.13.23
View file

@ -1,4 +0,0 @@
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)

4
prior-art/results/solution_II.13.34
View file

@ -1,4 +0,0 @@
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

4
prior-art/results/solution_II.15.4
View file

@ -1,4 +0,0 @@
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)))

3
prior-art/results/solution_II.2.42
View file

@ -1,3 +0,0 @@
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)

5
prior-art/results/solution_II.21.32
View file

@ -1,5 +0,0 @@
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))

3
prior-art/results/solution_II.34.11
View file

@ -1,3 +0,0 @@
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

2
prior-art/results/solution_II.34.29b
View file

@ -1,2 +0,0 @@
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

15
prior-art/results/solution_II.35.18
View file

@ -1,15 +0,0 @@
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

2
prior-art/results/solution_II.37.1
View file

@ -1,2 +0,0 @@
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

3
prior-art/results/solution_II.38.14
View file

@ -1,3 +0,0 @@
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)

4
prior-art/results/solution_II.38.3
View file

@ -1,4 +0,0 @@
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)

7
prior-art/results/solution_II.6.11
View file

@ -1,7 +0,0 @@
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)))

9
prior-art/results/solution_II.6.15a
View file

@ -1,9 +0,0 @@
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))

5
prior-art/results/solution_III.10.19
View file

@ -1,5 +0,0 @@
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

6
prior-art/results/solution_III.12.43
View file

@ -1,6 +0,0 @@
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)

2
prior-art/results/solution_III.13.18
View file

@ -1,2 +0,0 @@
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

5
prior-art/results/solution_III.14.14
View file

@ -1,5 +0,0 @@
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

6
prior-art/results/solution_III.15.27
View file

@ -1,6 +0,0 @@
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

4
prior-art/results/solution_III.17.37
View file

@ -1,4 +0,0 @@
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))

2
prior-art/results/solution_III.19.51
View file

@ -1,2 +0,0 @@
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

7
prior-art/results/solution_III.4.33
View file

@ -1,7 +0,0 @@
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

4
prior-art/results/solution_III.9.52
View file

@ -1,4 +0,0 @@
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)

2
results/solution_I.18.14
View file

@ -1,2 +0,0 @@
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)

3
results/solution_I.39.11
View file

@ -1,3 +0,0 @@
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)

9
results/solution_I.41.16
View file

@ -1,9 +0,0 @@
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)

3
results/solution_I.43.43
View file

@ -1,3 +0,0 @@
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)

5
results/solution_I.6.2a
View file

@ -1,5 +0,0 @@
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))

3
results/solution_II.10.9
View file

@ -1,3 +0,0 @@
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)

5
results/solution_II.11.27
View file

@ -1,5 +0,0 @@
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))

6
results/solution_II.11.28
View file

@ -1,6 +0,0 @@
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)

5
results/solution_II.11.3
View file

@ -1,5 +0,0 @@
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))

4
results/solution_II.13.23
View file

@ -1,4 +0,0 @@
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)

4
results/solution_II.13.34
View file

@ -1,4 +0,0 @@
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

4
results/solution_II.15.4
View file

@ -1,4 +0,0 @@
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)))

Some files were not shown because too many files have changed in this diff Show more