text:
"""LIANNtf_algorithmLIANN.py

# Author:
Richard Bruce Baxter - Copyright (c) 2020-2022 Baxter AI (baxterai.com)

# License:
MIT License

# Installation:
see LIANNtf_main.py

# Usage:
see LIANNtf_main.py

# Description:
LIANNtf algorithm LIANN - define local inhibition artificial neural network (force neural independence)

Emulates unsupervised singular value decomposition (SVD/factor analysis) learning for all hidden layers

"""

import tensorflow as tf
import numpy as np
from ANNtf2_operations import *	#generateParameterNameSeq, generateParameterName, defineNetworkParameters
import ANNtf2_operations
import ANNtf2_globalDefs
import LIANNtf_algorithmLIANN_math
import copy

#select learningAlgorithm (unsupervised learning algorithm for intermediate/hidden layers):
learningAlgorithmNone = False	#create a very large network (eg x10) neurons per layer, and perform final layer backprop only
learningAlgorithmCorrelationReset = False	#minimise correlation between layer neurons	#create a very large network (eg x10) neurons per layer, remove/reinitialise neurons that are highly correlated (redundant/not necessary to end performance), and perform final layer backprop only	#orig mode
learningAlgorithmPCA = False	#note layer construction is nonlinear (use ANNtf2_algorithmAEANN/autoencoder for nonlinear dimensionality reduction simulation)	#incomplete
learningAlgorithmCorrelationStocasticOptimise = False	#stochastic optimise weights based on objective function; minimise the correlation between layer neurons
learningAlgorithmIndependenceReset = False	#randomise neuron weights until layer output independence is detected
learningAlgorithmMaximiseAndEvenSignalStochasticOptimise = False	#stochastic optimise weights based on objective functions; #1: maximise the signal (ie successfully uninhibited) across multiple batches (entire dataset), #2: ensure that all layer neurons receive even activation across multiple batches (entire dataset)
learningAlgorithmUninhibitedImpermanenceReset = False	#increase the permanence of uninhibited neuron weights, and stocastically modify weights based on their impermanence
learningAlgorithmUninhibitedHebbianStrengthen = False	#strengthen weights of successfully activated neurons
learningAlgorithmPerformanceInhibitStocasticOptimise = True	#learn to inhibit neurons in net for a given task, maximising final layer performance	#neurons remain permanently inhibited, not just during training of weights; see Nactive
learningAlgorithmUnnormalisedActivityReset = False	#ensure that average layer activation lies between a minimum and maximum level	#regularise layer neural activity

#intialise network properties (configurable);
positiveExcitatoryWeights = False	#requires testing	#required for biological plausibility of most learningAlgorithms
positiveExcitatoryWeightsActivationFunctionOffsetDisable = False
supportSkipLayers = False #fully connected skip layer network	#TODO: add support for skip layers	#see ANNtf2_algorithmFBANN for template
supportMultipleNetworks = True

#intialise network properties;
largeBatchSize = False
generateLargeNetwork = False	#large number of layer neurons is required for learningAlgorithmUninhibitedHebbianStrengthen:useZAcoincidenceMatrix
generateNetworkStatic = False
generateDeepNetwork = False
generateVeryLargeNetwork = False

#debug parameters;
debugFastTrain = False
debugSmallBatchSize = False	#small batch size for debugging matrix output

#select learningAlgorithmFinalLayer (supervised learning algorithm for final layer/testing):
learningAlgorithmFinalLayerBackpropHebbian = True	#only apply backprop (effective hebbian) learning at final layer
if(learningAlgorithmFinalLayerBackpropHebbian):
	positiveExcitatoryWeightsFinalLayer = False	#allow negative weights on final layer to emulate standard backprop/hebbian learning

#default sparsity
estNetworkActivationSparsity = 0.5	#50% of neurons are expected to be active during standard propagation (no inhibition)

#intialise algorithm specific parameters;
inhibitionAlgorithmBinary = False	#simplified inhibition algorithm implementation - binary on/off
inhibitionAlgorithmArtificial = False	#simplified inhibition algorithm implementation
inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = False	#inhibit layer if more than x lateral neurons active
inhibitionAlgorithmArtificialSparsity = False	#inhibition signal increases with number of simultaneously active neurons
#inhibitionAlgorithmSimulation = True	#default inhibition algorithm (propagate signal through simulated inhibitory neurons)

#A thresholding;
Athreshold = False
AthresholdValue = 1.0	#do not allow output signal to exceed 1.0

#LIANN hidden layer vs final layer hebbian execution staging;
supportDimensionalityReductionLimitFrequency = False
supportDimensionalityReductionFirstPhaseOnly = True	#perform LIANN in first phase only (x epochs of training), then apply hebbian learning at final layer
if(supportDimensionalityReductionFirstPhaseOnly):
	supportDimensionalityReductionLimitFrequency = False
	supportDimensionalityReductionFirstPhaseOnlyNumEpochs = 1
else:
	supportDimensionalityReductionLimitFrequency = True
	if(supportDimensionalityReductionLimitFrequency):
		supportDimensionalityReductionLimitFrequencyStep = 1000

#intialise algorithm specific parameters;
enableInhibitionTrainAndInhibitSpecificLayerOnly = True
applyInhibitoryNetworkDuringTest = False
randomlyActivateWeightsDuringTrain = False

#enable shared stocastic optimisation parameters;
learningAlgorithmStochastic = False
if(learningAlgorithmCorrelationStocasticOptimise):
	learningAlgorithmStochastic = True
elif(learningAlgorithmMaximiseAndEvenSignalStochasticOptimise):
	learningAlgorithmStochastic = True

#learning algorithm customisation;
if(learningAlgorithmNone):
	#can pass different task datasets through a shared randomised net
	#note learningAlgorithmCorrelationReset requires supportSkipLayers - see LIANNtf_algorithmIndependentInput/AEANNtf_algorithmIndependentInput:learningAlgorithmLIANN for similar implementation
	#positiveExcitatoryWeights = True	#optional
	generateDeepNetwork = True	#optional	#used for algorithm testing
	generateVeryLargeNetwork = True
	generateNetworkStatic = True
elif(learningAlgorithmCorrelationReset):
	#note learningAlgorithmCorrelationReset requires supportSkipLayers - see LIANNtf_algorithmIndependentInput/AEANNtf_algorithmIndependentInput:learningAlgorithmLIANN for similar implementation
	#positiveExcitatoryWeights = True	#optional

	maxCorrelation = 0.95	#requires tuning
	supportDimensionalityReductionRandomise	= True	#randomise weights of highly correlated neurons, else zero them (effectively eliminating neuron from network, as its weights are no longer able to be trained)

	generateDeepNetwork = True	#optional	#used for algorithm testing
	generateVeryLargeNetwork = True
	generateNetworkStatic = True
elif(learningAlgorithmPCA):
	#positiveExcitatoryWeights = True	#optional
	largeBatchSize = True	#1 PCA is performed across entire dataset [per layer]
elif(learningAlgorithmIndependenceReset):
	#positiveExcitatoryWeights = True	#optional
	inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = True	#mandatory
	fractionIndependentInstancesAcrossBatchRequired = 0.3	#divide by number of neurons on layer	#if found x% of independent instances, then record neuron as independent (solidify weights)	#FUTURE: will depend on number of neurons on current layer and previous layer	#CHECKTHIS: requires calibration
	largeBatchSize = True
elif(learningAlgorithmStochastic):
	inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = True	#optional
	if(learningAlgorithmCorrelationStocasticOptimise):
		learningAlgorithmStochasticAlgorithm = "correlation"
		#positiveExcitatoryWeights = True	#optional
		#learning objective function: minimise the correlation between layer neurons
	elif(learningAlgorithmMaximiseAndEvenSignalStochasticOptimise):
		learningAlgorithmStochasticAlgorithm = "maximiseAndEvenSignal"
		#positiveExcitatoryWeights = True	#optional?
		#learning objective functions:
			#1: maximise the signal (ie successfully uninhibited) across multiple batches (entire dataset)
			#2: ensure that all layer neurons receive even activation across multiple batches (entire dataset)
		metric1Weighting = 1.0
		metric2Weighting = 1000.0	#normalise metric2Weighting relative to metric1Weighting; eg metric1 =  0.9575, metric2 =  0.000863842
	numberStochasticIterations = 10
	updateParameterSubsetSimultaneously = False	#current tests indiciate this is not required/beneficial with significantly high batch size
	if(updateParameterSubsetSimultaneously):
		numberOfSubsetsTrialledPerBaseParameter = 10	#decreases speed, but provides more robust parameter updates
		parameterUpdateSubsetSize = 5	#decreases speed, but provides more robust parameter updates
	else:
		numberOfSubsetsTrialledPerBaseParameter = 1
		parameterUpdateSubsetSize = 1
	NETWORK_PARAM_INDEX_TYPE = 0
	NETWORK_PARAM_INDEX_LAYER = 1
	NETWORK_PARAM_INDEX_H_CURRENT_LAYER = 2
	NETWORK_PARAM_INDEX_H_PREVIOUS_LAYER = 3
	NETWORK_PARAM_INDEX_VARIATION_DIRECTION = 4
elif(learningAlgorithmUninhibitedImpermanenceReset):
	inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = True	#optional
	#positiveExcitatoryWeights = True	#optional?
	enableInhibitionTrainAndInhibitSpecificLayerOnly = False	#always enable inhibition	#CHECKTHIS
	applyInhibitoryNetworkDuringTest = True	#CHECKTHIS (set False)
	Wpermanence = {}
	Bpermanence = {}
	WpermanenceInitial = 0.1
	BpermanenceInitial = 0.1
	WpermanenceUpdateRate = 0.1
	BpermanenceUpdateRate = 0.1
	permanenceNumberBatches = 10	#if permanenceUpdateRate=1, average number of batches to reset W to random values
	solidificationRate = 0.1
elif(learningAlgorithmUninhibitedHebbianStrengthen):
	tuneInhibitionNeurons = False	#optional
	useZAcoincidenceMatrix = True	#reduce connection weights for unassociated neurons
	positiveExcitatoryWeights = True	#mandatory (requires testing)
	positiveExcitatoryWeightsThresholds = True	#do not allow weights to exceed 1.0 / fall below 0.0 [CHECKTHIS]
	Athreshold = True	#prevents incremental increase in signal per layer
	alwaysApplyInhibition = False
	if(useZAcoincidenceMatrix):
		alwaysApplyInhibition = True	#inhibition is theoretically allowed at all times with useZAcoincidenceMatrix as it simply biases the network against a correlation between layer k neurons (inhibition is not set up to only allow X/1 neuron to fire)
	if(alwaysApplyInhibition):
		#TODO: note network sparsity/inhibition must be configured such that at least one neuron fires per layer
		positiveExcitatoryWeightsActivationFunctionOffsetDisable = True	#activation function will always be applied to Z signal comprising positive+negative components	#CHECKTHIS
		inhibitionAlgorithmArtificialSparsity = True
		generateLargeNetwork = True 	#large is required because it will be sparsely activated due to constant inhibition
		generateNetworkStatic = True	#equal number neurons per layer for unsupervised layers/testing
		enableInhibitionTrainAndInhibitSpecificLayerOnly = False	#always enable inhibition
		applyInhibitoryNetworkDuringTest = True
	else:
		inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = True	#optional
		enableInhibitionTrainAndInhibitSpecificLayerOnly = True
		applyInhibitoryNetworkDuringTest = False
	randomlyActivateWeightsDuringTrain = False	#randomly activate x weights (simulating input at simulataneous time interval t)
	if(randomlyActivateWeightsDuringTrain):
		randomlyActivateWeightsProbability = 1.0
	WinitialisationFactor = 1.0	#initialise network with relatively low weights	#network will be trained (weights will be increased) up until point where activation inhibited
	BinitialisationFactor = 1.0	#NOTUSED
	weightDecay = False
	if(useZAcoincidenceMatrix):
		useZAcoincidenceMatrix = True	#reduce connection weights for unassociated neurons
		if(useZAcoincidenceMatrix):
			#inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive = False	#useZAcoincidenceMatrix requires real negative weights
			normaliseWeightUpdates = False
		else:
			normaliseWeightUpdates = True	#unimplemented: strengthen/update weights up to some maxima as determined by input signal strength (prevent runaway increase in weight strength up to 1.0)
	else:
		weightDecay = True	#constant neural net weight decay, such that network can be continuously trained
		weightDecayRate = 0.0	#defined by defineTrainingParametersLIANN
		useZAcoincidenceMatrix = False
		normaliseWeightUpdates = False

	maxWeightUpdateThreshold = False	#max threshold weight updates to learningRate
	#TODO: ensure learning algorithm does not result in runnaway weight increases
elif(learningAlgorithmPerformanceInhibitStocasticOptimise):
	enableInhibitionTrainAndInhibitSpecificLayerOnly = False	#always enable inhibition
	inhibitionAlgorithmBinary = True
	inhibitionAlgorithmBinaryInitialiseRandom = True
	supportDimensionalityReductionRandomise	= True	#randomise weights of highly correlated neurons, else zero them (effectively eliminating neuron from network, as its weights are no longer able to be trained)
	generateDeepNetwork = True	#optional	#used for algorithm testing
	generateVeryLargeNetwork = True
	generateNetworkStatic = True
elif(learningAlgorithmUnnormalisedActivityReset):
	supportDimensionalityReductionRandomise	= True	#randomise weights of highly correlated neurons, else zero them (effectively eliminating neuron from network, as its weights are no longer able to be trained)
	generateDeepNetwork = True	#optional	#used for algorithm testing
	generateVeryLargeNetwork = True
	generateNetworkStatic = True
	supportDimensionalityReductionRandomise = True
	supportDimensionalityReductionRegulariseActivityMinAvg = 0.01	#requires tuning
	supportDimensionalityReductionRegulariseActivityMaxAvg = 0.99	#requires tuning
	supportDimensionalityReductionRandomise	= True

