Skip to content
Snippets Groups Projects
Commit b8fb10aa authored by FROGE Ewen's avatar FROGE Ewen
Browse files

Add new file

parent cd6bb1c5
No related branches found
No related tags found
No related merge requests found
# Import Libraries
using Plots
using Base.Sort
using Statistics
using StatsBase
using LinearAlgebra
using NPZ
# Define Functions
# Linear Regression Function
function LinearRegression(x::Matrix{T}, y::Vector{T}, weights::Matrix{T}) where T
X = hcat(ones(T, size(x, 1)), x)
coef = inv(X' * weights * X) * X' * weights * y
return coef
end
# K-Nearest Neighbors Prediction with Linear Regression
function KNN_pred_LinReg(PastArray, AnalogsArray, p, k, verbosity=false)
dists = [norm(PastArray - AnalogsArray[u-p:u-1]) for u = p+1:size(AnalogsArray)[1]]
function neighbours(p, k)
sort(collect(enumerate(dists)), alg=Sort.PartialQuickSort(k), by = x -> x[2])[1:k]
end
x_k = neighbours(p, k)
a = median([x_k[i][2] for i = 1:k])
w = diagm([exp(-x_k[i][2]/a) for i = 1:k])
x_linalg = zeros(k, p)
for i = 1:k, j = 1:p
x_linalg[i, j] = AnalogsArray[x_k[i][1]+j-1]
end
y_linalg = zeros(k)
for i = 1:k
y_linalg[i] = AnalogsArray[x_k[i][1]+p]
end
coef = LinearRegression(x_linalg, y_linalg, w)
if verbosity
@show [AnalogsArray[x_k[i][1]] for i = 1:k]
@show w
@show(coef)
display(x_linalg)
end
return coef[1] + dot(PastArray, coef[2:end])
end
# Predictions Function
function Preds(Data, AnalogsData, tmin, tmax, p, k)
[KNN_pred_LinReg(Data[t-p:t-1], AnalogsData, p, k) for t in tmin:tmax]
end
# Parameters
Nreal = 10 # Number of realisations
T = 2^16 # Size of realisations
p = 10 # Analog size
k = 70 # Number of analogs to use
BeginCutOff = 50 # Number to ensure stationarity
Data_path = "/path/to/data.dat" # Path to data
Saving_path = "/path/to/savefile.npz" # Saving path
# Import Modane Data
f = open(Data_path)
y = Vector{Float32}(undef, stat(f).size ÷ sizeof(Float32))
read!(f, y)
close(f)
y .= bswap.(y) #Change data to correct byte format
y = reshape(y, (T, Int(size(y)[1] / T))) #Reshape to an array of realisations of size T
y = permutedims(y)
y = y[BeginCutOff:end, :]
# Normalize Data
Meanflow = mean(y)
StdV = std(y)
y = (y .- mean(y, dims=2)) ./ std(y, dims=2) #Normalization over each realisations
# Compute Innovation
Tmin = p + 1 # First point for computing innovation
Analogs = Float64.(y[1, :])
KNNpreds = zeros(Nreal, T - Tmin + 1)
@time begin
Threads.@threads for i in 1:Nreal
real = 2 * i + 2 # Pick every two realisations, starting at the third
println(Threads.threadid())
Realisation = Float64.(y[real, :])
KNNpreds[i, :] = Preds(Realisation, Analogs, Tmin, T, p, k)
end
end
# Save associated velocity and time
V = zeros(Nreal, T - Tmin + 1)
t = zeros(Nreal, T - Tmin + 1)
for i = 1:Nreal
real = 2 * i + 2
Realisation = Float64.(y[real, :])
V[i, :] = Realisation[Tmin:T]
t[i, :] = range(1 + (real - 1) * T, real * T)[Tmin:T]
end
Innov = V - KNNpreds # Compute innovations
# Save Data
npzwrite(Saving_path,
t=t,V=V,Innov=Innov,Tmin=Tmin,p=p,k=k,Meanflow=Meanflow,StdV=StdV,Vtot=y)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment