diff --git a/.gitignore b/.gitignore
index eb7e1f7..058be9c 100755
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@
 /Manifest.toml
 /docs/build/
 /.CondaPkg/
+.DS_Store
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 1f4c43b..162c2a5 100755
--- a/Project.toml
+++ b/Project.toml
@@ -4,12 +4,16 @@ authors = ["Alexander Lulkin", "Benedikt V. Ehinger"]
 version = "0.1.0"
 
 [deps]
+AppleAccelerate = "13e28ba4-7ad8-5781-acae-3021b1ed3924"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 Diagonalizations = "9cd687f3-b62d-43f3-8fd3-ffcd9e581047"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GaussianMixtures = "cc18c42c-b769-54ff-9e2a-b28141a64aae"
+IntelVectorMath = "c8ce9da6-5d36-5c03-b118-5a70151be7bc"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
diff --git a/src/Amica.jl b/src/Amica.jl
index 7401a1a..39d5ec8 100755
--- a/src/Amica.jl
+++ b/src/Amica.jl
@@ -3,8 +3,12 @@ module Amica
     using LinearAlgebra
 	using GaussianMixtures
 	using Distributions
+    using IntelVectorMath
 	using SpecialFunctions
     using ProgressMeter
+    using LoopVectorization
+    using AppleAccelerate
+
     #using ComponentArrays
     using Diagonalizations
     using LogExpFunctions
diff --git a/src/helper.jl b/src/helper.jl
index bc2ee9e..e323293 100755
--- a/src/helper.jl
+++ b/src/helper.jl
@@ -3,17 +3,18 @@ function removeMean!(input)
 	mn = mean(input,dims=2)
 	(n,N) = size(input)
 	for i in 1:n
-		input[i,:] = input[i,:] .- mn[i]
+		input[i,:] .= input[i,:] .- mn[i]
 	end
 	return mn
 end
 
 #Returns sphered data x. todo:replace with function from lib
-function sphering(x)
-	(n,N) = size(x)
-	Us,Ss,Vs = svd(x*x'/N)
+function sphering!(x)
+	(_,N) = size(x)
+	Us,Ss = svd(x*x'/N)
 	S = Us * diagm(vec(1 ./sqrt.(Ss))) * Us'
-    return x = S*x
+    x .= S*x
+	return S
 end
 
 function bene_sphering(data)
@@ -33,11 +34,48 @@ end
 
 #taken from amica_a.m
 #L = det(A) * mult p(s|θ)
-function logpfun(x,rho)
-	return  (.-abs.(x).^rho .- log(2) .- loggamma.(1+1/rho))
+function logpfun(rho, y_rho)
+	return .- y_rho .- log(2) .- loggamma(1 + 1 / rho)
 end
 
+
 #taken from amica_a.m
-function ffun(x,rho)
-	return rho .* sign.(x) .* abs.(x) .^(rho.-1)
-end
\ No newline at end of file
+@views function ffun(x::AbstractArray{T, 1}, rho::T) where {T<:Real}
+	return @inbounds copysign.(optimized_pow(abs.(x), rho - 1), x) .* rho
+end
+
+# optimized power function for different cpu architectures
+function optimized_pow(lhs::AbstractArray{T, 1}, rhs::T) where {T<:Real}
+	optimized_pow(lhs, repeat([rhs], length(lhs)))
+end
+
+function optimized_pow(lhs::AbstractArray{T, 1}, rhs::AbstractArray{T, 1}) where {T<:Real}
+	if Sys.iswindows() || Sys.islinux()
+		return IVM.pow(lhs, rhs)
+	elseif Sys.isapple()
+		return AppleAccelerate.pow(lhs, rhs)
+	else 
+		return lhs .^ rhs
+	end
+end
+
+function optimized_log(val)
+	if Sys.iswindows() || Sys.islinux()
+		return IVM.log(val)
+	elseif Sys.isapple()
+		return AppleAccelerate.log(val)
+	else 
+		return log.(val)
+	end
+end
+
+
+function optimized_exp(val) 
+	if Sys.iswindows() || Sys.islinux()
+		return IVM.exp(val)
+	elseif Sys.isapple()
+		return AppleAccelerate.exp(val)
+	else 
+		return exp.(val)
+	end
+end
diff --git a/src/likelihood.jl b/src/likelihood.jl
index ce35c92..d133f37 100755
--- a/src/likelihood.jl
+++ b/src/likelihood.jl
@@ -24,96 +24,67 @@ function calculate_LL!(myAmica::MultiModelAmica)
 end
 
 #Update loop for Lt and u (which is saved in z). Todo: Rename
-function loopiloop(myAmica::SingleModelAmica)
-	(n,N) = size(myAmica.source_signals)
-	m = myAmica.m
-	Q = zeros(m,N)
-	@debug (:prop,myAmica.learnedParameters.proportions[1,1],:shape,myAmica.learnedParameters.shape[1],:y,myAmica.y[1])
+function loopiloop!(myAmica::SingleModelAmica, y_rho)
+	(n,_) = size(myAmica.source_signals)
+	
 	for i in 1:n
-		Q = calculate_Q(myAmica,Q,i) # myAmica.Q
+		Q = calculate_Q(myAmica,i, y_rho) # myAmica.Q
 		calculate_u!(myAmica,Q,i) # myAmica.z
-		@debug (:z,myAmica.z[1,1:5,1])
 		calculate_Lt!(myAmica,Q) # myAmica.Q
 	end
-	@debug (:Q,Q[1,1])
 end
 
-function loopiloop(myAmica::MultiModelAmica)
+function loopiloop!(myAmica::MultiModelAmica, y_rho)
 	M = size(myAmica.models,1)
-	(n,N) = size(myAmica.models[1].source_signals)
-	m = myAmica.m
-	Q = zeros(m,N)
+	(n,_) = size(myAmica.models[1].source_signals)
 
 	for h in 1:M #run along models
 		Threads.@threads for i in 1:n #run along components
-			#uses SingleModel versions of functions
-			Q = calculate_Q(myAmica.models[h],Q,i)
-			calculate_u!(myAmica.models[h],Q,i)
-			calculate_Lt!(myAmica.models[h],Q)
+			Q = calculate_Q(myAmica.models[h], i, y_rho)
+			calculate_u!(myAmica.models[h], @view(Q[:, n, :]),i)
+			calculate_Lt!(myAmica.models[h], @view(Q[:, n, :]))
 		end
 	end
 end
 
-#Applies location and scale parameter to source signals (per generalized Gaussian)
-function calculate_y!(myAmica::SingleModelAmica)
-	n = size(myAmica.A,1)
-	Threads.@threads for i in 1:n
-		for j in 1:myAmica.m
-			myAmica.y[i,:,j] = sqrt(myAmica.learnedParameters.scale[j,i]) * (myAmica.source_signals[i,:] .- myAmica.learnedParameters.location[j,i])
-		end
+@views function calculate_Q(myAmica::SingleModelAmica, i, y_rho)
+	(n,N) = size(myAmica.source_signals)
+	m = myAmica.m
+	Q = zeros(m,N)
+	
+	for j in 1:m
+		Q[j,:] .= log(myAmica.learnedParameters.proportions[j,i]) + 0.5 * log(myAmica.learnedParameters.scale[j,i]) .+ logpfun(myAmica.learnedParameters.shape[j,i], y_rho[i, :, j])
 	end
-end
 
-function calculate_y!(myAmica::MultiModelAmica)
-	n = size(myAmica.models[1].A,1)
-	for h in 1:size(myAmica.models,1)
-		#=Threads.@threads=# for i in 1:n
-			for j in 1:myAmica.m
-				myAmica.models[h].y[i,:,j] = sqrt(myAmica.models[h].learnedParameters.scale[j,i]) * (myAmica.models[h].source_signals[i,:] .- myAmica.models[h].learnedParameters.location[j,i])
-			end
-		end
-	end
+	return Q
 end
 
 #calculates u but saves it into z. MultiModel also uses the SingleModel version
-function calculate_u!(myAmica::SingleModelAmica, Q, i)
+@views function calculate_u!(myAmica::SingleModelAmica, Q, i)
 	m = size(myAmica.learnedParameters.scale,1)
 	if m > 1
 		for j in 1:m
-			Qj = ones(m,1) .* Q[j,:]'
-			myAmica.z[i,:,j] = 1 ./ sum(exp.(Q-Qj),dims = 1)
+            myAmica.z[i, :, j] .= (1 ./ sum(optimized_exp(Q .- Q[j, :]'), dims=1))[:]
 		end
 	end
 end
 
-#no longer in use, multimodel also uses the singlemodel version (see loopiloop)
-function calculate_u!(myAmica::MultiModelAmica,Q,i,h)
-	m = myAmica.m
-	
-	if m > 1 
-		for j in 1:m
-			Qj = ones(m,1) .* Q[j,:]'
-			myAmica.models[h].z[i,:,j] = 1 ./ sum(exp.(Q-Qj),dims = 1)
+#Applies location and scale parameter to source signals (per generalized Gaussian)
+@views function calculate_y!(myAmica::SingleModelAmica)
+	#for j in 1:myAmica.m
+	#	myAmica.y[:,:,j] .= sqrt.(myAmica.learnedParameters.scale[j,:]) .* (myAmica.source_signals .- myAmica.learnedParameters.location[j,:])
+	#end
+	for j in 1:myAmica.m
+		for k = 1: size(myAmica.m,1)
+			for i = 1:size(myAmica.m,2)
+				myAmica.y[k,i,j] = sqrt(myAmica.learnedParameters.scale[j,k]) * (myAmica.source_signals[k,i] - myAmica.learnedParameters.location[j,k])
+			end
 		end
 	end
 end
 
-#Calculates densities for each generalized Gaussian j. Currently used my MultiModelAmica too
-function calculate_Q(myAmica::SingleModelAmica, Q, i)
-	m = size(myAmica.learnedParameters.scale, 1) #m = number of GGs, can't use myAmica.m in case this gets used my MultiModelAmica
-
-	for j in 1:m
-		Q[j,:] = log(myAmica.learnedParameters.proportions[j,i]) + 0.5*log(myAmica.learnedParameters.scale[j,i]) .+ logpfun(myAmica.y[i,:,j],myAmica.learnedParameters.shape[j,i])
-	end
-	return Q
-end
-
-#no longer in use, Q is no longer saved in Amica object but instead in loopiloop
-function calculate_Q!(myAmica::MultiModelAmica, i, h)
-	m = myAmica.m
-	for j in 1:m
-		myAmica.Q[j,:] = log(myAmica.models[h].learnedParameters.proportions[j,i]) + 0.5*log(myAmica.models[h].learnedParameters.scale[j,i]) .+ logpfun(myAmica.models[h].y[i,:,j],myAmica.models[h].learnedParameters.shape[j,i])
-	end
+function calculate_y!(myAmica::MultiModelAmica)
+		calculate_y!.(myAmica.models[h])
 end
 
 #Calculates Likelihood for each time sample and for each ICA model
diff --git a/src/main.jl b/src/main.jl
index 82951d5..2adf9aa 100755
--- a/src/main.jl
+++ b/src/main.jl
@@ -33,32 +33,50 @@ function amica!(myAmica::AbstractAmica,
 	(n, N) = size(data)
 	m = myAmica.m
 
+	println("m: $(m), n: $(n), N: $(N)")
+
 	#Prepares data by removing means and/or sphering
 	if remove_mean
 		removed_mean = removeMean!(data)
 	end
 	if do_sphering
-		data = sphering(data)
+		S = sphering!(data)
+		myAmica.S = S
+		LLdetS = logabsdet(S)[1]
+	else
+		myAmica.S = I
+		LLdetS = 0
 	end
 	
 	dLL = zeros(1, maxiter)
-	fp = zeros(m ,N)
 
 	#todo put them into object
 	lambda = zeros(n, 1)
-	kappa = zeros(n, 1)
-	sigma2 = zeros(n, 1)
-
     prog = ProgressUnknown("Minimizing"; showspeed=true)
 
+	y_rho = zeros(size(myAmica.y))
+
 	for iter in 1:maxiter
 		#E-step
 		update_sources!(myAmica, data)
 		calculate_ldet!(myAmica)
 		initialize_Lt!(myAmica)
+		myAmica.Lt .+= LLdetS
 		calculate_y!(myAmica)
-		loopiloop(myAmica) #Updates y and Lt. Todo: Rename
+		
+		# pre-calculate abs(y)^rho
+		for j in 1:m
+			for i in 1:n
+				y_rho[i,:,j] .= optimized_pow(abs.(myAmica.y[i,:,j]), myAmica.learnedParameters.shape[j,i])
+			end
+		end
+
+
+
+		loopiloop!(myAmica, y_rho) #Updates y and Lt. Todo: Rename
 		calculate_LL!(myAmica)
+
+
 		@debug (:LL,myAmica.LL)
 		#Calculate difference in loglikelihood between iterations
 		if iter > 1
@@ -78,7 +96,7 @@ function amica!(myAmica::AbstractAmica,
 		#M-step
 		try
 			#Updates parameters and mixing matrix
-			update_loop!(myAmica, fp, lambda, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
+			update_loop!(myAmica, lambda, y_rho, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
 		catch e
 			#Terminates if NaNs are detected in parameters
 			if isa(e,AmicaNaNException)
@@ -96,6 +114,8 @@ function amica!(myAmica::AbstractAmica,
 	end
 	#If parameters contain NaNs, the algorithm skips the A update and terminates by jumping here
     @label escape_from_NaN
+
+
 	#If means were removed, they are added back
 	if remove_mean
 		add_means_back!(myAmica, removed_mean)
diff --git a/src/newton.jl b/src/newton.jl
index f48053a..ea8d1f5 100644
--- a/src/newton.jl
+++ b/src/newton.jl
@@ -7,8 +7,8 @@ function newton_method!(myAmica::SingleModelAmica, iter, g, kappa, do_newton, ne
 
 	sigma2 = sum(myAmica.source_signals.^2,dims=2) / N #is probably called sigma2 cause always squared
 
-	dA = Matrix{Float64}(I, n, n) - g * myAmica.source_signals[:,:]' 
-	bflag = 0
+	dA = Matrix{Float64}(I, n, n) - g * myAmica.source_signals' 
+	bflag = false
 	B = zeros(n,n)
 	
 	for i in 1:n
@@ -20,28 +20,28 @@ function newton_method!(myAmica::SingleModelAmica, iter, g, kappa, do_newton, ne
 				if denom > 0
 					B[i,k] = (-kappa[k] * sigma2[i] * dA[i,k] + dA[k,i]) / denom
 				else
-					bflag = 1
+					bflag = true
 				end
 			end
 		end		
 	end
-	if (bflag == 0) && (do_newton == 1) && (iter > newt_start_iter)
-		myAmica.A[:,:] = myAmica.A[:,:] + lrate * myAmica.A[:,:] * B
+	if (bflag == false) && (do_newton == 1) && (iter > newt_start_iter)
+		myAmica.A += lrate * myAmica.A * B
 	else
-		myAmica.A[:,:] = myAmica.A[:,:] - lnatrate * myAmica.A[:,:] * dA
+		myAmica.A -= lnatrate * myAmica.A * dA
 	end
 end
 
-function newton_method!(myAmica::MultiModelAmica, h, iter, g, kappa, do_newton, newt_start_iter, lrate::LearningRate, lambda)
+@views function newton_method!(myAmica::MultiModelAmica, h, iter, g, kappa, do_newton, newt_start_iter, lrate::LearningRate, lambda)
 	
 	lnatrate = lrate.natural_rate
 	lrate = lrate.lrate
 	M = size(myAmica.models)
 	(n,N) = size(myAmica.models[1].source_signals)
 
-	sigma2 = myAmica.models[h].source_signals[:,:].^2 * myAmica.ica_weights_per_sample[h,:] /myAmica.ica_weights[h]
+	sigma2 = myAmica.models[h].source_signals.^2 * myAmica.ica_weights_per_sample[h,:] /myAmica.ica_weights[h]
 
-	dA = Matrix{Float64}(I, n, n) - g * myAmica.models[h].source_signals[:,:]' 
+	dA = Matrix{Float64}(I, n, n) - g * myAmica.models[h].source_signals' 
 	bflag = 0
 	B = zeros(n,n)
 	
@@ -60,8 +60,8 @@ function newton_method!(myAmica::MultiModelAmica, h, iter, g, kappa, do_newton,
 		end		
 	end
 	if (bflag == 0) && (do_newton == 1) && (iter > newt_start_iter)
-		myAmica.models[h].A[:,:] = myAmica.models[h].A[:,:] + lrate * myAmica.models[h].A[:,:] * B
+		myAmica.models[h].A = myAmica.models[h].A + lrate * myAmica.models[h].A * B
 	else
-		myAmica.models[h].A[:,:] = myAmica.models[h].A[:,:] - lnatrate * myAmica.models[h].A[:,:] * dA
+		myAmica.models[h].A = myAmica.models[h].A - lnatrate * myAmica.models[h].A * dA
 	end
 end
\ No newline at end of file
diff --git a/src/parameters.jl b/src/parameters.jl
index b75a49a..ce89883 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -6,9 +6,9 @@ function reparameterize!(myAmica::SingleModelAmica, data)
 
 	for i in 1:n
 		tau = norm(myAmica.A[:,i])
-		myAmica.A[:,i] = myAmica.A[:,i] / tau
-		mu[:,i] = mu[:,i] * tau
-		beta[:,i] = beta[:,i] / tau^2
+		myAmica.A[:,i] ./= tau
+		mu[:,i] .*= tau
+		beta[:,i] ./= tau^2
 	end
 
 	myAmica.learnedParameters.location .= mu
@@ -50,16 +50,16 @@ function reparameterize!(myAmica::MultiModelAmica, data)
 end
 
 #Calculates sum of z. Returns N if there is just one generalized Gaussian
-function calculate_sumz(myAmica::SingleModelAmica,i,j)
+@views function calculate_sumz(myAmica::SingleModelAmica)
 	if myAmica.m == 1
-		return size(myAmica.source_signals,2)
+		return size(myAmica.source_signals, 2)
 	else
-		return sum(myAmica.z[i,:,j])
+		return sum(myAmica.z, dims=2)
 	end
 end
 
-function calculate_sumz(myAmica::MultiModelAmica,i,j,h)
-		return sum(myAmica.models[h].z[i,:,j])
+@views function calculate_sumz(myAmica::MultiModelAmica,h)
+	return sum(myAmica.models[h].z, dims = 2)
 end
 
 #Calculates densities for each sample per ICA model and per Gaussian mixture
@@ -68,7 +68,7 @@ function calculate_z!(myAmica::MultiModelAmica, i,j,h)
 	if myAmica.m > 1
 		myAmica.models[h].z[i,:,j] .= myAmica.ica_weights_per_sample[h,:] .* myAmica.models[h].z[i,:,j]
 	elseif myAmica.m == 1
-		myAmica.models[h].z[i,:,j] = myAmica.ica_weights_per_sample[h,:]
+		myAmica.models[h].z[i,:,j] .= myAmica.ica_weights_per_sample[h,:]
 	end
 end
 
@@ -145,9 +145,9 @@ function initialize_shape_parameter(myAmica::MultiModelAmica, shapelrate::Learni
 end
 
 #Updates Gaussian mixture Parameters and mixing matrix. todo: rename since its not a loop for single model
-function update_loop!(myAmica::SingleModelAmica, fp, lambda, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
+function update_loop!(myAmica::SingleModelAmica, lambda, y_rho, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
 		#Update parameters
-		g, kappa, lambda = update_parameters!(myAmica, fp, lambda, shapelrate, update_shape)
+		g, kappa = update_parameters!(myAmica, lambda, y_rho, shapelrate, update_shape)
 		
 		#Checks for NaN in parameters before updating the mixing matrix
 		if any(isnan, kappa) || any(isnan, myAmica.source_signals) || any(isnan, lambda) || any(isnan, g) || any(isnan, myAmica.learnedParameters.proportions)
@@ -158,7 +158,7 @@ function update_loop!(myAmica::SingleModelAmica, fp, lambda, shapelrate, update_
 end
 
 #Updates Gaussian mixture Parameters and mixing matrix.
-function update_loop!(myAmica::MultiModelAmica, fp, lambda, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
+function update_loop!(myAmica::MultiModelAmica, fp, lambda, y_rho, shapelrate, update_shape, iter, do_newton, newt_start_iter, lrate)
 	(n,N) = size(myAmica.models[1].source_signals)
 	M = size(myAmica.models,1)
 
@@ -178,7 +178,7 @@ function update_loop!(myAmica::MultiModelAmica, fp, lambda, shapelrate, update_s
 			continue
 		end
 
-		g, kappa, lambda = update_parameters!(myAmica, h, fp, lambda, shapelrate, update_shape)#todo: remove return
+		g, kappa, lambda = update_parameters!(myAmica, h, fp, lambda, y_rho, shapelrate, update_shape)#todo: remove return
 
 		#Checks for NaN in parameters before updating the mixing matrix
 		if any(isnan, kappa) || any(isnan, myAmica.models[h].source_signals) || any(isnan, lambda) || any(isnan, g) || any(isnan, myAmica.models[h].learnedParameters.proportions)
@@ -191,7 +191,7 @@ end
 
 #Updates Gaussian mixture parameters. It also returns g, kappa and lamda which are needed to apply the newton method.
 #Todo: Save g, kappa, lambda in structure, remove return
-function update_parameters!(myAmica::SingleModelAmica, fp, lambda, lrate_rho::LearningRate, update_shape)
+@views function update_parameters!(myAmica::SingleModelAmica, lambda, y_rho, lrate_rho::LearningRate, upd_shape)
 	alpha = myAmica.learnedParameters.proportions
 	beta = myAmica.learnedParameters.scale
 	mu = myAmica.learnedParameters.location
@@ -202,50 +202,85 @@ function update_parameters!(myAmica::SingleModelAmica, fp, lambda, lrate_rho::Le
 	kappa = zeros(n,1)
 	zfp = zeros(m, N)
 
-	
-	#Threads.@threads 
+	# update 
+	# - myAmica.learnedParameters.proportions 
+	# - myAmica.z
+
+	# this has to run because calculate_u updates z
+
+	# depends on 
+	# - myAmica.z
+	# - myAmica.source_signals
+	sumz = calculate_sumz(myAmica)
+	if m > 0
+		sumz[sumz .< 0] .= 1
+		myAmica.z ./= sumz
+	end
 	for i in 1:n
 		for j in 1:m
-			sumz = 0
-			sumz = calculate_sumz(myAmica,i,j)
-			update_mixture_proportions!(sumz,myAmica,j,i)
-			if sumz > 0
-				if (m > 1)
-					myAmica.z[i,:,j] .= myAmica.z[i,:,j] / sumz
-				end
-			else
-				continue
-			end
+			update_mixture_proportions!(sumz[i, 1, j],myAmica,j,i)
+		end
+	end
 	
-			fp[j,:] = ffun(myAmica.y[i,:,j], rho[j,i])
-			zfp[j,:] = myAmica.z[i,:,j] .* fp[j,:]
-			g[i,:] = g[i,:] .+ alpha[j,i] .* sqrt(beta[j,i]) .*zfp[j,:]
+	# update 
+	# - fp
+	# - zfp
+	# - g
+	# - kp
+	# - kappa
+	# - lambda
+	# - mu
+
+	# depends on 
+	# - myAmica.y
+	# - myAmica.z
+	# - rho
+	# - alpha
+	# - beta
+	# - mu
+	for i in 1:n
+		for j in 1:m
+			fp = ffun(myAmica.y[i,:,j], rho[j,i])
+			zfp[j,:] .= myAmica.z[i,:,j] .* fp
+			g[i,:] .+= alpha[j,i] .* sqrt(beta[j,i]) .*zfp[j,:]
 	
-			kp = beta[j,i] .* sum(zfp[j,:].*fp[j,:])
+			kp = beta[j,i] .* sum(zfp[j,:].*fp)
 	
-			kappa[i] = kappa[i]  + alpha[j,i] * kp
+			kappa[i] += alpha[j,i] * kp
 	
-			lambda[i] = lambda[i] .+ alpha[j,i] .* (sum(myAmica.z[i,:,j].*(fp[j,:].*myAmica.y[i,:,j] .-1).^2) .+ mu[j,i]^2 .* kp)
-			mu[j,i] =  update_location(myAmica,rho[j,i],zfp[j,:],myAmica.y[i,:,j],mu[j,i],beta[j,i],kp)
-			
-			
+			lambda[i] += alpha[j,i] * (sum(myAmica.z[i,:,j].*(fp.*myAmica.y[i,:,j] .-1).^2) + mu[j,i]^2 * kp)
+			mu[j,i] = update_location(myAmica,rho[j,i],zfp[j,:],myAmica.y[i,:,j],mu[j,i],beta[j,i],kp)
+
+
 			beta[j,i] = update_scale(zfp[j,:],myAmica.y[i,:,j],beta[j,i],myAmica.z[i,:,j],rho[j,i])
 
-			if update_shape == 1
-				update_shape!(myAmica, rho, j, i, lrate_rho)
+		end
+	end
+
+	# update rho
+	# depends on rho, zfp, myAmica.y, mu, beta
+	if upd_shape == 1
+		dr = sum(myAmica.z .* optimized_log(y_rho) .* y_rho, dims=2)
+
+		for i in 1:n
+			for j in 1:m
+				rho[j, i] = update_shape(rho[j, i], dr[i, 1, j], lrate_rho)
 			end
 		end
 	end
+
+
 	myAmica.learnedParameters.proportions = alpha
 	myAmica.learnedParameters.scale = beta
 	myAmica.learnedParameters.location = mu
 	myAmica.learnedParameters.shape = rho
-	return g, kappa, lambda
+
+	return g, kappa
 end
 
 #Updates Gaussian mixture parameters. It also returns g, kappa and lamda which are needed to apply the newton method.
 #Todo: Save g, kappa, lambda in structure, remove return
-function update_parameters!(myAmica::MultiModelAmica, h, fp, lambda, lrate_rho::LearningRate, update_shape)
+@views function update_parameters!(myAmica::MultiModelAmica, h, fp, y_rho, lambda, lrate_rho::LearningRate, upd_shape)
 	alpha = myAmica.models[h].learnedParameters.proportions #todo: move into loop and add h
 	beta = myAmica.models[h].learnedParameters.scale
 	mu = myAmica.models[h].learnedParameters.location
@@ -258,34 +293,42 @@ function update_parameters!(myAmica::MultiModelAmica, h, fp, lambda, lrate_rho::
 	zfp = zeros(m, N)
 
 	
+	sumz = calculate_sumz(myAmica,h)
+	if m > 0
+		sumz[sumz .< 0] .= 1
+		myAmica.models[h].z ./= sumz
+	end
+
+
 	#=Threads.@threads=# for i in 1:n
 		for j in 1:m
 			sumz = 0
 			calculate_z!(myAmica, i, j, h)
-			sumz = calculate_sumz(myAmica,i,j,h)
-			update_mixture_proportions!(sumz,myAmica,j,i,h)
-			if sumz > 0
-				myAmica.models[h].z[i,:,j] .= myAmica.models[h].z[i,:,j] / sumz
+			update_mixture_proportions!(sumz[i, 1, j],myAmica,j,i,h)
+			if sumz[i, 1, j] > 0
+				myAmica.models[h].z[i,:,j] ./= sumz[i, 1, j]
 			else
 				continue
 			end
 	
-			fp[j,:] = ffun(myAmica.models[h].y[i,:,j], rho[j,i])
-			zfp[j,:] = myAmica.models[h].z[i,:,j] .* fp[j,:]
-			g[i,:] = g[i,:] .+ alpha[j,i] .* sqrt(beta[j,i]) .*zfp[j,:]
+			fp[j,:] .= ffun(myAmica.models[h].y[i,:,j], rho[j,i])
+			zfp[j,:] .= myAmica.models[h].z[i,:,j] .* fp[j,:]
+			g[i,:] .+= alpha[j,i] .* sqrt(beta[j,i]) .*zfp[j,:]
 	
 			kp = beta[j,i] .* sum(zfp[j,:].*fp[j,:])
 	
-			kappa[i] = kappa[i]  + alpha[j,i] * kp
+			kappa[i] += alpha[j,i] * kp
 	
-			lambda[i] = lambda[i] .+ alpha[j,i] .* (sum(myAmica.models[h].z[i,:,j].*(fp[j,:].*myAmica.models[h].y[i,:,j] .-1).^2) .+ mu[j,i]^2 .* kp)
-			mu[j,i] =  update_location(myAmica,rho[j,i],zfp[j,:],myAmica.models[h].y[i,:,j],mu[j,i],beta[j,i],kp)
+			lambda[i] += alpha[j,i] .* (sum(myAmica.models[h].z[i,:,j].*(fp[j,:].*myAmica.models[h].y[i,:,j] .-1).^2) .+ mu[j,i]^2 .* kp)
+			mu[j,i] = update_location(myAmica,rho[j,i],zfp[j,:],myAmica.models[h].y[i,:,j],mu[j,i],beta[j,i],kp)
 			
 			
 			beta[j,i] = update_scale(zfp[j,:],myAmica.models[h].y[i,:,j],beta[j,i],myAmica.models[h].z[i,:,j],rho[j,i])
 
-			if update_shape == 1
-				update_shape!(myAmica.models[h], rho, j, i, lrate_rho) #uses SingleModel version on each model
+			if upd_shape == 1
+				log_y_rho = optimized_log(y_rho)
+				dr = sum(myAmica.models[h].z .* log_y_rho .* y_rho, dims=2)
+				rho[j, i] = update_shape(rho[j, i], dr[i, 1, j], lrate_rho)
 			end
 		end
 	end
@@ -293,49 +336,24 @@ function update_parameters!(myAmica::MultiModelAmica, h, fp, lambda, lrate_rho::
 	myAmica.models[h].learnedParameters.scale = beta
 	myAmica.models[h].learnedParameters.location = mu
 	myAmica.models[h].learnedParameters.shape = rho
-	return g, kappa, lambda
+	return g, kappa
 end
 
 #Updates the Gaussian mixture shape parameter
-function update_shape!(myAmica::SingleModelAmica, rho, j, i, lrate_rho::LearningRate)
-	rhomin, rhomax, shapelrate = lrate_rho.minimum, lrate_rho.maximum, lrate_rho.lrate
-	ytmp = abs.(myAmica.y[i,:,j]).^rho[j,i]
-	dr = sum(myAmica.z[i,:,j].*log.(ytmp).*ytmp)
-
-	if rho[j,i] > 2
-		dr2 = digamma(1+1/rho[j,i]) / rho[j,i] - dr
+@views function update_shape(rho, dr, lrate_rho::LearningRate)
+	if rho > 2
+		dr2 = digamma(1+1/rho) / rho - dr
 		if ~isnan(dr2)
-			rho[j,i] = rho[j,i] + 0.5 * dr2
+			rho += 0.5 * dr2
 		end
 	else
-		dr2 = 1 - rho[j,i] * dr / digamma(1+1/rho[j,i])
+		dr2 = 1 - rho * dr / digamma(1+1/rho)
 		if ~isnan(dr2)
-			rho[j,i] = rho[j,i] + shapelrate *dr2
+			rho += lrate_rho.lrate * dr2
 		end
 	end
-	rho[j,i] = min(rhomax, rho[j,i])
-	rho[j,i] = max(rhomin, rho[j,i])
-end
-
-#Updates the Gaussian mixture shape parameter. MultiModel version no longer in use
-function update_shape!(myAmica::MultiModelAmica, rho, j, i, h, lrate_rho::LearningRate)
-	rhomin, rhomax, shapelrate = lrate_rho.minimum, lrate_rho.maximum, lrate_rho.lrate
-	ytmp = abs.(myAmica.models[h].y[i,:,j]).^rho[j,i]
-	dr = sum(myAmica.models[h].z[i,:,j].*log.(ytmp).*ytmp)
 
-	if rho[j,i] > 2
-		dr2 = digamma(1+1/rho[j,i]) / rho[j,i] - dr
-		if ~isnan(dr2)
-			rho[j,i] = rho[j,i] + 0.5 * dr2
-		end
-	else
-		dr2 = 1 - rho[j,i] * dr / digamma(1+1/rho[j,i])
-		if ~isnan(dr2)
-			rho[j,i] = rho[j,i] + shapelrate *dr2
-		end
-	end
-	rho[j,i] = min(rhomax, rho[j,i])
-	rho[j,i] = max(rhomin, rho[j,i])
+	return clamp(rho, lrate_rho.minimum, lrate_rho.maximum)
 end
 
 #Calculates determinant of mixing Matrix A (with log). first log-likelihood part of L = |A| * p(sources)
@@ -353,8 +371,7 @@ end
 
 #Updates source singal estimations by unmixing the data
 function update_sources!(myAmica::SingleModelAmica, data)
-	b = pinv(myAmica.A) * data
-	myAmica.source_signals = b
+	myAmica.source_signals = pinv(myAmica.A) * data
 end
 
 function update_sources!(myAmica::MultiModelAmica, data)
diff --git a/src/types.jl b/src/types.jl
index 2459b4d..be17ff5 100755
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,34 +1,38 @@
-mutable struct GGParameters
-	proportions::AbstractArray{Float64} #source density mixture proportions
-	scale::AbstractArray{Float64} #source density inverse scale parameter
-	location::AbstractArray{Float64} #source density location parameter
-	shape::AbstractArray{Float64} #source density shape paramters
+mutable struct GGParameters{T}
+    proportions::Array{T, 2} #source density mixture proportions
+    scale::Array{T, 2} #source density inverse scale parameter
+    location::Array{T, 2} #source density location parameter
+    shape::Array{T, 2} #source density shape paramters
 end
 
+
 abstract type AbstractAmica end
 
-mutable struct SingleModelAmica <:AbstractAmica
-	source_signals					   #Unmixed signals
-	learnedParameters::GGParameters	   #Parameters of the Gaussian mixtures
-	m::Union{Integer, Nothing} 		   #Number of gaussians
-    A::AbstractArray 				   #Mixing matrix
-	z::AbstractArray				   #Densities for each sample per Gaussian (normalized)
-	y::AbstractArray				   #Source signals (scaled and shifted with scale and location parameter)
-	centers::AbstractArray 			   #Model centers
-	Lt::AbstractVector 				   #Log likelihood of time point for each model ( M x N )
-	LL::Union{AbstractVector, Nothing} #Log-Likelihood
-	ldet::Float64					   #log determinant of A
-	maxiter::Union{Int, Nothing} 	   #maximum number of iterations, can be nothing because it's not needed for multimodel
+mutable struct SingleModelAmica{T} <:AbstractAmica
+    source_signals::Array{T,2}
+    learnedParameters::GGParameters{T}
+	m::Int 		   #Number of gaussians
+    A::Array{T,2} # unmixing matrices for each model
+	S::Array{T,2} # sphering matrix
+    z::Array{T,3}
+    y::Array{T,3}
+    centers::Array{T,1} #model centers
+    Lt::Array{T,1} #log likelihood of time point for each model ( M x N )
+    LL::Array{T,1} #log likelihood over iterations todo: change to tuple 
+    ldet::T
+    maxiter::Int
 end
 
-mutable struct MultiModelAmica <:AbstractAmica
-	models::Array{SingleModelAmica} #Array of SingleModelAmicas
+
+
+mutable struct MultiModelAmica{T} <:AbstractAmica
+	models::Array{SingleModelAmica{T}} #Array of SingleModelAmicas
 	normalized_ica_weights 			#Model weights (normalized)
 	ica_weights_per_sample 			#Model weight for each sample
 	ica_weights						#Model weight for all samples
 	maxiter::Int					#Number of iterations
 	m::Int 							#Number of Gaussians
-	LL::AbstractVector				#Log-Likelihood
+	LL::Array{T,1}				#Log-Likelihood
 end
 
 #Structure for Learning Rate type with initial value, minumum, maximum etc. Used for learning rate and shape lrate
@@ -85,7 +89,7 @@ function SingleModelAmica(data::AbstractArray{T}; m=3, maxiter=500, A=nothing, l
 	ldet = 0.0
 	source_signals = zeros(n,N)
 
-	return SingleModelAmica(source_signals,GGParameters(proportions,scale,location,shape),m,A,z,y,#=Q,=#centers,Lt,LL,ldet,maxiter)
+	return SingleModelAmica{T}(source_signals,GGParameters{T}(proportions,scale,location,shape),m,A,I(size(A,1)), z,y,#=Q,=#centers,Lt,LL,ldet,maxiter)
 end
 
 #Data type for AMICA with multiple ICA models
diff --git a/test/compare_amica_implementations.jl b/test/compare_amica_implementations.jl
index 4b38a78..5f4c468 100644
--- a/test/compare_amica_implementations.jl
+++ b/test/compare_amica_implementations.jl
@@ -3,15 +3,17 @@ using Amica
 using SignalAnalysis
 using LinearAlgebra
 using Revise 
+#---
 using CairoMakie
 #---
-includet("src/simulate_data.jl")
-includet("src/fortran_tools.jl")
+includet("../../src/simulate_data.jl")
+includet("../../src/fortran_tools.jl")
 n_chan=3
-n_time=5000
+n_time=10000
 x,A,s = simulate_data(;T=50,n_chan,n_time,type=:gg6)
 i_scale,i_location,i_A = init_params(;n_chan,n_time)
 
+#---
 f  = Figure()
 series(f[1,1],s,axis=(;title="source"))
 #xlims!(1,20)
@@ -20,36 +22,49 @@ f
 
 #---
 
-maxiter = 500
+maxiter = 200
 
 # matlab run
 mat"""
+tic
 [mA,mc,mLL,mLt,mgm,malpha,mmu,mbeta,mrho] = amica_a($x,1,3,$maxiter,$i_location,$i_scale,$i_A,0);
+toc
 """
 mat"""
+tic
 [mAopt,mWopt,mSopt,mkhindsopt,mcopt,mLLopt,mLtopt,mgmopt,malphaopt,mmuopt,mbetaopt,mrhoopt] = amica_optimized($x,1,3,$maxiter,1,1,$i_location,$i_scale,$i_A);
+toc
 """
 mA = @mget(mA); # get matlab var
 mAopt = @mget(mAopt); # get matlab opt var
 
 # Fortran setup + ran
-fortran_setup(x;max_threads=1,max_iter=maxiter)
-run(`/scratch/projects/fapra_amica/fortran/amica15test julia.param`)
+fortran_setup(Float32.(x);max_threads=1,max_iter=maxiter)
+@time run(`/scratch/projects/fapra_amica/fortran/amica15test julia.param`)
 fA = reshape(reinterpret(Float64,(read("amicaout/A"))),n_chan,n_chan)
 
-# Julia run
-am = SingleModelAmica(x;maxiter=maxiter,A=i_A,location=i_location,scale=i_scale)
-fit!(am,x)
+fortran_setup(Float64.(x);max_threads=1,max_iter=maxiter,dble_data=1)
+@time run(`/scratch/projects/fapra_amica/fortran/amica15test julia.param`)
+
+#--- Julia run
+am32 = SingleModelAmica(Float32.(x);maxiter=maxiter,A=deepcopy(i_A),location=i_location,scale=i_scale)
+@time fit!(am32,Float32.(x))
+
+am64 = SingleModelAmica(Float64.(x);maxiter=maxiter,A=deepcopy(i_A),location=i_location,scale=i_scale)
+@time fit!(am64,Float64.(x))
+#am = SingleModelAmica(Float16.(x);maxiter=maxiter,A=i_A,location=i_location,scale=i_scale)
+#@time fit!(am,x)
 #vcat(@mget(mLL),am.LL')
 
 #---
 f2 = Figure(size=(800,800))
-series(f2[1,1],inv(am.A)*x, axis=(;title="unmixed julia"))
-series(f2[2,1],inv(mA)*x, axis=(;title="unmixed matlab"))
-series(f2[3,1],inv(mAopt)*x, axis=(;title="unmixed matlab_optimizd"))
-series(f2[4,1],.-inv(fA')*x, axis=(;title="unmixed fortran"))
-series(f2[5,1],s,axis=(;title="original source"))
-series(f2[6,1],x, axis=(;title="original mixed"))
+series(f2[1,1],inv(am64.A)*x, axis=(;title="unmixed julia64"))
+series(f2[2,1],inv(am32.A)*x, axis=(;title="unmixed julia32"))
+series(f2[3,1],inv(mA)*x, axis=(;title="unmixed matlab"))
+series(f2[4,1],inv(mAopt)*x, axis=(;title="unmixed matlab_optimizd"))
+series(f2[5,1],.-inv(fA')*x, axis=(;title="unmixed fortran"))
+series(f2[6,1],s,axis=(;title="original source"))
+series(f2[7,1],x, axis=(;title="original mixed"))
 hidedecorations!.(f2.content)
 
 linkxaxes!(f2.content...)
@@ -58,18 +73,17 @@ f2
 
 
 #--- compare LLs
-LL = am.LL
 mLL = @mget mLL
 mLLopt = @mget mLLopt
 fLL = reinterpret(Float64,(read("amicaout/LL")))
-f = Figure()
+f = Figure(size=(1024,1024))
 ax = f[1,1] = Axis(f)
-labels = ["julia","matlab", "matlab opt","fortran"]
-for (ix,d) = enumerate([LL,mLL[1,:],mLLopt[1,:],fLL])
+labels = ["julia64","julia32","matlab", "matlab opt","fortran"]
+for (ix,d) = enumerate([am64.LL, am32.LL,mLL[1,:],mLLopt[1,:],fLL])
     lines!(ax,d;label=labels[ix])
 end
 axislegend(ax)
-ylims!(ax,-1.3,-1)
+ylims!(ax,-2.3,-0.5)
 f
 
 #--- compare A
diff --git a/test/profile.jl b/test/profile.jl
new file mode 100644
index 0000000..14b582c
--- /dev/null
+++ b/test/profile.jl
@@ -0,0 +1,28 @@
+using Amica
+using MAT
+using Profile
+using Test
+
+file = matopen("test/testdata/eeg_data.mat")
+#file = matopen("test/testdata/pink_sinus_data.mat")
+
+x = read(file, "x")
+
+beta_init = read(file, "beta_init")
+mu_init = read(file, "mu_init")
+A_init = read(file, "A_init")
+
+close(file)
+# before: ca 1.2s / iter
+# last best: 0,53s / iter
+@profview myAmica = fit(SingleModelAmica, x; maxiter=30, do_sphering=true, remove_mean=true, m=3, scale=beta_init[:, :, 1], location=mu_init[:, :, 1], A=copy(A_init[:, :, 1]))
+
+x = [1 4; 4 1]*Float64.([1.0 2 3; 4 5 6])
+A_init = [1.0 0.003; -0.05 1.0]
+beta_init = [1.1 0.9; 1.0 0.9; 0.9 0.8]
+mu_init = [0.1 0.9; -0.01 0.0; 0.0 -0.02]
+am = fit(SingleModelAmica,x;maxiter=6, m=3, scale=beta_init, location=mu_init, A=copy(A_init), do_sphering = false)
+@test am.A ≈ [0.8761290481633254 0.7147631091971822; 0.48207664428431446 0.6993666404188701]
+@test am.LL[6] ≈ -1.701977346216155
+@test am.Lt[3] ≈ -3.4842563175935526
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 343ac21..1c86945 100755
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,9 +10,9 @@ using MAT
     beta_init = [1.1 0.9; 1.0 0.9; 0.9 0.8]
     mu_init = [0.1 0.9; -0.01 0.0; 0.0 -0.02]
     am = fit(SingleModelAmica,x;maxiter=6, m=3, scale=beta_init, location=mu_init, A=copy(A_init), do_sphering = false)
-    @test am.A == [0.8761290481633254 0.7147631091971822; 0.48207664428431446 0.6993666404188701]
-    @test am.LL[6] == -1.701977346216155
-    @test am.Lt[3] == -3.4842563175935526
+    @test am.A ≈ [0.8761290481633254 0.7147631091971822; 0.48207664428431446 0.6993666404188701]
+    @test am.LL[6] ≈ -1.701977346216155
+    @test am.Lt[3] ≈ -3.4842563175935526
 
     #Test multi model
     file = matopen("test/testdata/supersmall_data.mat")
@@ -23,8 +23,8 @@ using MAT
     close(file)
     file = matopen("test/testdata/supersmall_data_results.mat") #contains results for A and LL after 4 iterations (2 Models, 3 GGs)
     am = fit(MultiModelAmica,x;maxiter=4, m=3,M = 2,scale=beta_init, location=mu_init, A=copy(A_init), do_sphering = false)
-    @test am.LL == read(file,"LL_after_4iter")
-    @test am.models[1].A == read(file, "A1_after_4iter")
-    @test am.models[2].A == read(file, "A2_after_4iter")
+    @test am.LL ≈ read(file,"LL_after_4iter")
+    @test am.models[1].A ≈ read(file, "A1_after_4iter")
+    @test am.models[2].A ≈ read(file, "A2_after_4iter")
     close(file)
 end