learningRate = 0.0	#defined by defineTrainingParametersLIANN

#network/activation parameters;
#forward excitatory connections;
W = {}
B = {}
if(supportMultipleNetworks):
	WallNetworksFinalLayer = None
	BallNetworksFinalLayer = None
if(learningAlgorithmStochastic):
	Wbackup = {}
	Bbackup = {}
useBinaryWeights = False

if(generateVeryLargeNetwork):
	generateLargeNetworkRatio = 100	#100	#default: 10
else:
	if(generateLargeNetwork):
		generateLargeNetworkRatio = 3
	else:
		generateLargeNetworkRatio = 1

positiveExcitatoryWeightsActivationFunctionOffset = False
if(positiveExcitatoryWeights):
	if(positiveExcitatoryWeightsActivationFunctionOffsetDisable):
		positiveExcitatoryWeightsActivationFunctionOffset = False
	else:
		positiveExcitatoryWeightsActivationFunctionOffset = True
	normaliseInput = False	#TODO: verify that the normalisation operation will not disort the code's capacity to process a new data batch the same as an old data batch
	normalisedAverageInput = 1.0	#normalise input signal	#arbitrary
	if(positiveExcitatoryWeightsActivationFunctionOffset):
		positiveExcitatoryThreshold = 0.5	#1.0	#weights are centred around positiveExcitatoryThreshold, from 0.0 to positiveExcitatoryThreshold*2	#arbitrary
	Wmean = 0.5	#arbitrary
	WstdDev = 0.05	#stddev of weight initialisations	#CHECKTHIS
else:
	normaliseInput = False
	Wmean = 0.0
	WstdDev = 0.05	#stddev of weight initialisations

randomUniformMin = 0.0
randomUniformMax = 1.0
randomUniformMid = 0.5

if(inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive or inhibitionAlgorithmArtificialSparsity):
	inhibitionAlgorithmArtificial = True

if(inhibitionAlgorithmArtificial):
	if(inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive):
		inhibitionAlgorithmMoreThanXLateralNeuronActiveFraction = True
		if(inhibitionAlgorithmMoreThanXLateralNeuronActiveFraction):
			inhibitionAlgorithmMoreThanXLateralNeuronActiveFractionValue = 0.25	#fraction of the layer active allowed before inhibition
		else:
			inhibitionAlgorithmMoreThanXLateralNeuronActiveValue = 1	#maximum (X) number of neurons activate allowed before inhibition
else:
	inhibitionFactor1 = 1.0	#pass through signal	#positiveExcitatoryThreshold	#CHECKTHIS: requires recalibration for activationFunction:positiveExcitatoryWeights
	inhibitionFactor2 = estNetworkActivationSparsity	#-(WstdDev)	#inhibition contributes a significant (nullifying) effect on layer activation	#CHECKTHIS: requires calibration
	if(randomlyActivateWeightsDuringTrain):
		inhibitionFactor1 = inhibitionFactor1
		inhibitionFactor2 = (inhibitionFactor2*randomlyActivateWeightsProbability)	#the lower the average activation, the lower the inhibition
	#TODO: inhibitionFactor1/inhibitionFactor2 requires recalibration for activationFunction:positiveExcitatoryWeights
	singleInhibitoryNeuronPerLayer = False	#simplified inhibitory layer
	#lateral inhibitory connections (incoming/outgoing);
	IWi = {}
	IBi = {}
	IWo = {}
	IWiWeights = inhibitionFactor1	#need at least 1/IWiWeights active neurons per layer for the inhibitory neuron to become activated	#CHECKTHIS: requires calibration	#WstdDev*2	#0.5	#0.3
	IWoWeights = inhibitionFactor2	#will depend on activation sparsity of network (higher the positive activation, higher the inhibition required)	#requires calibration such that more than x (e.g. 1) active neuron on a layer will inhibit the layer
	In_h = []

if(inhibitionAlgorithmBinary):
  Nactive = {}  #effective bool [1.0 or 0.0]; whether neuron is active/inhibited

if(learningAlgorithmIndependenceReset):
	Bindependent = {}	#independent neurons previously identified	#effective boolean (0.0 or 1.0)	#FUTURE: consider making this a continuous variable, such that the higher the independence the less the variable is randomly shuffled per training iteration

#Network parameters
n_h = []
numberOfLayers = 0
numberOfNetworks = 0
datasetNumClasses = 0

#note high batchSize is required for learningAlgorithmStochastic algorithm objective functions (>= 100)
def defineTrainingParameters(dataset):
	global learningRate
	global weightDecayRate

	if(learningAlgorithmStochastic):
		learningRate = 0.001
	elif(learningAlgorithmUninhibitedHebbianStrengthen):
		learningRate = 0.001
		weightDecayRate = learningRate/10.0	#CHECKTHIS	#will depend on learningRate
	else:
		learningRate = 0.005

	if(debugSmallBatchSize):
		batchSize = 10
	else:
		if(largeBatchSize):
			batchSize = 1000	#current implementation: batch size should contain all examples in training set
		else:
			batchSize = 100	#3	#100
	if(generateDeepNetwork):
		numEpochs = 100	#higher num epochs required for convergence
	else:
		numEpochs = 10	#100 #10
	if(debugFastTrain):
		trainingSteps = batchSize
	else:
		trainingSteps = 10000	#1000

	displayStep = 100

	return learningRate, trainingSteps, batchSize, displayStep, numEpochs

def defineNetworkParameters(num_input_neurons, num_output_neurons, datasetNumFeatures, dataset, numberOfNetworksSet):

	global n_h
	global numberOfLayers
	global numberOfNetworks
	global datasetNumClasses
	if(not inhibitionAlgorithmArtificial):
		global In_h

	firstHiddenLayerNumberNeurons = num_input_neurons*generateLargeNetworkRatio
	if(generateDeepNetwork):
		numberOfLayers = 3
	else:
		numberOfLayers = 2
	n_h, numberOfLayers, numberOfNetworks, datasetNumClasses = defineNetworkParametersDynamic(num_input_neurons, num_output_neurons, datasetNumFeatures, dataset, numberOfNetworksSet, numberOfLayers, firstHiddenLayerNumberNeurons, generateNetworkStatic)
	#n_h, numberOfLayers, numberOfNetworks, datasetNumClasses = ANNtf2_operations.defineNetworkParameters(num_input_neurons, num_output_neurons, datasetNumFeatures, dataset, numberOfNetworksSet, generateLargeNetwork=generateLargeNetwork, generateNetworkStatic=generateNetworkStatic)

	if(not inhibitionAlgorithmArtificial):
		if(singleInhibitoryNeuronPerLayer):
			In_h = [1] * len(n_h)	#create one inhibitory neuron per layer
		else:
			In_h = copy.copy(n_h)	#create one inhibitory neuron for every excitatory neuron

	return numberOfLayers

def defineNeuralNetworkParameters():

	print("numberOfNetworks", numberOfNetworks)

	global randomNormal
	global randomUniformIndex
	randomNormal = tf.initializers.RandomNormal(mean=Wmean, stddev=WstdDev)
	#randomNormal = tf.initializers.RandomNormal()
	randomNormalFinalLayer = tf.initializers.RandomNormal()
	randomUniformIndex = tf.initializers.RandomUniform(minval=randomUniformMin, maxval=randomUniformMax)	#not available:	minval=0, maxval=numberOfSharedComputationalUnitsNeurons, dtype=tf.dtypes.int32;

	for networkIndex in range(1, numberOfNetworks+1):
		for l1 in range(1, numberOfLayers+1):
			#forward excitatory connections;
			EWlayer = randomNormal([n_h[l1-1], n_h[l1]])
			EBlayer = tf.zeros(n_h[l1])
			if(positiveExcitatoryWeights):
				EWlayer = tf.abs(EWlayer)	#ensure randomNormal generated weights are positive
				if((l1 == numberOfLayers) and not positiveExcitatoryWeightsFinalLayer):
					EWlayer = randomNormalFinalLayer([n_h[l1-1], n_h[l1]])
			if(learningAlgorithmUninhibitedHebbianStrengthen):
				EWlayer = tf.multiply(EWlayer, WinitialisationFactor)
				EBlayer = tf.multiply(EBlayer, BinitialisationFactor)
			W[generateParameterNameNetwork(networkIndex, l1, "W")] = tf.Variable(EWlayer)
			B[generateParameterNameNetwork(networkIndex, l1, "B")] = tf.Variable(EBlayer)

			if(learningAlgorithmIndependenceReset):
				Bindependent[generateParameterNameNetwork(networkIndex, l1, "Bindependent")] = tf.Variable(EBlayer)		#initialise all neurons to zero (false)
			elif(learningAlgorithmStochastic):
				Wbackup[generateParameterNameNetwork(networkIndex, l1, "W")] = tf.Variable(W[generateParameterNameNetwork(networkIndex, l1, "W")])
				Bbackup[generateParameterNameNetwork(networkIndex, l1, "B")] = tf.Variable(B[generateParameterNameNetwork(networkIndex, l1, "B")])
			elif(learningAlgorithmUninhibitedImpermanenceReset):
				EWlayerPermanence = tf.multiply(tf.ones([n_h[l1-1], n_h[l1]]), WpermanenceInitial)
				EBlayerPermanence = tf.multiply(tf.ones(n_h[l1]), BpermanenceInitial)
				Wpermanence[generateParameterNameNetwork(networkIndex, l1, "Wpermanence")] = tf.Variable(EWlayerPermanence)
				Bpermanence[generateParameterNameNetwork(networkIndex, l1, "Bpermanence")] = tf.Variable(EBlayerPermanence)

			if(not inhibitionAlgorithmArtificial):
				#lateral inhibitory connections (incoming/outgoing);
				#do not currently train inhibitory weights;
				IWilayer = tf.multiply(tf.ones([n_h[l1], In_h[l1]]), IWiWeights)		#CHECKTHIS: inhibitory neuron firing is a function of current (lateral) layer (not previous layer)
				IBilayer = tf.zeros(In_h[l1])
				if(singleInhibitoryNeuronPerLayer):
					IWoWeightsL = IWoWeights
				else:
					IWoWeightsL = IWoWeights/In_h[l1]	#normalise across number inhibitory neurons
				IWolayer = tf.multiply(tf.ones([In_h[l1], n_h[l1]]), IWoWeightsL)
				IWi[generateParameterNameNetwork(networkIndex, l1, "IWi")] = tf.Variable(IWilayer)
				IBi[generateParameterNameNetwork(networkIndex, l1, "IBi")] = tf.Variable(IBilayer)
				IWo[generateParameterNameNetwork(networkIndex, l1, "IWo")] = tf.Variable(IWolayer)

			if(inhibitionAlgorithmBinary):
				if(inhibitionAlgorithmBinaryInitialiseRandom):
					Nactivelayer = randomUniformIndex([n_h[l1]])	#tf.cast(), dtype=tf.dtypes.bool)
					Nactivelayer = tf.greater(Nactivelayer, randomUniformMid)
					Nactivelayer = tf.cast(Nactivelayer, dtype=tf.dtypes.float32)
				else:
					Nactivelayer = tf.ones(n_h[l1])
				Nactive[generateParameterNameNetwork(networkIndex, l1, "Nactive")] = tf.Variable(Nactivelayer)

	if(supportMultipleNetworks):
		if(numberOfNetworks > 1):
			global WallNetworksFinalLayer
			global BallNetworksFinalLayer
			WlayerF = randomNormal([n_h[numberOfLayers-1]*numberOfNetworks, n_h[numberOfLayers]])
			WallNetworksFinalLayer =  tf.Variable(WlayerF)
			Blayer = tf.zeros(n_h[numberOfLayers])
			BallNetworksFinalLayer	= tf.Variable(Blayer)	#not currently used

def neuralNetworkPropagation(x, networkIndex=1):
	return neuralNetworkPropagationLIANNtest(x, networkIndex)

def neuralNetworkPropagationLIANNtest(x, networkIndex=1, l=None):
	return neuralNetworkPropagationLIANNminimal(x, networkIndex, l)

def neuralNetworkPropagationLayer(x, networkIndex=1, l=None):
	return neuralNetworkPropagationLIANNminimal(x, networkIndex, l)
	#return neuralNetworkPropagationLIANN(x, None, networkIndex, trainWeights=False)

def neuralNetworkPropagationLIANNtrainIntro(x, y=None, networkIndex=1):
	if(enableInhibitionTrainAndInhibitSpecificLayerOnly):
		for l in range(1, numberOfLayers+1):
			if(l < numberOfLayers):
				return neuralNetworkPropagationLIANNtrain(x, y, networkIndex, layerToTrain=l)
	else:
		return neuralNetworkPropagationLIANNtrain(x, y, networkIndex, layerToTrain=None)

#if(supportMultipleNetworks):
def neuralNetworkPropagationAllNetworksFinalLayer(AprevLayer):
	Z = tf.add(tf.matmul(AprevLayer, WallNetworksFinalLayer), BallNetworksFinalLayer)
	#Z = tf.matmul(AprevLayer, WallNetworksFinalLayer)
	pred = tf.nn.softmax(Z)
	return pred

#minimal code extracted from neuralNetworkPropagationLIANN;
def neuralNetworkPropagationLIANNminimal(x, networkIndex=1, l=None):

	randomlyActivateWeights = False

	if(l == None):
		maxLayer = numberOfLayers
	else:
		maxLayer = l

	AprevLayer = x
	ZprevLayer = x
	for l in range(1, maxLayer+1):

		enableInhibition = False
		if(not enableInhibitionTrainAndInhibitSpecificLayerOnly):
			enableInhibition = True

		A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)

		if(learningAlgorithmFinalLayerBackpropHebbian):
			A = tf.stop_gradient(A)

		AprevLayer = A
		ZprevLayer = Z

	if(maxLayer == numberOfLayers):
		return tf.nn.softmax(Z)
	else:
		return A

def neuralNetworkPropagationLIANNtrain(x, y=None, networkIndex=1, layerToTrain=None):

	if(normaliseInput):
		#TODO: verify that the normalisation operation will not disort the code's capacity to process a new data batch the same as an old data batch
		averageTotalInput = tf.math.reduce_mean(x)
		#print("averageTotalInput = ", averageTotalInput)
		x = tf.multiply(x, normalisedAverageInput/averageTotalInput)	#normalise input wrt positiveExcitatoryThreshold
		#averageTotalInput = tf.math.reduce_mean(x)

	if(layerToTrain is None):
		maxLayer = numberOfLayers
	else:	#ie !enableInhibitionTrainAndInhibitSpecificLayerOnly
		maxLayer = layerToTrain

	AprevLayer = x
	ZprevLayer = x
	for l in range(1, maxLayer+1):

		trainLayer = False
		enableInhibition = False
		randomlyActivateWeights = False
		if(enableInhibitionTrainAndInhibitSpecificLayerOnly):
			if(l == layerToTrain):
				#enableInhibition = False
				enableInhibition = True
				trainLayer = True
		else:
			if(l < numberOfLayers):
				enableInhibition = True
				trainLayer = True
		if(randomlyActivateWeightsDuringTrain):
			randomlyActivateWeights = True

		if(trainLayer):
			#CHECKTHIS: verify learning algorithm (how to modify weights to maximise independence between neurons on each layer)
			if(learningAlgorithmNone):
				neuralNetworkPropagationLIANNlearningAlgorithmNone(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmCorrelationReset):
				neuralNetworkPropagationLIANNlearningAlgorithmCorrelationReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmPCA):
				neuralNetworkPropagationLIANNlearningAlgorithmPCA(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmIndependenceReset):
				neuralNetworkPropagationLIANNlearningAlgorithmIndependenceReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmStochastic):
				neuralNetworkPropagationLIANNlearningAlgorithmStochastic(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmUninhibitedImpermanenceReset):
				neuralNetworkPropagationLIANNlearningAlgorithmUninhibitedImpermanenceReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmUninhibitedHebbianStrengthen):
				neuralNetworkPropagationLIANNlearningAlgorithmUninhibitedHebbianStrengthen(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
			elif(learningAlgorithmPerformanceInhibitStocasticOptimise):
				neuralNetworkPropagationLIANNlearningAlgorithmPerformanceInhibitStocasticOptimise(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights, x, y)
			elif(learningAlgorithmUnnormalisedActivityReset):
				neuralNetworkPropagationLIANNlearningAlgorithmUnnormalisedActivityReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)

			A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition=(not enableInhibitionTrainAndInhibitSpecificLayerOnly), randomlyActivateWeights=False)
		else:
			A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights=False)

		AprevLayer = A
		ZprevLayer = Z

	return tf.nn.softmax(Z)

def calculatePropagationLoss(x, y, networkIndex=1):
	costCrossEntropyWithLogits = False
	pred = neuralNetworkPropagation(x, networkIndex)
	target = y
	lossCurrent = calculateLossCrossEntropy(pred, target, datasetNumClasses, costCrossEntropyWithLogits)
	#acc = calculateAccuracy(pred, target)	#only valid for softmax class targets
	return lossCurrent

def neuralNetworkPropagationLIANNlearningAlgorithmNone(networkIndex, AprevLayer, ZprevLayer, l1, enableInhibition, randomlyActivateWeights):
	pass

def neuralNetworkPropagationLIANNlearningAlgorithmCorrelationReset(networkIndex, AprevLayer, ZprevLayer, l1, enableInhibition, randomlyActivateWeights):
	A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l1)
	#measure and minimise correlation between layer neurons;
	neuronActivationCorrelationMinimisation(networkIndex, n_h, l1, A, randomNormal, Wf=W, Wfname="W", Wb=None, Wbname=None, updateAutoencoderBackwardsWeights=False, supportSkipLayers=supportSkipLayers, supportDimensionalityReductionRandomise=supportDimensionalityReductionRandomise, maxCorrelation=maxCorrelation)

def neuralNetworkPropagationLIANNlearningAlgorithmPCA(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights):
	#Afinal, Zfinal, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)	#batched
	SVDinputMatrix = LIANNtf_algorithmLIANN_math.generateSVDinputMatrix(l, n_h, AprevLayer)
	U, Sigma, VT = LIANNtf_algorithmLIANN_math.calculateSVD(M=SVDinputMatrix, k=n_h[l])
	AW = LIANNtf_algorithmLIANN_math.calculateWeights(l, n_h, SVDinputMatrix, U, Sigma, VT)
	W[generateParameterNameNetwork(networkIndex, l, "W")] = AW

	#weights = U -> Sigma -> VT	[linear]
	#M_reduced = reduce_to_k_dim(M=spikeCoincidenceMatrix, k=n_h[l])

def neuralNetworkPropagationLIANNlearningAlgorithmIndependenceReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights):
	layerHasDependentNeurons = True
	Bind = Bindependent[generateParameterNameNetwork(networkIndex, l, "Bindependent")]
	if(count_zero(Bind) > 0):	#more than 1 dependent neuron on layer
		layerHasDependentNeurons = True
	else:
		layerHasDependentNeurons = False

	while(layerHasDependentNeurons):
		Afinal, Zfinal, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)	#batched

		AnumActive = tf.math.count_nonzero(Afinal, axis=1)	#batched
		Aindependent = tf.equal(AnumActive, 1)	#batched
		Aindependent = tf.dtypes.cast(Aindependent, dtype=tf.dtypes.float32)	#batched
		Aindependent = tf.expand_dims(Aindependent, 1)	#batched
		#print("Afinal = ", Afinal)
		#print("AnumActive = ", AnumActive)
		#print("Aindependent = ", Aindependent)

		Aactive = tf.greater(Afinal, 0)	#2D: batched, for every k neuron
		Aactive = tf.dtypes.cast(Aactive, dtype=tf.dtypes.float32) 	#2D: batched, for every k neuron
		#print("Aactive = ", Aactive)
		#ex

		AactiveAndIndependent = tf.multiply(Aactive, Aindependent)	#2D: batched, for every k neuron
		AactiveAndIndependent = tf.reduce_sum(AactiveAndIndependent, axis=0) #for every k neuron

		AactiveAndIndependentPass = tf.greater(AactiveAndIndependent, fractionIndependentInstancesAcrossBatchRequired*n_h[l])	 #for every k neuron
		#print("AactiveAndIndependentPass = ", AactiveAndIndependentPass)

		BindBool = tf.dtypes.cast(Bind, dtype=tf.dtypes.bool)
		AactiveAndIndependentPassRequiresSolidifying = tf.logical_and(AactiveAndIndependentPass, tf.logical_not(BindBool))
		#print("AactiveAndIndependentPass = ", AactiveAndIndependentPass)
		#print("BindBool = ", BindBool)
		print("AactiveAndIndependentPassRequiresSolidifying = ", AactiveAndIndependentPassRequiresSolidifying)
		BindNew = tf.logical_or(BindBool, AactiveAndIndependentPassRequiresSolidifying)
		BdepNew = tf.logical_not(BindNew)

		#update layer weights (reinitialise weights for all dependent neurons);
		BindNew = tf.dtypes.cast(BindNew, dtype=tf.dtypes.float32)
		BdepNew = tf.dtypes.cast(BdepNew, dtype=tf.dtypes.float32)
		EWlayerDep = randomNormal([n_h[l-1], n_h[l]])
		if(positiveExcitatoryWeights):
			EWlayerDep = tf.abs(EWlayerDep)	#ensure randomNormal generated weights are positive
		EBlayerDep = tf.zeros(n_h[l])
		EWlayerDep = tf.multiply(EWlayerDep, BdepNew)	#requires broadcasting
		EBlayerDep = tf.multiply(EBlayerDep, BdepNew)
		EWlayerInd = W[generateParameterNameNetwork(networkIndex, l, "W")]
		EBlayerInd = B[generateParameterNameNetwork(networkIndex, l, "B")]
		EWlayerInd = tf.multiply(EWlayerInd, BindNew)	#requires broadcasting
		EBlayerInd = tf.multiply(EBlayerInd, BindNew)
		EWlayerNew = tf.add(EWlayerDep, EWlayerInd)
		EBlayerNew = tf.add(EBlayerDep, EBlayerInd)
		W[generateParameterNameNetwork(networkIndex, l, "W")] = EWlayerNew
		B[generateParameterNameNetwork(networkIndex, l, "B")] = EBlayerNew
		#print("EWlayerNew = ", EWlayerNew)

		#print("BdepNew = ", BdepNew)
		#print("BindNew = ", BindNew)

		Bindependent[generateParameterNameNetwork(networkIndex, l, "Bindependent")] = BindNew	#update independence record
		Bind = BindNew
		if(count_zero(Bind) > 0):	#more than 1 dependent neuron on layer
			layerHasDependentNeurons = True
			#print("layerHasDependentNeurons: count_zero(Bind) = ", count_zero(Bind))
		else:
			layerHasDependentNeurons = False
			#print("!layerHasDependentNeurons")

def neuralNetworkPropagationLIANNlearningAlgorithmStochastic(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights):

	if(learningAlgorithmStochastic):
		if(useBinaryWeights):
			variationDirections = 1
		else:
			variationDirections = 2

	#code from ANNtf2_algorithmLREANN_expSUANN;
	for s in range(numberStochasticIterations):
		for hIndexCurrentLayer in range(0, n_h[l]):
			for hIndexPreviousLayer in range(0, n_h[l-1]+1):
				if(hIndexPreviousLayer == n_h[l-1]):	#ensure that B parameter updates occur/tested less frequently than W parameter updates
					parameterTypeWorB = 0
				else:
					parameterTypeWorB = 1
				for variationDirectionInt in range(variationDirections):

					networkParameterIndexBase = (parameterTypeWorB, l, hIndexCurrentLayer, hIndexPreviousLayer, variationDirectionInt)

					metricBase = learningAlgorithmStochasticCalculateMetric(networkIndex, AprevLayer, ZprevLayer, l)

					for subsetTrialIndex in range(0, numberOfSubsetsTrialledPerBaseParameter):

						accuracyImprovementDetected = False

						currentSubsetOfParameters = []
						currentSubsetOfParameters.append(networkParameterIndexBase)

						for s in range(1, parameterUpdateSubsetSize):
							networkParameterIndex = getRandomNetworkParameter(networkIndex, currentSubsetOfParameters)
							currentSubsetOfParameters.append(networkParameterIndex)

						for s in range(0, parameterUpdateSubsetSize):
							networkParameterIndex = currentSubsetOfParameters[s]

							if(not useBinaryWeights):
								if(networkParameterIndex[NETWORK_PARAM_INDEX_VARIATION_DIRECTION] == 1):
									variationDiff = learningRate
								else:
									variationDiff = -learningRate

							if(networkParameterIndex[NETWORK_PARAM_INDEX_TYPE] == 1):
								#Wnp = W[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "W")].numpy()
								#currentVal = Wnp[networkParameterIndex[NETWORK_PARAM_INDEX_H_PREVIOUS_LAYER], networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]]
								currentVal = W[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "W")][networkParameterIndex[NETWORK_PARAM_INDEX_H_PREVIOUS_LAYER], networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]].numpy()

								#print("currentVal = ", currentVal)
								#print("W1 = ", W[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "W")])
								if(useBinaryWeights):
									if(useBinaryWeightsReduceMemoryWithBool):
										newVal = not currentVal
									else:
										newVal = float(not bool(currentVal))
										#print("newVal = ", newVal)
								else:
									newVal = currentVal + variationDiff

								if(positiveExcitatoryWeights):
									newVal = max(newVal, 0)	#do not allow weights fall below zero [CHECKTHIS]

								W[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "W")][networkParameterIndex[NETWORK_PARAM_INDEX_H_PREVIOUS_LAYER], networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]].assign(newVal)

								#print("W2 = ", W[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "W")])
							else:
								#Bnp = B[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "B")].numpy()
								#currentVal = Bnp[networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]]
								currentVal = B[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "B")][networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]].numpy()

								if(useBinaryWeights):
									if(useBinaryWeightsReduceMemoryWithBool):
										newVal = not currentVal
									else:
										newVal = float(not bool(currentVal))
								else:
									newVal = currentVal + variationDiff

								if(positiveExcitatoryWeights):
									newVal = max(newVal, 0)	#do not allow weights fall below zero [CHECKTHIS]

								B[generateParameterNameNetwork(networkIndex, networkParameterIndex[NETWORK_PARAM_INDEX_LAYER], "B")][networkParameterIndex[NETWORK_PARAM_INDEX_H_CURRENT_LAYER]].assign(newVal)

						metricAfterStochasticUpdate = learningAlgorithmStochasticCalculateMetric(networkIndex, AprevLayer, ZprevLayer, l)
						#print("metricBase = ", metricBase)
						#print("metricAfterStochasticUpdate = ", metricAfterStochasticUpdate)

						if(metricAfterStochasticUpdate > metricBase):
							#print("(metricAfterStochasticUpdate > metricBase)")
							accuracyImprovementDetected = True
							metricBase = metricAfterStochasticUpdate
						#else:
							#print("(metricAfterStochasticUpdate < metricBase)")

						if(accuracyImprovementDetected):
							#retain weight update
							Wbackup[generateParameterNameNetwork(networkIndex, l, "W")].assign(W[generateParameterNameNetwork(networkIndex, l, "W")])
							Bbackup[generateParameterNameNetwork(networkIndex, l, "B")].assign(B[generateParameterNameNetwork(networkIndex, l, "B")])
						else:
							#restore weights
							W[generateParameterNameNetwork(networkIndex, l, "W")].assign(Wbackup[generateParameterNameNetwork(networkIndex, l, "W")])
							B[generateParameterNameNetwork(networkIndex, l, "B")].assign(Bbackup[generateParameterNameNetwork(networkIndex, l, "B")])

def neuralNetworkPropagationLIANNlearningAlgorithmUninhibitedImpermanenceReset(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights):
	Afinal, Zfinal, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)

	#update W/B permanence;
	Afinal2D = tf.reduce_mean(Afinal, axis=0)	#average across batch
	Afinal2D = tf.expand_dims(Afinal2D, axis=0)	#make compatible shape to W
	WpermanenceUpdate = tf.multiply(Afinal2D, WpermanenceUpdateRate)	#verify that broadcasting works
	WpermanenceNew = tf.add(Wpermanence[generateParameterNameNetwork(networkIndex, l, "Wpermanence")], WpermanenceUpdate)	#increase the permanence of neuron weights that successfully fired
	Wpermanence[generateParameterNameNetwork(networkIndex, l, "Wpermanence")] = WpermanenceNew
	print("WpermanenceUpdate = ", WpermanenceUpdate)

	#stochastically modify weights based on permanence values:
	Wupdate = randomNormal([n_h[l-1], n_h[l]])
	Wupdate = tf.divide(Wupdate, Wpermanence[generateParameterNameNetwork(networkIndex, l, "Wpermanence")])
	Wupdate = tf.divide(Wupdate, permanenceNumberBatches)
	Wnew = tf.add(W[generateParameterNameNetwork(networkIndex, l, "W")], Wupdate)
	if(positiveExcitatoryWeights):
		Wnew = tf.maximum(Wnew, 0)	#do not allow weights fall below zero [CHECKTHIS]
	W[generateParameterNameNetwork(networkIndex, l, "W")] = Wnew
	#print("Wupdate = ", Wupdate)

def neuralNetworkPropagationLIANNlearningAlgorithmUninhibitedHebbianStrengthen(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights):
	AW = W[generateParameterNameNetwork(networkIndex, l, "W")]
	Afinal, Zfinal, EWactive = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
	#print("Zfinal = ", Zfinal)

	if(useZAcoincidenceMatrix):
		AWcontribution = tf.matmul(tf.transpose(ZprevLayer), Afinal)	#increase excitatory weights that contributed to the output signal	#hebbian
	else:
		AWcontribution = tf.matmul(tf.transpose(AprevLayer), Afinal)	#increase excitatory weights that contributed to the output signal	#hebbian

	if(randomlyActivateWeights):
		#do not apply weight updates to temporarily suppressed weights [CHECKTHIS];
		AWcontribution = tf.multiply(AWcontribution, EWactive)

	if(normaliseWeightUpdates):
		print("LIANNtf_algorithmLIANN:neuralNetworkPropagationLIANN error - normaliseWeightUpdates: normaliseWeightUpdatesReduceConnectionWeightsForUnassociatedNeurons unimplemented")
	else:
		if(maxWeightUpdateThreshold):
			AWcontribution = tf.minimum(AWcontribution, 1.0)

	AWupdate = tf.multiply(AWcontribution, learningRate)
	#print("AWupdate = ", AWupdate)

	AW = tf.add(AW, AWupdate)	#apply weight updates

	if(weightDecay):
		#apply decay to all weights;
		AWdecay = -weightDecayRate
		#print("AWdecay = ", AWdecay)
		AW = tf.add(AW, AWdecay)
		#print("AWdecay = ", AWdecay)

	if(positiveExcitatoryWeightsThresholds):
		AW = tf.minimum(AW, 1.0)	#do not allow weights to exceed 1.0 [CHECKTHIS]
		AW = tf.maximum(AW, 0)	#do not allow weights fall below zero [CHECKTHIS]

	W[generateParameterNameNetwork(networkIndex, l, "W")] = AW

def neuralNetworkPropagationLIANNlearningAlgorithmPerformanceInhibitStocasticOptimise(networkIndex, AprevLayer, ZprevLayer, l1, enableInhibition, randomlyActivateWeights, x=None, y=None):
	A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l1)
	#randomly select a neuron k on layer to trial inhibition performance;
	lossCurrent = calculatePropagationLoss(x, y, networkIndex)	#moved 15 Mar 2022
	Nactivelayer = Nactive[generateParameterNameNetwork(networkIndex, l1, "Nactive")]
	NactivelayerBackup = Nactivelayer #tf.Variable(Nactivelayer)
	layerInhibitionIndex = tf.cast(randomUniformIndex([1])*n_h[l1], tf.int32)[0].numpy()
	#print("layerInhibitionIndex = ", layerInhibitionIndex)
	if(inhibitionAlgorithmBinary):
		inhibitionValue = randomUniformIndex([1])
		inhibitionValue = tf.greater(inhibitionValue, randomUniformMid)
		inhibitionValue = tf.cast(inhibitionValue, dtype=tf.dtypes.float32)
		inhibitionValue = inhibitionValue[0].numpy()
		#print("inhibitionValue = ", inhibitionValue)
	else:
		inhibitionValue = 0.0
	Nactivelayer = tf.Variable(ANNtf2_operations.modifyTensorRowColumn(Nactivelayer, True, layerInhibitionIndex, inhibitionValue, False))	#tf.Variable added to retain formatting
	#print("NactivelayerBackup = ", NactivelayerBackup)
	#print("Nactivelayer = ", Nactivelayer)
	Nactive[generateParameterNameNetwork(networkIndex, l1, "Nactive")] = Nactivelayer
	loss = calculatePropagationLoss(x, y, networkIndex)
	#acc = calculateAccuracy(pred, target)	#only valid for softmax class targets
	if(loss < lossCurrent):
		lossCurrent = loss
		#print("loss < lossCurrent; loss = ", loss, ", lossCurrent = ", lossCurrent)
	else:
		Nactive[generateParameterNameNetwork(networkIndex, l1, "Nactive")] = NactivelayerBackup
		#print("loss !< lossCurrent; loss = ", loss, ", lossCurrent = ", lossCurrent)

def neuralNetworkPropagationLIANNlearningAlgorithmUnnormalisedActivityReset(networkIndex, AprevLayer, ZprevLayer, l1, enableInhibition, randomlyActivateWeights):
	A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l1)
	neuronActivationRegularisation(networkIndex, n_h, l1, A, randomNormal, Wf=W, Wfname="W", Wb=None, Wbname=None, updateAutoencoderBackwardsWeights=False, supportSkipLayers=supportSkipLayers, supportDimensionalityReductionRandomise=supportDimensionalityReductionRandomise, supportDimensionalityReductionRegulariseActivityMinAvg=supportDimensionalityReductionRegulariseActivityMinAvg, supportDimensionalityReductionRegulariseActivityMaxAvg=supportDimensionalityReductionRegulariseActivityMaxAvg)

def forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition=False, randomlyActivateWeights=False):
	#forward excitatory connections;
	EWactive = None
	EW = W[generateParameterNameNetwork(networkIndex, l, "W")]
	if(randomlyActivateWeights):
		#print("EW = ", EW)
		EWactive = tf.less(tf.random.uniform(shape=EW.shape), randomlyActivateWeightsProbability)	#initialised from 0.0 to 1.0
		EWactive = tf.dtypes.cast(EWactive, dtype=tf.dtypes.float32)
		#print("EWactive = ", EWactive)
		#EWactive = tf.dtypes.cast(tf.random.uniform(shape=EW.shape, minval=0, maxval=2, dtype=tf.dtypes.int32), dtype=tf.dtypes.float32)
		EW = tf.multiply(EW, EWactive)
	Z = tf.add(tf.matmul(AprevLayer, EW), B[generateParameterNameNetwork(networkIndex, l, "B")])
	A = activationFunction(Z, n_h[l-1])

	#lateral inhibitory connections (incoming/outgoing);
	if(enableInhibition):
		Afinal, Zfinal = forwardIterationInhibition(networkIndex, AprevLayer, ZprevLayer, l, A, Z)
	else:
		Zfinal = Z
		Afinal = A

	return Afinal, Zfinal, EWactive

def forwardIterationInhibition(networkIndex, AprevLayer, ZprevLayer, l, A, Z):
	if(inhibitionAlgorithmBinary):
		Afinal = tf.multiply(A, Nactive[generateParameterNameNetwork(networkIndex, l, "Nactive")])
		Zfinal = tf.multiply(Z, Nactive[generateParameterNameNetwork(networkIndex, l, "Nactive")])
	else:
		if(inhibitionAlgorithmArtificial):
			if(inhibitionAlgorithmArtificialSparsity):
				prevLayerSize = n_h[l-1]
				inhibitionResult = tf.math.reduce_mean(AprevLayer, axis=1)	#or ZprevLayer?	#batched
				#print("inhibitionResult = ", inhibitionResult)
				inhibitionResult = tf.multiply(inhibitionResult, prevLayerSize)	#normalise by prev layer size		#batched
				inhibitionResult = tf.multiply(inhibitionResult, Wmean)	#normalise by average weight
				inhibitionResult = tf.expand_dims(inhibitionResult, axis=1)	#batched
				Zfinal = tf.subtract(Z, inhibitionResult)	#requires broadcasting
				#print("Z = ", Z)
				#print("Zfinal = ", Zfinal)
				Afinal = activationFunction(Zfinal, prevLayerSize=prevLayerSize)
			elif(inhibitionAlgorithmArtificialMoreThanXLateralNeuronActive):
				layerSize = n_h[l]
				numActiveLateralNeurons = tf.math.count_nonzero(A, axis=1)
				if(inhibitionAlgorithmMoreThanXLateralNeuronActiveFraction):
					numberActiveNeuronsAllowed = inhibitionAlgorithmMoreThanXLateralNeuronActiveFractionValue*layerSize
				else:
					numberActiveNeuronsAllowed = inhibitionAlgorithmMoreThanXLateralNeuronActiveValue
				numberActiveNeuronsAllowed = int(numberActiveNeuronsAllowed)
				#print("numActiveLateralNeurons = ", numActiveLateralNeurons)
				#print("numberActiveNeuronsAllowed = ", numberActiveNeuronsAllowed)
				inhibitionResult = tf.greater(numActiveLateralNeurons, numberActiveNeuronsAllowed)

				inhibitionResult = tf.logical_not(inhibitionResult)
				inhibitionResult = tf.dtypes.cast(inhibitionResult, dtype=tf.dtypes.float32)
				inhibitionResult = tf.expand_dims(inhibitionResult, axis=1)
				#print("numActiveLateralNeurons = ", numActiveLateralNeurons)
				#print("inhibitionResult = ", inhibitionResult)
				Zfinal = tf.multiply(Z, inhibitionResult)	#requires broadcasting
				Afinal = tf.multiply(A, inhibitionResult)	#requires broadcasting
		else:
			#if((l < numberOfLayers) or positiveExcitatoryWeightsFinalLayer):

			#print("AprevLayer = ", AprevLayer)
			#print("Z = ", Z)

			IZi = tf.matmul(A, IWi[generateParameterNameNetwork(networkIndex, l, "IWi")])	#CHECKTHIS: inhibitory neuron firing is a function of current (lateral) layer (not previous layer)
			IAi = activationFunction(IZi, n_h[l-1])
			#print("IZi = ", IZi)
			#print("IAi = ", IAi)
			IZo = tf.matmul(IAi, IWo[generateParameterNameNetwork(networkIndex, l, "IWo")])
			#print("W = ", W[generateParameterNameNetwork(networkIndex, l, "W")])
			#print("IZo = ", IZo)

			#final activations;
			Zfinal = tf.add(Z, IZo)
			#print("Zfinal = ", Zfinal)
			Afinal = activationFunction(Zfinal, n_h[l-1])

	if(Athreshold):
		Afinal = tf.minimum(Afinal, AthresholdValue)
		#print("Afinal = ", Afinal)

	return Afinal, Zfinal

#LIANNlearningAlgorithmCorrelation metric:

