Homec4science
Diffusion LC-Sampling-Theory (master) Browse

R12761/src/test.py
e8a43e18a742master

Last Change
Raw File

test.py
View Options

from typing import Dict, List, Set
import numpy as np
import matplotlib.pyplot as plt
import numpy.fft as fft
from copy import copy
from tqdm import tqdm
class continuous_set:
start = None
stop = None
freq = None
def __init__(self,start = None,stop = None, freq = None):
assert start<stop
self.start = start
self.stop = stop
self.freq = freq
def __repr__(self) -> str:
return '[start: {}, stop {}, freq {}]'.format(self.start,self.stop,self.freq)
def val_is_in(self,value):
inside = False
if value>=self.start and value <= self.stop:
inside = True
return inside
def is_intersection_empty(self,set_to_compare):
empty = self.start>=set_to_compare.stop or self.stop<=set_to_compare.start
return empty
class Edge:
set:continuous_set = None
id:int = None
def __init__(self,set:continuous_set,id:int) -> None:
self.set = set
self.id = id
def __repr__(self) -> str:
return 'id: {}, set: {}'.format(self.id,self.set)
def is_connected(self,edge):
return not self.set.is_intersection_empty(edge.set)
class Graph: #Stupid graph, we don't need all the graph math for this task
edges:Dict[int,Edge] = {}
vertexes:Dict[int,Set[int]] = {}
def __init__(self) -> None:
pass
def __repr__(self) -> str:
rep = 'Number of edges: {}, number of vertexes {}'.format(len(self.edges), len(self.vertexes))
return rep
def add_edge(self,edge:Edge) -> None:
if edge.id not in self.vertexes:
self.vertexes[edge.id] = set()
self.edges[edge.id] = edge
def connect_edges(self,edge_id_0:int,edge_id_1:int) -> None:
assert edge_id_0 in self.vertexes
assert edge_id_1 in self.vertexes
self.vertexes[edge_id_0] = self.vertexes[edge_id_0].union([edge_id_1])
self.vertexes[edge_id_1] = self.vertexes[edge_id_1].union([edge_id_0])
def add_edge_and_connect(self,edge_to_add:Edge) -> None:
self.add_edge(edge_to_add)
for e_id in self.edges:
if e_id == edge_to_add.id:
continue
if self.edges[e_id].is_connected(edge_to_add):
self.connect_edges(edge_to_add.id,e_id)
def get_connected_subgraphs(self) -> List[List[Edge]]:
components = []
visited_edges = []
for edge in self.vertexes:
if edge not in visited_edges:
#Simple flooding algo
edges_this_component = [edge]
neighbours = list(self.vertexes[edge])
for neighbour in neighbours:
if neighbour not in edges_this_component:
edges_this_component.append(neighbour)
neighbours.extend(list(self.vertexes[neighbour]))
visited_edges.extend(edges_this_component)
edges_this_component_all_info = [self.edges[id] for id in edges_this_component]
components.append(edges_this_component_all_info)
return components
def create_signal(coefs, ws, dt, t_span):
t = np.linspace(0,t_span,int(t_span/dt))
x = np.zeros(len(t))
for w,c in tqdm(zip(ws, coefs)):
x += c*np.sin(w*t)
return t,x
def get_max_freq(signal, dt, t_base, stopping_prec = 100, v = False):
signal -= np.mean(signal)
f_sig = fft.fft(signal)[0:len(signal)//2]
freq_sig = np.fft.fftfreq(t_base.size, d=dt)[0:len(signal)//2]
threshold = max(abs(f_sig))/stopping_prec
last_significant_coef_pos = 0
for i,f_c in enumerate(f_sig[-1::-1]):
magn_f_c = abs(f_c)
if magn_f_c > threshold:
last_significant_coef_pos = len(f_sig)-1-i
break
max_freq = freq_sig[last_significant_coef_pos]
if v:
print("\tf max: {} (value: {}, threshold: {} --> delta: {})"
.format(freq_sig[last_significant_coef_pos],abs(f_sig[last_significant_coef_pos]), threshold, abs(f_sig[last_significant_coef_pos])-threshold
)
)
plt.figure()
plt.plot(freq_sig,abs(f_sig))
plt.show()
if max_freq > 10000:
print("\tf max: {} (value: {}, threshold: {} --> delta: {})"
.format(freq_sig[last_significant_coef_pos],abs(f_sig[last_significant_coef_pos]), threshold, abs(f_sig[last_significant_coef_pos])-threshold
)
)
plt.figure()
plt.plot(freq_sig,abs(f_sig),"d-")
plt.figure()
plt.plot(t_base,signal,"d-")
plt.show()
return max_freq
def invert_fun(t,v):
if v[0] == max(v):
v = v[::-1]
t = t[::-1]
v_base = np.linspace(min(v),max(v),len(t)*2)
t_resamp = np.interp(v_base,v,t)
dv = abs(v_base[1]-v_base[0])
return v_base,t_resamp,dv
def injective_subdiv_signal(ts,vs):
sections = []
section = ([],[])
grad = np.gradient(vs)
g_prev = grad[0]
for t,v,g in zip(ts,vs,grad):
if (g_prev<0 and g<0) or (g_prev>=0 and g>=0):
section[0].append(t)
section[1].append(v)
else:
if len(section[0]) > 1:
sections.append(np.array(section))
section = ([t],[v])
g_prev = g
return sections
def sample_signal_at_freq(signal, support, freq):
if support[0] == max(support):
support = support[::-1]
signal = signal[::-1]
support_s = np.arange(support[0],support[-1],1/freq)
signal_s = np.interp(support_s,support,signal)
return support_s,signal_s
def main():
rng = np.random.default_rng(31415926514)
w_orig = rng.choice(100, size=10, replace=False)
coefs = rng.choice(3000, size=10)
BW = max(w_orig)
f_m = BW/(np.pi*2)
t_span = 100
dt_orig = 1/(150*2*f_m)
print("Creating original signal")
t_orig,x_orig = create_signal(coefs, w_orig,dt_orig,t_span)
print("Analyzing original signal")
f__max_numerical = get_max_freq(x_orig,dt_orig,t_orig, v = False)
print("Max freq: ",sorted(w_orig/(2*np.pi))[-1], "Max numerical freq: ",f__max_numerical)
print("Sampling original signal")
dt_sampl = 1/(2.1*f__max_numerical)
t_sampled,x_sampled = create_signal(coefs, w_orig,dt_sampl,t_span)
print("Computing Injective subdivision")
injective_sections = injective_subdiv_signal(t_orig,x_orig)
print("Analyzing Injective subdivision")
freqs_invert = []
id = 0
graph = Graph()
for section in tqdm(injective_sections):
v,t,dv = invert_fun(section[0],section[1])
f = get_max_freq(t,dv,v,stopping_prec = 2)
freqs_invert.append(f)
graph_edge = Edge(continuous_set(min(v),max(v),f),id)
graph.add_edge_and_connect(copy(graph_edge))
id += 1
print(graph)
components = graph.get_connected_subgraphs()
print("-------------------------------------------------------------------\n")
print("-------------------------------- CONNECTED COMPONENTS: {} ---------------------------\n".format(len(components)))
for c in components:
print(len(c))
print("-------------------------------------------------------------------")
''' OLD TILING, WRONG MATH BEHIND
#This should work, oth, I could always hang mslf :|
#Tiling of the values set so to have tiles with highest freqeuncies
val_edges.sort()
tiling_sets = []
tile = continuous_set(start = -1e6, freq = -1e6)
for val in val_edges:
for set in sets:
if set.get_start() == val and set.get_freq() > tile.get_freq(): #START OF A NEW SET WITH HIGHER FREQ
tile.set_stop(val)
tiling_sets.append(copy(tile))
tile = continuous_set(start = val, freq = set.get_freq())
break
elif set.get_stop() == val and set.get_freq() == tile.get_freq(): #END OF CURRENT SET
tile.set_stop(val)
tiling_sets.append(copy(tile))
max_f = 0
for set_IDGAF in sets: #FIND NEW SET WITH THIS VALUE INSIDE AND HIGHEST FREQEUNCEY
if set_IDGAF.val_is_in(val) and set_IDGAF.get_freq() > max_f and tile.get_freq() != set_IDGAF.get_freq():
max_f = set_IDGAF.get_freq()
tile = continuous_set(start = val, freq = max_f)
break
'''
print("Sampling injective sections")
print(min(freqs_invert),max(freqs_invert),2.1*max(freqs_invert))
print(np.average(freqs_invert),np.std(freqs_invert),np.median(freqs_invert))
plt.figure()
plt.hist(freqs_invert,10000)
plt.show(block = True)
lvl_sampling_freq = 2.1*max(freqs_invert)
max_f_pos = np.argmax(freqs_invert)
sections_s = []
pts = 0
for section in tqdm(injective_sections):
sampled = sample_signal_at_freq(section[0],section[1],lvl_sampling_freq)
sections_s.append(sampled)
pts += len(sampled[0])
plt.figure()
plt.plot(injective_sections[max_f_pos][0],injective_sections[max_f_pos][1], "-^")
plt.plot(sections_s[max_f_pos][1],sections_s[max_f_pos][0], "-d")
plt.show(block = False)
plt.figure()
for i in range(max_f_pos-2,max_f_pos+2):
plt.plot(sections_s[i][1],sections_s[i][0],"-d")
print(pts, len(t_sampled))
plt.plot(t_orig,x_orig)
plt.plot(t_sampled,x_sampled,"o-")
plt.show()
'''
plt.figure()
for i in range(5):
v,t = invert_fun(injective_sections[i][0],injective_sections[i][1])
plt.plot(v,t,"-o")
plt.plot(injective_sections[i][1],injective_sections[i][0],"-d")
plt.show()
'''
def test_graph():
set_1 = continuous_set(0,10,1)
set_2 = continuous_set(10,20,1)
set_3 = continuous_set(19,20,1)
set_4 = continuous_set(8,10,1)
set_5 = continuous_set(30,35,1)
edge_1 = Edge(set_1,1)
edge_2 = Edge(set_2,2)
edge_3 = Edge(set_3,3)
edge_4 = Edge(set_4,4)
edge_5 = Edge(set_5,5)
graph = Graph()
graph.add_edge_and_connect(edge_1)
graph.add_edge_and_connect(edge_2)
graph.add_edge_and_connect(edge_3)
graph.add_edge_and_connect(edge_4)
graph.add_edge_and_connect(edge_5)
components = graph.get_connected_subgraphs()
for c in components:
print("component --> ",c)
def main2():
import numpy as np
import numpy.fft as fft
import jax.numpy as jnp
from jax import lax,vmap
from functools import partial
def multi_sinc_interpolation(samples,sinc_freqs,ampls,signal_freq,time_base):
dt_samples = 1/signal_freq
tensor = jnp.zeros((len(samples),len(time_base),len(ampls)))
for i in range(len(ampls)):
base_sinc = lambda n:ampls[i]*jnp.sinc((time_base-n*dt_samples)*sinc_freqs[i])
tensor = tensor.at[:,:,i].set(vmap(base_sinc)(jnp.array(range(len(samples)))))
sincs_projections = jnp.tensordot(samples,tensor,axes = 1)# axes = 1 -> stupid nnumpy notation for saying "tensor dot product", why is this not the default? fuck me I guess
return jnp.sum(sincs_projections,axis = 1)
def multi_sinc_interpolation_2(samples,sinc_freqs,ampls,signal_freq,time_base):
dt_samples = 1/signal_freq
tensor = jnp.zeros((len(samples),len(time_base),len(ampls)))
for i in range(len(ampls)):
base_sinc = lambda n:ampls[i]*jnp.sinc((time_base-n*dt_samples)*sinc_freqs[i])
tensor = tensor.at[:,:,i].set(vmap(base_sinc)(jnp.array(range(len(samples)))))
sincs_projections = jnp.tensordot(samples,tensor,axes = 1)# axes = 1 -> stupid numpy notation for saying "tensor dot product", why is this not the default? fuck me I guess
return jnp.sum(sincs_projections,axis = 1)
s = jnp.array([1,2,3,4,5,6])
fs = jnp.array([1,2])
a = jnp.array([1,1])
f = 10
t = jnp.array([1,2,3,4])
r = multi_sinc_interpolation(s,fs,a,f,t)
print(r)
if __name__ == "__main__":
#main()
#test_graph()
main2()
c4science · Help