def neuronActivationCorrelationMinimisation(networkIndex, n_h, l1, A, randomNormal, Wf, Wfname="W", Wb=None, Wbname=None, updateAutoencoderBackwardsWeights=False, supportSkipLayers=False, supportDimensionalityReductionRandomise=True, maxCorrelation=0.95):

	resetNeuronIfSameValueAcrossBatch = False #reset neuron if all values of a neuron k being the same value across the batch
	randomlySelectCorrelatedNeuronToReset = False	#randomly select one of each correlated neuron to reset

	useCorrelationMatrix = True	#only implementation currently available

	Atransposed = tf.transpose(A)
	if(useCorrelationMatrix):
		correlationMatrix = LIANNtf_algorithmLIANN_math.calculateOffDiagonalCorrelationMatrix(A, nanReplacementValue=0.0, getOffDiagonalCorrelationMatrix=True)	#off diagonal correlation matrix is required so that do not duplicate k1->k2 and k2->k1 correlations	#CHECKTHIS: nanReplacementValue
		#nanReplacementValue=0.0; will set the correlation as 0 if all values of a neuron k being the same value across the batch
		#print("correlationMatrix = ", correlationMatrix)
		#print("correlationMatrix.shape = ", correlationMatrix.shape)

	if(useCorrelationMatrix):
		if(randomlySelectCorrelatedNeuronToReset):
			correlationMatrixRotated = np.transpose(correlationMatrix)
			k1MaxCorrelation = correlationMatrix.max(axis=0)
			k2MaxCorrelation = correlationMatrixRotated.max(axis=0)
			#print("k1MaxCorrelation = ", k1MaxCorrelation)
			#print("k2MaxCorrelation = ", k2MaxCorrelation)
			kSelect = np.random.randint(0, 2, size=k1MaxCorrelation.shape)
			mask1 = kSelect.astype(bool)
			mask2 = np.logical_not(mask1)
			mask1 = mask1.astype(float)
			mask2 = mask2.astype(float)
			k1MaxCorrelation = np.multiply(k1MaxCorrelation, mask1)
			k2MaxCorrelation = np.multiply(k2MaxCorrelation, mask2)
			kMaxCorrelation = np.add(k1MaxCorrelation, k2MaxCorrelation)
			#print("correlationMatrix = ", correlationMatrix)
			#print("correlationMatrixRotated = ", correlationMatrixRotated)
			#print("k1MaxCorrelation = ", k1MaxCorrelation)
			#print("k2MaxCorrelation = ", k2MaxCorrelation)
			#print("mask1 = ", mask1)
			#print("mask2 = ", mask2)
			#print("kMaxCorrelation = ", kMaxCorrelation)
		else:
			k1MaxCorrelation = correlationMatrix.max(axis=0)
			k2MaxCorrelation = correlationMatrix.max(axis=1)
			#k1MaxCorrelation = np.amax(correlationMatrix, axis=0)	#reduce max
			#k2MaxCorrelation = np.amax(correlationMatrix, axis=1)	#reduce max
			kMaxCorrelation = np.maximum(k1MaxCorrelation, k2MaxCorrelation)
		#kMaxCorrelationIndex = correlationMatrix.argmax(axis=0)	#or axis=1
		kMaxCorrelation = tf.convert_to_tensor(kMaxCorrelation, dtype=tf.dtypes.float32)	#make sure same type as A
		#print("kMaxCorrelation;", kMaxCorrelation)

		if(resetNeuronIfSameValueAcrossBatch):
			AbatchAllZero = tf.reduce_sum(A, axis=0)
			AbatchAllZero = tf.equal(AbatchAllZero, 0.0)
			AbatchAllZero = tf.cast(AbatchAllZero, tf.float32)
			kMaxCorrelation = tf.add(kMaxCorrelation, AbatchAllZero)	#set kMaxCorrelation[k]=1.0 if AbatchAllZero[k]=True
			#print("AbatchAllZero;", AbatchAllZero)

	else:
		#incomplete;
		for k1 in range(n_h[l1]):
			#calculate maximum correlation;
			k1MaxCorrelation = 0.0
			for k2 in range(n_h[l1]):
				if(k1 != k2):
					Ak1 = Atransposed[k1]	#Ak: 1d vector of batchsize
					Ak2 = Atransposed[k2]	#Ak: 1d vector of batchsize
					k1k2correlation = calculateCorrelation(Ak1, Ak2)	#undefined

	#generate masks (based on highly correlated k/neurons);
	#print("kMaxCorrelation = ", kMaxCorrelation)
	kPassArray = tf.less(kMaxCorrelation, maxCorrelation)
	randomiseLayerNeurons(networkIndex, n_h, l1, kPassArray, randomNormal, Wf, Wfname, Wb, Wbname, updateAutoencoderBackwardsWeights, supportSkipLayers, supportDimensionalityReductionRandomise)

def neuronActivationRegularisation(networkIndex, n_h, l1, A, randomNormal, Wf, Wfname="W", Wb=None, Wbname=None, updateAutoencoderBackwardsWeights=False, supportSkipLayers=False, supportDimensionalityReductionRandomise=True, supportDimensionalityReductionRegulariseActivityMinAvg=0.1, supportDimensionalityReductionRegulariseActivityMaxAvg=0.9):
	#CHECKTHIS: treat any level/intensity of activation the same
	Aactive = tf.cast(A, tf.bool)
	AactiveFloat = tf.cast(Aactive, tf.float32)
	neuronActivationFrequency = tf.reduce_mean(AactiveFloat, axis=0)
	#print("neuronActivationFrequency = ", neuronActivationFrequency)
	kPassArray = tf.logical_and(tf.greater(neuronActivationFrequency, supportDimensionalityReductionRegulariseActivityMinAvg), tf.less(neuronActivationFrequency, supportDimensionalityReductionRegulariseActivityMaxAvg))
	#print("kPassArray = ", kPassArray)
	randomiseLayerNeurons(networkIndex, n_h, l1, kPassArray, randomNormal, Wf, Wfname, Wb, Wbname, updateAutoencoderBackwardsWeights, supportSkipLayers, supportDimensionalityReductionRandomise)

def randomiseLayerNeurons(networkIndex, n_h, l1, kPassArray, randomNormal, Wf, Wfname="W", Wb=None, Wbname=None, updateAutoencoderBackwardsWeights=False, supportSkipLayers=False, supportDimensionalityReductionRandomise=True):
	kFailArray = tf.logical_not(kPassArray)
	#print("kPassArray = ", kPassArray)
	#print("kFailArray = ", kFailArray)
	kPassArrayF = tf.expand_dims(kPassArray, axis=0)
	kFailArrayF = tf.expand_dims(kFailArray, axis=0)
	kPassArrayF = tf.cast(kPassArrayF, tf.float32)
	kFailArrayF = tf.cast(kFailArrayF, tf.float32)
	if(updateAutoencoderBackwardsWeights):
		kPassArrayB = tf.expand_dims(kPassArray, axis=1)
		kFailArrayB = tf.expand_dims(kFailArray, axis=1)
		kPassArrayB = tf.cast(kPassArrayB, tf.float32)
		kFailArrayB = tf.cast(kFailArrayB, tf.float32)

	#apply masks to weights (randomise specific k/neurons);
	if(supportSkipLayers):
		for l2 in range(0, l1):
			if(l2 < l1):
				#randomize or zero
				if(supportDimensionalityReductionRandomise):
					WlayerFrand = randomNormal([n_h[l2], n_h[l1]])
				else:
					WlayerFrand = tf.zeros([n_h[l2], n_h[l1]], dtype=tf.dtypes.float32)
				Wf[generateParameterNameNetworkSkipLayers(networkIndex, l2, l1, Wfname)] = applyMaskToWeights(Wf[generateParameterNameNetworkSkipLayers(networkIndex, l2, l1, Wfname)], WlayerFrand, kPassArrayF, kFailArrayF)
				if(updateAutoencoderBackwardsWeights):
					if(supportDimensionalityReductionRandomise):
						WlayerBrand = randomNormal([n_h[l1], n_h[l2]])
					else:
						WlayerBrand = tf.zeros([n_h[l1], n_h[l2]], dtype=tf.dtypes.float32)
					Wb[generateParameterNameNetworkSkipLayers(networkIndex, l2, l1, Wbname)] = applyMaskToWeights(Wb[generateParameterNameNetworkSkipLayers(networkIndex, l2, l1, Wbname)], WlayerBrand, kPassArrayB, kFailArrayB)
	else:
		if(supportDimensionalityReductionRandomise):
			WlayerFrand = randomNormal([n_h[l1-1], n_h[l1]])
		else:
			WlayerFrand = tf.zeros([n_h[l1-1], n_h[l1]], dtype=tf.dtypes.float32)
		Wf[generateParameterNameNetwork(networkIndex, l1, Wfname)] = applyMaskToWeights(Wf[generateParameterNameNetwork(networkIndex, l1, Wfname)], WlayerFrand, kPassArrayF, kFailArrayF)
		if(updateAutoencoderBackwardsWeights):
			if(supportDimensionalityReductionRandomise):
				WlayerBrand = randomNormal([n_h[l1], n_h[l1-1]])
			else:
				WlayerBrand = tf.zeros([n_h[l1], n_h[l1-1]], dtype=tf.dtypes.float32)
			Wb[generateParameterNameNetwork(networkIndex, l1, Wbname)] = applyMaskToWeights(Wb[generateParameterNameNetwork(networkIndex, l1, Wbname)], WlayerBrand, kPassArrayB, kFailArrayB)

def applyMaskToWeights(Wlayer, WlayerRand, kPassArray, kFailArray):
	WlayerFail = tf.multiply(WlayerRand, kFailArray)
	#print("WlayerFail = ", WlayerFail)
	WlayerPass = tf.multiply(Wlayer, kPassArray)
	#print("WlayerPass = ", WlayerPass)
	Wlayer = tf.add(WlayerPass, WlayerFail)
	return Wlayer

#LIANNlearningAlgorithmStochastic metric:

def learningAlgorithmStochasticCalculateMetric(networkIndex, AprevLayer, ZprevLayer, l):
	randomlyActivateWeights = False
	if(randomlyActivateWeightsDuringTrain):
		randomlyActivateWeights = true

	if(learningAlgorithmCorrelationStocasticOptimise):
		enableInhibition = False
		A, Z, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
		metric = learningAlgorithmStochasticCalculateMetricCorrelation(A)
	elif(learningAlgorithmMaximiseAndEvenSignalStochasticOptimise):
		enableInhibition = True
		Afinal, Zfinal, _ = forwardIteration(networkIndex, AprevLayer, ZprevLayer, l, enableInhibition, randomlyActivateWeights)
		metric = learningAlgorithmStochasticCalculateMetricMaximiseAndEvenSignal(Afinal, metric1Weighting, metric2Weighting)
	return metric

def learningAlgorithmStochasticCalculateMetricCorrelation(A):
	#print("A = ", A)
	meanCorrelation = LIANNtf_algorithmLIANN_math.calculateCorrelationMean(A)
	print("meanCorrelation = ", meanCorrelation)
	metric = 1 - meanCorrelation
	#print("metric = ", metric)
	return metric

def learningAlgorithmStochasticCalculateMetricMaximiseAndEvenSignal(Afinal, metric1Weighting, metric2Weighting):
	#learning objective functions:
	#1: maximise the signal (ie successfully uninhibited) across multiple batches (entire dataset)
	#2: ensure that all layer neurons receive even activation across multiple batches (entire dataset)

	#print("Afinal = ", Afinal)

	AfinalThresholded = tf.greater(Afinal, 0.0)	#threshold signal such that higher average weights are not preferenced
	AfinalThresholded = tf.dtypes.cast(AfinalThresholded, dtype=tf.dtypes.float32)
	#print("Afinal = ", Afinal)
	#print("AfinalThresholded = ", AfinalThresholded)

	metric1 = tf.reduce_mean(AfinalThresholded)	#average output across batch, across layer

	#stdDevAcrossLayer = tf.math.reduce_std(Afinal, axis=1)	#stddev calculated across layer [1 result per batch index]
	#metric2 = tf.reduce_mean(stdDevAcrossLayer)	#average output across batch

	stdDevAcrossBatches = tf.math.reduce_mean(Afinal, axis=0)	 #for each dimension (k neuron in layer); calculate the mean across all batch indices
	metric2 = tf.math.reduce_std(stdDevAcrossBatches)	#then calculate the std dev across these values

	metric1 = metric1.numpy()
	metric2 = metric2.numpy()
	#print("metric1 = ", metric1)
	#print("metric2 = ", metric2)

	metric1 = metric1*metric1Weighting
	metric2 = metric2*metric2Weighting
	#print("metric1 = ", metric1)
	#print("metric2 = ", metric2)
	if(metric2 != 0):
		metric = metric1/metric2
	else:
		metric = 0.0

	return metric

def getRandomNetworkParameter(networkIndex, currentSubsetOfParameters):

	foundNewParameter = False
	while not foundNewParameter:

		variationDirection = random.randint(2)
		layer = random.randint(1, len(n_h))
		parameterTypeWorBtemp = random.randint(n_h[layer-1]+1)	#ensure that B parameter updates occur/tested less frequently than W parameter updates	#OLD: random.randint(2)
		if(parameterTypeWorBtemp == n_h[layer-1]):
			parameterTypeWorB = 0
		else:
			parameterTypeWorB = 1
		hIndexCurrentLayer = random.randint(n_h[layer])	#randomNormal(n_h[l])
		hIndexPreviousLayer = random.randint(n_h[layer-1]) #randomNormal(n_h[l-1])
		networkParameterIndex = (parameterTypeWorB, layer, hIndexCurrentLayer, hIndexPreviousLayer, variationDirection)

		matches = [item for item in currentSubsetOfParameters if item == networkParameterIndex]
		if len(matches) == 0:
			foundNewParameter = True

	return networkParameterIndex

def activationFunction(Z, prevLayerSize=None):
	return reluCustomPositive(Z, prevLayerSize)

def reluCustomPositive(Z, prevLayerSize=None):
	if(positiveExcitatoryWeightsActivationFunctionOffset):
		#CHECKTHIS: consider sigmoid instead of relu
		#offset required because negative weights are not used:

		#Zoffset = tf.ones(Z.shape)
		#Zoffset = tf.multiply(Zoffset, normalisedAverageInput)
		#Zoffset = tf.multiply(Zoffset, Wmean)
		#Zoffset = tf.multiply(Zoffset, prevLayerSize)
		Zpred = prevLayerSize*normalisedAverageInput*Wmean
		Zoffset = Zpred
		#print("Zoffset = ", Zoffset)

		Z = tf.subtract(Z, Zoffset)
		A = tf.nn.relu(Z)
		A = tf.multiply(A, 2.0)	#double the slope of A to normalise the input:output signal
		#print("A = ", A)
	else:
		A = tf.nn.relu(Z)
		#A = tf.nn.sigmoid(Z)
	return A

def count_zero(M, axis=None):	#emulates count_nonzero behaviour
	if axis is not None:
		nonZeroElements = tf.math.count_nonzero(M, axis=axis)
		totalElements = tf.shape(M)[axis]
		zeroElements = tf.subtract(totalElements, nonZeroElements)
	else:
		totalElements = tf.size(M)
		nonZeroElements = tf.math.count_nonzero(M).numpy()
		zeroElements = tf.subtract(totalElements, nonZeroElements)
		zeroElements = zeroElements.numpy()
	return zeroElements

max_stars_repo_name: baxterai/LIANNtf

max_issues_repo_name: baxterai/LIANNtf

max_forks_repo_name: baxterai/LIANNtf


text:
import numpy as np
from scipy import linalg
from copy import copy

def distance(delta,N):
	''' assume periodic boundary conditions'''
	signs = np.sign(delta)
	# return np.where(np.abs(delta) > 0.5, -signs*(1 - np.abs(delta)), delta)
	return np.where(np.abs(delta) > (0.5 * N),\
	 		-signs*(N - np.abs(delta)), delta)

def gaussian(x,y,s):
	return np.exp(-(x**2+y**2)/2./s**2)#1./2.*np.pi/s**2*

class Network:
	def __init__(self, from_size, to_size, Version):
		self.from_size = from_size
		self.to_size = to_size
		self.rng = np.random.RandomState(20210311 + Version)

	def create_matrix(self, conn_params, profile, **kwargs):
		grid = np.linspace(0,self.from_size[0],self.from_size[0],endpoint=False)
		xfrom,yfrom = np.meshgrid(grid,grid)
		grid = np.linspace(0,self.to_size[0],self.to_size[0],endpoint=False)
		xto,yto = np.meshgrid(grid,grid)
		xdelta = distance(xto[None,None,:,:]-xfrom[:,:,None,None],self.to_size[0])
		ydelta = distance(yto[None,None,:,:]-yfrom[:,:,None,None],self.to_size[0])

		if profile in ("Gaussian",):
			sigma = conn_params["sigma"]
			ampl = conn_params["ampl"]
			disc_gaussian = gaussian(xdelta,ydelta,sigma) #/ np.prod(self.from_size)
			# norm_factor = np.sum(disc_gaussian,axis=(2,3))[:,:,None,None]
			conn_matrix = ampl * disc_gaussian #/ norm_factor
		elif profile=="Gabor":
			sigma = conn_params["sigma"] * np.nanmax(xdelta)
			ampl = conn_params["ampl"]
			theta = conn_params["theta"]
			psi = conn_params["psi"]
			disc_gaussian = ampl * gaussian(xdelta,ydelta,sigma)
			Lambda = 2*np.pi/conn_params["freq"] * np.nanmax(xdelta)
			gamma = 1.
			sigma_x = sigma
			sigma_y = float(sigma) / gamma

			# Bounding box
			nstds = 3  # Number of standard deviation sigma
			xmax = np.max(np.concatenate([np.abs(nstds * sigma_x * np.cos(theta)),\
									   np.abs(nstds * sigma_y * np.sin(theta))]))
			xmax = np.ceil(max(1, xmax))
			ymax = np.max(np.concatenate([np.abs(nstds * sigma_x * np.sin(theta)),\
					 				   np.abs(nstds * sigma_y * np.cos(theta))]))
			ymax = np.ceil(max(1, ymax))
			print("xd",np.nanmax(xdelta),xmax,ymax,sigma,Lambda)
			x = xdelta * 2 * xmax / np.nanmax(xdelta)
			y = ydelta * 2 * ymax / np.nanmax(ydelta)
			# Rotation
			x_theta = x * np.cos(theta[None,None,:,:]) - y * np.sin(theta[None,None,:,:])
			y_theta = x * np.sin(theta[None,None,:,:]) + y * np.cos(theta[None,None,:,:])
			if isinstance(psi,np.ndarray):
				conn_matrix =\
					 np.exp(-.5 * (x_theta ** 2 / sigma_x ** 2 + y_theta ** 2 / sigma_y ** 2)) *\
				 	 np.cos(2 * np.pi / Lambda * x_theta + psi[:,:,None,None])
			else:
				conn_matrix =\
					 np.exp(-.5 * (x_theta ** 2 / sigma_x ** 2 + y_theta ** 2 / sigma_y ** 2)) *\
				 	 np.cos(2 * np.pi / Lambda * x_theta + psi)
		elif profile=="Mexican-hat":
			sigma1 = conn_params["sigma1"]
			sigma2 = conn_params["sigma2"]
			ampl1 = 1./2/np.pi/sigma1**2
			ampl2 = 1./2/np.pi/sigma2**2
			disc_gaussian1 = gaussian(xdelta,ydelta,sigma1) #/ np.prod(self.from_size) / self.Nvert[0]
			# disc_gaussian1 /= np.sum(disc_gaussian1,axis=(0,1))[None,None,:,:]
			disc_gaussian2 = gaussian(xdelta,ydelta,sigma2) #/ np.prod(self.from_size) / self.Nvert[0]
			# disc_gaussian2 /= np.sum(disc_gaussian2,axis=(0,1))[None,None,:,:]
			conn_matrix = ampl1 * disc_gaussian1 - ampl2 * disc_gaussian2
			conn_matrix /= np.nanmax(conn_matrix)
		else:
			print("Specified connectivity profile ({}) not found.".format(profile))

		self.conn_matrix = conn_matrix.reshape(np.prod(self.from_size),\
							np.prod(self.to_size))

		if "r_A" in kwargs:
			profile_A = kwargs.get("profile_A", "heaviside")
			arbor = self.create_arbor(radius=kwargs["r_A"],profile=profile_A)
			self.conn_matrix[np.logical_not(arbor)] = 0.

		return self.conn_matrix

	def create_arbor(self, radius, profile="heaviside"):
		"""
		arbor function: crude model of overlap of geniculocortical terminal
		arbors with cortical dendritic arbors
		"""
		grid = np.linspace(0,self.from_size[0],self.from_size[0],endpoint=False)
		xfrom,yfrom = np.meshgrid(grid,grid)
		grid = np.linspace(0,self.to_size[0],self.to_size[0],endpoint=False)
		xto,yto = np.meshgrid(grid,grid)
		xdelta = distance(xto[None,None,:,:]-xfrom[:,:,None,None],self.to_size[0])
		ydelta = distance(yto[None,None,:,:]-yfrom[:,:,None,None],self.to_size[0])

		d = np.sqrt(xdelta**2 + ydelta**2)		## absolute distance |vec(x)-vec(y)|
		d = d.reshape(np.prod(self.to_size),np.prod(self.from_size))
		if profile=="heaviside":
			arbor = (np.around(d,2)<=radius).astype(float)

		elif profile=="gaussian":
			arbor = gaussian(d,0,radius)
			arbor[np.around(d,2)>radius] = 0.0

		elif profile=="overlap":
			d_unique = np.unique(d[np.around(d,2)<=radius])
			cA = 0.5
			N = radius*(cA+2)

			delta_bin = 0.05
			i,j = np.meshgrid(np.arange(-(self.to_size[0]//2),self.to_size[0]//2+delta_bin,\
				delta_bin),\
				np.arange(-(self.from_size[0]//2),self.from_size[0]//2+delta_bin,delta_bin))
			delta1 = np.sqrt(i**2+j**2)
			delta2 = np.sqrt((d_unique[:,None,None]-i[None,:,:])**2+(j[None,:,:])**2)
			circles = (delta1 < radius)[None,:,:]*(delta2 < (radius*cA))
			norm = np.sum(delta2 < (radius*cA),axis=(1,2))
			overlap = 1.*np.sum(circles,axis=(1,2))/norm

			arbor = np.zeros_like(d,dtype=float)
			for iduni in range(d_unique.size):
				arbor[d==d_unique[iduni]] = overlap[iduni]

		elif profile=="linear_falloff":
			arbor = linear_arbor_falloff(np.around(d,2),radius=radius,\
										 inner_radius=radius*0.5)
			arbor[np.around(d,2)>radius] = 0.0

		else:
			print("Specified arbor profile ({}) not found.".format(profile))

		return arbor

ext: py

reason: import numpy,from scipy


text:
'''
/************************************************************************
 MIT License

 Copyright (c) 2018 Harsh Kakashaniya,Koyal Bhartia

 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.
 *************************************************************************/

/**
 *  @file    final.py
 *  @author  Harsh Kakashaniya and Koyal Bhartia
 *  @date    27/2/2019
 *  @version 1.0
 *
 *  @brief Project 1,AR tag detection
 *
 *  @section DESCRIPTION
 *
 *  This is code has 3 parts,
 1. tag id Tag_id_detection
 2. Superimposing image on AR tag with orientation
 3. Superimposing of cube on AR tag with orientation
 *
 */
'''
# Libraries
import argparse
import numpy as np
import os, sys
from numpy import linalg as LA
from numpy import linalg as la
import math
from PIL import Image
import random

try:
    sys.path.remove('/opt/ros/kinetic/lib/python2.7/dist-packages')
except:
    pass
import cv2
#-------------------------------------------------------------------------------

# @brief Function for converting gray scale to binary
#
#  @param Matrix
#
#  @return Matrix
#
def binary(mat):
    for row in range(0,len(mat)):
        for col in range(0,len(mat[0])):
            if (mat[row,col]>150):
                mat[row,col]=1
            else:
                mat[row,col]=0
    return mat

# @brief To give 4 points in Edgedetection
#
#  @param Image and old_ctr points
#
#  @return corners
#
def Edgedetection(image,old_ctr):

    gray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
    blurred = cv2.medianBlur(gray,3)
    (T, thresh) = cv2.threshold(blurred, 180, 255, cv2.THRESH_BINARY)
    contours, hierarchy=cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    ctr=[]
    for j, cnt in zip(hierarchy[0], contours):
        cnt_len = cv2.arcLength(cnt,True)
        cnt = cv2.approxPolyDP(cnt, 0.02*cnt_len,True)
        if cv2.contourArea(cnt) > 1000 and cv2.isContourConvex(cnt) and len(cnt) == 4  :
            cnt=cnt.reshape(-1,2)
            if j[0] == -1 and j[1] == -1 and j[3] != -1:
                ctr.append(cnt)
        #print(np.shape(ctr))
        old_ctr=ctr
    return ctr

# @brief Calculate homography
#
#  @param source and destination
#
#  @return homography
#
def homography_calc(src,des):
    #print(des)
    #print(np.shape(des))
    c1 = des[0,0]
    c2 = des[0,1]
    c3 = des[0,2]
    c4 = des[0,3]

    w1 = src[0]
    w2 = src[1]
    w3 = src[2]
    w4 = src[3]

    A=np.array([[w1[0],w1[1],1,0,0,0,-c1[0]*w1[0],-c1[0]*w1[1],-c1[0]],
                [0,0,0,w1[0], w1[1],1,-c1[1]*w1[0],-c1[1]*w1[1],-c1[1]],
                [w2[0],w2[1],1,0,0,0,-c2[0]*w2[0],-c2[0]*w2[1],-c2[0]],
                [0,0,0,w2[0], w2[1],1,-c2[1]*w2[0],-c2[1]*w2[1],-c2[1]],
                [w3[0],w3[1],1,0,0,0,-c3[0]*w3[0],-c3[0]*w3[1],-c3[0]],
                [0,0,0,w3[0], w3[1],1,-c3[1]*w3[0],-c3[1]*w3[1],-c3[1]],
                [w4[0],w4[1],1,0,0,0,-c4[0]*w4[0],-c4[0]*w4[1],-c4[0]],
                [0,0,0,w4[0], w4[1],1,-c4[1]*w4[0],-c4[1]*w4[1],-c4[1]]])

    #Performing SVD
    u, s, vt = la.svd(A)

            # normalizing by last element of v
            #v =np.transpose(v_col)
    v = vt[8:,]/vt[8][8]

    req_v = np.reshape(v,(3,3))

    return req_v
# @brief To fix image on top of AR tag
#
#  @param image , corners , source
#
#  @return superimposed image
#

def Superimposing(ctr,image,src):
    pts_dst = np.array(ctr,dtype=float)
    pts_src = np.array([[0,0],[511, 0],[511, 511],[0,511]],dtype=float)
    h = homography_calc(pts_src, pts_dst)

    temp = cv2.warpPerspective(src, h,(image.shape[1],image.shape[0]));
    cv2.fillConvexPoly(image, pts_dst.astype(int), 0, 16);

    image = image + temp;

    return image,h
# @brief To find get grid of tag image
#
#  @param image , corners
#
#  @return tag image and gray scale image
#
def perspective_for_tag(ctr,image):
    dst1 = np.array([
        [0, 0],
        [100, 0],
        [100, 100],
        [0, 100]], dtype = "float32")

    M1,status = cv2.findHomography(ctr[0], dst1)
    warp1 = cv2.warpPerspective(image.copy(), M1, (100,100))
    warp2=cv2.medianBlur(warp1,3)
    #warp2= warp1-warp1_5

    tag_image=cv2.resize(warp2, dsize=None, fx=0.08, fy=0.08)
    return tag_image,warp2
# @brief To find tag ID and tag image
#
#  @param image , corners
#
#  @return tag and tag_id
#
def Tag_id_detection(ctr,tag_image):
    gray = cv2.cvtColor(tag_image,cv2.COLOR_BGR2GRAY)
    pixel_value=binary(gray)
    print(pixel_value,'Tag Value')
    status=0
    A_ctr=ctr[0][0]
    #print(A_ctr,'ctr A')
    B_ctr=ctr[0][1]
    #print(B_ctr,'ctr B')
    C_ctr=ctr[0][2]
    #print(C_ctr,'ctr B')
    D_ctr=ctr[0][3]
    #print(D_ctr,'ctr C')
    if (pixel_value[2,2] == 1):
        L1=A_ctr
        L2=B_ctr
        L3=C_ctr
        L4=D_ctr
        status=0
        one = pixel_value[4,4]
        two = pixel_value[4,3]
        three = pixel_value[3,3]
        four = pixel_value[3,4]

    elif pixel_value[5,2]==1:
        L1=D_ctr
        L2=A_ctr
        L3=B_ctr
        L4=C_ctr
        status=1
        one = pixel_value[3,4]
        two = pixel_value[4,4]
        three = pixel_value[4,3]
        four = pixel_value[3,3]

    elif pixel_value[5,5] == 1:
        L1=C_ctr
        L2=D_ctr
        L3=A_ctr
        L4=B_ctr
        status=2
        one = pixel_value[3,3]
        two = pixel_value[3,4]
        three = pixel_value[4,4]
        four = pixel_value[4,3]

    elif pixel_value[2,5] == 1:
        L1=B_ctr
        L2=C_ctr
        L3=D_ctr
        L4=A_ctr
        status=3
        one = pixel_value[4,3]
        two = pixel_value[3,3]
        three = pixel_value[3,4]
        four = pixel_value[4,4]

    else:
        L1=A_ctr
        L2=B_ctr
        L3=C_ctr
        L4=D_ctr
        one = pixel_value[4,4]
        two = pixel_value[4,3]
        three = pixel_value[3,3]
        four = pixel_value[3,4]

    new_ctr=np.array([[L1,L2,L3,L4]])

    #print(new_ctr,'new_ctr')

    tag_id = four*8 + three*4 + two*2 + one*1
    print('Tag id value will be',tag_id)
    return new_ctr,tag_id
# @brief To draw 3d cube
#
#  @param image , perspective points in camera frame
#
#  @return image of cube
#
def draw(img, imgpts):
    imgpts = np.int32(imgpts).reshape(-1,2)
    # draw ground floor in green
    img = cv2.drawContours(img, [imgpts[:4]],-1,(255,0,0),-3)
    # draw pillars in blue color
    for i,j in zip(range(4),range(4,8)):
        img = cv2.line(img, tuple(imgpts[i]), tuple(imgpts[j]),(0,150,150),3)
        # draw top layer in red color
        img = cv2.drawContours(img, [imgpts[4:]],-1,(0,0,255),3)
    return img

# @brief Calculates projection matrix
#
#  @param Homography
#
#  @return Projection matrix for 3D transformation
#
def Projection_mat(homography):
   # homography = homography*(-1)
   # Calling the projective matrix function
   K =np.array([[1406.08415449821,0,0],
       [ 2.20679787308599, 1417.99930662800,0],
       [ 1014.13643417416, 566.347754321696,1]])

   K=K.T
   rot_trans = np.dot(la.inv(K), homography)
   col_1 = rot_trans[:, 0]
   col_2 = rot_trans[:, 1]
   col_3 = rot_trans[:, 2]
   l = math.sqrt(la.norm(col_1, 2) * la.norm(col_2, 2))
   rot_1 = col_1 / l
   rot_2 = col_2 / l
   translation = col_3 / l
   c = rot_1 + rot_2
   p = np.cross(rot_1, rot_2)
   d = np.cross(c, p)
   rot_1 = np.dot(c / np.linalg.norm(c, 2) + d / np.linalg.norm(d, 2), 1 / math.sqrt(2))
   rot_2 = np.dot(c / np.linalg.norm(c, 2) - d / np.linalg.norm(d, 2), 1 / math.sqrt(2))
   rot_3 = np.cross(rot_1, rot_2)

   projection = np.stack((rot_1, rot_2, rot_3, translation)).T
   return np.dot(K, projection)

    #def Three_d_cube(K, homography):
    #    return 0
# @brief it draws cube on image
#
#  @param Projection matrix and iamge
#
#  @return None
#
def Cube3D(proj_mat,image):
    axis = np.float32([[0,0,0,1],[0,512,0,1],[512,512,0,1],[512,0,0,1],[0,0,-512,1],[0,512,-512,1],[512,512,-512,1],[512,0,-512,1]])
    Proj= np.matmul(axis,proj_mat.T)
    # Normalize the matrix
    Norm1 = np.divide(Proj[0],Proj[0][2])
    Norm2 = np.divide(Proj[1],Proj[1][2])
    Norm3 = np.divide(Proj[2],Proj[2][2])
    Norm4 = np.divide(Proj[3],Proj[3][2])
    Norm5 = np.divide(Proj[4],Proj[4][2])
    Norm6 = np.divide(Proj[5],Proj[5][2])
    Norm7 = np.divide(Proj[6],Proj[6][2])
    Norm8 = np.divide(Proj[7],Proj[7][2])

    points = np.vstack((Norm1,Norm2,Norm3,Norm4,Norm5,Norm6,Norm7,Norm8))
    final_2d=np.delete(points,2, axis=1)
    draw(image,final_2d)
# @brief Main loop to run the code and video breaking
#
#  @param path of video, source of lena
#
#  @return image_array and size of image
#
def Imageprocessor(path,src):

    vidObj = cv2.VideoCapture(path)
    count = 0
    success = 1
    img_array=[]

    while (success):
        if (count==0):
            success, image = vidObj.read()
        height,width,layers=image.shape

        size = (width,height)
        if (count==0):
            old_corners=0
        corners=Edgedetection(image,old_corners)
        #img = cv2.drawContours(image, corners,0,(0,255,0),1)
        if(len(corners)==0):
            corners=old_corners

        tag_image,Tag=perspective_for_tag(corners,image)
        new_corners,tag_id=Tag_id_detection(corners,tag_image)

        image,h=Superimposing(new_corners,image,src)
        proj_mat=Projection_mat(h)
        #Cube3D(proj_mat,image)
        old_corners=corners
        count += 1
        print('Number of frames is',count)
        #cv2.imwrite('%d.jpg' %count,image)
        img_array.append(image)
        success, image = vidObj.read()

    return img_array,size
#--------------------------------------------------------------
# @brief Loop to run video from image array
#
#  @param final image array and size
#
#  @return None
#
def video(img_array,size):
    video=cv2.VideoWriter('video_Tag.avi',cv2.VideoWriter_fourcc(*'DIVX'), 16.0,size)
    #print(np.shape(img_array))
    for i in range(len(img_array)):
        video.write(img_array[i])
    video.release()
#---------------------------------------------------------------
# main
if __name__ == '__main__':

    # Calling the function
    src=cv2.imread('lena.jpg')
    #print(np.size(src))
    Image,size=Imageprocessor('Tag0.mp4',src)
    video(Image,size)

max_stars_repo_name: harshkakashaniya/AR-tag-detection

max_issues_repo_name: harshkakashaniya/AR-tag-detection

max_forks_repo_name: harshkakashaniya/AR-tag-detection

reason: import numpy,from numpy


text:
"""

Monte Carlo sampling with UQpy
==================================
"""

# %% md
#
# We'll be using UQpy's Monte Carlo sampling functionalities. We also use Matplotlib to display results graphically.
#
# Additionally, this demonstration opts to use Numpy's random state management to ensure that results are reproducible
# between notebook runs.

# %%

from UQpy.sampling import MonteCarloSampling
import matplotlib.pyplot as plt

from numpy.random import RandomState

# %% md
#
# Step-by-step: continuous univariate distribution
# -------------------------------------------------
# First, we import UQpy's normal distribution class.

# %%
from UQpy.distributions import Normal

# %% md
#
# We'll start by constructing two identical standard normal distributions :code:`normal1` and :code:`normal2`

# %%

normal1 = normal2 = Normal()

# %% md
#
# Next, we'll construct a :code:`MonteCarloSampling` object :code:`mc` to generate random samples following those
# distributions. Here, we specify an optional initial number of samples, :code:`nsamples` to be generated at the
# object's construction. For teh purposes of this demonstration, we also supply a random seed :code:`random_state`.
#
# We access the generated samples via the :code:`samples` attribute.

# %%

mc = MonteCarloSampling(distributions=[normal1, normal2],
                        nsamples=5,
                        random_state=RandomState(123))

mc.samples

# %% md
#
# To generate more samples on :code:`mc` after construction, we call :code:`mc.run` and once again specify
# :code:`nsamples`.

# %%

mc.run(nsamples=2, random_state=RandomState(23))

mc.samples

# %% md
#
# We can transform the samples onto the unit hypercube via applying the probability integral transformation on the
# samples to yield similar samples from the uniform distribution. We call :code:`mc.transform_u01`, from which results
# are stored in the :code:`samplesU01` attribute.

# %%

mc.transform_u01()

mc.samplesU01

# %% md
#
# We can visualize the (untransformed) samples by plotting them on axes of each distribution's range.

# %%

fig, ax = plt.subplots()
plt.title('Samples')

plt.scatter(x=mc.samples[:, 0],
            y=mc.samples[:, 1],
            marker='o')

plt.setp(ax, xlim=(-1.7, 1.7), ylim=(-2.6, 2.6))
ax.yaxis.grid(True)
ax.xaxis.grid(True)

# %% md
#
# As well, we can visualize each distribution's sample densities via histograms.

# %%

fig, ax = plt.subplots(1, 2)

for i in (0, 1):
    ax[i].set_title('Distribution ' + str(i + 1))
    ax[i].hist(mc.samples[:, i])
    ax[i].yaxis.grid(True)
    ax[i].xaxis.grid(True)

plt.setp(ax, xlim=(-3, 3), ylim=(0, 2));

# %% md
#
# Additional Examples
# -------------------------------------------------
# Continuous multivariate distribution
# """"""""""""""""""""""""""""""""""""""""
# We'll use UQpy's multivariate normal class.

# %%

from UQpy.distributions import MultivariateNormal

# %% md
#
# And we construct a multivariate normal distribution :code:`mvnormal` specifying parameters :code:`mean` with a vector
# of mean values and :code:`cov` with a covariance matrix.

# %%

mvnormal = MultivariateNormal(mean=[1, 2],
                              cov=[[4, -0.9],
                                   [-0.9, 1]])

# %% md
#
# With this distribution, we construct a :code:`MonteCarloSampling` object :code:`mvmc` and generate five samples on
# construction.

# %%

mvmc = MonteCarloSampling(distributions=mvnormal,
                          nsamples=5,
                          random_state=RandomState(456))

mvmc.samples

# %% md
#
# Mixing a multivariate and a univariate continuous distribution
# """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
# Here, we use one of our normal distributions and our multivariate normal distribution. Notice how each distribution
# has its own bundle (array) of samples per run of sampling—even when that bundle contains a single value.

# %%

mixedmc = MonteCarloSampling(distributions=[normal1, mvnormal],
                             nsamples=5,
                             random_state=RandomState(789))

mixedmc.samples

# %% md
#
# Mixing a continuous and a discrete distribution
# """""""""""""""""""""""""""""""""""""""""""""""""
# We'll use UQpy's binomial distribution class for our discrete distribution.

# %%

from UQpy.distributions import Binomial

# %% md
#
# With that, we'll construct a :code:`binomial` distribution binomial with five trials :code:`n` and a 40% probability
# :code:`p` of success per trial.

# %%

binomial = Binomial(n=5, p=0.4)

# %% md
#
# And we construct a :code:`MonteCarloSampling` object :code:`cdmv` with five initial samples using our binomial
# distribution and one of our normal distributions.

# %%

cdmv = MonteCarloSampling(distributions=[binomial, normal1],
                          nsamples=5,
                          random_state=RandomState(333))

cdmv.samples

max_stars_repo_name: SURGroup/UncertaintyQuantification

max_issues_repo_name: SURGroup/UncertaintyQuantification

max_forks_repo_name: SURGroup/UncertaintyQuantification

reason: from numpy


text:
from __future__ import division

import math

import matplotlib
import matplotlib.pyplot as plt
import numpy as np

class Wedge(object):
    def __init__(self, a, b, c, alpha):
        self.a = a
        self.b = b
        self.c = c
        self.alpha = alpha

    def __repr__(self):
        return "W["+str(self.a)+", \t"+str(self.b)+", \t"+str(self.c)+"]"

def getWedges(vertex_list):
    '''Constructs all inner Angles as Wedges of the Vertices A-B-C so that there exists exactly one Wedge B-C-X1.'''
    from operator import itemgetter, attrgetter
    allWedges = []

    for vertex in vertex_list:
        orderedconnections = []
        number_of_neighbours = 0
        for neighbour in vertex.neighbours:

            alpha = np.arctan2(
                neighbour.coords[1]-vertex.coords[1], neighbour.coords[0]-vertex.coords[0])
            if alpha < 0:
                alpha += 2*np.pi
            orderedconnections.append([neighbour.selfindex, alpha])
            number_of_neighbours += 1

        orderedconnections.sort(key=itemgetter(1))

        for i in range(number_of_neighbours):
            this = orderedconnections[i-1]
            afterthis = orderedconnections[i]
            if i == 0:
                winkel = this[1] - afterthis[1] + 2*np.pi
            else:
                winkel = afterthis[1]-this[1]
            winkel %= 2*np.pi
            neueswedge = Wedge(this[0], vertex.selfindex, afterthis[0], winkel)
            allWedges.append(neueswedge)

    allWedges.sort(key=attrgetter('a', 'b'))
    return allWedges

def getPolygon2Ds(vertex_list):
    '''Finds all closed Polygon2Ds. The algorithm starts with Wedge A-B-C and looks for Wedge B-C-X1, C-X1-X2... A polygon is found when Xn == A'''
    wedgeliste = getWedges(vertex_list)

    from bisect import bisect_left as search
    allpolygons = []

    def search(x):
        s1 = x.b
        s2 = x.c
        for wedge in wedgeliste:
            if wedge.a == s1 and wedge.b == s2:
                return wedge

    while len(wedgeliste) > 0:
        start = x = wedgeliste[0]
        currentpolygon = []
        while True:
            if not (np.pi - 0.001 < x.alpha < np.pi + 0.001):
                currentpolygon.append(x)
            x = search(x)
            if x is not None:
                wedgeliste.remove(x)
                if x is start:

                    allpolygons.append(currentpolygon)
                    break
    return allpolygons

def main(vertex_list=None):
    '''Input: list of vertices representing the Roadmap
    Output: List of all Polygon2Ds representing Lots,
    List of all Polygon2Ds representing Blocks
    List of all Polygon2Ds which are too large to be Lots
    Polygon2D representing the road-network'''
    import sys
    sys.path.append(
        '/home/lenny/Documents/Stadtprojekt/procedural_city_generation')
    if vertex_list is None:
        from procedural_city_generation.additional_stuff import pickletools

        vertex_list = pickletools.reconstruct()
        print("Reconstructing of data structure finished")

    import os
    import procedural_city_generation
    path = os.path.dirname(procedural_city_generation.__file__)

    with open(path+"/temp/border.txt", "r") as f:
        border = f.read()
    border = [int(x) for x in border.split(" ") if x is not '']

    print("Extracting Polygon2Ds")
    from procedural_city_generation.polygons import construct_polygons
    polylist = construct_polygons.getPolygon2Ds(vertex_list)
    print("Polygon2Ds extracted")
    return polylist, vertex_list

if __name__ == '__main__':
    main()

reason: import numpy


text:
import torch
from torch import nn
from torch.nn import functional as F
import numpy as np

import models.extractors as extractors

class PSPModule(nn.Module):
    def __init__(self, features, out_features=1024, sizes=(1, 2, 3, 6)):
        super(PSPModule, self).__init__()
        self.stages = []
        self.stages = nn.ModuleList([self._make_stage(features, size) for size in sizes])
        self.bottleneck = nn.Conv2d(features * (len(sizes) + 1), out_features, kernel_size=1)
        self.relu = nn.ReLU()

    def _make_stage(self, features, size):
        prior = nn.AdaptiveAvgPool2d(output_size=(size, size))
        conv = nn.Conv2d(features, features, kernel_size=1, bias=False)
        return nn.Sequential(prior, conv)

    def forward(self, feats):
        h, w = feats.size(2), feats.size(3)
        priors = [F.upsample(input=stage(feats), size=(h, w), mode='bilinear') for stage in self.stages] + [feats]
        bottle = self.bottleneck(torch.cat(priors, 1))
        return self.relu(bottle)

class PSPUpsample(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(PSPUpsample, self).__init__()
        self.conv = nn.Sequential(
            nn.Upsample(scale_factor=2, mode='bilinear', align_corners=True),
            nn.Conv2d(in_channels, out_channels, 3, padding=1),
            nn.PReLU()
        )

    def forward(self, x):
        return self.conv(x)

class UNetResNet34(nn.Module):
# class PSPNet(nn.Module):
    def __init__(self, num_classes=20, p = 0.0,sizes = (1, 2, 3, 6), psp_size = 512, deep_features_size = 256, backend = 'resnet34',
                 pretrained=False):
        # super(PSPNet, self).__init__()
        super(UNetResNet34, self).__init__()
        self.num_classes = num_classes
        self.feats = getattr(extractors, backend)(pretrained)
        self.psp = PSPModule(psp_size, 1024, sizes)
        self.drop_1 = nn.Dropout2d(p=0.3)

        self.up_1 = PSPUpsample(1024, 256)
        self.up_2 = PSPUpsample(256, 64)
        self.up_3 = PSPUpsample(64, 64)

        self.drop_2 = nn.Dropout2d(p=0.15)
        self.final = nn.Sequential(
            nn.Conv2d(64, num_classes, kernel_size=1),
            # nn.LogSoftmax()
        )

        self.classifier = nn.Sequential(
            nn.Linear(deep_features_size, 256),
            nn.ReLU(),
            nn.Linear(256, num_classes)
        )

    def forward(self, data_dict):
        x = data_dict['image']
        h, w = x.shape[2], x.shape[3]
        # padding
        min_size = 16
        pad_h = int((h + min_size - 1) / min_size) * min_size - h
        pad_w = int((w + min_size - 1) / min_size) * min_size - w
        if pad_h > 0 or pad_w > 0:
            # Pad 0 here. Not sure whether has a large effect
            x = F.pad(x, [0, pad_w, 0, pad_h])
        # assert h % 16 == 0 and w % 16 == 0

        preds = dict()

        f, class_f = self.feats(x)
        p = self.psp(f)
        p = self.drop_1(p)

        p = self.up_1(p)
        p = self.drop_2(p)

        p = self.up_2(p)
        p = self.drop_2(p)

        p = self.up_3(p)

        # crop
        if pad_h > 0 or pad_w > 0:
            p = p[:, :, 0:h, 0:w]

        seg_logit = self.final(p)

        preds['seg_logit'] = seg_logit

        preds['feature'] = p

        return preds

    def get_loss(self, cfg):
        from safnet.models.loss import SegLoss
        if cfg.TRAIN.LABEL_WEIGHTS_PATH:
            weights = np.loadtxt(cfg.TRAIN.LABEL_WEIGHTS_PATH, dtype=np.float32)
            weights = torch.from_numpy(weights).cuda()
        else:
            weights = None
        return SegLoss(weight=weights)

    def get_metric(self, cfg):
        from safnet.models.metric import SegAccuracy, SegIoU
        metric_fn = lambda: [SegAccuracy(), SegIoU(self.num_classes)]
        return metric_fn(), metric_fn()

if __name__ == '__main__':
    b, c, h, w = 2, 20, 120, 160
    image = torch.randn(b, 3, h, w).cuda()
    net = PSPNet(c, pretrained=False)
    net.cuda()
    preds = net({'image': image})
    for k, v in preds.items():
        print(k, v.shape)

max_stars_repo_name: lqzhao/SAFNet

max_issues_repo_name: lqzhao/SAFNet

max_forks_repo_name: lqzhao/SAFNet


text:
#!/usr/bin/env python3

"""A python script to perform watermark embedding/detection
   in the wavelet domain."""

# Copyright (C) 2020 by Akira TAMAMORI

# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

import numpy as np
from scipy.io import wavfile
from scipy.signal import windows
import pywt

HOST_SIGNAL_FILE = "bass_half.wav"  # 透かし埋め込み先のファイル
WATERMARK_SIGNAL_FILE = "wmed_signal.wav"           # 透かしを埋め込んだファイル
PSEUDO_RAND_FILE = 'pseudo_rand.dat'                # 疑似乱数列のファイル
WATERMARK_ORIGINAL_FILE = 'watermark_ori.dat'       # オリジナルの透かし信号

REP_CODE = True                 # 繰り返し埋め込みを使う
FRAME_LENGTH = 2048             # フレーム長
CONTROL_STRENGTH = 1000         # 埋め込み強度
OVERLAP = 0.5                   # フレーム分析のオーバーラップ率 (固定)
NUM_REPS = 3                    # 埋め込みの繰り返し数

WAVELET_BASIS = 'db4'
WAVELET_LEVEL = 3
WAVELET_MODE = 'symmetric'
THRESHOLD = 0.0

def fix(xs):
    """
    A emuration of MATLAB 'fix' function.
    borrowed from https://ideone.com/YjJwOh
    """

    # res = [np.floor(e) if e >= 0 else np.ceil(e) for e in xs]
    if xs >= 0:
        res = np.floor(xs)
    else:
        res = np.ceil(xs)
    return res

def embed():
    """
    perform embedding.
    """

    # ホスト信号
    sr, host_signal = wavfile.read(HOST_SIGNAL_FILE)
    signal_len = len(host_signal)

    # フレームの移動量
    frame_shift = int(FRAME_LENGTH * (1 - OVERLAP))

    # 隣接フレームとのオーバーラップ長
    overlap_length = int(FRAME_LENGTH * OVERLAP)

    # 埋め込みの総ビット数
    embed_nbit = fix((signal_len - overlap_length) / frame_shift)

    if REP_CODE:
        # 実効的な埋め込み可能ビット数
        effective_nbit = np.floor(embed_nbit / NUM_REPS)

        embed_nbit = effective_nbit * NUM_REPS
    else:
        effective_nbit = embed_nbit

    # 整数化
    frame_shift = int(frame_shift)
    effective_nbit = int(effective_nbit)
    embed_nbit = int(embed_nbit)

    # オリジナルの透かし信号を作成（0と1のビット列）
    wmark_original = np.random.randint(2, size=int(effective_nbit))

    # オリジナルの透かし信号を保存
    with open(WATERMARK_ORIGINAL_FILE, 'w') as f:
        for d in wmark_original:
            f.write("%d\n" % d)

    # 透かし信号を拡張する
    if REP_CODE:
        wmark_extended = np.repeat(wmark_original, NUM_REPS)
    else:
        wmark_extended = wmark_original

    # 透かしの埋め込み強度
    alpha = CONTROL_STRENGTH

    # 透かしが埋め込まれた信号(watermarked signal)を生成 in wavelet domain
    pointer = 0
    count = 0
    wmed_signal = np.zeros((frame_shift * embed_nbit))  # watermarked signal
    prev = np.zeros((FRAME_LENGTH))
    for i in range(embed_nbit):
        frame = host_signal[pointer: pointer + FRAME_LENGTH]

        # Wavelet係数を計算
        coeffs = pywt.wavedec(data=frame, wavelet=WAVELET_BASIS,
                              level=WAVELET_LEVEL, mode=WAVELET_MODE)

        # 透かしの埋め込み強度を平均と同じオーダーに設定する（adaptive）
        # coef_size = int(np.log10(np.abs(np.mean(coeffs[0])))) + 1
        # alpha = 10 ** coef_size

        # 透かしの埋め込み
        if wmark_extended[count] == 1:
            coeffs[0] = coeffs[0] - np.mean(coeffs[0]) + alpha
        else:
            coeffs[0] = coeffs[0] - np.mean(coeffs[0]) - alpha

        # 再構成
        wmarked_frame = pywt.waverec(coeffs=coeffs, wavelet=WAVELET_BASIS,
                                     mode=WAVELET_MODE)

        # 窓をかける (Hann window)
        wmarked_frame = wmarked_frame * windows.hann(FRAME_LENGTH)

        wmed_signal[frame_shift * i: frame_shift * (i+1)] = \
            np.concatenate((prev[frame_shift:FRAME_LENGTH] +
                            wmarked_frame[0:overlap_length],
                            wmarked_frame[overlap_length:frame_shift]))

        prev = wmarked_frame
        count = count + 1
        pointer = pointer + frame_shift

    # ホスト信号の残りと結合
    wmed_signal = np.concatenate(
        (wmed_signal, host_signal[len(wmed_signal): signal_len]))

    # 透かしが埋め込まれた信号をwavとして保存
    wmed_signal = wmed_signal.astype(np.int16)  # convert float into integer
    wavfile.write(WATERMARK_SIGNAL_FILE, sr, wmed_signal)

def detect():
    """
    perform detecton.
    """

    sr, host_signal = wavfile.read(HOST_SIGNAL_FILE)

    # 埋め込み済みの音声ファイルを開く
    _, eval_signal = wavfile.read(WATERMARK_SIGNAL_FILE)
    signal_len = len(eval_signal)

    # オリジナルの透かし信号をロード
    with open(WATERMARK_ORIGINAL_FILE, 'r') as f:
        wmark_original = f.readlines()
    wmark_original = np.array([float(w.rstrip()) for w in wmark_original])

    # フレームの移動量
    frame_shift = int(FRAME_LENGTH * (1 - OVERLAP))

    # 埋め込みビット数
    embed_nbit = fix((signal_len - FRAME_LENGTH * OVERLAP) / frame_shift)

    if REP_CODE:
        # 実質的な埋め込み可能ビット数
        effective_nbit = np.floor(embed_nbit / NUM_REPS)

        embed_nbit = effective_nbit * NUM_REPS
    else:
        effective_nbit = embed_nbit

    frame_shift = int(frame_shift)
    effective_nbit = int(effective_nbit)
    embed_nbit = int(embed_nbit)

    # オリジナルの透かし信号をロード
    with open(WATERMARK_ORIGINAL_FILE, 'r') as f:
        wmark_original = f.readlines()
    wmark_original = np.array([int(w.rstrip()) for w in wmark_original])

    # 透かし情報の検出
    pointer = 0
    detected_bit = np.zeros(embed_nbit)
    for i in range(embed_nbit):
        wmarked_frame = eval_signal[pointer: pointer + FRAME_LENGTH]

        # wavelet decomposition
        wmarked_coeffs = pywt.wavedec(
            data=wmarked_frame, wavelet=WAVELET_BASIS, level=WAVELET_LEVEL,
            mode=WAVELET_MODE)

        thres = np.sum(wmarked_coeffs[0])

        if thres >= THRESHOLD:
            detected_bit[i] = 1
        else:
            detected_bit[i] = 0

        pointer = pointer + frame_shift

    if REP_CODE:
        count = 0
        wmark_recovered = np.zeros(effective_nbit)

        for i in range(effective_nbit):

            # ビットを集計（平均値）
            ave = np.sum(detected_bit[count:count + NUM_REPS]) / NUM_REPS

            if ave >= 0.5:      # 過半数
                wmark_recovered[i] = 1
            else:
                wmark_recovered[i] = 0

            count = count + NUM_REPS
    else:
        wmark_recovered = detected_bit

    # ビット誤り率を表示
    denom = np.int(np.sum(np.abs(wmark_recovered - wmark_original)))
    BER = np.sum(np.abs(wmark_recovered - wmark_original)) / \
        effective_nbit * 100
    print(f'BER = {BER} % ({denom} / {effective_nbit})')

    # SNRを表示
    SNR = 10 * np.log10(
        np.sum(np.square(host_signal.astype(np.float32)))
        / np.sum(np.square(host_signal.astype(np.float32)
                           - eval_signal.astype(np.float32))))
    print(f'SNR = {SNR:.2f} dB')

    # bpsを表示
    print('BPS = {:.2f} bps'.format(embed_nbit / (len(host_signal) / sr)))

def main():
    """Main routine. """
    embed()
    detect()

if __name__ in '__main__':
    main()